Member, IEEE, Sami Muhaidat, Senior Member, IEEE Lina Mohjazi
email: l.mohjazi@ieee.org

Senior Member, IEEE Qammer H Abbasi
email: qammer.abbasi@glasgow.ac.uk

Senior Member, IEEE Muhammad Ali Imran
email: muhammad.imran@glasgow.ac.uk

Fellow, IEEE Octavia A Dobre
email: odobre@mun.ca

Fellow, IEEE Marco Di Renzo
email: marco.direnzo@centralesupelec.fr

M A Imran

Battery Recharging Time Models for Reconfigurable Intelligent Surface-Assisted Wireless Power Transfer Systems

Keywords: Reconfigurable intelligent surfaces (RIS), statistical models, wireless power transfer, recharging time. I

In this paper, we develop an analytical framework for the statistical analysis of the battery recharging time (BRT) in reconfigurable intelligent surfaces (RISs) aided wireless power transfer (WPT) systems. Specifically, we derive novel closedform expressions for the probability density function (PDF), cumulative distribution function, and moments of the BRT of the radio frequency energy harvesting wireless nodes. Moreover, closed-form expressions of the the PDF of the BRT is obtained for two special cases: i) when the RIS is equipped with one reflecting element (RE), ii) when the RIS consists of a large number of REs. Capitalizing on the derived expressions, we offer a comprehensive treatment for the statistical characterization of the BRT and study the impact of the system and battery parameters on its performance. Our results reveal that the proposed statistical models are analytically tractable, accurate, and efficient in assessing the sustainability of RIS-assisted WPT networks and in providing key design insights for large-scale future wireless applications. For example, we demonstrate that a 4-fold reduction in the mean time of the BRT can be achieved by doubling the number of RIS elements. Monte Carlo simulation results corroborate the accuracy of the proposed theoretical framework.

systems. As a result, the number of connected Internet-of-Everything (IoE) devices (e.g. sensors, wearables, implantables, tablets) is anticipated to witness a phenomenal growth in the next few years, reaching up to tens of billions [START_REF] Tariq | A speculative study on 6G[END_REF]. This poses a fundamental challenge on provisioning a ubiquitous seamless connectivity, while concurrently prolonging the lifetime of a massive number of energy-constrained low-power low-cost devices.

Wireless power transfer (WPT) has been highly recognized in both academia and industry as a promising technology to address the energy sustainability problem of wireless nodes, and has rapidly gained a growing interest in the research of B5G communication networks [START_REF] Varshney | Transporting information and energy simultaneously[END_REF], [START_REF] Grover | Shannon meets tesla: Wireless information and power transfer[END_REF]. This is mainly due to its capability to deliver on-demand wireless energy to a large number of wireless devices in a controllable and low-cost manner, and thus, eliminating the need for battery replacement. In this framework, radio frequency (RF) signals, that are received from dedicated wireless power transmitters, are leveraged either to recharge the batteries of wireless nodes or to directly power the wireless transmissions of battery-less devices [START_REF] Bi | Wireless powered communication: opportunities and challenges[END_REF], [START_REF] Mohjazi | RF-powered cognitive radio networks: technical challenges and limitations[END_REF].

However, it is demonstrated in recent studies [START_REF] Mohjazi | Performance analysis of SWIPT relaying systems in the presence of impulsive noise[END_REF]- [START_REF] Nasir | Relaying protocols for wireless energy harvesting and information processing[END_REF] that the distance between the RF transmitter and the corresponding RF energy harvesting (RFEH) receiver creates a performance bottleneck for practical wirelessly-powered wireless networks. This stems from the fact that the efficiency of WPT is inversely proportional to the distance, and hence, conventional relayaided wireless communications were proposed to realize WPT range expansion [START_REF] Mohjazi | Performance analysis of SWIPT relaying systems in the presence of impulsive noise[END_REF]- [START_REF] Nasir | Relaying protocols for wireless energy harvesting and information processing[END_REF]. More recently, advanced technologies, such as massive multiple-input multiple-output (MIMO), are studied as potential candidates to achieve considerable WPT efficiency gains through the exploitation of beamforming with or without relaying techniques [START_REF] Yang | Digital beamforming-based massive MIMO transceiver for 5G millimeter-wave communications[END_REF], [START_REF] Larsson | Joint beamforming and broadcasting in massive MIMO[END_REF]. Nonetheless, this comes at the price of severe energy consumption, higher computational complexity, and increased hardware cost, all of which are more pronounced at higher RF frequencies, such as millimeter-wave frequencies for future wireless systems [START_REF] Hao | Beamforming design in SWIPT-based joint multicastunicast mmWave massive MIMO with lens-antenna array[END_REF]. In addition to the distance factor, the efficiency of WPT, achieved through these conventional mechanisms, is degraded by the RF signal attenuation resulting from high penetration loss, multi-path fading, molecular absorption, and Doppler shift [START_REF] Liaskos | A new wireless communication paradigm through software-controlled metasurfaces[END_REF]. This effect is even more noticeable in ultra-dense network deployments, which feature highly dynamic radio environments.

To address the afore-mentioned challenges, reconfigurable intelligent surfaces (RISs) have recently emerged as a promising technology that can potentially offer fundamental performance improvements in wireless systems, in terms of spectrum and energy efficiencies, in a cost-effective manner [START_REF] Imran | Smart radio environments empowered by reconfigurable AI meta-surfaces: an idea whose time has come[END_REF]- [START_REF] Liaskos | A new wireless communication paradigm through software-controlled metasurfaces[END_REF], [START_REF] Mohjazi | An outlook on the interplay of AI and softwaredefined metasurfaces[END_REF]- [START_REF] Huang | Holographic MIMO surfaces for 6G wireless networks: Opportunities, challenges, and trends[END_REF]. Practically, RISs can be realized by different approaches. This includes (i) surfaces equipped with large arrays of discrete, inexpensive, and tiny antenna elements, called unit cells, (ii) implementations based on conformal large surfaces or metamaterial-based planar with scattering elements spaced apart at distances much smaller than the wavelength [START_REF] Renzo | Smart radio environments empowered by reconfigurable intelligent surfaces: How it works, state of research, and road ahead[END_REF]. Unlike conventional approaches that lack full control over the propagation environment, RISs allow a transformative control of electromagnetic (EM) waves [START_REF] Imran | Smart radio environments empowered by reconfigurable AI meta-surfaces: an idea whose time has come[END_REF], [START_REF] Wu | Towards smart and reconfigurable environment: Intelligent reflecting surface aided wireless network[END_REF], [START_REF] Renzo | Smart radio environments empowered by reconfigurable intelligent surfaces: How it works, state of research, and road ahead[END_REF]. For example, in RISs made of large arrays of inexpensive antennas, each element is individually configured and optimized to manipulate the impinging EM waves in arbitrary ways [START_REF] Dunna | ScatterMIMO: enabling virtual MIMO with smart surfaces[END_REF], [START_REF] Arun | RFocus: Beamforming using thousands of passive antennas[END_REF]. This is achieved by jointly manipulating the reflected signal amplitude and/or phase shift at each of the RIS elements in real time according to the dynamic and implicit randomness of wireless channels. For instance, the signal component arriving from an RF source node, and reflected by the RIS elements, can be steered towards an intended destination node to enhance its received signal power. Traditional transmission techniques, such as phased arrays, multi-antenna transmitters, and relays, involve active components of complex hardware that exhibit high power consumption. On the other hands, RISs require a large number of scattering elements, each of them is supported by the lowest number of small-sized, lowpower, and inexpensive components [START_REF] Dunna | ScatterMIMO: enabling virtual MIMO with smart surfaces[END_REF], [START_REF] Arun | RFocus: Beamforming using thousands of passive antennas[END_REF]. Details of the key similarities and differences between RIS and relays are provided in [START_REF] Ntontin | Reconfigurable intelligent surfaces vs. relaying: Differences, similarities, and performance comparison[END_REF].

A. Related Work

The opportunities opened by RIS have spurred, in a short span of time, research in many areas related to wireless communication systems. This includes multi-user resource allocation, beamforming optimization, design of efficient enabling mechanisms, and performance analysis of RIS-assisted wireless networks. For example, in the area of resource allocation, the authors in [START_REF] Huang | Reconfigurable intelligent surfaces for energy efficiency in wireless communication[END_REF] developed energy-efficient power allocation approaches subject to individual link budget guarantees for multiple mobile users. The study in [START_REF] Yang | IRS-enhanced OFDM: Power ation and passive array optimization[END_REF] proposed an achievable rate optimization framework for orthogonal frequency division multiplexing (OFDM) that jointly identifies the transmit power at the base station (BS) and the reflection coefficients at the RIS. Furthermore, several research works designed enabling mechanisms, such as channel estimation schemes in an effort to achieve the passive beamforming gains of RIS [START_REF] Nadeem | reflecting surface assisted multi-user MISO communication: Channel estimation and beamforming design[END_REF], [START_REF] Zappone | Overhead-aware design of reconfigurable intelligent surfaces in smart radio environments[END_REF]. Also, to reduce the overhead in channel training, the authors in [START_REF] Yang | Intelligent reflecting surface meets OFDM: Protocol design and rate maximization[END_REF] proposed a practical transmission protocol that involves estimating the combined channel of a group of RIS elements. Furthermore, the authors in [START_REF] Basar | Reconfigurable intelligent surface-based index modulation: A new beyond MIMO paradigm for 6G[END_REF] and [START_REF] Zhang | Space-time-coding digital metasurfaces[END_REF] investigated the realization of index modulation and space-time modulated digital coding, respectively, in RIS-assisted transmissions to improve the spectral efficiency. From the performance analysis point of view, recent research studies provided a theoretical framework to characterize the performance of RIS-assisted wireless systems in terms of outage probability [START_REF] Jung | Reliability analysis of large intelligent surfaces (LISs): Rate distribution and outage probability[END_REF], reflection probability [START_REF] Di Renzo | Reflection probability in wireless networks with metasurface-coated environmental objects: an approach based on random spatial processes[END_REF], spectral efficiency [START_REF] Jung | Spectral efficiency in large intelligent surfaces: Asymptotic analysis under pilot contamination[END_REF], [START_REF] Zhou | Spectral and energy efficiency of irs-assisted MISO communication with hardware impairments[END_REF], and capacity [START_REF] Zhang | Capacity characterization for intelligent reflecting surface aided MIMO communication[END_REF]- [START_REF] Perović | Channel capacity optimization using reconfigurable intelligent surfaces in indoor mmWave environments[END_REF]. Additionally, the fundamental limits of the error probability performance of RIS-aided backscatter and non-orthogonal multiple access (NOMA) were, respectively, examined in [START_REF] Zhao | Performance analysis of large intelligent surface aided backscatter communication systems[END_REF] and [START_REF] Thirumavalavan | BER analysis of reconfigurable intelligent surface assisted downlink power domain NOMA system[END_REF]. Meanwhile, numerous research studies focused on proposing active and passive beamforming strategies to achieve secrecy [START_REF] Chu | Intelligent reflecting surface aided multi-antenna secure transmission[END_REF] and sum-rate [START_REF] Zhang | Sum rate optimization for two way communications with intelligent reflecting surface[END_REF] enhancements and transmit power reduction [START_REF] Wu | Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming[END_REF]- [START_REF] Zhou | Robust beamforming design for intelligent reflecting surface aided MISO communication systems[END_REF] for multi-antenna and/or multi-user RIS-assisted networks.

While the preceding research works focus on leveraging the benefits of RIS to improve the performance of information transfer, RIS is shown to offer significant enhancements in the power transfer efficiency in wireless systems powered by near-field WPT [START_REF] Lang | Optimization of wireless power transfer systems enhanced by passive elements and metasurfaces[END_REF]. Apart from this, the deployment of RIS is demonstrated to enhance the far-field WPT and establish effective RFEH zones through compensating RF signals over long distances [START_REF] Wu | Weighted sum power maximization for intelligent reflecting surface aided SWIPT[END_REF]- [START_REF] Wu | Joint active and passive beamforming optimization for intelligent reflecting surface assisted SWIPT under QoS constraints[END_REF]. Specifically, the work in [START_REF] Wu | Weighted sum power maximization for intelligent reflecting surface aided SWIPT[END_REF] and [START_REF] Pan | Intelligent reflecting surface aided MIMO broadcasting for simultaneous wireless information and power transfer[END_REF] considered the weighted sum-power and sum-rate optimization problems, respectively, for RIS-aided simultaneous wireless information and power transfer (SWIPT) systems, where an RIS is deployed in the vicinity of two separate groups of energy and information receivers. Moreover, joint active and passive beamforming design for RIS-aided SWIPT is proposed in [START_REF] Wu | Joint active and passive beamforming optimization for intelligent reflecting surface assisted SWIPT under QoS constraints[END_REF] taking into consideration the signal-tointerference-plus-noise ratio constraints imposed by the information receivers (IRs) and RFEH constraints imposed by the energy receivers (ERs).

B. Motivation and Contribution

The afore-mentioned studies in [START_REF] Wu | Weighted sum power maximization for intelligent reflecting surface aided SWIPT[END_REF]- [START_REF] Wu | Joint active and passive beamforming optimization for intelligent reflecting surface assisted SWIPT under QoS constraints[END_REF] assume batteryfree RFEH ERs, whose harvested energy is directly used for future transmissions. In this case, the amount of the received RF signals, and consequently, the amount of harvested energy, is considered to be sufficient and predictable over a certain period of time. However, in scenarios where RFEH nodes are equipped with batteries [START_REF] Altinel | Statistical models for battery recharging time in rf energy harvesting systems[END_REF], the harvested energy is stored first in the battery before being used for future transmissions. Since the power of the received RF signal depends on the distribution of the probabilistic wireless fading channel between the transmitter and the receiver, the RFEH process and similarly, the time required to recharge the battery of an RFEH node, called the battery recharging time (BRT), become stochastic processes. The statistical characterization and modeling of BRT were developed in [START_REF] Altinel | Statistical models for battery recharging time in rf energy harvesting systems[END_REF]- [START_REF] Salahat | Statistical models for battery recharge time from RF energy scavengers in generalized wireless fading channels[END_REF] for conventional SWIPT systems operating over multi-path fading channels. Their results demonstrated that the RFEH process is significantly impacted by the system, fading, and battery parameters, including the capacity and discharging depth, and the BRT. Although the previously mentioned research work, carried out in [START_REF] Wu | Weighted sum power maximization for intelligent reflecting surface aided SWIPT[END_REF]- [START_REF] Wu | Joint active and passive beamforming optimization for intelligent reflecting surface assisted SWIPT under QoS constraints[END_REF], provide useful results to improve the RFEH process in RIS-assisted SWIPT systems, their approach focuses on developing transmission protocols and their insights are limited to battery-free ERs.

Despite being a fundamental figure of merit in designing and quantifying the sustainability of RIS-assisted WPT networks in various operational setups in B5G systems, to the best of our knowledge, the statistical characterization of BRT is not yet studied in the open literature. Motivated by this, the prime focus of this work is to develop a novel theoretical framework to characterize the statistical properties of BRT for RIS-assisted WPT systems, consisting of ERs with limited battery capacity. In our work, we consider that RIS comprises passive REs spaced half of the wavelength apart, and that each element is individually configured for realizing optimal WPT. The main contributions of this paper are listed in more details as follows:

• First, we derive novel accurate closed-form approximations for the probability density function (PDF) and cumulative distribution function (CDF) of the instantaneous total received power of the energy harvesting node. The obtained expressions take into account the number of RIS REs and the distances of the source (S)→RIS and RIS→ER links, and assume that all wireless channels are subject to Rayleigh fading. • Building upon the received power analytical expressions, we obtain closed-form expressions for the PDF and CDF of the BRT. The derived results are shown to accurately capture the impact of the battery parameters (i.e., discharge depth, battery capacity, and charging voltage), system parameters, and the number of RIS elements on the BRT statistical properties. To the best of our knowledge, these expressions are novel in literature. • Next, analytical closed-form expressions for the PDF of the BRT are derived for two special cases, namely 1) when the RIS is equipped with only one RE, which serves as a benchmark to quantify the gains obtained by increasing the number of RIS elements and 2) when the number of REs grows large. For the later case, we exploit the central limit theorem (CLT) to demonstrate that the PDF of the BRT converges to an impulse response, revealing that the deployment of RIS is highly promising in practically realizing WPT in future large-scale IoE networks. • To further investigate the properties of the BRT, we derive a computationally simple closed-form expression for its moments. We employ this result to obtain statistical tools to evaluate the mean value, variance, skewness and kurtosis of the BRT. • The derived moments expression is also exploited to examine the amount of fading (AoF) of the BRT as a function of RIS elements. Our study unveils that employing RISs of large size can effectively boost the efficiency of WPT over fading channels. • Finally, we present Monte Carlo simulation and numerical results to validate the accuracy of the developed theoretical framework.

C. Organization

The rest of the paper is organized as follows: Section II provides the RIS-assisted WPT system model, as well as the statistical characterization of the corresponding power received by the energy harvesting node. In Section III, the analytical expressions of the statistics of the BRT for RIS-assisted WPT systems are derived. Simulation and numerical results are presented in Section IV, while concluding remarks are given in Section V.

II. SYSTEM AND CHANNEL MODEL

In this paper, we consider a single-antenna RF source node, S, and a single-antenna energy-constrained ER, as depicted in Fig. 1. The ER could be a low-power sensor node equipped with a battery with a finite capacity. In order to extend the operational range of the ER, while ensuring that its harvested energy is sufficient for real-life operation, we propose WPT assisted by an RIS.

The end-to-end (E2E) channel gain between S and ER characterizes the power received at the ER, and accordingly, defines the behavior of the overall RFEH process, including the instantaneous BRT at ER. Therefore, to quantify the impact of RIS-assisted WPT on the required time to charge the battery of the ER node, in this section, we analytically present the distribution of the instantaneous received power at ER, which will be exploited next to develop the statistical characterization of the underlying instantaneous BRT.

We further assume that a direct link does not exist between S and ER. This is motivated by the fact that this link is subject to strong attenuation, due to deep fading, or shadowing effects, due to surrounding physical obstacles, or both and WPT can be achieved only via the RIS. It is worth mentioning that such an assumption is widely adopted in research studies related to WPT communication systems [START_REF] Mohjazi | Performance analysis of SWIPT relay networks with noncoherent modulation[END_REF], [START_REF] Pan | Intelligent reflecting surface aided MIMO broadcasting for simultaneous wireless information and power transfer[END_REF].

In our setup, we consider an RIS that is made of N REs. Each of the elements can be reconfigured by a communication oriented software through a controller, as illustrated in Fig. 1. The power transmitted from S, being either a BS or an RF source, which is reflected by the RIS towards the ER, is harvested and stored in a battery with a limited capacity before being used in future signal transmissions.

As shown in the block diagram of Fig. 1, h i and g i denote the small scale complex channel fading coefficients of the S →RIS and RIS→ER links, respectively, where i ∈ {1, 2, ..., N } denotes the i-th element of RIS. It is assumed that the envelopes of the two wireless links are modeled as independent and identically distributed (i.i.d) Rayleigh fading channels with scale parameter, σ, being equal to 1, i.e., |h i |, |g i |∼CN (0, 2σ 2 ) for i∈ {1, 2, ..., N }, where CN (0, κ) stands for a zero-mean complex Gaussian distribution with κ variance. The assumption of Rayleigh fading channels is representative of scenarios in which line-of-sight (LOS) propagation cannot be established due to random RIS deployments, e.g., if the RISs are deployed on spatial blockages. In this case, the system designer has no control over optimizing RISs locations [START_REF] Qian | Beamforming through reconfigurable intelligent surfaces in single-user MIMO systems: SNR distribution and scaling laws in the presence of channel fading and phase noise[END_REF].

Let P s denote the transmit wireless power of the source node. The instantaneous total power received at ER through the i-th element of RIS is expressed as where d 1 and d 2 represent the distance between S and the center of the RIS and between the center of the RIS and ER, respectively, and δ is the path loss exponent. Furthermore, θ i specifies the adjustable phase induced by the i-th RE of the RIS, respectively [START_REF] Basar | Wireless communications through reconfigurable intelligent surfaces[END_REF]. The total received power expression given in ( 1) is applicable in the far-field regime, as defined in [START_REF] Tang | Wireless communications with reconfigurable intelligent surface: Path loss modeling and experimental measurement[END_REF], [START_REF] Renzo | Analytical modeling of the path-loss for reconfigurable intelligent surfaces-anomalous mirror or scatterer?[END_REF]. Accordingly, N can be large but it needs to be finite [20, Sec. IV-D].

P r = N i=1 |h i ||g i |e jθi 2 d δ 1 d δ 2 P s , (1)
It is assumed that the channel phases of h i and g i , denoted respectively as φ hi and φ gi , are perfectly known to the RIS and accordingly, it is able to provide optimal phase shifting, i.e., θ i = -(φ hi + φ gi ). This idealized scenario sets a system operation performance benchmark for practical operations [START_REF] Basar | Wireless communications through reconfigurable intelligent surfaces[END_REF]. Therefore, the results obtained represent a lower bound on the BRT required to sustain the operation of RIS-assisted WPT systems. In this paper, we assume that the amplitude of the reflection coefficient of each RE is equal to 1, 1 for all i ∈ N . Consequently, the instantaneous received power is maximized and (1) can be formulated as

P r = P s d δ 1 d δ 2 B 2 (2)
where

B = N i=1 |h i ||g i | (3)
is the E2E channel gain.

A. Statistical Characterization of the E2E Channel

As previously mentioned, the BRT, T r , is determined by the amount of power received and then harvested at ER. Therefore, it is necessary to have in hand the statistical characterization of the E2E channel fading coefficient, B, in order to derive the distributions of P r and T r . Note that B presents a sum of N double Rayleigh random variables (RVs). An accurate approximation of its PDF is delivered in the following proposition.

Proposition 1: The PDF of the E2E channel coefficient of RIS-assisted WPT is accurately approximated in a closed-form 1 Recent studies proposed efficient designs achieving a reflection coefficient value as high as 1 [START_REF] Wu | Towards smart and reconfigurable environment: Intelligent reflecting surface aided wireless network[END_REF].

as f B (x) ≈ a 1 G 2,0 1,2 x a 2 -; a 3
a 3 = 4ϕ 4 -9ϕ 3 + 6ϕ 2 -µ 1 -ϕ 4 + 3ϕ 3 -3ϕ 2 + µ 1 , (5)
a 2 = a 3 2 (ϕ 4 -2ϕ 3 + ϕ 2 ) + 2ϕ 4 -3ϕ 3 + ϕ 2 , (6)
a 4 = a 6 + a 7 2 , ( (7)
) 8
a 5 = a 6 -a 7 2 , (9)
with

a 6 = a 3 (ϕ 2 -µ 1 ) + 2ϕ 2 -µ 1 a 2 -3, (10)
a 7 = a 3 (ϕ 2 -µ 1 ) + 2ϕ 2 -µ 1 a 2 -1 2 -4 µ 1 (a 3 + 1) a 2 , (11)
and ϕ i = µ j µ j-1 , j > 1. (12)
Above, µ j is the j-th moment of B and Γ(.) is the Gamma function defined in [56, Eq. (6.1.1)].

Proof: The RV B, given in (3), is written as a sum of the RVs |h i ||g i |. Therefore, its PDF can be formulated in a closed-form approximation using the moment-based density approximants method presented in [START_REF] Bouanani | Accurate closed-form approximations for the sum of correlated Weibull random variables[END_REF]. The evaluation of the first four moments, ensuring an accurate approximation for the PDF of B, are obtained in [START_REF] Bariah | Large intelligent surfaces-based non-orthogonal multiple access: Performance analysis[END_REF] and are presented here for clarity and completeness of the work. Specifically, the first, second, third, and fourth moments are expressed as ( 13), ( 14), [START_REF] Larsson | Joint beamforming and broadcasting in massive MIMO[END_REF], and [START_REF] Hao | Beamforming design in SWIPT-based joint multicastunicast mmWave massive MIMO with lens-antenna array[END_REF], respectively. This completes the proof.

µ 1 = N π 2 , (13)
µ 2 = 4 + (N -1) π 2 4 N, (14)
µ 3 =      N π 9 2 + 6(N -1) + (N -1)(N -2) π 2 8 , N ≥ 3 9π + 3 × 2 × π 2 × 4 = 21π, N = 2 9π 2 , N = 1. (15)
Based on the PDF obtained in (4), an accurate approximation of the CDF of B can be analytically computed as [START_REF] Bouanani | Accurate closed-form approximations for the sum of correlated Weibull random variables[END_REF]

F B (x) ≈ a 1 a 2 G 2,1 2,3 x a 2 1; a 3 + 1 a 4 + 1, a 5 + 1; 0 , x ≥ 0. ( 17
)
µ 4 =            64N + 48N (N -1) + 9N (N -1)π 2 + 6N (N -1)(N -2)π 2 + N (N -1)(N -2)(N -3)π 4 16 , N ≥ 4 480 + 90π 2 , N = 3 224 + 18π 2 , N = 2 64, N = 1. (16)

B. Statistical Characterization of the Received Power

Capitalizing on the statistical model of the E2E channel gain presented earlier, we derive the distribution of the instantaneous total received power at ER in the following proposition. Proposition 2: For an RIS-assisted WPT system, the CDF of the instantaneous power received at the ER node can be expressed as

F Pr (x) ≈ a 1 a 2 G 2,1 2,3 1 a 2
x Pr 1; a 3 + 1

a 4 + 1, a 5 + 1; 0 , x ≥ 0. (18)
Proof: The CDF of P s is given as

F Pr (x) = Pr(P r ≤ x). (19)
Substituting ( 2) in ( 19) yields

F Pr (x) = Pr B ≤ x Pr (20)
or equivalently

F Pr (x) = F B x Pr , (21)
where

Pr = P r d δ 1 d δ 2 ( 22
)
denotes the average power received at ER. By invoking the expression given in ( 17), (18) can be obtained, which concludes the proof.

Accordingly, the PDF of P s can be obtained by applying [59, Eq. (07.34.20.0001.01)] to differentiate [START_REF] Yang | A programmable metasurface with dynamic polarization, scattering and focusing control[END_REF], i.e.,

f Pr (x) = dF Pr (x) dx , (23)
yielding

f Pr (x) ≈ a 1 a 2 2x G 2,0 1,2 1 a 2
x Pr -; a 3 + 1

a 5 + 1, a 4 + 1; - , x ≥ 0. (24)

III. BATTERY RECHARGING TIME STATISTICAL MODELS

In this section, the statistical models developed in Section II are employed to derive analytical expressions for the statistical characterization of the instantaneous BRT, T r at the ER node.

Based on [START_REF] Altinel | Energy harvesting from multiple RF sources in wireless fading channels[END_REF], the BRT of the ER node is inversely proportional to its received power, P r and is defined as

T r = α P r , ( 25
)
where P r is given in [START_REF] Basar | Wireless communications through reconfigurable intelligent surfaces[END_REF]. Moreover, α denotes the conversion coefficient, which is a function of the battery and the RFEH circuit parameters; this can be expressed as follows

α = C b D d V b η , (26)
where C b is the battery capacity, D d is the battery discharge depth, V b is the charging voltage, and η is the RF to direct currect efficiency. In this section, we will exploit the PDF and CDF of the instantaneous received power, derived in ( 24) and ( 18), respectively, to obtain the statistical distribution of the BRT considering RIS-assisted WPT.

A. Probability Density Function (PDF) of BRT

The following proposition returns a closed-form expression for the PDF of the BRT when WPT is completed through an RIS.

Proposition 3: For RIS-assisted WPT systems, the PDF of the battery recharging time at ER node is given as

f Tr (τ ) ≈ a 1 a 2 2τ G 2,0 1,2 1 a 2 α Pr τ -; a 3 + 1 a 5 + 1, a 4 + 1; - , τ > 0. (27)
Proof: Using (25) and with the help of the Jacobian transformation method [START_REF] Papoulis | Probability, Random Variables, and Stochastic Processes[END_REF], the PDF of T r can be given as

f Tr (τ ) = α τ 2 f Pr α τ , (28)
which, with the aid of (24), can be expressed in a closed-form as [START_REF] Zappone | Overhead-aware design of reconfigurable intelligent surfaces in smart radio environments[END_REF]. This completes the proof.

It is worth noting that ( 27) is simple and incorporates the Meijer G-function, which is a standard built-in function in most of the well-known mathematical software packages, such as MATLAB, MAPLE, and MATHEMATICA, and can therefore, be efficiently evaluated.

B. Cumulative Distribution Function (CDF) of BRT

The CDF of the BRT is defined as the probability that the instantaneous BRT falls below a predetermined threshold, τ th , i.e., F Tr (τ th ) = P (τ ≤ τ th ). Taking into account [START_REF] Yang | IRS-enhanced OFDM: Power ation and passive array optimization[END_REF], it is straightforward to note that the relation between the CDFs of the received power and the BRT becomes

F Tr (τ th ) ≈ 1 -F Pr α τ th . (29)
Therefore, by substituting ( 18) into (29), we obtain the CDF of the BRT in a closed-form as

F Tr (τ th ) = 1 -a 1 a 2 G 2,1 2,3 1 a 2 α Pr τ th 1; a 3 + 1 a 4 + 1, a 5 + 1; 0 , τ th > 0. ( 30
)
f τr (τ ) = 2α τ 2 N (16 -π 2 ) Pr τ N 2 π 2 Pr 4α 1 4 exp - (4α/τ ) + N 2 π 2 Pr 2N (16 -π 2 ) Pr I -1 2 2π (16 -π 2 ) Pr Pr α τ , τ > 0. (40)

C. BRT Mean Value, Variance, Skewness, Kurtosis, and AoF

The n-th order moment of the BRT, denoted by µ T r (n), is a very useful statistical tool, as it enables the characterization of the mean value of the BRT, in addition to other underlying useful properties such as its skewness and kurtosis. Moreover, it can be employed to quantify the AoF, as will be elaborated next.

Having [START_REF] Zappone | Overhead-aware design of reconfigurable intelligent surfaces in smart radio environments[END_REF] in hand, one can derive the n-th moment of the BRT through the n-th order statistical expectation, as presented in the following proposition.

Proposition 4: The n-th moment of the BRT of an RISassisted WPT system can be expressed in a simple closed-form as

µ Tr (n) ≈ a 1 a (1-2n) 2 α Pr n Γ(a 4 + 1 -2n)Γ(a 5 + 1 -2n) Γ(a 3 + 1 -2n) . (31)
Proof: The n-th moment of T r can be evaluated in a straightforward manner by taking the statistical expectation

µ Tr (n) = ∞ 0 τ n f Tr (τ )dτ. (32)
Then, by substituting ( 27) in [START_REF] Di Renzo | Reflection probability in wireless networks with metasurface-coated environmental objects: an approach based on random spatial processes[END_REF], we obtain

µ Tr (n) = a 1 a 2 2 J ( 33
)
where

J = ∞ 0 τ n-1 G 2,0 1,2 1 a 2 α Pr τ -; a 3 + 1 a 5 + 1, a 4 + 1; - dτ.
(34) By applying [59, Eq. (07.34.21.0009.01)] to solve the integral, J , in a closed-form, the expression in (31) can be obtained. This completes the proof.

The expression in [START_REF] Jung | Reliability analysis of large intelligent surfaces (LISs): Rate distribution and outage probability[END_REF] can be used to obtain the mean value of the BRT, Tr by setting n = 1, i.e., Tr = µ Tr [START_REF] Imran | Smart radio environments empowered by reconfigurable AI meta-surfaces: an idea whose time has come[END_REF]. Also, the variance, σ 2 τr , can be computed as

σ 2 τr = µ Tr (2) -T 2 r . ( 35
)
The AoF parameter is viewed as a unified statistical measure of the severity of fading. Therefore, it is useful in quantifying the robustness of the RIS-assisted WPT technology against channel fading. The AoF is defined in [START_REF] Simon | Digital Communications over Fading Channels[END_REF]Eq.(1.27)] as the ratio of the variance to the square mean of the instantaneous received power. Subsequently, by employing [START_REF] Jung | Reliability analysis of large intelligent surfaces (LISs): Rate distribution and outage probability[END_REF], AoF = σ 2 τr / T 2 r . Remark: By examining (31), when n is set to 1, we note that the scaling law of the mean value of the BRT as a function of the RIS distances from the source and the ER nodes, respectively denoted as d 1 and d 2 , and which are implicitly defined through Pr in [START_REF] Arun | RFocus: Beamforming using thousands of passive antennas[END_REF], dictates that Tr increases with the square of the product of the distances d 1 and d 2 . This suggests that the minimum value of Tr is achieved when RIS is located either closer to the source or the ER node.

In addition to the mean and variance of the BRT, statistical properties such as skewness, denoted by , and kurtosis, denoted by Ψ, can be evaluated from [START_REF] Jung | Reliability analysis of large intelligent surfaces (LISs): Rate distribution and outage probability[END_REF] to provide a deeper insight of the distribution of the BRT. More specifically, the skewness measures the asymmetry of the PDF of the BRT about its mean value, while the kurtosis is an indicator of its peakedness or flatness and the heaviness of its tail. The skewness can be expressed as [START_REF] Papoulis | Probability, Random Variables, and Stochastic Processes[END_REF] = µ Tr (3)

µ 3/2 Tr (2) , ( 36
)
while the kurtosis can be given as

Ψ = µ Tr (4) µ 2 Tr (2) -3. ( 37
)
Special Case: In the special case when RIS consists of one RE (N = 1), B = |h| |g| is modeled as a double Rayleigh distribution. Accordingly, the PDF of the total received power at ER is given as [START_REF] Simon | Digital Communications over Fading Channels[END_REF] f

Pr (x) = 2 Pr K 0 2 x Pr , x ≥ 0. (38)
By substituting ( 38) in ( 28), the PDF of the BRT can be expressed as

f τr (τ ) = 2α Pr τ 2 K 0 2 α Pr τ , τ > 0. ( 39
)
Special Case: In the special case when the RIS consists of an asymptotically large number of REs (N 1), the following lemma returns a closed-form expression for the PDF of the BRT for WPT systems.

Lemma 1: For a sufficiently large number of RIS REs, N 1, the PDF of the BRT for RIS-assisted WPT system can be expressed in a closed-form as [START_REF] Chu | Intelligent reflecting surface aided multi-antenna secure transmission[END_REF], shown at the top of this page. Proof: Recalling that |h i | and |g i | are independently Rayleigh distributed RVs, then the mean and variance of B i are, respectively, given as

E[B i ] = E[|h i ||g i |] = π 2 (41)
and

VAR[B i ] = 16 -π 2 4 . ( 42
)
As the number of RIS elements becomes sufficiently large, then according to the CLT B converges to a Gaussian distribution with mean

µ N = E[B] = N π 2 ( 43
)
and variance

σ 2 N = N 16 -π 2 4 , ( 44
)
f Pr (x) = 2 N (16 -π 2 ) Pr N 2 π 2 Pr 4x 1 4 exp - 4x + N 2 π 2 Pr 2N (16 -π 2 ) Pr I -1 2 2π (16 -π 2 ) Pr Pr x , x ≥ 0 ( 45
)
where µ N and σ 2 N can be deduced by using ( 13) and ( 14). As a result, the instantaneous received power at ER, P s , defined in (2), presents a non-central chi-square distribution with one degree of freedom [START_REF] Proakis | Digital Communications[END_REF], and its PDF can be expressed as [START_REF] Lang | Optimization of wireless power transfer systems enhanced by passive elements and metasurfaces[END_REF], shown at the top of this page, where I v (.) is the modified Bessel function of order v defined in [START_REF] Abramowitz | Handbook of Mathematical Functions, With Formulas, Graphs, and Mathematical Tables[END_REF]Eq. (9.6.20)]. Finally, The expression in [START_REF] Chu | Intelligent reflecting surface aided multi-antenna secure transmission[END_REF] can be obtained by employing [START_REF] Lang | Optimization of wireless power transfer systems enhanced by passive elements and metasurfaces[END_REF] into [START_REF] Yang | IRS-enhanced OFDM: Power ation and passive array optimization[END_REF]. This completes the proof.

Based on Lemma 1 and using [START_REF] Chu | Intelligent reflecting surface aided multi-antenna secure transmission[END_REF], the mean value of the BRT when N 1 converges to

Tr = E[τ r ] = 4α (N 2 π 2 + N (16 -π 2 )) Pr . ( 46
)
By direct inspection of ( 46), we note that the mean value of the BRT in RIS-assisted WPT systems is shown to be inversely proportional to the square of the total number of REs, N , i.e., E[τ r ] ∝ 1 N 2 Pr . This result is in agreement with the BRT definition given in [START_REF] Yang | IRS-enhanced OFDM: Power ation and passive array optimization[END_REF].

IV. NUMERICAL AND SIMULATION RESULTS

In this section, numerical and Monte Carlo simulation results are presented to validate the accuracy of the proposed theoretical framework. This section also focuses on characterizing the properties of the BRT in RIS-assisted WPT wireless systems. The term Monte Carlo simulations refers to the use of actual fading channel variates with a number of repetitions of 10 6 trials.

Unless otherwise stated, the RFEH efficiency factor η = 0.5, as a worst case scenario, capturing the effects of low-cost hardware, and the total distance, d tot , between the source node, S, and the ER node is set to 5 m. In order to ensure far-field WPT, we assume that the size of RIS is relatively smaller than the transmission distance. It is assumed that the RIS is located mid-way between S and ER, i.e., d 1 = d 2 = d tot /2, and the path-loss exponent, δ = 2.7 [START_REF] Mohjazi | Performance analysis of SWIPT relaying systems in the presence of impulsive noise[END_REF]. Also, it is recalled that P s defines the total transmit power of the system. All simulation parameters, including the ER battery parameters, are summarized in Table 1. In Fig. 2, the PDF of the BRT for RIS-assisted WPT systems is illustrated for N RIS elements, using the analytical expression given in [START_REF] Zappone | Overhead-aware design of reconfigurable intelligent surfaces in smart radio environments[END_REF], where N = 4, 6, 8, 10, 50. Additionally, the PDF of the BRT obtained through the CLT, as given in [START_REF] Chu | Intelligent reflecting surface aided multi-antenna secure transmission[END_REF], is also shown for N = 10 and 50. It is observed that the simulation results are in full agreement with the analytical PDF curves, reflecting the accuracy of our proposed mathematical model and its effectiveness in capturing the statistical properties of the BRT. Additionally, we notice that as N increases, the BRT value of the ER node decreases. In more details, the highest probable value of the BRT drops from 0.64 to 0.2 hr (38.4 to 12 mins) when the number of REs increases from 4 to 8. Finally, it is evident that as the number of RIS elements increases through the CLT, the PDF converges to an impulse function, indicating that increasing N causes a significant reduction in the BRT of WPT systems, due to the enhanced spatial diversity gain of RIS-assisted systems, and thus, making them particularly attractive for large-scale RFEH applications. It is worth mentioning that the effect of increasing the number of RF sources on the BRT was studied in [START_REF] Altinel | Energy harvesting from multiple RF sources in wireless fading channels[END_REF]. Although their results demonstrated that the BRT decreases notably by increasing the number of RF sources, this comes at the cost of extra transmission power needed for each additional RF source. In Fig. 3, we investigate the effect of varying the value of the transmit power, P s , on the BRT performance of RIS-assisted WPT. We plot the PDF of the BRT when N = 1, 2, and 4, using the analytical expressions obtained in [START_REF] Zappone | Overhead-aware design of reconfigurable intelligent surfaces in smart radio environments[END_REF] and [START_REF] Thirumavalavan | BER analysis of reconfigurable intelligent surface assisted downlink power domain NOMA system[END_REF]. Figs. (3a), (3b), and (3c) present, respectively, low, moderate, and high transmit power regimes, i.e., P s =7, 15, and 40 dBm. By closely inspecting these figures, it is clearly observed that when higher transmit power values are encountered, the variance of the distribution of the BRT decreases while its kurtosis increases (i.e., sharper PDF peak), indicating the improvement in the BRT predictability. Interestingly, it can be also noted that even in the low transmit power scenario, depicted in Fig. 3a, the mean value of the BRT exhibits significantly promising in enhancing the reliability of the WPT process, and may be attained by adding low-cost passive RIS elements (featuring no transmit power consumption) instead of increasing the source transmission power. Further details about this aspect will be provided in subsequent discussions.

To gain more insights about the effect of varying the transmit power on the statistical distribution of the BRT in RIS-assisted WPT systems, we provide in Fig. 4 the CDF as a function of the BRT threshold, τ th , of the RIS-assisted system using the expression in [START_REF] Basar | Reconfigurable intelligent surface-based index modulation: A new beyond MIMO paradigm for 6G[END_REF]. The examination is carried out for different N values and assuming two transmit power scenarios, namely low-(P s = 7 dBm) and high-(P s = 30 dBm) transmit power. The excellent fit between the simulation and the analytical results verify the accuracy of our developed theoretical framework. As expected, for a fixed N , as the τ th increases, the CDF value increases. For example, for N = 4 and P s = 7 dBm, as τ th changes from 5 to 10 hrs, the CDF value approximately doubles. Additionally, for a given τ th value, as N increases, the CDF value improves. This indicates that the efficiency of the RFEH process is remarkably improved in an RIS-assisted WPT system by increasing the number of REs, N . For instance, at P s = 30 dBm, τ th can be reduced by about 13 hrs by employing an RIS with 2 elements instead of 1 in order to achieve a targeted probability of charging of 0.9. Fig. 4 also demonstrates that the steepness of the CDF curve increases as the value of the transmit power shifts from low to high. This indicates that adding more RIS elements is more rewarding in the low transmit power than the high transmit power case, where the reduction in the BRT value is more significant. Fig. 5 depicts the mean value of the BRT, Tr , as a function of the total transmit power, P s , for different values of N . We also examine the convergence of the mean value, obtained analytically in [START_REF] Jung | Reliability analysis of large intelligent surfaces (LISs): Rate distribution and outage probability[END_REF], towards that obtained through the CLT, given in [START_REF] Wu | Weighted sum power maximization for intelligent reflecting surface aided SWIPT[END_REF]. It can be deduced from Fig. 5 that Tr converges to that computed via the CLT even at relatively small values of N ; for example, when N > 8. This suggests that in such setups, the mean value of the BRT in RIS-assisted WPT systems can be mathematically evaluated using the simpler expressions in (46) instead of [START_REF] Jung | Reliability analysis of large intelligent surfaces (LISs): Rate distribution and outage probability[END_REF], which is computationally more demanding, since it involves the evaluation of the parameters a 1 through a 5 . We further notice that as N and P s increase, the mean value of the BRT linearly decreases. This finding is in agreement with the theoretical result obtained in [START_REF] Wu | Weighted sum power maximization for intelligent reflecting surface aided SWIPT[END_REF]. In more details, for a given value of P s , doubling the number of deployed RIS elements, N , yields a reduction in the mean value of the BRT by about 4 times. In summary, the choice of N depends on whether the application scenario operates in the low or high transmit power scenarios.

Remark: We emphasize that the analytical expressions derived in [START_REF] Jung | Reliability analysis of large intelligent surfaces (LISs): Rate distribution and outage probability[END_REF] and (46) lend themselves as effective tools in determining the minimum number of RIS elements that should be deployed in order to achieve a feasible BRT value for a targeted energy efficiency in WPT systems, while avoiding the need of unnecessary phase adjustments to account for the extra deployed RIS elements.

To quantify the impact of the RIS location on the BRT performance, in Fig. 6, we inspect the behavior of the mean value of the BRT, Tr , as a function of the normalized S →RIS distance, d 1 /d tot , for different values of N and fixed P s =20 dBm. We set the normalized RIS→ER distance to d 2 /d tot = 1 -d 1 /d tot . It is noted that for a given RIS location, the energy efficiency of a WPT system can be enhanced by equipping the RIS with more REs. This key BRT performance insight is indeed beneficial when there is no flexibility in choosing the location of the RIS due to the layout or geometry of the application environment. Additionally, it is straightforwardly observed that the minimum mean value of BRT is attained when the RIS is located either closer to S or the ER node. This verifies our theoretically proven result, as discussed in Sec. III-C following [START_REF] Jung | Reliability analysis of large intelligent surfaces (LISs): Rate distribution and outage probability[END_REF]. Therefore, it is not surprising to note that adding RIS elements is more rewarding, with respect to Tr , when RIS is located mid-way between the S and ER nodes.

Remark: The investigation carried out above would enable system design engineers to extract fundamental performance insights to choose the optimal placement of the RIS and/or IoE nodes in large-scale WPT systems.

To address the effect of varying the battery parameters, in Fig. 7 we plot the mean value of the BRT, Tr , with respect to the battery capacity, C b , of the ER node when the value of the transmit power is fixed to 20 dBm. As expected, for a given N value, it is observed that as C b increases, the mean value of the BRT increases as well. However, one can take the advantage of adding RIS elements to dramatically boost the sustainability of WPT systems without changing the battery parameters of the RFEH device. In other words, this figure reveals that, independent of C b , as N doubles in value, Tr decreases by about 5 times. For example, for a given C b = 10 mAh, as N is varied from 5 to 10 to 20, the mean value of the BRT is reduced from 1 hr to 12.65 mins to 2.86 mins. Therefore, it turns out that our proposed mathematical tools are useful in identifying the design parameters of IoE devices for future RIS-assisted WPT applications. This can be of considerable advantage for applications where the battery capacity implies constraints on the shape and size of the RFEH nodes, such as wearables or implantables.

Finally, in Fig. 8, we examine the AoF for RIS-assisted WPT as a function of the number of REs, N . We also plot the AoF for single-input single-output (SISO) WPT as a benchmark. Fig. 8 clearly shows that our analytical framework is useful in quantifying the robustness of RIS-assisted WPT to fading channels. Particularly, Fig. 8 reveals that increasing the number of REs reduces the AoF of the end-to-end channel, thereby improving the efficiency of WPT undergoing fading.

V. CONCLUSION

In this paper, we developed a theoretical framework to investigate the energy sustainability of RIS-assisted WPT systems, from the BRT perspective of an RFEH node. In particular, assuming that WPT is completed over Rayleigh fading channels, we provided the statistical characterization of the instantaneous received power of the system, including its PDF and CDF. Based on that, we derived novel lowcomplexity tight closed-form approximations for the PDF, CDF, and moments of the BRT as functions of the received power, battery parameters, and number of RIS REs. Additionally, using the CLT, we derived closed-form expressions for the PDF and mean value of the BRT considering that RIS is equipped with a large number of elements. Besides being accurate and mathematically tractable, our results reveal that the proposed statistical tools provide an efficient means to evaluate RIS-assisted WPT systems and extract useful design insights. For example, our results show that doubling the number of RIS elements improves the predictability of the BRT of the RFEH nodes in the network and offers a 4-fold reduction in its mean value. Moreover, it is reported that the characteristics of the BRT are significantly impacted not only by the system parameters, such as the distance between the nodes, but also by the battery parameters of the RFEH node, such as the battery capacity. Finally, our results illustrated that significant performance gains in the BRT have been observed

by locating the RIS close to the source or to the RFEH node.

Figure 1 :

 1 Figure 1: RIS-assisted WPT system model.

a 4 ,

 4 a 5 ; -, x ≥ 0 (4) where G .,. .,. [.|.] denotes the Meijer G-function defined in [55, Eq. (8.2.1.1)] and a 1 = Γ(a 3 + 1) a 2 Γ(a 4 + 1)Γ(a 5 + 1) ,

Figure 2 :

 2 Figure 2: The PDF of the BRT of an RIS-assisted WPT system, for different values of N , where the transmit power, Ps = 15 dBm. The PDFs obtained through the CLT are also plotted for N = 10 and 50.

Figure 3 :Figure 4 :

 34 Figure 3: The PDF of the BRT in RIS-assisted WPT systems, in (a) low transmit power, Ps = 7 dBm, (b) moderate transmit power, Ps = 15 dBm, and (c) high transmit power, Ps = 40 dBm, regimes, and for different values of N .

Figure 5 :Figure 6 :

 56 Figure 5: Mean value of the BRT as a function of the transmit power in RIS-assisted WPT systems, for different values of N .

Figure 7 :

 7 Figure 7: Mean value of the BRT as a function of the battery capacity of the ER node in RIS-assisted WPT systems, for different values of N and when Ps = 20 dBm.

Figure 8 :

 8 Figure 8: AoF as a function of N when Ps = 20 dBm.

Table I :

 I Simulation Parameters.

	Name	Symbol	Value
	RFEH efficiency of the ER circuit	η	0.5
	Normalized S →RIS distance	d 1 /dtot	0.5
	Normalized RIS→ER distance	d 2 /dtot	0.5
	Path-loss exponent	δ	2.7
	Battery capacity	C b	10 mAh
	Discharge depth	D d	0.4
	Battery charging voltage	V b	1.2 V

The work of M. Di Renzo was supported by the European Commission through the H2020 ARIADNE Project under Grant 871464.


A Valdenaire

Alaa E Giba

M Stoffel

X Devaux

J M Poumirol

C Bonafos

S Guehairia

R Demoulin

E Talbot

M Vergnat

H Rinnert

Heavily doped Si-nanocrystals formed in P-(SiO/SiO2) multilayers: a novel route to infrared Si-based plasmonics

As building blocks of novel multifunctional materials involving coupling at the nanoscale, highly doped semiconductor nanocrystals are of great interest for potential applications in nanophotonics. In this work, we investigate the plasmonic properties of highly doped Si nanocrystals embedded in a silica matrix. These materials are obtained by evaporation of heavily Phosphorus-doped SiO/SiO2 multilayers in a ultrahigh vacuum chamber followed by rapid thermal annealing. For P contents between 0.7 and 1.9 at%, structural investigations at the nanoscale give clear evidence that P atoms are mainly located in the core of Si nanocrystals with concentrations reaching up to 10 at%, i.e. well beyond the solid solubility limit of P in bulk Si. Alloying and formation of SiP nanoparticles is observed for P contents exceeding 4 at% in the multilayer. Infrared absorption measurements give evidence of a localized surface plasmon resonance located in the 3 to 6 µm range. A core-shell structure was used to model Si nanocrystals embedded in a silica matrix. Based on the Mie theory and the Drude model, both the mobility and the free charge carrier density were extracted from the simulation, with values reaching 27 cm 2 V -1 s -1 and 2.3×10 20 cm -3 , respectively. This results in a dopant activation rate of about 8 %.

Localized surface plasmon resonance (LSPR) corresponds to the collective oscillations of free charge carriers confined in a nanoparticle driven by an electromagnetic field. Many different plasmonic materials have been considered till now in the literature. Among them, noble metals such as Au or Ag nanostructures have first gained a significant interest. Indeed, their large free charge carrier concentration (10 22 -10 23 cm -3 ) is known to lead to a plasmon resonance in the visible spectral range. The coupling between the plasmon and the electric field of the incoming light leads to a significant enhancement of the absorption or fluorescence which has been exploited for various applications such as detection of single molecules [START_REF] Hoffman | Infrared reflection absorption spectroscopy of adsorbed molecules[END_REF] and more generally label-free biosensing [START_REF] Unser | Localized surface plasmon resonance biosensing: current challenges and approaches[END_REF] . However, due to the fact that Au and Ag-based nanostructures are not compatible with the Complementary Metal Oxide Semiconductor (CMOS) technology, and with the aim to extent the spectral range in the infrared range with nano-sized structures, alternative materials like heavily doped quantum dots and metal oxide nanocrystals have been explored [START_REF] Comin | New materials for tunable plasmonic colloidal nanocrystals[END_REF][START_REF] Liu | Heavily doped colloidal semiconductor and metal oxide nanocrystals: an emerging new class of plasmonic nanomaterials[END_REF] . By reducing the charge carrier concentration, the plasmon resonance can be possibly shifted to the near and mid-infrared which may further extend the range of potential applications [START_REF] Zhong | Review of mid-infrared plasmonic materials[END_REF] . Many different materials such as silicides [START_REF] Maex | Silicides for integrated circuits: TiSi2 and CoSi2[END_REF] , germanides [START_REF] Krusin-Elbaum | Unusually low resistivity of copper germanide thin films formed at low temperatures[END_REF] but also conducting oxides and nitrides as well as transition metal nitrides [START_REF] Naik | Oxides and nitrides as alternative plasmonic materials in the optical range[END_REF][START_REF] Naik | Demonstration of Al:ZnO as a plasmonic component for near infrared metamaterials[END_REF][START_REF] Naik | Titanium nitride as a plasmonic material for visible and near infrared wavelengths[END_REF] were shown to exhibit a strong potential for near-infrared plasmonics. Alternatively, doped semiconductor nanocrystals may also be of great interest [START_REF] Agrawal | Localized surface plasmon resonance in semiconductor nanocrystals[END_REF] . Several studies have concerned II-VI semiconductor nanocrystals such as Cu2-xS, Cu2-xSe, Cu2-xTe [START_REF] Luther | Localized surface plasmon resonances arising from free carriers in doped quantum dots[END_REF][START_REF] Kriegel | Tuning the excitonic and plasmonic properties of copper chalcogenide nanocrystals[END_REF][START_REF] Kriegel | Shedding light on vacancy doped copper chalcogenides: shape-controlled synthesis, optical properties and modelling of copper telluride nanocrystals with near-infrared plasmon resonances[END_REF] but also ZnO [START_REF] Schimpf | Charge tunable quantum plasmons in colloidal semiconductor nanocrystals[END_REF] and ITO [START_REF] Li | Infrared plasmonics with indium tin oxide nanorod arrays[END_REF] . Due to their low cost, their compatibility with the CMOS technology, their ease of fabrication and their non-toxicity, Si nanocrystals (Si-NCs) appear particularly promising. Although some theoretical works suggested that insertion of dopants at the nanoscale is not favored by thermodynamics [START_REF] Dalpian | Self-purification in semiconductor nanocrystals[END_REF] , several authors have nevertheless shown that Si-NCs can contain P atoms or B atoms [START_REF] Khelifi | Efficient n-type doping of Si nanocrystals embedded in SiO2 by ion beam synthesis[END_REF][START_REF] Gnaser | Phosphorus doping of Si nanocrystals embedded in silicon oxynitride determined by atom probe tomography[END_REF][START_REF] Trad | Influence of phosphorus on the growth and the photoluminescence properties of Si-NCs formed in P-doped SiO/SiO2 multilayers[END_REF] . Moreover, hyperdoping has emerged as a possible route to generate LSPR in Si-NCs. Using tight binding calculations, Pi and Delerue [START_REF] Pi | Tight binding calculations of the optical response of optimally P-doped Si nanocrystals: a model for localized surface plasmon resonance[END_REF] have investigated the optical response of hyperdoped Si-NCs.

For Si-NCs having sizes of about 4 nm, a LSPR can be possibly obtained when the nanocrystals contain more than 10 P atoms. Rowe et al. [START_REF] Rowe | Phosphorus doped Si nanocrystals exhibiting mid-infrared localized surface plasmon resonance[END_REF] gave the first experimental evidence of a tunable mid-infrared LSPR in free-standing phosphorus doped Si-NCs. The latter were obtained by a nonthermal plasma method and P was supplied by PH3 gas. The plasmonic properties of both heavily phosphorus and boron doped Si-NCs prepared by non-thermal plasma were further studied, reporting different behaviors depending on whether the Si-NCs are doped with boron or phosphorus [START_REF] Kramer | Plasmonic properties of Si nanocrystals doped with boron and phosphorus[END_REF][START_REF] Zhou | Comparative study of the localized surface plasmon resonance of Boron-and Phosphorus doped Si nanocrystals[END_REF][START_REF] Zhang | Doped silicon nanocrystal plasmonics[END_REF] . Poumirol et al. [START_REF] Poumirol | Hyper-doped silicon nanoantennas and metasurfaces for tunable infrared plasmonics[END_REF] used a combination of e-beam lithography and reactive ion etching to fabricate ordered arrays of P-hyperdoped Si nanodisks on a Silicon-On-Insulator substrate. By adjusting the free carrier concentration, the tunability of the LSPR was demonstrated in a large spectral window ranging from 2 μm to 5 μm. In addition, numerical simulations based on the Green Dyadic Method (GDM) were used to calculate the extinction efficiency of a single Si-nanodisk, which is found to be in good agreement with the experimental LSPR. The reflectance was further studied as a function of the spacing between nanodisks. Majorel at al. [START_REF] Majorel | Theory of plasmonic properties of hyper-doped silicon nanostructures[END_REF] used the same kind of simulations to investigate doped silicon nanostructures having an arbitrary shape. Both the extinction coefficient as well as the electromagnetic field in and around any kind of nanostructure can be calculated precisely thus enabling the control of light in and around tiny structures. Very recently, Zhang et al. [START_REF] Zhang | Infrared nanoplasmonic properties of hyperdoped embedded Si nanocrystals in the few electrons regime[END_REF] reported the observation of LSPR in a single layer of hyperdoped Si-NCs, with characteristic sizes lower than 6 nm, prepared by low energy ion implantation in a silica matrix and rapid thermal annealing. They demonstrated than only few electrons in a nanocrystal can give rise to LSPR. Till now, only the plasmonic properties of heavily doped, free-standing Si-NCs as well as single layers of hyperdoped Si-NCs embedded in a silica matrix have been considered. In the case of free-standing Si-NCs, the plasmonic properties turned out to be sensitive to oxidation [START_REF] Kramer | Plasmonic properties of Si nanocrystals doped with boron and phosphorus[END_REF] and more generally to the modification of the chemical environment [START_REF] Zhou | Comparative study of the localized surface plasmon resonance of Boron-and Phosphorus doped Si nanocrystals[END_REF] , which could be detrimental to applications. Moreover, in the case of a single layer of hyperdoped Si-NCs embedded in a silica matrix, the absorption remains weak [START_REF] Zhang | Infrared nanoplasmonic properties of hyperdoped embedded Si nanocrystals in the few electrons regime[END_REF] .

In this paper, we report on a tunable LSPR from heavily P-doped Si-NCs formed in P-(SiO/SiO2) multilayers. The multilayer approach has been successfully applied to achieve highly luminescent and size-controlled Si-NCs [START_REF] Zacharias | Size-controlled highly luminescent Si nanocrystals: a SiO/SiO2 superlattice approach[END_REF][START_REF] Jambois | Photoluminescence and electroluminescence of size-controlled Si nanocrystallites embedded in SiO2 thin films[END_REF] . To date, there is no report in the literature concerning heavily doped Si-NCs showing LSPR obtained using the well-known multilayer architecture. This approach offers several advantages: firstly, thanks the SiO2 barrier layers, the size of the Si-NCs can be controlled by the SiO layer thickness, and secondly, several layers of heavily doped Si-NCs can contribute to the LSPR thereby improving the signal. In a first part, the structure and the local chemical composition of the multilayers are determined by combining atom probe tomography (APT), advanced scanning transmission electron microscopy (STEM) and electron energy loss spectroscopy (EELS) mapping at the atomic scale.

In a second part, the absorption properties measured by Fourier-transform infrared absorption spectroscopy (FTIR) are presented. Simulations based on the theory of Mie and assuming a core-shell model were able to reproduce correctly both the vibrational and the electronic contributions to the absorption spectra. They allow us to determine both the mobility and the density of the charge carriers from which the activation rate of the dopants inserted into Si-NCs were deduced. P-doped SiO/SiO2 multilayers were grown by alternating evaporation of both SiO and SiO2

provided by e-beam guns on Si(100) substrates held at room temperature. The deposition rates of SiO and SiO2 were controlled by a quartz microbalance and fixed at 0.1 nm/s. The SiO/SiO2 layer thicknesses were kept at 10 nm and 6 nm, respectively. The multilayers consist of 20 doublets of SiO/SiO2 with a total nominal thickness of 320 nm. P atoms were supplied by a GaP decomposition source [START_REF] Lippert | Heavy phosphorus doping in molecular beam epitaxial grown silicon by a GaP decomposition source[END_REF] held at a constant temperature during the whole deposition process thus leading to incorporation of P atoms in both SiO and SiO2 layers. As-deposited multilayers were amorphous. The samples were then post-grown annealed in a rapid thermal annealing furnace at various temperatures for 5 or 15 minutes under N2 flow (200 sccm) to destabilize the SiO layer and to form Si-NCs. The microstructural properties were investigated by energy filtered transmission electron microscopy (EFTEM), STEM, STEM-EELS and STEM energy dispersive spectroscopy (EDS). Electron microscopy and the associated spectroscopies were carried out on a probe corrected JEOL ARM 200 F equipped with a cold-field emission gun.

EELS spectrum images (SI) were recorded with a GATAN GIF quantum ER filter. To avoid sample destruction caused by the electron beam, all experiments were carried out at 80 kV with a maximal probe current of about 50 pA for EELS spectrum images and EDS mapping. Two EELS-SI images were simultaneously recorded for the low loss part containing the zero loss and for the core-loss which allows advanced data post-processing like energy drift and multiple scattering corrections. The spectra were recorded with an electron beam having a half convergence angle of 24 mrad and a half collection angle of 56 mrad with an energy dispersion of 0.05 eV per channel and a pixel time of 0.005 s. SI were recorded with a spatial sampling of about 0.2 x 0.2 nm 3 per pixel. A multivariate statistical analysis software was used to improve the quality of the STEM-EELS data by de-noising the core-loss SI. Low loss spectra were processed by zero-loss extraction and Fourier-log deconvolution for multiple scattering corrections. To obtain chemical maps, multiple linear least square fitting was applied on SiL and PL signals. Due to the overlapping of both SiL and PL edges and to separate different chemical states of Si (Si 0 and Si IV ), reference spectra were recorded on pure Si, on thermal silica grown on pure Si and on SiP2. The reference spectra were extracted from SI recorded on lamellae with the same thickness and with the same STEM and spectrometer parameters than that used for the sample under study. More detailed information on this procedure was previously published [START_REF] Geiskopf | Formation of SiP2 nanocrystals embedded in SiO2 from phosphorus-rich SiO1,5 thin films[END_REF][START_REF] Geiskopf | Thin films of SiP lamellar alloys: a first step toward 2D SiP[END_REF] . For STEM observations, lamellas thinned up to t/λ <0.8 (ratio between the thickness and the mean free path of primary electrons) were prepared by focused ion beam (FIB) using a FEI Helios Nanolab 600i dual beam microscope. During FIB operation, the multilayers were capped by a pure carbon layer prior to deposition of a Pt-C layer to protect the sample during the first steps of the lamella extraction. After FIB milling, the lamellae were stored in vacuum to prevent oxidation. Atom probe tomography (APT) was further used to investigate the structural and chemical properties at the nanoscale. The experiments were carried out using a laser assisted wide angle tomographic atom probe (LAWATAP, Cameca) using a femtosecond UV pulsed laser (λ=343 nm, 350 fs, 33 nJ). The analysis was performed in an ultra-high vacuum chamber (10 -10 mbar) at T = 80 K. Tip shaped samples were prepared using a Thermoscientific plasma FIB Helios G4 XCe. The data treatments and the 3D reconstructions were performed by using GPM3D software. The infrared absorption spectra were obtained using a Fourier transform infrared (FTIR) absorption spectrometer. The spectra of the thin films were obtained after subtraction of the reference spectrum measured on a Si substrate. The steady state photoluminescence (PL) was excited by an arc source xenon lamp.

A monochromator was used to select an excitation centered at 325 nm with a full width at half maximum of 5 nm. The signal was collected by a monochromator equipped with a 150 grooves per mm grating and measured by a cooled CCD camera operating in the visible-near infrared range. All spectra were corrected from the response of the system.

Figure 1 shows the evolution of the chemical composition for P-doped SiO/SiO2 multilayers grown with different temperatures of the GaP decomposition source ranging between 780°C and 860°C. The elemental concentrations were determined by energy dispersive spectroscopy measurements performed on the whole multilayer structure. The obtained values were further cross-checked by atom probe tomography. The P content in the multilayers can be simply adjusted by varying the temperature of the GaP decomposition source. The P content was thus controlled from about 0.3 at% up to 4.1 at% when the GaP decomposition source temperature is increased from 780°C up to 860°C. Figure 2(a) shows an energy filtered transmission electron microscopy (EFTEM) image of a P-doped SiO/SiO2 multilayer containing 1.1 at% of P annealed at 1100°C for 15 minutes. The image is obtained by selecting the silica plasmon peak at 25 eV for imaging (slit width 1 eV). The silica matrix appears in light grey while Si and PtC layer are darker. One can clearly identify a rather high density of nanoparticles almost randomly dispersed in the silica matrix. A high resolution TEM micrograph of a single nanoparticle is shown in Fig. 2(b) with a fast Fourier transform in the inset. They show that the nanoparticles are mainly single crystalline with the Si diamond-like cubic structure. The expected multilayer structure is no longer fully preserved as a result of the P-induced softening of the silica matrix, as already mentioned in a previous work [START_REF] Trad | Influence of phosphorus on the growth and the photoluminescence properties of Si-NCs formed in P-doped SiO/SiO2 multilayers[END_REF] or by the increase of atomic diffusion induced by dopants [START_REF] Hao | Effects of Phosphorus Doping on Structural and Optical Properties of Silicon Nanocrystals in a SiO2 Matrix[END_REF] . The Si-NCs size remains however limited and almost controlled by the presence of the SiO2 barrier layers, as shown by the histograms of the size distribution (Fig. 2(f)). For the undoped layer, the corresponding histogram shows that the average Si-nanocrystal size is close to 7 nm, which is lower than the nominal thickness of the SiO layers (i.e. 10 nm). This is due to the weak density of the SiO layer that is obtained by evaporation in ultrahigh vacuum with a substrate maintained at room temperature. The densification of the SiO layer induced by the annealing treatment leads to the formation of Si-NCs having diameters slightly lower than the nominal SiO layer thickness. For the doped multilayers, an increase of the Si-NCs size is observed up to a mean size close to 15 nm. To get further information on the local chemistry, we used STEM-EELS to map the chemical distribution of the elements at the nanoscale. Figures 2(c)-(e) show the chemical maps relative to Si 0 , P and Si 4+ , respectively. Si 0 atoms are localized in the substrate and in the nanocrystals formed in the sample. Si 4+ atoms (corresponding to Si atoms present in SiO2) are located in the matrix surrounding the nanocrystals. Interestingly, P atoms appear localized mainly in the nanocrystals. However, since EELS signals originate from the volume irradiated by the electron probe, P atoms cannot be localized in the thickness of the TEM lamella and, one cannot exclude the presence of P atoms at the interface between Si-NCs and the surrounding silica matrix. To further clarify the location of P atoms, APT measurements have been performed on P-doped SiO/SiO2 multilayers containing 0.7 and 1.1 at% of P. Figure 3 (a)-(b) shows typical slice views (46 × 46 × 10 nm 3 ) drawn from the 3D reconstructed volume and representing both Si atoms (in red) and P atoms (in black) for the multilayer containing 1.1 at% of P. One can clearly identify Si-NCs represented by a concentration of red spots (Fig. 3(a)). Their average size is about (6.5 ± 0.1) nm, in good agreement with the results obtained by STEM. Fig. 3(b) shows that P atoms can be found mainly inside Si-NCs to an amount of 60 % and in the surrounding silica matrix to an amount of 40 % thus confirming the preferential localization of P in Si-NCs evidenced by STEM-EELS measurements. A careful analysis of the P content inside Si-NCs indicates that this content is almost independent on the Si-NC diameter (Fig. 3(c)). For the multilayer containing 1.1 at% of P, the mean P content reaches a value of (10.3 ± 1.1) at%, i. e. well beyond the solid solubility limit of P in bulk Si (i.e. 6 x 10 20 cm -3 at 1100°C). As the P content in the multilayers increases from 0.7 at% to 1.1 at%, the mean Si-NC size increases from about (5.4 ± 0.1) nm to (6.5 ± 0.1) nm, which is well explained by the enhanced Si diffusion in presence of P [START_REF] Trad | Influence of phosphorus on the growth and the photoluminescence properties of Si-NCs formed in P-doped SiO/SiO2 multilayers[END_REF][START_REF] Demoulin | Atomic scale characterization of N-doped Si nanocrystals embedded in SiO2 by atom probe tomography[END_REF] . Meanwhile the P content in the Si-NCs also increases from about (5.6 ± 1.4) at% to (10.3 ± 1.1) at% at constant Si-NC size. However, despite the large P contents in the Si-NCs, alloying and formation of SiP nanoparticles is not observed.

Figure 4(a) shows a STEM-HAADF image of the P-doped SiO/SiO2 multilayer with the largest P content considered in this study, i.e. 4.1 at%. Spherical shaped nanoparticles with typical sizes ranging between 10 nm and 50 nm can be identified. Typical EDS spectra measured either in a nanoparticle or in the surrounding matrix are shown in Fig. 4(b). In the matrix, both Si and O are detected but no phosphorus. In the particles, one can identify a significant PK peak originating from phosphorus atoms. A quantitative analysis yields a P content close to 50 %, thus suggesting that the large sized nanoparticles consist of SiP. For such high P content in the multilayer, Si-P alloying was already observed in annealed Si:P thin films [START_REF] Geiskopf | Thin films of SiP lamellar alloys: a first step toward 2D SiP[END_REF] .

Figure 5(a) shows the evolution of the PL spectra measured at room temperature on an undoped SiO/SiO2 multilayer prior to and after annealing for 5 minutes at various temperatures between 950°C and 1100°C. For the as-deposited multilayer, no PL can be observed. This can be easily understood since phase separation and nanocrystal formation has not yet started. For annealing temperatures larger than 950°C, a rather broad PL band is observed centered at about 880 nm. This band, which is due to radiative recombinations between electrons and holes confined in Si-NCs, is found to increase and to redshift as the annealing temperature further increases up to 1100°C. The observed redshift can be well explained by the increase of the mean Si-NCs size and the induced decrease of the bandgap. When the annealing time increases from 5 to 15 minutes, the intensity of the PL band further increases and the position of the peak is slightly redshifted. The PL intensity increase can be related to an improved crystalline quality of the Si-NCs. In the following, we will investigate the evolution of the PL from P-doped SiO/SiO2 multilayers. Figure 5(b) shows the PL spectra measured on both undoped and P-doped SiO/SiO2 multilayers with the P content varying between 0.3 at% to 4.1 at% after annealing at 1100°C for 15 minutes. For the undoped multilayer, a strong PL band is observed centered at about 920 nm. For increasing P contents, the PL band first redshifts and its intensity strongly decreases before quenching occurs for P concentrations larger than 0.7 at%. The observed redshift is well explained by the increase of the mean Si-NC size as shown by APT and STEM measurements.

As the P content in the multilayer increases, the fraction of P atoms located in the Si-NCs also increases as evidenced by APT measurements (see fig. 3(c)). We can thus tentatively explain the observed PL quenching by the Auger effect and by the increased mean size of the Si-NCs in which the quantum confinement effect does not occur. In addition, when the solubility limit is exceeded, P atoms may further occupy interstitial sites thus generating defects that may further contribute to quench the observed PL. Finally, for the largest P content investigated in this study (i.e. 4.1 at%), Si-NCs are no longer observed thus explaining the disappearance of the associated PL band.

Figure 6 (a) shows the Fourier transform infrared absorption spectrum measured on a multilayer containing 1.1 at% of P and annealed at 1100 °C during 15 min. We can identify several sharp features at 461 cm -1 , 816 cm -1 and 1085 cm -1 which are due to the rocking, bending and stretching mode of Si-O-Si bonds in silica, respectively [START_REF] Pai | Infrared spectroscopy study of SiOx thin films produced by plasma enhanced chemical vapor deposition[END_REF] . In addition, we observe the presence of oscillations on the baseline which originate from interferences in the multilayers.

Interferences are due to multiple reflections occurring at the multilayer/substrate and air/multilayer interfaces. Finally, a rather broad band is observed centered at about 2000 cm -1 which cannot be attributed to a vibrational contribution. According to its position, this band could possibly be the signature of a LSPR in P-doped Si-NCs. To clarify the origin of this band and to eliminate the contribution arising from interferences in the multilayers, we have firstly performed a simulation of the interferences based on the model developed by Swanepoel [START_REF] Swanepoel | Determination of the thickness and optical constants of amorphous silicon[END_REF] . The result of the simulation is shown in Fig. 6(a) as a red solid line. To carry out the simulation of the interferences, the thickness of the layer was taken from the TEM analysis and the refractive index was modeled by the Cauchy relation. Figure 6(b) displays the Fourier-transform infrared absorption spectra for different P-doped multilayers with P contents varying between 0.3 at% and 4.1 at% after subtracting the contribution related to interferences. For the multilayer with the lowest P content (i.e. 0.3 at%), we can only identify the vibrational features related to the Si-O-Si bonds. No signal can be detected in the spectral range 1500-2000 cm -1 . However, for multilayers containing 0.7 at% and 1.1 at% of P, a rather broad band extending from 1500 cm - 1 till 3700 cm -1 can be identified. Interestingly, this band shifts to larger wavenumbers as the P content increases from 0.7 at% up to 1.1 at%. Considering the position of this broad band, which is in good agreement with previous results [START_REF] Rowe | Phosphorus doped Si nanocrystals exhibiting mid-infrared localized surface plasmon resonance[END_REF][START_REF] Poumirol | Hyper-doped silicon nanoantennas and metasurfaces for tunable infrared plasmonics[END_REF][START_REF] Zhang | Infrared nanoplasmonic properties of hyperdoped embedded Si nanocrystals in the few electrons regime[END_REF] , and the observed blueshift induced by the increased P concentration, we can attribute this band to a LSPR in P-doped Si-NCs. Moreover, as expected in presence of free carriers in the nanocrystals, the emergence of a LSPR related band occurs concomitantly with the PL quenching. This occurs for a minimal P content in the Si-NCs slightly larger than 5 at%. When the P content further increases, the fraction of P atoms located in the Si-NCs increases thus leading to an increased free charge carrier density. As a result, the LSPR mode blueshifts as expected from the Drude's model. However, for a P content of about 1.9 at%, the observed trend does not continue but instead a redshift is observed. For such large P concentrations, an increased amount of P atoms either occupy interstitial sites or form small clusters, which does not lead to additional free charge carriers. Finally, for the highest P content considered in this study (i.e. 4.1 at%), we can only observe vibrational features. In this case, the STEM-EDS measurements clearly showed the formation of SiP nanoparticles (see fig. 4), which then explain the absence of any Si-NC related LSPR band.

To determine the free charge carrier density as well as the mobility of the free charge carriers, simulations of the infrared absorption spectra were carried out assuming a core-shell model (Fig. 7(a)). Spherical shaped Si-NCs at the core being surrounded by a SiO2 shell to simulate the effect of the silica matrix are embedded in an effective medium characterized by a refractive index neff that is allowed to vary between the refractive indices of SiO2 (1.42) and Si (3.47). The absorption cross section was calculated assuming the theory of Mie for a spherical particle within the quasi-static approximation: with the dielectric function of heavily doped Si that follows the Drude model :

= - + = 8π ω - ω + 2 !
where " is the bulk plasma frequency, is the free carrier damping constant (the inverse of the charge carrier scattering time) and is the high frequency dielectric constant of Si-NCs.

" is proportional to the square root of the density of free charge carriers while the free carrier damping constant is inversely proportional to the charge carrier mobility. Figure 7(b) shows the experimental infrared absorption spectra (red lines) superimposed with the simulated spectra (blue dotted lines calculated using PyMieScatt package) for P-doped multilayers with P concentrations ranging between 0.7 and 1.9 at%. For all the samples, the diameter of the core is maintained equal to 10 nm. In the diameter range explored in this work, i.e. between 5 and 12 nm, no significant modification is observed in the simulation. However in our model, which is able to simulate both vibrational and electronic contributions, the diameter of the shell dshell is a parameter that needs to be adjusted to obtain the correct intensity ratio between both vibrational and electronic features. Moreover, the refractive index of the effective medium neff is varied between the refractive indices of SiO2 and Si to simulate correctly the position of the vibrational peaks. For all the samples, a good agreement is obtained for neff =3, which suggests that the core-shell nanoparticles feel a surrounding medium that is strongly influenced by the presence of other Si-NCs. The simulated spectra are found to correctly match the experimental spectra for the parameters listed in Table I. The free charge carrier density is found to increase from 2×10 20 cm -3 to 2.3×10 20 cm -3 when the P content increases from 0.7 at% to 1.1 at%, i.e. a P concentration in the Si-NCs equal to 5.6 at% and 10.3 at%, respectively. Assuming an atomic concentration of Si of 5×10 22 cm -3 , and based on the P atomic concentration given in Fig. 3(c), this corresponds to a P activation rate of 8 and 4.6 %, respectively. In bulk silicon, this value is highly dependent on the doping level and on the annealing treatments carried out to activate dopants. The values obtained in this study are lower than those reported for highly P-doped Si nanodisks, which can reach values as high as 70 % [START_REF] Poumirol | Hyper-doped silicon nanoantennas and metasurfaces for tunable infrared plasmonics[END_REF] . In the latter, dopants were introduced by low energy implantation of P into bulk Si and the dopant activation was promoted by laser thermal annealing allowing strong out of equilibrium conditions. The activation efficiency reported in this work is, however, significantly larger than the P activation efficiencies reported for free-standing P-doped Si-NCs which is lower than 1 % 24 or for P-doped Si-NCs obtained by diffusion doping of freestanding Si/SiOx core/shell Si-NCs via high temperature annealing in phosphorus vapor for which the activation efficiency is in the range of 0.8-7 % [START_REF] Bubenov | Diffusion doping route to plasmonic Si/SiOx nanoparticles[END_REF] . In the latter work, this value is underestimated by the presence of inactive P atoms in the oxide shell which give rise to activation efficiency of 34 % after chemical etching of the oxide shell. The value reported here is further in good agreement with activation efficiency values between 3 and 5 % determined for heavily P-doped Si-NCs embedded in SiO2 [START_REF] Zhang | Infrared nanoplasmonic properties of hyperdoped embedded Si nanocrystals in the few electrons regime[END_REF] . Moreover, the spectral width of the LSPR peak is rather narrow as compared to the values reported in the literature [START_REF] Kramer | Plasmonic properties of Si nanocrystals doped with boron and phosphorus[END_REF][START_REF] Bubenov | Diffusion doping route to plasmonic Si/SiOx nanoparticles[END_REF] . One reason is probably the narrow P concentration distribution of doped Si-NCs, as indicated in Fig.

3(c)

. Even though we see a heterogeneous broadening of the LSPR peak due to different free carrier concentrations within Si-NCs from one Si-NC to another one, the width of the peak is mainly controlled by the scattering time of free carriers. Based on the simulation, it is then meaningful to extract the charge carrier mobility which is estimated to be around 25 cm 2 V -1 s -1 .

This value is three times lower than values reported for P-doped bulk Si (i.e. ~ 60 cm 2 V -1 s -1 ) [START_REF] Masetti | Modeling of carrier mobility against carrier concentration in Arsenic, Phosphorus and Boron doped silicon[END_REF] .

The decreased charge carrier mobility could have different origins. According to the recent work of Zhang et al. [START_REF] Zhang | Infrared nanoplasmonic properties of hyperdoped embedded Si nanocrystals in the few electrons regime[END_REF] , the size-induced effect on the scattering time may occur for Si-NCs having diameters lower than 5 nm, which is however not the case in our study. It is then rather expected that the weak free charge carrier mobility is induced by the presence of non-active dopants which may generate defects or to small clusters that may act as scattering centers.

To summarize, using P-doped (SiO/SiO2) multilayers prepared by evaporation, we have shown that P atoms can be inserted in the core of Si-NCs embedded in a silica matrix with concentrations well beyond the solid solubility limit of P in bulk Si. Structural and chemical analysis allowed us to demonstrate that most of P dopants are localized in Si-NCs and the P content in Si-NCs can reach values as high as 10 at%. Infrared absorption measurements gave evidence of a tunable LSPR band induced by the activation of dopants that give rise to free electrons. After subtraction of optical interferences in the thin layer and thanks to a core-shell modeling of Si-NCs embedded in a silica matrix, simulations based on the theory of Mie and the Drude model allowed us to reproduce correctly both the vibrational and the electronic contributions to the absorption spectra. Both the charge carrier mobility and their density were determined. A free carrier concentration as high as 2×10 20 cm -3 corresponding to a doping efficiency of 8 % was obtained. The charge carrier mobility is estimated to 25 cm 2 V -1 s -1 , i.e.

between 2 and 3 times lower than the charge carrier mobility in bulk Si. This low carrier mobility can be explained by P atoms that are not electronically active. Our work is the first report on a tunable LSPR in highly doped Si-NCs obtained using the multilayer architecture. It is an interesting and promising method to implement such plasmonic structures in future nanophotonic devices operating in the infrared range. Si IV present in SiO2 (e), (f) Histograms of the Si-NC size distributions for both undoped and doped SiO/SiO2 multilayers with the P content varying between 0.7 at% and 1.9 at% . Table I

Figure 1 :

 1 Figure 1: Chemical composition of the P-doped SiO/SiO2 multilayers as a function of the

Figure 2 :

 2 Figure 2: Transmission electron microscopy analysis of a P-doped SiO/SiO2 multilayer with a

Figure 3 :

 3 Figure 3: (a)-(b) 2D slice views (46 x 46 x 10 nm 3 ) taken from the 3D reconstructed volume

Figure 4 :

 4 Figure 4: (a) STEM-HAADF image of a P-doped SiO/SiO2 multilayer containing 4,1 at% of P.

Figure 5 :

 5 Figure 5: (a) Evolution of the PL spectra measured at RT on an undoped SiO/SiO2 multilayer

Figure 6 :

 6 Figure 6: (a) Infrared absorption spectrum measured for the P-doped SiO/SiO2 multilayer

Figure 7 :

 7 Figure 7: (a) Scheme of the core-shell model used for the simulations, (b) Experimental infrared

  Figure 1

  Figure 3