Jad Oueis

Razvan Stanica

Fabrice Valois

Energy Harvesting Wireless Sensor Networks: From Characterization to Duty Cycle Dimensioning

Energy harvesting capabilities are challenging our understanding of wireless sensor networks by adding recharging capacity to sensor nodes. This has a significant impact on the communication paradigm, as networking mechanisms can benefit from these potentially infinite renewable energy sources. In this work, we study the consequences of implementing photovoltaic energy harvesting on the duty cycle of a wireless sensor node, in both outdoor and indoor scenarios. We show that for the static duty cycle approach in outdoor scenarios, very high duty cycles, in the order of tens of percents, are achieved. This further eliminates the need for additional energy conservation schemes. In the indoor case, our analysis shows that the dynamic duty cycle approach based solely on the battery residual energy does not necessarily achieve better results than the static approach. We identify the main reasons behind this behavior, and test new design considerations by adding information on the battery level variation to the duty cycle computation. We demonstrate that this approach always outperforms static solutions when perfect knowledge of the harvestable energy is assumed, as well as in realistic deployments, where this information is not available.

I. INTRODUCTION

During the last decade, the energy-efficiency of wireless sensor networks (WSN) received a lot of attention from academic and industrial communities, with the extension of network lifetime under limited battery constraints as the main goal [START_REF] Anastasi | Energy Conservation in Wireless Sensor Networks: A Survey[END_REF]. Recently, energy harvesting has emerged as an alternative paradigm for powering a node. Energy harvesting is the act of collecting energy from the environment (e.g. solar) or other ambient energy sources (e.g. vibrations), and converting it into electrical energy [START_REF] Sudevalayam | Energy Harvesting Sensor Nodes: Survey and Implications[END_REF]. By implementing energy harvesting, potentially infinite energy becomes available to the sensor node. Energy management is still essential in this context. However, the focus shifts from the recurrent performance-lifetime trade-off towards exploiting harvesting opportunities in order to improve the WSN performance, e.g. by extending the nodes duty cycles [START_REF] Vigorito | Adaptive Control of Duty Cycling in Energy-Harvesting Wireless Sensor Networks[END_REF]. To adapt the wireless node behavior to the harvested energy, a multitude of harvesting-aware mechanisms have been proposed, mainly focused on dynamic duty cycle scheduling. Such mechanisms aim at achieving longer node activity time, which further leads to lower communication delays and reduced network latency.

In this paper, we address the photovoltaic (PV) energy harvesting paradigm in WSNs, by studying its impact on the node duty cycle, in both outdoor and indoor scenarios. This allows us to outline three important contributions: i) Operation with relatively high duty cycles (more than 30%) is feasible in outdoor scenarios, even when a static duty cycle approach is used. This questions the need for sophisticated energy management techniques in the case of outdoor WSN.

ii) In indoor applications, where PV energy is scarce, a dynamic approach based solely on the battery residual energy is not always better than a static duty cycle solution. This is unexpected considering that adaptive solutions are commonly considered to be more efficient.

iii) By adding supplementary information regarding the battery level variation to the duty cycle computation, the dynamic duty cycle approach becomes more performant than the static one. This is verified when assuming a perfect knowledge of the environment, with optimized parameters of the dynamic mechanism, as well as in practical cases, in which parameters are fixed to default values.

To study the relation between energy harvesting and duty cycle dimensioning, we propose a detailed energy harvesting analytical model for a node, taking into account energy collection, energy management and energy consumption. We provide as input to this model real-world energy datasets representative of several indoor and outdoor environments, and consider different power storage technologies.

The remainder of the paper is organized as follows. Sec. II presents related work on dynamic duty cycle scheduling algorithms for non-harvesting and harvesting-based WSNs. We introduce our complete energy harvesting model in Sec. III, detailing each of its components and the relation between energy production and consumption with duty cycle dimensioning. In Sec. IV, we define our use cases and the input parameters of the model. In Sec. V, we consider a static duty cycle approach and quantify the gain achieved by implementing energy harvesting. We address dynamic duty cycling and underline the limitations of battery residual energy-based solutions in Sec. VI, while the addition of information on the battery level variation is evaluated in Sec. VII. Finally, we conclude the paper and discuss future work in Sec. VIII.

II. WSN DUTY CYCLE BACKGROUND

In WSNs, dynamic duty cycles are usually applied as energy conservation techniques [START_REF] Anastasi | Energy Conservation in Wireless Sensor Networks: A Survey[END_REF]. Several approaches can be used to implement duty cycle mechanisms [START_REF] Maheswar | A Survey on Duty Cycling Schemes for Wireless Sensor Networks[END_REF]: i) in on-demand techniques [START_REF] Gu | Radio-Triggered Wake-Up for Wireless Sensor Networks[END_REF], nodes wake-up only when they are involved in communication. The wake-up procedure requires a dedicated hardware, with a separate, low-energy radio interface; ii) scheduled rendezvous solutions [START_REF] Ye | Medium Access Control with Coordinated Adaptive Sleeping for Wireless Sensor Networks[END_REF] determine the duty cycle schedule of all the nodes, building on a tight synchronization within the network, which leads to additional communication overhead; iii) in asynchronous protocols [START_REF] Buettner | X-MAC: A Short Preamble MAC Protocol for duty cycled Wireless Sensor Networks[END_REF], each node wakes up independently of the others, which requires more frequent wake-up periods in order to reduce latency.

In WSNs, where energy is the main challenge, a higher duty cycle increases the energy consumption. Hence, duty cycle dimensioning is generally a compromise between the desired network lifetime and performance. However, in energy harvesting-based networks, taking profit from the harvested energy allows improving the network performance by increasing the node duty cycle, while extending the network lifetime by recharging depleted batteries. As a matter of fact, experimental studies [START_REF] Raghunathan | Design Considerations for Solar Energy Harvesting Wireless Embedded Systems[END_REF] showed that a harvesting-capable node, placed outdoors for a week, can operate continuously with a 20% static duty cycle, which is much higher than the 1%-5% envisioned for most WSN applications. In this paper, we further investigate the theoretical limits of static duty cycles that guarantee continuous operation for longer deployment periods, and under several indoor and outdoor conditions.

As the amount of available energy is highly variable depending on the harvesting opportunities, adaptive duty cycle solutions appear to be a viable solution. The adaptive approaches proposed in the literature are generally based on the harvested energy rate (using historic [START_REF] Romaniello | Sustainable Traffic Aware Duty-Cycle Adaptation in Harvested Multi-Hop Wireless Sensor Networks[END_REF] or forecast [START_REF] Kansal | Power Management in Energy Harvesting Sensor Networks[END_REF] models), the battery level [START_REF] Vigorito | Adaptive Control of Duty Cycling in Energy-Harvesting Wireless Sensor Networks[END_REF], [START_REF] Yoo | Dynamic Duty Cycle Scheduling Schemes for Energy-Harvesting Wireless Sensor Networks[END_REF], [START_REF] Le | Power Manager with PID Controller in Energy Harvesting Wireless Sensor Networks[END_REF], or a combination of both [START_REF] Yoo | Dynamic Duty Cycle Scheduling Schemes for Energy-Harvesting Wireless Sensor Networks[END_REF]. Using predictions of energy availability can be an efficient solution in outdoor homogeneous environments, where the solar energy follows daily and yearly periodicities [START_REF] Kansal | Power Management in Energy Harvesting Sensor Networks[END_REF]. However, as shown by Gorlatova et al. [START_REF] Gorlatova | Networking Low-Power Energy Harvesting Devices: Measurements and Algorithms[END_REF], the availability of indoor PV energy is not necessarily periodic, and hardly predictable due to the heterogeneity of most indoor locations.

On the other hand, prediction-free algorithms were proposed by Vigorito et al. [START_REF] Vigorito | Adaptive Control of Duty Cycling in Energy-Harvesting Wireless Sensor Networks[END_REF], and Yoo et al. [START_REF] Yoo | Dynamic Duty Cycle Scheduling Schemes for Energy-Harvesting Wireless Sensor Networks[END_REF]. Nevertheless, these algorithms are based solely on the battery residual energy. As detailed in Sec. VI-A, we show that an approach using only information on battery residual energy suffers from several insufficiencies, such as aggressive initial phase battery depletion and the lack of short-term duty cycle variation. Starting from these observations, we evaluate a prediction-free dynamic duty cycle approach, which adds information about the node battery level variation to the duty cycle dimensioning process. To the best of our knowledge, this is the first study evaluating the integration of battery level variation in the duty cycle computation, showing that this approach solves the main problems discussed above.

III. ENERGY HARVESTING MODEL: FROM THE LIGHT TO

THE SENSOR Focused on a single wireless sensor node, our PV energy harvesting model links energy collection, energy consumption, and energy management. In a previous work [START_REF] Oueis | Linking the Environment, the Battery, and the Application in Energy Harvesting Wireless Sensor Networks[END_REF], we used this model to determine the operational limits of its three main components: the environment, the battery and the application.

In the following, we consider that time is divided into timeslots of duration T seconds each, and denote by T the set of all timeslots of the observation period. We suppose the duty cycle DC(t) to be constant for the duration of a timeslot t. We denote by E h (t) the energy harvested by the node during a timeslot t, and by E r (t) the battery residual energy at the beginning of the timeslot. We define E h = {E h (t); ∀t ∈ T } as a vector representing the energy harvested by the node over the observation period. Similarly, DC = {DC(t); ∀t ∈ T } gives the node duty cycle values and E r = {E r (t); ∀t ∈ T } the energy stored by the node. The proposed model takes E h and DC as input, and returns E r as output.

A. Photovoltaic energy collection

PV energy derives from solar or artificial light, and is harnessed by a solar panel formed of PV cells. The energy harvested by the node during a timeslot t, E h (t), depends on several parameters. First of all, the global horizontal irradiance I gh (t), which is the total amount of shortwave radiation received by the PV panel, during timeslot t, measured in W/m 2 . The solar panel itself represents the second major actor in the energy harvesting process, with two important parameters: S P V , the solar panel illuminated area, in m 2 , and η P V , the PV cell efficiency, defined as the ratio of the power actually delivered by the solar cell to the incident power it received. The PV cell efficiency may present some variations over time, therefore we consider the mean cell efficiency over the deployment period. Thus, considering timeslots of constant duration T , the energy harvested during timeslot t is computed as:

E h (t) = I gh (t) • S P V • E[η P V ] • T (1)

B. Power storage

Rechargeable batteries are generally used as a power storage component in harvesting-based sensor nodes [START_REF] Sudevalayam | Energy Harvesting Sensor Nodes: Survey and Implications[END_REF]. Practically, nodes are equipped with rechargeable batteries, providing an initial reserve of energy and capable of storing the harvested energy delivered by the solar panel. The most common rechargeable technologies are Nickel Metal Hydride (NiMH) and Lithium Ion (Li).

The batteries are mainly characterized by their capacity E max , and their charge/discharge efficiency η bat . E max represents the maximum amount of energy a battery is capable of storing. η bat is a ratio strictly lower than one, meaning that there is always energy loss when charging or discharging the battery [START_REF] Basagni | Wireless Sensor Networks with Energy Harvesting[END_REF]. Moreover, batteries suffer from self-discharge, which is the loss of energy even when the battery is not being used. The battery lifetime can also be limited, even when recharge is possible, since its capacity is reduced at each deep recharge cycle. It should be noted that under different operating temperatures, changes in battery characteristics are expected [START_REF] Buchmann | Batteries in a Portable World: A Handbook on Rechargeable Batteries for Non-engineers[END_REF]. However, without loss of generality, we consider in this study a constant battery-specific discharge efficiency representing the average of the charge losses a battery may undergo under varying operating conditions.

C. Energy consumption

The main power consumption sources in a wireless sensor node are the sensor, the microcontroller, and the radio transceiver [START_REF] Sudevalayam | Energy Harvesting Sensor Nodes: Survey and Implications[END_REF]. The CPU power consumption is several orders of magnitude lower than the others, and it can be seen as negligible [START_REF] Prayati | A Modeling Approach on the TelosB WSN Platform Power Consumption[END_REF]. The power consumption induced by the sensors is application-dependent. For example, a GPS-equipped sensor will result in a significant energy consumption. However, sensors used in general (e.g. temperature, humidity) have low energy requirements. Therefore, in the following, the energy consumption model only considers the radio transceiver activity, but including sensing or CPU information in the model is straightforward.

The activity of the node radio transceiver is dictated by the duty cycle: during given timeslots, the transceiver will be ON for a time DC(t) • T , and OFF for the rest of the time. During the OFF period, several states are possible, depending on whether functions such as voltage regulator, crystal oscillator or RAM access are enabled [START_REF]Chipcon CC2420 Datasheet[END_REF]. The energy consumption in these states shows slight variations, but the intensity of the drawn current in these states remains in the order of µA. During the ON period, the transceiver can also be in different states: receive (Rx), transmit (Tx) or listen (Lx). When not transmitting, the radio module is always sampling the channel (Lx mode) in order to find the preamble of a frame, in which case it switches to Rx mode.

Both transceiver datasheets [START_REF]Chipcon CC2420 Datasheet[END_REF] and experiments [START_REF] Prayati | A Modeling Approach on the TelosB WSN Platform Power Consumption[END_REF] show that the current drawn in these three states is quite similar, in the order of mA. Consequently, the exact time spent by a node in each of the three states during an ON period has a very small impact on the energy consumption: for instance, during a given timeslot, if the transceiver is only transmitting, its energy consumption would be approximately the same had it been transmitting during half of that time, and listening during the other half. This allows us to simplify the consumption model by considering only the impact of the node duty cycle. We therefore denote by I avg the average of the current intensity required in Tx, Rx and Lx, and V op the operating voltage. Thus, during an active time period T , the transceiver consumes:

E 0 = I avg • V op • T (2)
However, when considering duty cycled activity, during a time period T , the node radio module switches between ON and OFF states. Thus, the energy consumed by the node during timeslot t, denoted by E c (t) is:

E c (t) = E 0 • DC(t) (3)

D. Power manager

The power manager is in charge of delivering power to the node and includes an output regulator that matches the operating ranges of the batteries and the node [START_REF] Taneja | Design, Modeling, and Capacity Planning for Micro-solar Power Sensor Networks[END_REF]. An output regulator has an output efficiency η out , with a value strictly less than one. This affects the supply/demand ratio: more energy than initially needed by the node must be provided in order to satisfy the node demand. Thus, the node energy requirement during a timeslot t is equal to Ec(t) ηout . An input regulator can be required in order to condition the output of the panel to meet the operational constraints of the battery [START_REF] Taneja | Design, Modeling, and Capacity Planning for Micro-solar Power Sensor Networks[END_REF]. Because this parameter can be easily integrated in η out , we do not consider it specifically in our model.

Regarding the power manager, we adopt the energy management model presented in [START_REF] Taneja | Design, Modeling, and Capacity Planning for Micro-solar Power Sensor Networks[END_REF]. We extend this model by adding the battery charge/discharge efficiency to the battery model, as discussed in Sec. III-B. We assume that the battery is initially charged, i.e. E r (0) = E max . For the following timeslots, the node consumes energy as described in Eq. 3, and this energy is provided by the power manager. Several possibilities can be distinguished, based on the difference between the energy harvested by the node E h (t) and the total consumed energy,

Ec(t)
ηout . We define E o (t) as the overall energy gain during timeslot t. More precisely:

E o (t) = E h (t) - E c (t) η out (4)
Based on the value of E o (t), two possibilities open up for the power manager:

• E o (t) ≥ 0: in this case, the node is powered solely by the solar panel. Moreover, if the solar panel output exceeds the node requirements, the harvested energy surplus is stored in the battery, i.e. the node is charging. The node battery level is therefore written in an iterative manner:

E r (t + 1) = min (E r (t) + η bat • E o (t), E max ) (5)
• E o (t) < 0: the solar panel output is not sufficient to meet the node energy requirements, and the power manager needs to use power from the battery as well. In this case the battery is discharging, possibly until complete depletion. The battery level evolution is written as:

E r (t + 1) = max E r (t) + E o (t) η bat , 0 (6)
Using Eq. 5 and Eq. 6, we compute E r and use it in dimensioning the node duty cycle, both statically and dynamically, as discussed in Sec. V and Sec. VI, respectively. For this, we focus on the continuous operation of the node, i.e. the case where the node can sustain its activity without suffering from a battery depletion that forces it to stop and wait to recharge, thus becoming non-operational. Continuous operation (c.o.) is formally defined as:

c.o. ⇐⇒ E h (t) + η bat • E r (t) > E c (t) η out , ∀t ∈ T (7)
In our analysis we focus on the overall activity (OA) of the node under the constraint of continuous operation. We define the overall activity as the percentage of time during which the node is active out of the whole observation period | T |:

OA =    ∀t∈T DC(t)•T |T | • 100, if c.o. 0, otherwise (8)
This metric is based on the idea that a duty cycle mechanism rendering the node non-operational is not well dimensioned for applications requiring continuous operation. However, we note that this constraint can be easily removed if applications tolerating non-operational time are considered [START_REF] Oueis | Linking the Environment, the Battery, and the Application in Energy Harvesting Wireless Sensor Networks[END_REF]. In all the cases, a higher OA means the node can be kept ON for longer periods, which is directly linked to an improvement in network metrics, such as throughput or delay.

IV. USE CASE AND INPUT PARAMETERS

We use the model previously described to evaluate the overall node activity in several indoor and outdoor conditions. This section discusses the parameters of the model and the numerical values we use in our analysis.

A. PV energy

As input for the energy harvesting model, we determine the amount of harvested energy each timeslot E h (t), by using realworld datasets providing light availability for both outdoor and indoor locations. As we focus on the recharge of small wireless nodes, we consider a solar panel with an illuminated area S P V = 10 cm 2 , having an efficiency E[η P V ] = 0.25 [START_REF] Green | Solar Cell Efficiency Tables (v. 45)[END_REF]. We consider observation periods | T |≈ 1 year, divided in hourly timeslots, i.e. T = 3600 s.

1) Outdoor: To model light availability outdoor, we use datasets provided by the U.S. Department of Energy, presenting hourly measurements of I gh (t). For our study, we choose datasets collected from Los Angeles, during the whole year of 2014 [START_REF] Andreas | Solar Resource & Meteorological Assessment Project (SOLRMAP)[END_REF]. Using these datasets, we compute the amount of electrical energy delivered to the system by the solar panel E h (t). While consecutive days show very similar patterns, significant differences are mainly observed on a monthly level (Fig. 1a). 2) Indoor: For indoor scenarios, we use datasets provided by Gorlatova et al. [START_REF] Gorlatova | CRAWDAD Dataset Columbia/Enhants[END_REF], who conducted long term light energy measurements in office buildings in New York City [START_REF] Gorlatova | Networking Low-Power Energy Harvesting Devices: Measurements and Algorithms[END_REF]. We consider four particular locations indexed by A, B, C, and D. Tab. I shows the characteristics of each indoor location with its respective deployment period and average daily irradiance. Indoor data presented in Fig. 1b shows that, even in the same location, two consecutive days can have significantly different values of harvested energy during the same period of the day. This is due to several factors that make indoor energy availability harder to predict, such as the location of the sensor, the presence or not of people in the room, and their luminosity preferences.

B. Power storage and manager

Tab. II summarizes the main characteristics of the different rechargeable batteries we consider as power storage component. Regarding the output regulator, an efficiency η out = 0.8 is considered [START_REF] Taneja | Design, Modeling, and Capacity Planning for Micro-solar Power Sensor Networks[END_REF].

C. Energy consumption in the wireless sensor node

For the remainder of this work, we base our numerical application upon the Chipcon CC2420 RF transceiver as the node radio component [START_REF]Chipcon CC2420 Datasheet[END_REF]. Based on the datasheet, we consider the average current drawn by the sensor during active periods I avg = 16.7 mA, and V op = 3 V, resulting in an average hourly energy consumption E 0 = 180 J when no duty cycle is considered.

V. STATIC DUTY CYCLE ANALYSIS

In harvesting-based WSNs, the node performance is enhanced by adjusting parameters such as duty cycle, transmission power, and sensing reliability to name a few [START_REF] Sudevalayam | Energy Harvesting Sensor Nodes: Survey and Implications[END_REF]. In traditional networks this may lead to more energy consumption and, subsequently, network lifetime limitation. However, with energy harvesting, the recharge opportunities compensate for the energy consumption.

In this section, we highlight the gains achieved regarding the node duty cycle when implementing energy harvesting. Using the model proposed in Sec. III and the input data presented in Sec. IV, we compute the maximum achievable static duty cycle guaranteeing continuous operation for the node, for both outdoor and indoor locations, and for different rechargeable battery models. We first focus on a static duty cycle, meaning that DC(t) = DC static , ∀t ∈ T . In this case, we define the maximum achievable static duty cycle that is capable of guaranteeing continuous operation for the node:

DC max = max c.o. (DC static ) (9)
A. Outdoor location Fig. 2 shows DC max in the outdoor location, for a one year deployment, starting January. We compare results for a node with no harvesting capability and for a harvestingcapable node, for different battery models. We notice that, [START_REF] Sudevalayam | Energy Harvesting Sensor Nodes: Survey and Implications[END_REF].

0 20 40 60 80 100 A A -N i M H A A A -N i M H A A -L i U l t r a t h i n -2 0 0 U l t r a t h i n -1 0 0 U l t
without energy harvesting, the maximum achievable duty cycle is negligible in comparison with the harvesting-capable case. Introducing energy harvesting allows for much higher duty cycles, even for batteries with small capacities. For instance, the two smallest ultrathin batteries, with capacities of 43 mAh and 10 mAh, can achieve duty cycles of 18% and 4%, respectively, when the node has harvesting capabilities. However, for the same batteries, continuous operation is not even possible over a whole year without energy harvesting. While duty cycles of 1 to 5% are usually envisioned for sensor network applications [START_REF] Raghunathan | Design Considerations for Solar Energy Harvesting Wireless Embedded Systems[END_REF], implementing energy harvesting for outdoor applications allows duty cycles as high as 50% for an AAsized Li-Ion battery for example. This eliminates the urgent need for energy conservation schemes and energy consumption minimization protocols. It is noteworthy that these results are strongly dependent on the specific location we chose. In fact, harvested energy is affected by several conditions that vary from one location to another, such as weather conditions, and the presence of shading which limits the node exposure to direct light. While we consider a scenario particularly suitable for PV energy harvesting (i.e. sunny Los Angeles, no shadowing), the order of magnitude of outdoor achievable duty cycles remains generally the same in other cases as well.

B. Indoor locations

Fig. 3 shows DC max in the four indoor locations (Tab. I), for their corresponding deployment periods. Due to the significant difference in the amount of harvestable energy between outdoor and indoor locations, the maximum achievable duty cycles indoor are much lower than those achieved outdoor. The highest duty cycle of 2.9% is achieved at location C with AA-sized NiMH batteries as a power storage. For the ultrathin batteries of capacities 43 mAh and 10 mAh, continuous operation is not even guaranteed for locations A and B, where the amount of harvestable energy is too low.

Therefore, although implementing energy harvesting guar- antees continuous operation with higher values of the duty cycle in most of the cases, the low and random indoor harvesting opportunities limit the gain. We note that battery depletion is still faced for smaller batteries in particular shadowed locations.

0 0.5 1 1.5 2 2.5 3 3.5 4 A A -N i M H A A A -N i M H A A -L i U l t r a t h i n -2 0 0 U l t r a t h i n -1 0 0 U l t

VI. DYNAMIC DUTY CYCLE ANALYSIS

After studying the performance of a static duty cycle, we now focus on its dynamic adaptation. Considering the results in Sec. V, any comparison with a dynamic approach is superfluous in outdoor scenarios, as the performance of the static solution is largely sufficient for any realistic WSN application. We therefore limit our study in this section to the indoor scenario, considering the four different indoor locations described in Sec. IV. We conduct an analysis of two dynamic duty cycle approaches: the first based on Battery Residual Energy (BRE), and the second based on Battery Residual Energy and battery level Variation (BREV).

A. Insufficiencies in duty cycle scheduling based on battery residual energy (BRE)

As an example of a dynamic duty cycle approach based on the battery residual energy, we consider in the following a mechanism proposed by Yoo et al. [START_REF] Yoo | Dynamic Duty Cycle Scheduling Schemes for Energy-Harvesting Wireless Sensor Networks[END_REF]. This BRE mechanism uses two parameters, E th and DC min , determined by the application requirements: DC min is the minimum duty cycle acceptable by the application for correct functioning; E th is the minimum battery level required to guarantee the activity of the node with a fixed duty cycle DC min during a predefined period of time. Therefore, when the battery level is below the threshold E th , the duty cycle is set to this minimal value DC min . The functioning of BRE is based on the relative residual battery level, defined as:

r l (t) = E r (t) -E th E max -E th (10)
At each timeslot t, the node computes its duty cycle as follows:

DC(t) = max(DC min +(100-DC min )•r l (t), DC min ) (11)
Using the proposed model, we investigate the behavior of the BRE approach, detecting two major problems:

Aggressive battery depletion at the beginning of the deployment. When the battery is fully charged (E r (t) = E max ), the node functions with the maximum allowed value of the duty cycle DC(t) = 100%. At the beginning of the deployment, when the battery is full, the node sets high duty cycles accordingly, which leads to increased energy consumption. However, the node does not necessarily have significant recharge opportunities to compensate for this energy loss. This leads to a faster battery depletion at the beginning of the deployment, which can have a negative impact on the node overall activity, possibly rendering the node non-operational. We depict this problem in Fig. 4, in which we represent the battery level variation function of time at the beginning of the deployment, for a node implementing the BRE mechanism at indoor location C. Lack of short-term duty cycle adaptation. Since harvestable energy is usually low at indoor locations, the battery does not always have significant recharge opportunities, and is mostly discharging. With a mechanism based solely on the residual energy, the duty cycle variation will exactly follow the decreasing battery level variation. Higher duty cycles are only set at the beginning of the deployment when the battery is full, while lower duty cycles are set afterwards when the battery is close to depletion, as shown in Fig. 5. Furthermore, there is practically no short-term duty cycle variation, meaning no adaptation to the harvested energy pattern during the day. This is noticeable in Fig. 5, where we zoom in on the duty cycle variation for five consecutive days during the deployment, at indoor location C. We notice that there is only a negligible variation of the duty cycle, practically defeating the purpose of a dynamic scheduling algorithm, and suggesting that energy harvested during the day is not optimally used. These problems can heavily degrade the performance of BRE solutions. To verify this, we compare the performance of a static duty cycle approach with the BRE approach. For both cases, we compute the overall activity using Eq. 8, denoted by OA static , and OA BRE , respectively. In the case of the static approach, we consider DC(t) = DC max , ∀t ∈ T , where DC max is given by Eq. 9. For BRE, we iteratively compute DC(t) using Eq. 11.

Fig. 6 shows the ratio OA BRE /OA static for all indoor locations. Results show that, at indoor locations C and D, the ratio is lower than 1, meaning that the BRE approach has a lower overall activity than the static duty cycle. Our analysis therefore proves that inadequate design choices in the dimensioning of an adaptive duty cycle lead, in some indoor locations, to lower performance than a static approach.

B. Duty cycle scheduling based on battery residual energy and battery level variation (BREV)

In this section, we consider the use of supplementary information in the duty cycle dimensioning process, namely the battery level variation (i.e. the charge/discharge behavior of the battery). This allows us to propose a new class of dynamic duty cycle mechanisms, which we denote as BREV.

For comparison purposes, we focus in this section on a specific mechanism representative of the BREV approach. However, the good properties we observe, such as the better duty cycle adaptation to the available energy, are a result of the general design guidelines we propose for BREV, and not a by-product of the specific details of the selected mechanism.

Our proposal is based on the idea that an increase (resp. decrease) in the battery level indicates that the battery is charging (resp. discharging). This means that the node duty cycle is too small (resp. high) compared with the level of available harvestable energy, and an increase (resp. decrease) of the node duty cycle is due. In order to introduce this factor, we consider the battery level variation between two consecutive timeslots, denoted by ∆E r (t

) = E r (t) -E r (t -1). As discussed below, ∆E r (t) ∈ [-max d ; max c ],
where max d (Eq. 12) is the maximum possible discharge between two consecutive slots, and max c (Eq. 14) is the maximum charging value.

Similarly to the BRE mechanism, BREV uses a predefined battery threshold E th , and a minimum acceptable duty cycle DC min . Moreover, we consider that a maximum duty cycle value exists, DC M AV , above which the benefit of increasing the node activity is negligible. In real deployments, DC M AV depends on the application requirements and on the network topology. As shown below, our model can easily integrate different values for this parameter. Since we impose that DC(t) ≤ DC M AV , ∀t ∈ T , the maximum discharge value the battery can undergo between two consecutive timeslots happens when the duty cycle is the maximum allowed, i.e. DC M AV , and no energy is harvested:

max d = E 0 • DC M AV η out • η bat (12)
On the contrary, the maximum charging value appears when the energy consumption is minimum, and the harvested energy is maximum. The minimal energy consumption during a timeslot, E min , is written as:

E min = E 0 • DC min η out (13)
In this case, we have:

max c = η bat • (H max -E min ) (14)
where H max is the maximum value of the energy harvested by the node over the observation period. We also use a parameter r to combine information regarding the battery level variation and the battery residual energy, with different possible weights. More precisely, the node computes at each timeslot a duty cycle variation, ∆DC(t), as follows:

∆DC(t) = (DC M AV -DC min ) • r • ∆E r (t) + max d max c + max d + (1 -r) • E r (t) E max (15)
If the battery level reaches a threshold E r (t) ≤ E th , then the duty cycle is set to its minimum value. In other cases, the duty cycle is adapted based on ∆DC(t):

DC(t) = DC min if E r (t) ≤ E th DC min + ∆DC(t) if E r (t) > E th (16)
The parameter r in Eq. 15 serves us in order to investigate the contribution of the two factors: the battery level variation and the battery residual energy. The BRE solutions can be seen in this case as a particular case of BREV. More precisely, for r = 0, only residual energy is taken into consideration when computing the duty cycle, resulting in a BRE mechanism.

We notice that when information on the charge/discharge behavior of the battery is added to the node duty cycle computation, the problems affecting the BRE approach are alleviated. Fig. 4 shows that for BREV, contrarily to the BRE case, slower battery depletion happens at the beginning of the deployment, since the assigned values of the duty cycle are adapted to the charge/discharge behavior. As expected, the higher the values of r, the higher the impact on this initial phase. Moreover, BREV introduces a short-term duty cycle variation which follows the harvested energy pattern, as shown in Fig. 5 for the extreme case when r= 1.

VII. ANALYSIS OF DUTY CYCLE SCHEDULING BASED ON BATTERY RESIDUAL ENERGY AND BATTERY LEVEL VARIATION (BREV)

In this section, we propose an extensive evaluation study of the BREV approach. We show that, by considering battery level variation, a dynamic duty cycle strategy always performs better than static solutions, a condition not always verified by the BRE approach. We verify that this is true when the different parameters of the dynamic mechanism are optimized, as well as when these parameters are set to default values.

A. Oracle evaluation

Both the static duty cycle solution discussed in Sec. V and BREV use parameters whose optimal values can be computed only by an oracle with a priori knowledge of future harvestable energy availability. This is the case for DC static in the static approach, and for DC M AV and r in BREV. In the following, we make the assumption that this perfect information exists, in order to compare the two solutions under their optimal settings. Tab. III summarizes the numerical values we use in our analysis, with the optimal values for those parameters impacted by the location of the sensor.

The application-dependent DC min parameter is set to 0.01%, a value assumed to guarantee a minimum functionality for the node in critical low energy situations. The battery level threshold E th is set to 10 J, allowing the node to maintain a minimal activity for up to 23 days, even when no harvesting opportunities exist. In the case of the static duty cycle, we set DC static at each indoor location to DC max , the maximum achievable duty cycle capable of guaranteeing continuous operation, given in Eq. 9. The parameter DC M AV , which limits the possible values of the duty cycle, has a significant effect on the output of the BREV mechanism. In fact, with some values of DC M AV , we notice that continuous operation of the node cannot be guaranteed. This is because higher values of DC M AV lead to higher values of assigned duty cycles, which can cause early battery depletion in some cases. Thus, we determine, for each indoor location, the maximum value of DC M AV that is capable of guaranteeing continuous operation. We observe, in Fig. 7, that this value depends on the weight r associated to the battery level variation. We consider the best value of DC M AV , shown in Tab. III for each location, to be the one capable of guaranteeing continuous operation for all values of r. Even when using the best DC M AV value, the weight r has an impact on the performance of BREV. We compute the overall activity at each location for the different values of parameter r in order to investigate its effect.

As shown in Fig. 8, in locations A and B, the impact of the weight r is practically unnoticeable. This is due to the rare harvesting opportunities in these locations, which result in a steady decrease of the battery residual energy and very small battery level variation. Nevertheless, a small gain in the overall activity (in the order of 0.01%) exists for r ∈ [0.1, 0.4], in the case of location A, and for r ∈ [0.1, 0.2], in the case of location B. For locations C and D, where the available energy is more heterogeneously distributed, higher values of r achieve longer node activity time, which means that the battery state information has a significant impact. In fact, for both locations, the best overall activity is obtained for r = 1, when we completely disregard the residual energy and consider solely the battery level variation. Using these optimal values for the weight r, we compute, in Fig. 6, the ratio OA BREV /OA static . Unlike the BRE solutions, BREV always performs better than the static approach, showing the significant impact of the additional information we use in the dimensioning of the node duty cycle.

B. Practical case evaluation

Finding optimal values for parameters such as DC M AV , DC static and r is not trivial, as these values are location dependent and require knowledge of future events. In order to determine optimal values in a particular indoor location, the availability of harvestable energy for the entire deployment period would need to be accurately predicted.

Therefore, we consider in this section a more realistic approach, in which we fix the different parameters of the static and BREV mechanisms to default values in all indoor locations. We compute the corresponding overall activity, and prove that a dynamic approach that takes into consideration the battery state information always performs better than a fixed duty cycle.

Based on extensive results for the overall activity of the BREV approach with different parameter settings, we set, for all locations, DC M AV = 8% and r = 0.3. For the static duty cycle, we arbitrarily fix DC static = {0.7, 1, 1.5, 2}. We note that DC static = 0.7 is the maximum static duty cycle value that allows continuous operation in all the studied locations.

Fig. 9 shows a comparison between the obtainable overall activity when implementing BREV, and different values of the static duty cycle. We remind that the overall activity is set to zero when the static duty cycle is not capable of guaranteeing continuous operation (Eq. 8). We notice that, for all locations, the BREV dynamic duty cycle always achieves a better overall activity than the static duty cycle, when the latter is capable of guaranteeing continuous operation.

We argue that, unlike outdoor deployments, the heterogeneity of indoor environments is not suitable for static duty cycle solutions, as no single DC static value performs well in all locations. For example, a static duty cycle of 2% can be sustained in locations with high light availability (e.g. location C), but completely depletes the battery in more shadowed locations (e.g. location A). Smaller duty cycles enable continuous operation in more constrained environments, but underperform in better energy conditions. Well designed dynamic duty cycle solutions are therefore required, and the results in this section confirm that, unlike BRE, the BREV approach is a flexible solution, adapting to all environments, and increasing the node overall activity in both oracle and practical scenarios.

VIII. CONCLUSION

In this paper, we focus on duty cycle dimensioning in PV energy harvesting-based WSNs, for both outdoor and indoor scenarios. We propose a detailed energy harvesting model for a wireless sensor node that links energy collection, energy consumption, and energy management. By inputting realworld datasets to this model, we are able to dimension both static and dynamic duty cycle. Results show that for outdoor locations, solar energy is able to provide a sustained energy supply to enable WSNs continuously, with relatively high static duty cycles. For indoor scenarios, we show that state of the art mechanisms, basing dynamic duty cycle computation solely on battery residual energy, are not necessarily more performant than a static approach. Thus, we propose an approach that takes into account the battery level variation for dynamic duty cycle computation. Our analysis shows that, with this approach, inefficiencies of residual energy-based algorithms are resolved, and performance is better than the static approach.

For future work, it is of utmost importance to further evaluate and compare the performance of the adaptive duty cycle mechanisms in the context of the whole WSN, in order to evaluate and compare several performance criteria such as end-to-end delay and packet delivery ratio. Eventually, energy harvesting raises new questions since "infinite" lifetime does not necessarily mean "continuous". A node may cease its activity following a battery depletion, then return to the network once recharged by the harvested energy. Tolerance and feasibility of discontinuity in the network lifetime must be further investigated.

Fig. 1 .

 1 Fig. 1. Harvestable PV energy for a typical day of each month outdoor, and over 5 days in different indoor locations.

Fig. 2 .

 2 Fig. 2. Maximum achievable duty cycle for continuous operation -Outdoor.

Fig. 4 .

 4 Fig. 4. Battery discharge at the beginning of deployment -Indoor location C.

Fig. 5 .

 5 Fig. 5. Duty cycle variation during the complete deployment period -Indoor location C.

Fig. 6 .

 6 Fig.6. Ratio between overall activity with dynamic approaches and overall activity with a static duty cycle.

Fig. 7 .

 7 Fig. 7. Maximum value of DC M AV capable of guaranteeing continuous operation function of r at different indoor locations.

Fig. 8 .

 8 Fig. 8. Overall activity of BREV function of r at different indoor locations.

Fig. 9 .

 9 Fig.9. Overall activity comparison between the BREV dynamic approach with DC M AV = 8% and r = 0.3, and the static approach with fixed values of DC static , at indoor locations.

TABLE III PARAMETERS

 III FOR INDOOR LOCATIONS A, B, C AND D.


Stéphanie Buil

Julien Laverdant

Bruno Berini

Pierre Maso

Jean-Pierre Hermier

Xavier Quélin
email: quelin@physique.uvsq.fr

FDTD simulations of localization and enhancements on fractal plasmonics nanostructures

A parallelized 3D FDTD (Finite-Difference Time-Domain) solver has been used to study the near-field electromagnetic intensity upon plasmonics nanostructures. The studied structures are obtained from AFM (Atomic Force Microscopy) topography measured on real disordered gold layers deposited by thermal evaporation under ultra-high vacuum. The simulation results obtained with these 3D metallic nanostructures are in good agreement with previous experimental results: the localization of the electromagnetic intensity in subwavelength areas ("hot spots") is demonstrated; the spectral and polarization dependences of the position of these "hot spots" are also satisfactory; the enhancement factors obtained are realistic compared to the experimental ones. These results could be useful to further our understanding of the electromagnetic behavior of random metal layers.

The understanding of the interaction between light and disordered matter is of great interest for many research fields. For example, the scattering through random media have attracted considerable attraction in a double perspective: on account of fundamental issues in physics [START_REF] Sheng | Introduction to wave scattering, localization, and mesoscopic phenomena[END_REF] and for imaging problems in complex media [START_REF]Waves and imaging through complex media[END_REF]. In particular, due to their unique optical properties, the metallic disordered rough layers on dielectric surfaces have been widely studied both theoretically and experimentally over the two last decades [START_REF] Grésillon | Experimental observation of localized optical excitations in random metal-dielectric films[END_REF][START_REF] Sarychev | Electromagnetic field fluctuations and optical nonlinearities in metaldielectric composites[END_REF][START_REF] Ducourtieux | Near-field optical studies of semicontinuous metal films[END_REF][START_REF] Bozhevolnyi | Statistics of local field intensity enhancements at nanostructured surfaces investigated with a near-field optical microscope[END_REF][START_REF] Stockman | Localization versus delocalization of surface plasmons in nanosystems: can one state have both characteristics?[END_REF]. These random films are made of clusters of metal nanoparticles. Depending on the metallic filling fraction, three different regimes have been identified: (i) for low concentration, metal nano-islands are distributed on the dielectric substrate. (ii) For intermediate concentration at the vicinity of the percolation transition, growing clusters have coalesced and they exhibit self-similarity with fractal dimension. (iii) For high metal filling fraction, the deposited film is now an almost continuous metallic layer with dielectric voids. In the intermediate regime (ii) of interest, two main properties have been identified: a large absorption over the red and near-infrared spectrum [START_REF] Buil | Local field intensity enhancements on gold semicontinuous films investigated with an aperture near-field optical microscope in collection mode[END_REF] and strong spatial fluctuations of the near-field intensity [START_REF] Seal | Nearfield intensity correlations in semi-continuous metal-dielectric films[END_REF]. Theoretical studies give a quite simple description of these behaviors [START_REF] Shalaev | Optical properties of nanostructured random media[END_REF]. In this specific metallic filling fraction, the multi-scale clusters exhibit a self-similarity character, meaning that parts resemble the whole in a statistical sense. In the case of these disordered metal films, the selfsimilarity corresponds to a wide distribution of the sizes and shapes of the clusters. Due to surface plasmon resonances [START_REF] Raether | Surface plasmons on smooth and rough surfaces and on gratings[END_REF], an anomalous spectral absorption and a spectral broadening are observed. The random character of the self-organized metal clusters is also responsible of the multiple scattering, leading to light confinement and enhancement at a sub-wavelength scale. This spatial localization of the electromagnetic field leads to "hot spots" [START_REF] Grésillon | Experimental observation of localized optical excitations in random metal-dielectric films[END_REF] and huge fluctuations of the near-field intensity [START_REF] Buil | Local field intensity enhancements on gold semicontinuous films investigated with an aperture near-field optical microscope in collection mode[END_REF].

Fluctuations and enhancements of the near-field intensity above such random metal films have been reported by different groups using Near-Field Optical Microscopy (NSOM) [START_REF] Grésillon | Experimental observation of localized optical excitations in random metal-dielectric films[END_REF][START_REF] Ducourtieux | Near-field optical studies of semicontinuous metal films[END_REF][START_REF] Tsai | Photon scanning tunneling microscopy images of optical excitations of fractal metal colloid clusters[END_REF][START_REF] Bozhevolnyi | Direct observation of localized dipolar excitations on rough nanostructured surfaces[END_REF]. More recently, hot spots have been studied using the PEEM (Photo-Emission Electron Microscopy) technique [START_REF] Word | Photoelectron emission control with polarized light in plasmonics metal random structure[END_REF][START_REF] Awada | Experimental study of hot spots in gold/glass nanocomposite films by photoemission electron microscopy[END_REF]. The fluctuations of the near-field electromagnetic energy, evidenced by the hot spots, have also been related to the local density of states (LDOS) through the statistical distribution of fluorescent-molecule lifetimes [START_REF] Krachmalnicoff | Fluctuations of the local density of states probe localized surface plasmons on disordered metal films[END_REF].

Numerical simulations, based on effective-medium theory [START_REF] Sheng | Theory of the dielectric function of granular composite media[END_REF], quasi-static calculations [START_REF] Stockman | Localization versus delocalization of surface plasmons in nanosystems: can one state have both characteristics?[END_REF] or scaling theory [START_REF] Shalaev | Nonlinear optics of random metal-dielectric films[END_REF][START_REF] Shalaev | Nonlinear optics of random media[END_REF], have also predicted anomalous absorption and localized "hot spots". More recently, full wave numerical simulations (FDTD for Finite-Difference Time Domain) have been used to describe the near and the far-field optical responses of disordered metal-dielectric composites [START_REF] Chan | Electromagnetic scattering of waves by random rough surface: a finite-difference time-domain approach[END_REF][START_REF] Chettiar | FDTD modeling of realistic semicontinuous metal films[END_REF][START_REF] Thoreson | Fabrication and realistic modeling of three-dimensional metal-dielectric composites[END_REF]. But in all these numerical studies, the metal-dielectric structure was either: a sample computer-generated random film [START_REF] Markel | Small-particle composites: linear optical properties[END_REF]; or generated with a Kinetic Monte-Carlo algorithm [START_REF] Aubineau | Modélisation de couches minces métalliques fractales et calcul d'exaltations de champs électromagnétiques[END_REF]; or obtained from SEM (Scanning Electron Microscopy) or TEM (Transmission Electron Microscopy) images of a real structure [START_REF] Chettiar | FDTD modeling of realistic semicontinuous metal films[END_REF]. Despite the diversity of the techniques used, the structures under study were in all cases 2D films: the roughness of the metal layer is never taken into account. Some unresolved discrepancies between experimental and numerical results still remain: for example, the measured enhancement factors are much lower than the theoretical or numerical predicted ones; the effect of the incident polarization is not described by theoretical predictions.

In this article, we show that the combination between FDTD simulations and the roughness of random metal samples (real 3D nanostructures) gives results that are in very good agreement with experimental data. We study a disordered gold layer at the vicinity of the percolation threshold. We verify some previously confirmed experimental results: the spatial distribution of the near-field intensity strongly depends on the incident photon energy; the incident polarization state plays a major role on the local energy distributions; the nearfield intensity is dominated by the non-propragative contribution; the sizes of the hot spots is in good agreement with experimental measurements. We also show that the Probability Density Function (PDF) of the enhancement factors is well fitted by a log-normal function [START_REF] Awada | Experimental study of hot spots in gold/glass nanocomposite films by photoemission electron microscopy[END_REF][START_REF] Stockman | Giant fluctuations of local optical fields in fractal clusters[END_REF][START_REF] Sanchez-Gil | Calculations of the direct electromagnetic enhancement in surface enhanced Raman scattering on random self-affine fractal metal surfaces[END_REF].

The metallic structures studied are random gold films near the percolation threshold, just below the transition's concentration. Gold is evaporated on a glass cover-slip, under ultra-high vacuum conditions (10 -9 Torr). The mass thicknesses of the samples are between 6.1 nm and 6.4 nm and the typical rms roughness is about 4.3 nm. An AFM image [START_REF] Horcas | WSXM : a software for scanning probe microscopy and a tool for nanotechnology[END_REF] of a standard layer, 512 x 512 points, is shown in Fig. 1(a). The filling factor is about 64%. This image and the corresponding typical profile shown in Fig. 1(b) were obtained with an ultra-sharp probe (Veeco TESP-SS), in order to get the most accurate topography. At the vicinity of the percolation threshold, it is well-known that Euclidian and fractal clusters co-exist. The perimeter (P) versus the surface (S) of the clusters for a given sample is shown in Fig. 1(c). From the well-known relation P S D/2 , it is possible to determine the dimension D of the largest non-Euclidian clusters. We obtain D ~ 1.9 which is in very good agreement with the theoretical value 1.88.

The topographic data z(x,y) describing the surface of the layer are then transformed in an array of complex ( n n ik ) refractive index [START_REF] Johnson | Optical constants of the noble metals[END_REF], for the dielectric (air) and for the metal (gold) on each point. Then the (n, k) array is implemented in the parallelized 3D FDTD solver (Lumerical, FDTD Solutions) in order to define the 3D structure studied, supported by a glass substrate. The ensemble (substrate + gold layer) is then illuminated with various incident wavelengths under normal incidence through the substrate. The FDTD method or Yee algorithm [START_REF] Yee | Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media[END_REF][START_REF] Taflove | Computational Electrodynamics: The Finite-Difference Time-Domain Method[END_REF] [START_REF]Waves and imaging through complex media[END_REF] where media , are considered to be isotropic, linear and non-magnetic. It is possible to substitute the space and time derivatives by central differences using a staggered mesh where electric E and magnetic H field components are located at different points. For the time derivative, the electric and magnetic fields are disposed according to the leapfrog scheme. So the time and space evolutions of every field component can be written as, for example the x component of the electric field,

1/ 2 1/ 2 1/ 2 1/ 2 , , 1 , , , , ,, 1/ 2
nn xx ï j k ï j k i j k nn nn yy n zz i j k i j k i j k i j k , 1, , , , , 1 , ,,
x source i j k t EE HH HH J yz [START_REF] Grésillon | Experimental observation of localized optical excitations in random metal-dielectric films[END_REF] where i, j, k, n are integers that are running on the space and time mesh, the terms 1 / 2 are due to the staggered field components, y and z (and non appearing x) are the space mesh steps and t is the time mesh step.

Due to the limitation of the computer memory, the simulation cell has to be restricted to a finite region of space (500x500x2500 nm 3 ). To avoid non-physical reflections on the boundary of the simulation cell and due to the strong scattering behavior of the metallic structure which leads to a wide range of incidence of the scattered fields on the boundaries, we have used the Perfectly Matched Layers (PML) as boundary conditions with at least 216 layers. This ensures a maximum reflection of 10 -7 over all incidences in the half-space. Due to the cubic nature of the FDTD mesh, it is not possible to avoid stairs effects on the typical random metallic surfaces described in Fig. 1. To minimize this issue as much as possible, we have defined the value of the maximum mesh step equal to 0.8 nm on the volume of the layer that gives at least 11 x 10 9 simulation grid points to describe the metallic layer. The substrate which is a homogeneous glass coverslip, is also discretized on half of its thickness, just below the metallic nanostructure. Before performing the simulations whose results are presented in this article, we have checked that a double-sized meshing (1,6 nm) did not fundamentally modify the localization of the hot spots but, of course, with a poorest resolution, wider hot spots and weaker intensities. In Fig. 2, the results of our FDTD simulations are shown for a single geometry in order to compare the localization and the values of local intensities. All the calculated intensities above the gold layer are normalized with respect to the incident intensity: in other words, rather than intensity maps, these are enhancement factor's maps. For presenting the simulation results in the same manner as the experimental ones previously obtained with SNOM technique [START_REF] Bozhevolnyi | Statistics of local field intensity enhancements at nanostructured surfaces investigated with a near-field optical microscope[END_REF][START_REF] Buil | Local field intensity enhancements on gold semicontinuous films investigated with an aperture near-field optical microscope in collection mode[END_REF], these maps were calculated for each (x, y) point, at a distance 10 nm [START_REF] Grésillon | Experimental observation of localized optical excitations in random metal-dielectric films[END_REF] on top of the corresponding z(x, y) value. This leads to representations similar to the ones obtained in "constant distance mode" for SNOM technique. To improve clarity, in the plots of Fig. 2(b) and (c), the maximum value of the color bar was set at the mean value of the respective maximum. One has to underline that the polarization states of the normal incident plane waves are identical. From the enhancement factors ( ) maps, large local-field enhancements ( 14 at maximum for 650 nm and 27 for 1 µm) are revealed. And the positions of the highest enhancements are clearly wavelength-dependent.

Two points have to be noticed. First, the values are rather low (one or two orders of magnitude less) compared to the previously simulated and published ones [START_REF] Sarychev | Electromagnetic field fluctuations and optical nonlinearities in metaldielectric composites[END_REF][START_REF] Chettiar | FDTD modeling of realistic semicontinuous metal films[END_REF], but they are in a quite good agreement with the experimental results measured with aperture-SNOM technique [START_REF] Bozhevolnyi | Statistics of local field intensity enhancements at nanostructured surfaces investigated with a near-field optical microscope[END_REF][START_REF] Buil | Local field intensity enhancements on gold semicontinuous films investigated with an aperture near-field optical microscope in collection mode[END_REF] or with the PEEM technique [START_REF] Awada | Experimental study of hot spots in gold/glass nanocomposite films by photoemission electron microscopy[END_REF]. And the maximum values increase with the wavelength.

Secondly, we can analyze the Probability Density Functions (PDF) of the distribution. The PDFs, shown in Fig. 2(d) and (e), are well fitted by the log-normal distribution given by 2 ln ln exp 22

K P D F (4)
where K is a constant, is the standard deviation of the distribution ln and ln is the mean value of the logarithm of the enhancement factor . This behavior was previously identified in [START_REF] Genov | Plasmon localization and local-field distribution in metal-dielectric films[END_REF]. It is clearly seen that, for the same topography, larger values of are obtained for wavelengths in the near infrared. And, as theoretically predicted [START_REF] Shalaev | Optical properties of nanostructured random media[END_REF], the probability of occurrence for the largest values is very low which means that the enhancement factor distribution displays a long tail. Our simulation results also demonstrate the evidence of a distribution's broadening as it was recently observed in [START_REF] Awada | Experimental study of hot spots in gold/glass nanocomposite films by photoemission electron microscopy[END_REF].

Another result of our FDTD simulations is that the locations of the hot spots are depending on the incident polarization state. In Fig. 3, we present, for the same topography presented in Fig. 2(a), the spatial fluctuations of the enhancement factor for two incident wavelengths and for two linear orthogonal polarization states: = 1 µm and (a) and (b) for linear X and Y polarizations respectively; = 1.3 µm and (c) and (d) for linear X and Y polarizations respectively. This behavior, previously observed in experimental results [START_REF] Buil | Local field intensity enhancements on gold semicontinuous films investigated with an aperture near-field optical microscope in collection mode[END_REF][START_REF] Word | Photoelectron emission control with polarized light in plasmonics metal random structure[END_REF], clearly means that, in addition to a dependence of the incident photon's energy, the hot spots positions are also due to a complex interaction between the topographic configurations (nanogaps, tips,…) and the incident polarization of the electromagnetic field.

One can observe in Fig. 2 and 3 that the sizes and the shapes of the hot spots are very different, depending both on the incident wavelength and on the local structure. Considering the hot spots with intensities larger than 90% of the maximum intensity, we can deduce that the sizes (Half-Width at Half Maximum) of the hot spots are ranging between 7 nm to 32 nm which is in very good agreement with the values determined in [START_REF] Grésillon | Experimental observation of localized optical excitations in random metal-dielectric films[END_REF][START_REF] Word | Photoelectron emission control with polarized light in plasmonics metal random structure[END_REF][START_REF] Cang | Probing the electromagnetic field of a 15-nanometre hotspot by single molecule imaging[END_REF]. The shapes of the hot spots in the (x, y) plane vary from slightly to strongly stretched, depending on the gaps of coupled metal junction and/or the polarization of the incident electromagnetic field. The relationship between the polarization direction of the incident field and the orientation of the inter-particles is not as clear as previously demonstrated in [START_REF] Xu | Polarization-dependent surface-enhanced Raman spectroscopy of isolated silver nanoaggregates[END_REF][START_REF] Wei | Polarization dependence of surfaceenhanced Raman scattering in gold nanoparticles-nanowire systems[END_REF], but, in our case, the metallic structure studied is much more complex than 2 or 3 metal spheres in interaction.

Extracting the different components (x, y, z) of the calculated intensity is possible with the FDTD method. If we consider the near-field intensity distribution at constant height, for example z = 12 nm in Fig. 4(a), (b), (c), it is clear that the main contribution to the intensity is due to 2 zz IE , where E z is the component of the local electromagnetic field parallel to the incident direction. It is no more true at some distance (z = 50 nm)) from the metal layer (Fig. 4(d), (e) and (f)) and the main contribution is parallel to the incident field. So it is confirmed that, with our FDTD simulations, the major contribution to the intensity in the very near-field zone is due to the non-propagative component, parallel to the incident wave-vector.

In conclusion, the introduction of surface roughness in the 3D parallel FDTD simulations performed on fractal gold films gives results in good agreement with the experimental nearfield measurements. The local fluctuations in the near-field intensity are quite well reproduced and their statistical behavior is very satisfying compared to the experimentally observed ones. These results could be useful to further the understanding of the local fields. These local field fluctuations are a signature of the inhomogeneous character of the disordered metallic films, as for the fluctuations of the LDOS. So, by modifying the nature of the excitation source, the FDTD simulations could be a useful tool to understand the interactions between such random metallic layers and single nano-emitters deposited upon it.

Fig. 1 .

 1 Fig. 1. Structural characteristics of the gold disordered layer: (a) the 3D AFM topography; (b) a typical profile of the topography; (c) perimeter versus surface of the gold clusters.

Fig. 2 .

 2 Fig. 2. (a) Topography of the gold film used in the FDTD simulation. (b) Enhancement factor ( ) map for an incident wavelength = 650 nm with n = 0.175 and k = 3.13. (c) Same as (b) for = 1000 nm with n = 0.25 and k = 6.82. (d) Black circles: Probability Density Function of corresponding to image (b); blue line: the Log-Normal data fit with parameters values indicated. (e) Black circles: PDF of for image (c); blue line: Log-Normal data fit.

Fig. 3 .

 3 Fig. 3. Spatial distribution of the enhancement factor calculated for = 1 µm with n = 0.22769 and k = 6.4731, (a) incident X linear polarization; (b) incident Y linear polarization; and for = 1.3 µm with n = 0.38797 and k = 8.7971, (c) X polarization; (d) Y polarization. The maximum value of each color bar was set at the mean value of the respective maximum (a and b) and (c and d).

Fig. 4 .

 4 Fig. 4. Maps of the calculated intensities for = 650 nm for two different constant heights: toprow (a, b, c) for z = 12 nm and bottom-row (d, e, f) for z = 50 nm and 3 different directions: (a, d) X component, (b, e) Y component, (c, f) Z component. The incident polarization is along the X direction for all. Same color bar is used for the 3 top-row and the 3 bottom-row images respectively.

  is a well-known technique based on space and time discretizations of the Maxwell curl equations

		H	1	E	(1)
		t
	E	11 HJ	source
	t

Acknowledgments

X. Q. wishes to acknowledge support from the ANR PNANO "NANOEC".