Moufida Maimour
email: moufida.maimour@univ-lorraine.fr

The Impact of Interference aware Metrics on the Performance of Incremental Multipath Routing Protocols in WSN

Keywords: Wireless Sensor Networks (WSN), Interference Aware Metrics, Multipath Routing Protocols

Multipath routing has been considered in the literature to provide bandwidth aggregation for high data rate applications. Because of the phenomena of inter-path interference in wireless networks, the overall throughput is far from being the summation of the provided ones by the different paths individually. One promising solution to this problem is to make use of interference aware metrics to build multiple paths such that the inter-path interference is minimised. In this work, we are interested in assessing the impact of a set of interference aware metrics on the performance of incremental multipath routing in Wireless Sensor Networks (WSN). Incremental multipath routing builds one path per request/reply session and hence is more suitable to be combined to interference aware metrics. We review, implement and evaluate both active and passive monitoring metrics and show that for constrained networks such as WSNs, the latter are more adapted. Additionally, we show that within the passive monitoring metrics, the best one is not necessarily the most expensive to measure.

Introduction

Wireless Sensor Networks (WSN) are composed of low-cost sensor nodes characterised by their scarce resources in terms of energy, bandwidth, processing and storage means. As a consequence, WSN applications were restricted to those that generate low rate data. More recently, high data rate applications have emerged. They range from monitoring and surveillance to wild-life tracking and telemedicine services [START_REF] Borges | Survey on the characterization and classification of wireless sensor network applications[END_REF]. Due to the scarce resources of sensor nodes mainly and in particular to their low-rate radio modules, the available bandwidth is insufficient with respect to the requirements of high data rate applications. Multipath routing has been considered in the literature to provide bandwidth aggregation in the network. However, research results showed that the obtained bandwidth when multiple routes are used is far from being he summation of those provided individually by each path. This is due to the shared nature of the transmission medium where simultaneous use of adjacent routes especially for high rates results in significant inter-path interference. This is referred to as the route coupling problem in the literature and it affects seriously the capacity of wireless networks [START_REF] Gupta | The Capacity of Wireless Networks[END_REF]. Efficient multipath routing protocols have to be carefully designed where both inter-path and intra-path interference have to be minimised in order to increase network throughput.

Solutions to the problem of interference based on special hardware support have been proposed by the research community. They mainly include the use of, if available, directional antennas [START_REF] Noguchi | Adaptive location-aware routing with directional antennas in mobile adhoc networks[END_REF][START_REF] Tarter | Directional Antennas for Convergecast in Wireless Sensor Networks: Are They a Good Idea ?[END_REF] or multiple channels [START_REF] Liao | Path Selection of Reliable Data Delivery in Wireless Sensor Networks[END_REF][START_REF] Da Silva | A distributed algorithm to schedule TSCH links under the SINR interference model[END_REF]. Location information can also help building physically separated paths with low interference as proposed in [START_REF] Zijian | Energy Efficient Collision Aware Multipath Routing for Wireless Sensor Networks[END_REF][START_REF] Bennis | Carrier sense aware multipath geographic routing protocol[END_REF]. Even though these hardware-dependent solutions are efficient in producing non-interfering paths, they still require special hardware to be available. Moreover, localisation algorithms may incur significant overhead in terms of computational and communication complexity. As for multichannel approach, the required channel scheduling activities may introduce significant overhead that may exceed 200 µs on CC2420 chipcon [9] for instance. Due to these issues, routing protocols that require special hardware support may not be cost-effective especially in dense wireless sensor networks.

One promising approach to reduce interference without any special hardware support is the design of routing protocols based on new metrics that consider the amount of interference that could be experienced by the routes to build. Interference aware routing metrics have received a keen interest in the literature of wireless mesh networks. ETX (Expected Transmission Count) [START_REF] De Couto | A High-Throughput Path Metric for Multi-Hop Wireless Routing[END_REF], a popular one, characterises a link loss ratio by estimating the expected number of transmissions, including retransmissions, required to successfully deliver a unicast packet. Its computation is based on the delivery ratio of probes sent on the forward and the reverse directions of a link. In order to improve ETX metric, Expected Transmission Time (ETT) [START_REF] Draves | Routing in Multi-radio, Multi-hop Wireless Mesh Networks[END_REF] considers the link data rate. Still, both ETX and ETT do not consider explicitly the effects of intra-path and inter-path interference. Based on ETX and ETT, other metrics have been proposed by the research community such as interference AWARE (iAWARE) and Interferer Neighbours Count (INX) [START_REF] Langar | Mobility-aware clustering algorithms with interference constraints in wireless mesh networks[END_REF]. Periodic probes exchange in ETX based metrics introduce significant overhead that limits the achievable throughput. As a result, passive monitoring metrics were developed such as Contention Aware Transmission Time (CATT) [START_REF] Genetzakis | A contention-aware routing metric for multi-rate multi-radio mesh networks[END_REF] and Metric for INterference and channel Diversity (MIND) [START_REF] Borges | Routing metric for interference and channel diversity in multi-radio wireless mesh networks[END_REF].

Regarding multipath routing, Wu et al. [START_REF] Wu | Performance Study of a Multipath Routing Method for Wireless Mobile Ad Hoc Networks[END_REF] introduced the correlation factor metric to measure the relative degree of inter-path Independence. The correlation factor of two paths is defined as the number of links connecting two paths to each other. A DSR [START_REF] Johnson | DSR: The Dynamic Source Routing Protocol for Multi-Hop Wireless Ad Hoc Networks[END_REF] based interference aware multipath routing is proposed using this metric. Node-disjoint paths with a small correlation factor are built while taking care of obtaining small length difference between them. DCHT (Delay-Constrained High Throughput) [START_REF] Li | Delay-constrained high throughput protocol for multi-path transmission over wireless multimedia sensor networks[END_REF] is another multipath routing protocol based on EDGE [START_REF] Li | Edge: A routing algorithm for maximizing throughput and minimizing delay in wireless sensor networks[END_REF] with application to video transport. In DCHT, paths construction is based on an improved method for estimating aggregate ETX for a route such that the end-to-end throughput is increased and the delay minimised.

Incremental multipath approach had also been proposed in the literature to deal with the problem of inter-path interference. An iterative path discovery is adopted where only one path is built at once [START_REF] Maimour | Maximally Radio-Disjoint Multipath Routing for Wireless Multimedia Sensor Networks[END_REF][START_REF] Teo | Interference-Minimized Multipath Routing with Congestion Control in Wireless Sensor Network for High-Rate Streaming[END_REF]. Any other path is built by avoiding nodes surrounding already built ones in order to minimise their correlation. The iterative path discovery can be performed along with an interference aware metric. In [START_REF] Radi | Interference-Aware Multipath Routing Protocol for QoS Improvement in Event-Driven Wireless Sensor Networks[END_REF], an ETX based metric is employed. Both proposals in [START_REF] Radi | IM2PR: interference-minimized multipath routing protocol for wireless sensor networks[END_REF] and [START_REF] Bidai | Interference-aware Multipath Routing Protocol for Video Transmission over ZigBee Wireless Sensor Networks[END_REF] make use of a recovery tree combined to a path discovery using metrics that we consider as derived from CATT. This is also the case of the proposed protocol in [START_REF] Hossain | Neighbour-Disjoint Multipath for Low-Power and Lossy Networks[END_REF] which however does not rely on a recovery tree.

In this paper, we assess the impact of using interference aware metrics on the performance of multipath routing protocols. We target the iterative path discovery mechanism to build multiple paths that we believe more relevant to WSN. On the one hand, it does not require neither special hardware support nor a pre-established recovery tree. On the other hand, it allows to achieve low overhead in terms of control messages required to build paths as will be shown in Section 5. In our previous work [START_REF] Maimour | Interference-aware multipath routing for WSNs: Overview and performance evaluation[END_REF], we made a performance evaluation where only one interference metric is considered. In this paper we focus on the impact of five other metrics on the performance of multipath routing that adopts an incremental paths discovery. For the best of our knowledge, this is the first work that assesses the impact of interference aware metrics of the performances of multipath routing protocols. Most of the cited work compares either to a single path or to another multipath routing protocol that does not account for interference. For instance, in [START_REF] Li | Delay-constrained high throughput protocol for multi-path transmission over wireless multimedia sensor networks[END_REF], DCHT is compared to EDGE [START_REF] Li | Edge: A routing algorithm for maximizing throughput and minimizing delay in wireless sensor networks[END_REF] and Directed Diffusion [START_REF] Intanagonwiwat | Directed diffusion for wireless sensor networking[END_REF]. In [START_REF] Hossain | Neighbour-Disjoint Multipath for Low-Power and Lossy Networks[END_REF], the proposed protocol is compared to single path RPL [START_REF]RPL: IPv6 routing protocol for low-power and lossy networks[END_REF] and LOADng [START_REF] Clausen | The lightweight on-demand ad hoc distance-vector routing protocol-next generation (LOADng)[END_REF]. Finally, IM2PR [START_REF] Radi | IM2PR: interference-minimized multipath routing protocol for wireless sensor networks[END_REF] is compared to EERP [START_REF] Boulfekhar | A novel energy efficient and lifetime maximization routing protocol in wireless sensor networks[END_REF] and MSMRP [START_REF] Gao | Multi-path routing protocol with unavailable areas identification in wireless sensor networks[END_REF] that are not interference aware.

The paper is organised as follows. The related work is reviewed in Section 2. Afterwards, in Section 3, we overview the interference aware metrics we chose to consider in our performance evaluation. Sections 4 describes two generic multipath routing protocols, object of our study. Simulation results are discussed in Section 5 before concluding the paper with some future directions. DCHT [START_REF] Li | Delay-constrained high throughput protocol for multi-path transmission over wireless multimedia sensor networks[END_REF] Zeghilet et al. [START_REF] Zeghilet | The Effects of Interference on Video Quality over Wireless S ensor Networks[END_REF] CFMPR-M [START_REF] Ghrab | Correlation-Free MultiPath Routing for Multimedia Traffic in Wireless Sensor Networks, in: Ad-hoc, Mobile, and Wireless Networks[END_REF] Iterative paths discovery No metric used MR2 [START_REF] Maimour | Maximally Radio-Disjoint Multipath Routing for Wireless Multimedia Sensor Networks[END_REF] I2MR [START_REF] Teo | Interference-Minimized Multipath Routing with Congestion Control in Wireless Sensor Network for High-Rate Streaming[END_REF] IAMR [START_REF] Wang | Interference aware multipath routing protocol for wireless sensor networks[END_REF] Metric based LIEMRO [START_REF] Radi | Interference-Aware Multipath Routing Protocol for QoS Improvement in Event-Driven Wireless Sensor Networks[END_REF] IMMR [START_REF] Bidai | Interference-aware Multipath Routing Protocol for Video Transmission over ZigBee Wireless Sensor Networks[END_REF] IM2PR [START_REF] Radi | IM2PR: interference-minimized multipath routing protocol for wireless sensor networks[END_REF] NDM [START_REF] Hossain | Neighbour-Disjoint Multipath for Low-Power and Lossy Networks[END_REF] Figure 1. Summary of the overviewed interference aware multipath protocols

Related Work

Most of multipath protocols in the literature are on-demand reactive routing protocols. When a node (source) has to send data but does not have route information to the sink, it broadcasts a route request message (RREQ) toward the sink. Selected paths are sent by the sink back to the source using route reply messages (RREP). Interference aware multipath routing has generated a great deal of interest in the WSN research community. Minimum interference can be achieved through the use of an interference aware metric. The adopted metric depends on paths discovery mechanism which can be either concurrent or iterative (incremental). In concurrent multipath routing, multiple paths are built using one request/reply session while in the incremental approach, only one path is built at a time. Figure 1 summarises the interference aware multipath routing protocols detailed in what follows.

Interference Aware Concurrent Multipath Routing

One of the first proposed metrics to characterise inter-path interference in multipath routing is the route coupling metric [START_REF] Pearlman | On the impact of alternate path routing for load balancing in mobile ad hoc networks[END_REF]. It describes the average number of nodes in a path that are blocked for receiving packets while a node on another path is transmitting. However, to obtain paths with minimum interference, the network connectivity graph has to be established which implies significant additional calculations especially in dense WSNs.

Wu et al. [START_REF] Wu | Performance Study of a Multipath Routing Method for Wireless Mobile Ad Hoc Networks[END_REF] introduced the correlation factor (CF) metric to measure the relative degree of inter-path Independence. The correlation factor of two node-disjoint paths is defined as the number of the links connecting the two paths. If there is no link between two node-disjoint paths, we say the two node-disjoint paths are unrelated. A link (u, v) is said to be connecting two paths P 1 and P 2 if (u ∈ P 1 and v ∈ P 2 ) or (u ∈ P 2 and v ∈ P 1 ). The total correlation factor of a set of multiple paths is defined as the sum of the correlation factor of each pair of paths. A DSR [START_REF] Johnson | DSR: The Dynamic Source Routing Protocol for Multi-Hop Wireless Ad Hoc Networks[END_REF] based interference aware multipath routing is proposed based on this metric. Node-disjoint paths with a small correlation factor are built while taking care of getting small length difference between them. To compute the correlation factor, each RREP message brings additionally neighbourhood information of the advertised route.

Piggybacking neighbourhood information is also adopted by Zeghilet et al. [START_REF] Zeghilet | The Effects of Interference on Video Quality over Wireless S ensor Networks[END_REF] where RREQ instead of RREP messages are used. Based on neighbourhood information, the sink picks routes with the least common neighbours. This allows avoiding nodes in the neighbourhood of a chosen path. The same mechanism is used in [START_REF] Ghrab | Correlation-Free MultiPath Routing for Multimedia Traffic in Wireless Sensor Networks, in: Ad-hoc, Mobile, and Wireless Networks[END_REF] while the correlation issue is considered in a duty-cycled environment. An issue of this technique is that additional data conveyed by control messages may produce large control overhead in dense networks. This is even more exacerbated in [START_REF] Wu | Performance Study of a Multipath Routing Method for Wireless Mobile Ad Hoc Networks[END_REF] where path selection is done by the source. DCHT (Delay-Constrained High Throughput) [START_REF] Li | Delay-constrained high throughput protocol for multi-path transmission over wireless multimedia sensor networks[END_REF] follows the basics of Directed Diffusion [START_REF] Intanagonwiwat | Directed diffusion for wireless sensor networking[END_REF] with multipath extension. Instead of pure delay, it introduces a path p cost metric Cost p , the product of expected transmission count (ETX) and delay to maximise the throughput and improve the delay performance :

Cost p = ET X α p × Delay β α and β are non negative integers and ET X p is computed using :

ET X p = Max N-3 i=0 ( i+2 j=i ET X j )
where N is the number of hops in the path and ET X j is the ETX value of the j th hop. In [START_REF] Li | Delay-constrained high throughput protocol for multi-path transmission over wireless multimedia sensor networks[END_REF], it is suggested to compute ETX based on SNR (Signal to Noise Ratio) read from each exploratory packet at intermediate nodes. The ETX of the three previous hops is considered to ensure intra-path interference. It is worth noting that the authors do not provide the method to obtain SNR measures. Results show that DCHT applied to video streaming over WSNs obtains higher throughput than its single path counterpart (EDGE [START_REF] Li | Edge: A routing algorithm for maximizing throughput and minimizing delay in wireless sensor networks[END_REF]). However, there was no comparison with another interference aware multipath routing protocol. One limitation of DCHT is that it considers more link quality and focuses on intra-path rather than inter-path interference in multipath routing.

Interference Aware Incremental Multipath Routing

Iterative paths discovery allows obtaining low inter-path interference either by explicitly preventing nodes that may interfere with an already built path from being considered in subsequent paths construction or by using an adequate interference aware metric. In MR2 (Maximally Radio-disjoint multipath Routing) [START_REF] Maimour | Maximally Radio-Disjoint Multipath Routing for Wireless Multimedia Sensor Networks[END_REF], neighbour nodes of an already built path are put in a passive mode where the radio can be switched off. This is twofold : energy is saved and interfering nodes are prevented from taking part in subsequent routes. Regarding interference awareness, I2MR (Interference-minimised multipath routing) [START_REF] Teo | Interference-Minimized Multipath Routing with Congestion Control in Wireless Sensor Network for High-Rate Streaming[END_REF] adopts the same idea while assisted by location information on the source and the destination nodes.

In order to consider the fact that interference range can reach twice the communication range [START_REF] Zhou | Impact of radio irregularity on wireless sensor networks[END_REF], interferencezone marking can be done up to two-hop neighbours of intermediate nodes of an already built route. IAMR (Interference aware multipath routing) [START_REF] Wang | Interference aware multipath routing protocol for wireless sensor networks[END_REF] marks the shortest path nodes as well as its one-hop neighbours. Afterwards, it builds two paths surrounding the marked zone. Note that two-hop neighbours marking may increase routes hop count which may result in more losses and large buffer at the destination to reorder packets.

Based on ETX metric, LIEMRO (Low-Interference Energy-Efficient Multipath Routing Protocol) [START_REF] Radi | Interference-Aware Multipath Routing Protocol for QoS Improvement in Event-Driven Wireless Sensor Networks[END_REF] makes use of the iterative approach to build multiple routes. A source node starts to establish the first path by transmitting an RREQ message towards the sink node. Whenever a node receives an RREQ packet, it computes the transmission cost for every neighbouring node not included in any path from the source to the sink. The newly estimated accumulated ETX cost is broadcast if it is lower than the current one. The cost update process is performed at the initialisation phase and whenever a node obtains a new transmission cost to the sink. The active monitoring required by ETX based metrics impacts the overall achieved throughput as will be shown in Section 5. Moreover, the numerous updates may consume significant computation resources that are limited in sensor nodes.

Flowing packets overhearing mechanism is also adopted by IMMR (Interference-Minimised Multipath Routing) [START_REF] Bidai | Interference-aware Multipath Routing Protocol for Video Transmission over ZigBee Wireless Sensor Networks[END_REF]. IMMR is an interference aware extension of ZMTR [START_REF] Bidai | Node disjoint multi-path routing for ZigBee cluster-tree wireless sensor networks[END_REF], a multipath routing protocol targeted to ZigBee lowrate networks. Routes are built incrementally based on a spanning tree rooted at the sink. The first one is the tree path that follows the child-parent links toward the sink. Subsequent ones are selected based on a metric called Interference Level (IL) estimated based on the overhearing of data packets that already circulate in the vicinity as well as potential path explore messages. The source selects paths with the least IL when more than the required number of paths are discovered. Simulation results show that IMMR achieves better performances with respect to single and multipath routing without interference awareness. Following a similar approach, authors of [START_REF] Radi | IM2PR: interference-minimized multipath routing protocol for wireless sensor networks[END_REF] build a minimum cost recovery tree and propose Interference-Minimised Multipath Routing Protocol (IM2PR). The recovery tree construction is based on a quality metric estimated using beacon messages exchange at the initialisation phase. The first path construction is triggered by an RREQ issued by the source toward the sink. A node is chosen to be in this path using a metric with two components : a transmission cost to this node based on link quality and an energy cost based on the remaining energy with respect to initial one. Subsequent paths are built by considering the interference level of the next node. The interference level of node v is measured by how many neighbours of v are in the other active paths. These neighbours are identified using similar mechanism used in IMMR and LIEMRO i.e. packets overhearing exploitation. Paths in IM2PR are built one by one even if the discovery of a subsequent path can be done before the previous path discovery is fully completed. Considering the fact that sensors have limited computation resources, the numerous required updates of the different measures to estimate links quality as well as paths interference level in IM2PR are resource consuming.

More recently, the work in [START_REF] Hossain | Neighbour-Disjoint Multipath for Low-Power and Lossy Networks[END_REF] also adopted the iterative approach and proposed NDM (Neighbour-Disjoint Multipath). Less interfering backup paths with respect to a primary path are built based on a combined metric that includes hop count and correlation factor. The primary path that is the shortest (hop count) one between a sensor and the sink is selected. Then a set of backup paths with nodes that are not a neighbour of any node of the primary path except source or sink, are selected. NDM makes use of the RREQ/RREP messages of LOADng [START_REF] Clausen | The lightweight on-demand ad hoc distance-vector routing protocol-next generation (LOADng)[END_REF], a lightweight AODV reactive routing protocol. Unicast and broadcast RREP messages are transmitted by the sink and are sent one after another. Unicast RREP are processed as it would have been for a traditional AODV protocol. When a broadcast RREP is received, the node checks if it is inside the primary path and then it it is a neighbour to any node of the primary path. Accordingly, it updates the correlation and hop-count metrics. Only RREPs with a better correlation metric or a lower hop-count in case of equal correlation parameter are broadcast.

At this stage, it is worth noting that all the above multipath routing protocols were evaluated and compared to a single path and/or multipath routing without interference awareness.

interference aware Routing Metrics

In this section, we overview some of the interference aware routing metrics proposed in the literature. We suggest to integrate these selected metrics to an incremental multipath protocol that will be presented in the following Section. Building subsequent paths while others are in service allows to capture the potential interference level of the former with respect to the latter by the interference aware metrics. The considered metrics are chosen based on the following criteria. First, we consider a sensor network where each node is equipped with a radio that operates using only one channel. In addition to the multi-channel cost, its integration, usually results in non isotonic metrics [START_REF] Borges | Routing metric for interference and channel diversity in multi-radio wireless mesh networks[END_REF][START_REF] Subramanian | Interference aware routing in multi-radio wireless mesh networks[END_REF] which require employing a virtual network with isotonicity verified in order to be able to use Bellman-Ford or Dijkstra's algorithms. Therefore, if a metric includes a multi-channel component, this component will not be considered in our study. Second, a selected metric has to be of a reasonable complexity to be suitable to WSN. Finally, we selected metrics that are the most representative of those of the literature.

Overall, interference aware routing can be classified into two categories, namely active and passive monitoring metrics. As opposed to the latter, the former require periodic probing messages to assess the amount of the potential interference a link may experience.

Active Monitoring Metrics

In this section, we describe ETX metric in addition to two other ETX-based metrics iAWARE and INX.

ETX [10]

ETX (Expected Transmission Count) measures the link loss ratio using the expected number of transmissions, counting retransmissions, needed to successfully receive a unicast packet through a link. The path metric is the summation of each link ETX in the path and the one with the minimum metric value is selected. A node u computes the ETX of the link to a node v using the delivery ratio of probes sent on the forward (p f ) and reverse (p r ) directions. These delivery ratios are, respectively, the fraction of successfully received probes from u announced by v and the fraction of successfully received probes from v, at the same period. The ETX of link l(u, v) is given by :

ET X l = 1 p f × p r
ETX is bidirectional and appropriately handles asymmetry by incorporating loss ratios in both direction. Measurements on a wireless test-bed show that ETX penalises routes with more hops which have lower throughput due to interferences between different hops of the same path [START_REF] De Couto | A High-Throughput Path Metric for Multi-Hop Wireless Routing[END_REF]. Although ETX does very well in homogeneous environments, it does not perform well in environments with different data rates. This is due to the small size of probes compared to typical data packets in addition to the fact that it does not consider the data rate at which the packets are transmitted. Moreover, ETX does not capture the interference experienced by the links completely since the periodic probes are sent at a slow interval (usually 1 sec) which does not reflect how busy a link is [START_REF] Subramanian | Interference aware routing in multi-radio wireless mesh networks[END_REF].

iAWARE [37]

The interference AWARE (iAWARE) routing metric has two components. The first component is based on the physical interference model. It uses SNR (Signal to Noise Ratio) and SINR (Signal to Interference and Noise Ratio) to continuously reproduce neighbouring interference variations. In the physical model, a communication between nodes u and v on the link l(u, v) is successful if the SINR at the receiver v is above a certain threshold. The second component exploits the channel diversity to find paths with least intra-flow interference. In this work, we only consider the first component since in our network model, we assumed that only one channel is available.

iAWARE is also based on the Expected Transmission Time (ETT) metric [START_REF] Draves | Routing in Multi-radio, Multi-hop Wireless Mesh Networks[END_REF] which estimates the time a data packet requires to be transmitted successfully to each neighbour. It improves upon ETX by capturing the data rate used by each link, that is :

ET T l = ET X l × S B
where the number of expected transmissions is multiplied by the link bandwidth giving the time required to transmit the packet. S is the size of the used packet and B is the bandwidth (raw data rate) of the link.

Being based on the physical interference model. iAWARE employs a correlation between the ratio S INR/S NR and ETT, thus, capturing the effects of variation in link loss-ratio, differences in transmission rate, as well as inter-path interference in a dynamic way. iAWARE metric of link l is given as follows :

iAWARE l = ET T l IR l
where IR l = min(IR l (u), IR l (v)) is the interference ratio of link l(u, v). IR l (u), IR l (v) ∈]0, 1] are the interference ratio of nodes u and v respectively. For a node u, it is given by :

IR l (u) = S INR l (u) S NR l (u)
S NR l (u) and S INR l (u) are given by :

S NR l (u) = P u (v) N S INR l (u) = P u (v) N + w∈η(u)-v τ(w)P u (w)
where P u (v) denotes the signal strength of a packet from node u to node v. N is the background noise. η(u) denotes the set of nodes from which u can hear a packet and τ(w) is the normalised rate at which node w generates traffic It is worth noting that when there is no interference IR l = 1 which implies that iAWARE is simply ET T l . This latter captures the link loss ratio as well as its packet transmission rate. ET T l is weighted with IR l to capture the interference experienced by the link from its neighbours. A link with low ETT and high IR will have a low iAWARE value. The lower the iAWARE of a link, the better is the link. The cumulative path metric iAWARE(p) of a path p is defined as iAWARE(p) = l∈P iAWARE l iAWARE gives more importance to ETT than interference. Moreover, iAWARE only considers interference at the physical layer and does not capture the interference at the MAC layer. The required measurements of the SNR, the signal strength and the background noise are not easy to obtain unless a cross-layer method is employed. This results also in considering all heard packets by the routing layer of a node even if they were not intended to this node.

INX [12]

Interferer Neighbours Count (INX) routing metric is improved on the basis of ETX. It takes into account interference through measuring the sum of the transmission rates of the links that can interfere with the transmission on link l as :

INX l = ET X l × j∈S l r j
where r j is the data rate of link j; it is set to the same value for all links in the network. ET X l is the expected transmission time of link l. S l is the set of interfering links. For a given path p, INX(p) = l∈p INX l /N where N is the number of links in the WSN.

Passive Monitoring Metrics

All active monitoring metrics overviewed in the previous section are based on ETX. Therefore, probing packets are exchanged to acquire delivery ratio of each link in both directions. Since probe messages are usually shorter than data packets, a node may experience higher losses for actual packets [START_REF] De Couto | A High-Throughput Path Metric for Multi-Hop Wireless Routing[END_REF]. Moreover, the number of broadcast probes in an n-node network is O(n) which causes additional overhead especially in dense networks. Due to the restricted resources of WSN, passive monitoring is desirable to help acquiring measurements without introducing extra overhead disturbing the normal operation in the network. Providing sufficiently accurate interference metric without active monitoring is challenging.

CATT [13]

Contention aware transmission time (CATT) takes into account the effect of interferers on the transmission time of packets over the interfered link. It depicts the intra-path and inter-path interference as well as the traffic load in a uniform way by adding up the delays of the interfering neighbour links that are one and two hops away. CATT of link l is defined as :

CAT T l = j∈Nl L j R j
where L j is the packet size of link j; it is set to the same value for all packets in the network. R j is the data rate of link j. N l is the set of links (including l) whose transmission can interfere with the transmission on link l. The CATT associated to a path p is CAT T (p) = l∈p CAT T l It is worth mentioning that CAT T LD , an extension of CATT metric, was also proposed in [START_REF] Genetzakis | A contention-aware routing metric for multi-rate multi-radio mesh networks[END_REF]. CAT T LD as opposed to CATT makes use of ETX, and thus is an active monitoring metric, to include loss-dependence as follows :

CAT T LD l = ET X l × j∈Nl L j R j

MIND [14]

MIND (Metric for INterference and channel Diversity) is an interference aware routing metric that relies on passive measurements to depict interference and traffic load. Like iAWARE, it is composed of two components. One that concerns inter-path interference and the other exploits channel diversity to limit intra-path interference. When considering one channel, MIND employs the interference ratio as done in iAWARE (Section 3.1.2) to capture the experienced interference based on the physical interference model. Additionally, as the Channel Busy Time (CBT) provides a more precise way of measuring traffic load [START_REF] Wu | Passive and accurate traffic load estimation for infrastructure-mode wireless lan[END_REF], CBT is also considered. That is :

MIND(p) = l∈p ((1 -IR l ) * τ) * CBT l
where IR l is the interference ratio of link l (Section 3.1.2) and τ is a configurable parameter used to provide a higher weight to interference. CBT is defined by : CBT = T otalT ime -IdleT ime T otalT ime

Multipath Routing Protocols Description

The aim of this paper is not to propose yet another multipath routing protocol, only to assess the advantage the incremental discovery approach can make when integrating interference aware metrics. Instead of giving a detailed description of a specific multipath protocol, we consider two generic protocols that will be the object of the performance evaluation conducted in Section 5. The first one corresponds to the concurrent multipath routing, we call MP, in which only one request/reply session is required to build multiple paths. The second one referred to as IM2R, consists in an interference aware iterative (incremental) paths discovery where only one path is built at once. Different interference metrics (considered as path cost) are integrated to this latter in order to assess their impact on its performance.

Both MP and IM2R are on-demand reactive routing protocols. The source transmits an RREQ packet with its ID. Each intermediate node that receives an RREQ, appends its own ID before re-broadcasting the packet. This allows avoiding loops as well as selecting disjoint paths since the list of nodes that belong to each traversed path is available. RREP messages also carry the whole path information between the source and the sink. They are broadcast throughout the advertised route instead of unicast to take into account non symmetric links. In fact, the reverse path does not necessarily exist. Based on a received RREP, each intermediate node maintains a route table that indicates the path to the sink. The source maintains a similar table with multiple entries, one per built path. Finally, each node in the network has to maintain an up-to-date neighbours table . In what follows, particular behaviour of the two considered multipath protocols are described.

Concurrent Paths Discovery (MP)

In this approach, all routes are built based on the flooding of one RREQ. An important point to consider is RREQ strategy suppression at intermediate nodes since dropping all duplicate RREQs may reduce drastically the probability to find multiple disjoint paths. Thus, RREQ caching and rebroadcast are required in multipath routing. When an intermediate node receives a new RREQ, it records its sequence number, the source ID and the cost of the corresponding path cost from the source to this node in a path record. An RREQ is considered as a duplicate if both its sequence number and source ID have already been recorded. Only duplicate RREQs that bring higher route cost are dropped.

Upon the reception of an RREQ by the sink, a timer is set and the advertised path nodes along with the corresponding cost are recorded. On timeout, the sink selects disjoint paths with the least cost and sends back to the source an RREP per chosen path. An intermediate node that belongs to the advertised path by a received RREP, updates its routing table accordingly. Finally, upon the reception of the RREP, the source adds an entry to its routing table for the advertised path.

interference aware Metric-based Incremental Multipath Routing (IM2R)

In IM2R, the source sends a first RREQ to build one route and transmits subsequent ones each time it receives an RREP with a new disjoint path whenever the needed number of paths is not achieved. As opposed to the concurrent approach, this requires that RREQ and RREP messages contain an additional field that gives the rank of the built path. Additionally, only one path record with the list of its nodes is maintained at the intermediate nodes.

When an intermediate node, that does not belong to another path, receives an RREQ, it checks if it is the first time it receives such RREQ. If so, the newly received path is saved with its cost and a timer is set to allow the reception of other RREQs with paths having better cost. During the timer period, subsequent RREQs are considered by the intermediate node if they convey a lower cost and the path record is updated. On timeout, an RREQ is sent using broadcast. It conveys information on the best path among those advertised by the received RREQs during the timer period.

Upon the reception of an RREQ, the sink adds the conveyed path to its paths table and sets a timer for a given period to allow considering subsequent RREQs. On timeout, the sink selects one path with the best cost metric value. An RREP that contains the elected path record is sent back to the source. When an RREP is received by an intermediate node that belongs to the advertised path, an entry is set in the routing table using the information provided by this RREP. Otherwise, the RREP is ignored. Finally, upon the reception of the RREP, the source adds an entry to its routing table. If the number of required paths is not achieved then a new RREQ is sent by the source.

In addition to the traditional hop count (HC) metric, to be interference aware, metrics defined in Section 3, namely ETX, iAWARE, INX, CATT and MIND are implemented along with the incremental multipath routing. Moreover, we also considered the metric used in DCHT presented in Section 2.1 that we refer to as ETX3 since it considers three-hop ETX value. All ETX based metrics require that each entry in the neighbours table is extended to maintain for each neighbour, the number of received probes from this neighbour and those received by the neighbour from this node. In turn, iAWARE and MIND require to record additionally, the receive power from each neighbour, an estimation of its neighbours's transmission rate as well as the interference ration (IR) of its links. Finally, ETX3 requires three cost fields in RREQs to keep ETX measures for the three last hops.

Simulation Results

Using Castalia [39], a wireless sensor network Omnet++ based simulator, we implemented both concurrent multipath variant (MP) as well as the interference aware incremental variant (IM2R). Metrics ETX, iAWARE, INX, ETX3, CATT and MIND were implemented along with IM2R as well as hop count (HC) metric in order to assess their impact on its performance. In active monitoring metrics, probes are sent every one second and ETX value is computed every 10 seconds.

Sensor nodes are equipped with a CC2420 radio operating at 250 kbps transmitting at a power of 0 dBm which results in a transmission range of about 46.42 meters when an ideal radio mode is assumed. Nodes are deployed in a randomised grid and spaced in average 21 meters from each other which results in a mean node degree of 12. The sink and the source are placed respectively at the lower left and the upper right corners of the sensors area.

To model interference, we adopted the Additive Interference model provided by Castalia. At the MAC layer, we made use of a contention based CSMA. Each sensor node has a 32-message MAC buffer and a 16-frame radio one. In our simulations, we varied the data packet rate from 20 to 50 packets per second with a data packet payload of 64 bytes giving a useful data rate ranging from 10.24 Kb/s to 25.6 Kb/s. The transmission rate is equally assigned on the built routes. In all the performed simulations, there was no buffer overflow at the source even with the highest data rate. At least 20 simulations using different seeds are performed for each scenario.

We limit our simulations to the case of two paths even if more paths can be discovered. This is sufficient to emphasise the performances of each metric with respect to each other. Obtained results can be generalised to more paths as well as for the case of multiple sources [START_REF] Maimour | Interference-aware multipath routing for WSNs: Overview and performance evaluation[END_REF]. Simulations we performed using up to three sources showed that the relative performance of the different multipath protocols is still the same regardless the number of involved 1 where default values are highlighted in bold.

Routing Overhead Versus Reliability

One important issue in routing protocols and in particular multipath ones is the amount of overheard caused by the different control messages required to build and maintain routes. These messages include mainly the RREQ and RREP messages in addition to the probe messages needed to maintain neighbourhood information. In our simulations, the number of control messages is normalised by dividing their number by the number of sensors. Figure 2a plots the boxes that summarise the normalised number of RREQs received by the sensors in the network during the path discovery phase for MP and IM2R combined with the different implemented metrics that occur at the x-axis. Each box plot displays the distribution of the number of RREQs obtained by each protocol variant. The central box spans the the interquartile range with a segment giving the median. The mean value is represented with a black square inside each box. Minimum and maximum achieved values are represented by the whiskers above and below each box.

As opposed to the incremental paths discovery, we observe that there is an implosion in the number of the exchanged RREQs in MP (concurrent discovery). There is an average of more than 250 RREQs reception per node in MP while its maximum does not exceed 50 in IM2R whatever the used metric. The very big number of RREQs in MP is due to the RREQ suppression strategy used in which a minimum number of RREQs have to be rebroadcast to allow multiple disjoint paths discovery. This problem does not exist in IM2R whatever the used metric since paths are discovered sequentially.

Figure 2b depicts the normalised number of received RREPs. Even though, one RREP is sent by the sink for each chosen path in all the evaluated variants, we note that their overall number varies. Regardless the used metric, IM2R produces more RREPs than MP. This is still the case when the number of built paths is increased from 2 to 6 as shown by Figure 3a that plots the amount of exchanged RREPs for MP (concurrent) and IM2R combined with HC (incremental) as a function of the number of built paths. The observed higher overhead in terms of RREPs in IM2R is due to the fact that paths built by the incremental multipath routing are likely to be longer as depicted in Figure 3b. This results in more RREPs to be exchanged along each path from the sink to the source. We note however that with respect to the concurrent paths discovery in MP, this overhead in the incremental approach, drops from 39% of additional RREPs for two built paths to 14% when the number of paths is rised to 6.

Figure 4a brings the obtained path lengths for MP and IM2R with the seven considered metrics. Path length may give an indication on the average delay required for a packet to go from the source to a destination since delay is in part related to the path length. We can see that regardless of the used metric, IM2R produces longer paths compared to MP. This is due to the suppression strategy used to limit the number of exchanged RREQs in MP. Compared to IM2R, the sink selects the paths based on a larger number of RREQs, thus increasing the probability of selecting the shortest ones. Our simulation results show that the first selected route in MP is in average shorter than those selected first in IM2R. We can also see that interference aware metrics in IM2R produce longer paths than the hop count metric. In the former, subsequent paths are built so they are the least correlated to previous ones which results in longer paths. Nodes in the vicinity of already built paths are likely to be penalised by the interference aware metrics.

Reliability is one important aspect of any WSN application. Even if reliability is mainly a matter of the transport layer, the routing layer must still do its best. In our performance evaluation, data reliability is assessed using the percentage of data packet successfully received by the sink with respect to the total number of sent messages by the source node. Figures 4b and4c show respectively the success ratio and the throughput achieved by MP and IM2R using the different implemented metrics. The throughput is computed using the amount of useful data (payload) received by the sink without considering the different headers of the transmitted packets. Overall, from Figures 4b and4c, two main observations can be made. First, while generating the least amount of routing overhead in terms of the number of RREQs and RREPs, passive monitoring metrics (except ETX3) and even HC allow better performances in IM2R when compared to MP. Second, we see that in terms of both success ratio and throughput, IM2R combined with active monitoring metrics exhibits poor performances when compared to MP. This is mainly due to the periodic probes used in active monitoring which consumes bandwidth. All active monitoring metrics are based on ETX that measures link quality rather than its experienced interference. The poor performance of ETX3 when compared to other passive monitoring metrics can be explained by the fact that ETX3 deals with intra-path instead of inter-path interference. In fact, ETX3 builds the paths with the highest correlation factor. Figure 4d plots the correlation factor we computed off-line based on the network connectivity graph and the built paths by the different evaluated protocols. The interference range is set to a slightly higher value than the transmission range. We note that ETX3 obtains an average correlation factor that is greater than the one obtained in MP and HC i.e. IM2R without an interference aware metric. Despite that when compared to MP, they achieve less performances in terms of throughput and success ratio, active monitoring metrics build less interfering paths. Among all metrics, CATT allows to build the least correlated paths which confirms the highest performances obtained in terms of success ratio and throughput. Finally, Figure 4d confirms that less correlated paths are likely to be longer (Figure 4a) and hence justifying the amount of additional RREPs exchanged.

Active Monitoring Metrics

To get more insight on the performance of active monitoring metrics when compared to each other, we performed simulations for different data transmission rates. Figure 5 plots the mean throughput achieved by IM2R using iAWARE, ETX and INX as a function of the transmission rate. We can see that ETX achieves higher throughput than iAWARE when the data transmission rate is increased and that INX outperforms both ETX and iAWARE metrics. This is mainly due to the fact that INX allows to build the least correlated paths as shown by Figure 4d. We can also note that iAWARE paths are less correlated than those produced with ETX. Nevertheless, the achieved throughput by iAWARE is not superior. The reason behind that is the additional overhead introduced by the technique of packet pairs [START_REF] Keshav | A control-theoretic approach to flow control[END_REF] to estimate ETT values required by iAWARE. Furthermore, employing ETT may also lead to links quality to be overestimated [START_REF] Borges | The impact of interference-aware routing metrics on video streaming in Wireless Mesh Networks[END_REF]. Finally, iAWARE paths are more correlated than INX ones and hence iAWARE is less interference aware. This is due to the fact that iAWARE only considers interference at the physical level and does not capture interference at the MAC layer [START_REF] Genetzakis | A contention-aware routing metric for multi-rate multi-radio mesh networks[END_REF].

Passive Monitoring Metrics

In what follows, we are interested in the performances of passive monitoring metrics. Figure 6a depicts the evolution of the success ratio for IM2R when passive monitoring metrics are employed as a function of the number of nodes ranging from 225 to 729. As can be expected, we can see that the success delivery ratio decreases with the network size. The sink and the source become more distant from each other since they are located at opposite corners. Thus, the number of transmission and failures increases since traversed routes are longer. Regardless the network size, as opposed to ETX3, we can see that MIND and CATT allow higher success ratio than the traditional HC metric. The same observations can be done following Figure 6b that plots the achieved throughput by the different passive metrics.

The main reason for the above observations is that paths produced by CATT and MIND are less correlated than those of HC, that is in turn less correlated when ETX3 is employed as a metric. This is confirmed by Figure 6c that plots the average value of the correlation factor of the built paths as a function of the network size. IM2R routing using HC takes advantage from the overhearing property. When RREQs are broadcast to establish the following path, data packets are already flowing through the first path. Due to interference with this latter, second path RREQs are most likely to be lost in the vicinity of the first path. As a result, these RREQs have lower probability to reach the sink and hence the second path is more likely to be less interfering with the first one.

When the data transmission rate varies, we obtain almost the same behaviour for the different evaluated passive metrics. As shown by Figure 7a, that plots the throughput as a function of the transmission rate, CATT achieves the highest throughput. CATT is able to estimate the interference degree in the neighbourhood (up to two links away) of a path in a more accurate manner than the other metrics. MIND allows to build less interfering paths than HC since it considers how busy is the channel in the vicinity. The use of CBT instead of ETT allows MIND to build less interfering paths than iAWARE (Figure 4d). ETX3 achieves the worst throughput since it considers more link quality and intra-path rather than inter-path interference. Once again, the mean correlation factor depicted in Figure 7b for different loads confirms that the achieved throughput depends on the degree of correlation of built paths. In fact, CATT builds the least interfering paths while ETX3 paths are more correlated even when compared to the HC metric. (c) Mean correlation factor for passive metrics as the network size is varied.

Conclusion

In this paper, we were interested in evaluating the impact of a set of interference aware metrics on the performance of incremental multipath routing in WSN. We first, considered traditional (concurrent) multipath routing where multiple paths are built using one request/reply session. We showed that the concurrent multipath routing introduce considerable amount of routing overhead when compared to the incremental approach. Moreover, even without interference awareness (using hop count as a metric), incremental multipath achieves higher performances in terms of both throughput and success ratio.

Naturally, incremental multipath is obviously more suitable to integrate and take benefit of interference aware metrics. We qualified these metrics either as active or passive monitoring metrics. The former require periodic probing messages to assess the amount of the potential interference a link may experience. We showed that even if active monitoring metrics are able to build less interfering paths, they still obtain poor performances in terms of throughput and success ratio. This is due to the bandwidth overhead introduced by the probe messages required to estimate ETX and/or ETT component in the corresponding metrics.

Passive monitoring metrics, except ETX3, allow to achieve better performances when compared to active ones. ETX3 considers more link quality and intra-path rather than inter-path interference. The best performance can be obtained using CATT followed by MIND. We also showed that thanks to packets overhearing, HC combined with an increment multipath is somehow interference aware and can obtain better performances when compared to traditional concurrent multipath. With respect to metric estimation complexity, we argue that CATT is once again the best candidate and is best suited for WSN in which nodes are characterised by their scarce resources. As a future work, we expect to implement the passive metrics for experimentation on a real WSN test-bed.

  period of time. τ(w) is 1 when node w sends out packets at the full data rate supported. τ(w) is used to weight the signal strength from an interfering node w as τ(w) gives the fraction of time node w occupies the channel.

Figure 2 .

 2 Figure 2. Control messages

Figure 3 .

 3 Figure 3. Concurrent (MP) vs. incremental (IM2R) multipath routing

Figure 4 .

 4 Figure 4. Performances in a 400-node network and a transmission rate of 25 packets per second

Figure 5 .

 5 Figure 5. Mean throughput for active metrics as the data rate is varied.

Figure 6 .

 6 Figure 6. Performances of passive metrics as the network size is varied.

Figure 7 .

 7 Figure 7. Performances of passive metrics as the data rate is varied.

Table 1

 1 In what follows, unless stated otherwise, we present our simulation results based on a 400-node sensor network where the source sends 25 data packets (64-byte) per second. Simulation parameters are summarised in Table

	Simulation Parameters
	Number of sensors	225, 324, 400, 529, 625, 729.
	Transmission rate	20, 25, 30, 35, 40, 45, 50 pps.
	Data packet payload	64 bytes.
	Node average degree	12.
	Transmission power	0 dBm
	Transmission range	46.42 m (average)
	Collision model	Additive interference
	Active monitoring metrics	ETX, iAWARE, INX
	Passive monitoring metrics	HC, ETX3, CATT, MIND
	sources.


domain: sde spi

text:
Yannick Abautret
email: yannick.abautret@fresnel.fr

Dominique Coquillat

Myriam Zerrad

Xavier Buet

Ryad Bendoula

Gabriel Soriano

Nicolas Brouilly

Daphné Heran

Bruno Grezes-Besset

Frédéric Chazallet

Daphné Héran

Bruno Grèzes-Besset

Claude Amra

Terahertz probing of sunflower leaf multilayer organization

come

Introduction

Over the past decades, climate change has challenged agronomist to find new techniques and plant varieties genetically adapted to water stress. Indeed, lack of water and drought episodes that threaten several regions of the world in the close future have increased the number of scientific investigations on new phenotyping tools and techniques. Until now, phenotyping has required statistical analysis with a huge number of data to measure and process. In this context, the availability of massive characterization vectors, like phenotyping platforms, offers new opportunities for the development of optical investigation tools. For example, multi-spectral and infrared imaging are already operational at the level of the crop [START_REF] Suárez | Assessing canopy PRI for water stress detection with diurnal airborne imagery[END_REF][START_REF] Gitelson | Remote estimation of canopy chlorophyll content in crops[END_REF].

In this new context, there is also a growing demand to develop non-contact (optics, acoustics, mechanics... ) techniques to analyze leaves microstructure. Among them are optical techniques that have been used for decades to probe elastic and non-elastic properties of the plants [START_REF] Onoda | Global patterns of leaf mechanical properties: Global patterns of leaf mechanical properties[END_REF]. Luminescence properties were extensively explored and mainly provide chemical information (chlorophyll-related) about the samples and photosynthesis [START_REF] Baker | Applications of chlorophyll fluorescence can improve crop production strategies: an examination of future possibilities[END_REF][START_REF] Genty | The relationship between the quantum yield of photosynthetic electron transport and quenching of chlorophyll fluorescence[END_REF][START_REF] Murchie | Chlorophyll fluorescence analysis: a guide to good practice and understanding some new applications[END_REF], while few information concerns the opto-geometrical properties, that is, the thicknesses and indices of the layers which constitute the leaf. Optical coherence tomography (OCT) is another modern technique which allows to probe the thickness of the sample and may allow to emphasize the multilayer structure [START_REF] Sapozhnikova | Visualization of plant tissues by optical coherence tomography[END_REF]. However, the depth probe is often limited around 50 m, a value which may be decreased depending on the leaf heterogeneities responsible for high scattering that is plant-dependent.

More classical are the spectro-photometric techniques though they have today led to the development of hyperspectral imaging systems. Relevant information [START_REF] Jacquemoud | PROSPECT: A model of leaf optical properties spectra[END_REF][START_REF] Féret | PROSPECT-D: Towards modeling leaf optical properties through a complete lifecycle[END_REF] is collected with these systems but they do not provide accurate data about the opto-geometrical parameters of the leaf. Indeed, the leaf structure gathers rough interfaces and inhomogeneous bulks which totally scatter the incident light and loose the coherence properties of light; the resulting speckle pattern is complex (see [START_REF] Ghabbach | Depolarization and enpolarization DOP histograms measured for surface and bulk speckle patterns[END_REF]) and does not allow to retrieve the information about the multilayer structure, even when exact electromagnetic theories are used.

In order to solve this last point (related to the complexity of the speckle pattern created by leaf heterogeneities), one basic idea consists in increasing the wavelength ( ) of light. Indeed we know that the weight of the inhomogeneities decreases at large wavelengths, and that the optical resolution is given by the probing wavelength. Surface light scattering for instance would be strongly reduced for low roughness-to-wavelength ratios ( / ), with the root-mean-square of the leaf roughness. Since roughness measurements obtained with white light interferometry microscopy indicate that is of the order of a few micro-meters, a wavelength in the mm range would guarantee that most of the reflected or the transmitted energies is carried by the specular beams, that is, in the directions of Snell/Descartes law. This is the reason why this work is focused on a THz investigation of leaves. As the result of negligible scattering in the THz frequency-range, the leaf structure is seen by the THz system as a quasi-planar homogeneous multilayer, and this allows to use most techniques of the optical thin film [START_REF] Amra | Electromagnetic Optics of Thin-Film Coatings: Light Scattering, Giant Field Enhancement, and Planar Microcavities[END_REF] (optical coatings) community to retrieve the opto-geometrical parameters (the organization) of the leaf from spectral data. Roughly speaking, for single layers these techniques take full advantage of the location of extrema that can be seen in the wavelength variations of reflection and transmission, and of the contrast between these extrema. For more complex organizations, a specific RET is required.

To summarize, thin film RET is used in this paper to analyse the leaf multilayer through THz pulsed spectroscopy data measured by reflection and transmission. By comparison with optical thin film techniques, we must stress on a number of key differences (in addition to the different wavelength ranges and associated technologies for emitters, detectors and components in the optics and THz regimes). The first difference is related to the high absorption level of plants in specific spectral regions (which are mainly due to water and pigments), which does not exist in high precision multi-dielectric coatings; such difference will give more complexity to the inverse problem to solve in leaves. The second difference is connected with the nature (temporal or spectral) of the source; indeed, while broad-band CW sources are currently used in optics, this is not the case for THz sources, for which reason broad-band THz sources most often work in a temporal (pulsed) regime. Though the theory is the same if we use the time-frequency equivalence, this will require to discuss and adjust a number of numerical procedures. This paper presents in section 2 the THz technics and main caracteristics used for this study. A validation campaign on samples is performed and depicted in section 3. Once the validation step is done, the study focuses on the sunflower leaf sample in section 4. Finally, results and perspectives of this work are summarized in section 5.

THz Time-Domain Spectroscopy (TDS)

THz spectroscopy is today largely used to probe films of inert materials whose thickness is of the order of 1mm [START_REF] Duvillaret | A reliable method for extraction of material parameters in terahertz time-domain spectroscopy[END_REF][START_REF] Duvillaret | Highly precise determination of optical constants and sample thickness in terahertz time-domain spectroscopy[END_REF]. The measurement techniques are similar to those of optics, that is, the thickness and index of the film are derived from the reflection or transmission data.

In the field of agronomy, TDS was used to analyze and monitor the water content status of plants by measuring transmission. Water is a very absorbing material in THz-domain and this property has been used to evaluate water content and the physiologic behavior of plants during water stress situation [START_REF] Nie | Detection of Water Content in Rapeseed Leaves Using Terahertz Spectroscopy[END_REF][START_REF] Borovkova | Terahertz time-domain spectroscopy for non-invasive assessment of water content in biological samples[END_REF][START_REF] Santesteban | Terahertz time domain spectroscopy allows contactless monitoring of grapevine water status[END_REF]. Correlation between water content and transmission level has been established and these results offer a new way to develop a non-invasive and faster tool to measure water plant status.

Experimental set-up

In this study we performed the measurements campaign on a Terapulse 4000 spectrometer with external sample chamber compatible with a specular reflectance and transmission module (TeraView Ltd., Cambridge, UK). THz radiation is produced by ultra-short pulse fibre laser (wavelength 790 nm, pulse duration <100 fs). The pulse is separated into a pump and probe beam, where the pump beam generates THz radiation and the probe beam detects THz pulses using laser-gated photoconductive emitters and detectors. The incident THz pulse was sampled with a spacing of the data points in the time-domain Δ of 0.02 ps over an optical delay extend of 1000 ps, enabling a high spectral resolution of 1.0 GHz. Spectra were the average of 100 scans (approximately 10-min measuring time). Parabolic mirrors focus the THz radiation onto the sample and the set-up can be used in transmission mode or in reflection mode. The sample compartment can be nitrogen-purged to make sure that no residual water vapour peaks are visible in the spectra.

Notice that, as a non-destructive and non-invasive technique, THz TDS does not require sample preparation. Thus the samples under investigation were directly placed on the sample holder which is an aluminum disk with a 2-cm-diameter-hole at the center. The interested region of the leaf, placed in the middle of the disk, is maintained steady in an average plane. Concerning the leaves, the upper epidermis was facing the incoming THz pulse. More details about THz technology can be found in numerous papers such as [START_REF] Jepsen | Terahertz spectroscopy and imaging -Modern techniques and applications[END_REF].

Calibration

Let us denote by ( ) the THz pulse incident on the sample, and by ( ) the signal delivered by the leaf under illumination. Assuming a linear interaction, the input/output relationship can be written as:

( ) = ( ) * ℎ( ) (1)
Where ( * ) designates a convolution product over time , and ℎ( ) the impulse response (or transfer function) of the sample under study. This impulse response is characteristic of a linear filter and takes into account of the inertia of the leaf tissues. It carries the information which is leaf-related. Note that inertia is absorption-related and takes account of the dispersion of the complex indices of the leaf. Reflection and transmission measurements require different set-up geometries. For that reason it is necessary to measure ( ) in both configurations (reflection and transmission modes) and we will note the two incident pulses as ( ) and ( ) before leaf interaction. We precise here that reference in reflection is taken with a sample of a gold sheet which is equivalent to a mirror in the optical range. Hence, according to the mode under study, the leaf signal is written as:

( ) = ( ) * ℎ ( ) (2) ( ) = ( ) * ℎ ( ) (3)
where ℎ ( ) and ℎ ( ) are considered for reflection and transmission data. Actually they are connected with the inverse Fourier transforms of the amplitude reflection ( ) and transmission ( ) coefficients currently measured in the Fourier space in optics. The two measured (input) signals ( ) and ( ) are plotted in Fig. 1(a) and 1(b) in a 15 ps time interval. These curves are the signals delivered by the detector in response to the incident THz pulse in reflection and transmission modes in the absence of the leaf. The difference in these curves is attributed to the difference in the configuration modes. Their characterization is therefore necessary in the calibration process.

Fourier spectrum of the reference signals

Since RET will be performed in the Fourier space, it is interesting at this step to analyse the Fourier spectrum of the input signals. These spectra are given in Fig. 2 at frequencies below 5 THz. We observe the presence of noise at high frequencies ( > 3 THz), for which reason RET will be limited to the frequency window ( < 3 THz) of meaningful data.

Etalon samples

Before proceeding to the leaf analysis, it is prudent to test the procedure on etalon samples. These 3. In these plots, we only consider the frequency range < 1.5 THz due the weak signal-to-noise ratio beyond this limit. Notice here that the Fourier spectra are those of the samples, that is, the transfer functions given by the ouput-to-input ratios of the signals:

ĥ( ) = ˆ ( ) ˆ ( ) (4)

Menzel

The input signal is also plotted in dotted line for each time-domain curve. This allows to emphasize the different echoes by reflection and transmission, that are related to the thickness and index dispersion of the samples. The higher the thickness, the later the echoes. Depending on the sample parameters, such echoes are easy to observe or not. The sample parameters are also responsible for the quasi-periodic oscillations in the spectrum of the signal. The higher the thickness, the smaller the period. The oscillations should contain absorption information due to finger-prints of the chemical substances. However, these finger-prints are negligable in our case for the frequency range < 1.5 THz. This oscillation vanishes at high frequencies for the first two samples, due to absorption. Note in the Fourier space that we plotted the square modulus of the transfer function, that is: ĥ ( ) 2 and ĥ ( )

2 (5)
Also, the absorption is taken as:

= 1 -( ĥ ( ) 2 + ĥ ( ) 2 ) (6)
We notice that the high-resistivity silicon wafer sample has a negligible absorption in the whole frequency range f < 1.5 THz. This can be seen by the repeated zero values taken by reflection at specific frequencies. Indeed at these frequencies the sample is a multiple half-wave layer between 2 identical surrounding media, so that reflection cancels in the absence of absorption (the transparent half-wave layer is known to be absentee). Hence, it would be useless to measure transmission since we can deduce it by T = 1-R. We also notice the presence of noise in the vicinity of the zero frequency, whose origins are not clearly identified at this step (the values approach the mean value of the time signal and this could require a lower time sampling).

Corning Wafer

Transfer function, reflection and transmission coefficients

Here we connect the transfer function to the amplitude reflection and transmission coefficients of the sample. The input pulsed signal is collimated at quasi-normal illumination on the sample. Hence the THz pulse ( ) can be written as a collimated frequency wave packet:

( ) = ∫ ( ) 0 ( ) ( -0 ) - (7)
with 0 ( 0, or 0, ) is related to a phase origin, and where the width Δ of the amplitude distribution ( ) is the inverse of the pulse duration. The parameter follows 0 ( ) = 0 / = 2 0 / , with the speed of light in vacuum, the illumination wavelength in vacuum and 0 the refractive index of the incident medium. Most often, except in vegetal domain where water is majority, this index dispersion is negligible, which leads to:

( ) = ˇ [ -( ) ( -0 )] = ˇ ( ) * [ -( ) ( -0 )] (8)
with ˇ ( ) the inverse Fourier transform of ( ) and the Dirac distribution.

The transmitted field is given by a formula similar to (8) but in which the amplitude distribution is weighted at each frequency by the transmission coefficient ( ) of the leaf, that is:

( ) = ∫ ( ) ( ) 0 ( ) ( -0, ) - (9)
This yields:

( ) = ˇ ( ) * ˇ ( ) * [ -( -0, )] (10)
with ˇ ( ) the inverse Fourier transform of the transmission coefficient ( ). The result is analogous for the reflected field, that is:

( ) = ˇ ( ) * ˇ ( ) * [ -( -0, )] (11)
Relations [START_REF] Ghabbach | Depolarization and enpolarization DOP histograms measured for surface and bulk speckle patterns[END_REF] and [START_REF] Amra | Electromagnetic Optics of Thin-Film Coatings: Light Scattering, Giant Field Enhancement, and Planar Microcavities[END_REF] show that the transfer functions are identified to the inverse Fourier transforms of the reflection and transmission coefficients:

ℎ ( ) = ȟ -( 0 -) ⇔ ℎ ( ) = ( ) [- ( 0 -)] (12)
ℎ ( ) = ȟ -( 0 -) ⇔ ℎ ( ) = ( ) [- ( 0 -)] (13)
and finally:

|ℎ ( )| 2 = | ( )| 2 = ( ) (14) |ℎ ( )| 2 = | ( )| 2 = ( ) (15)
These last intensity quantities ( 14) and ( 15) can be easily measured whatever the position of detectors by reflection and transmission, and on which reverse engineering techniques can be applied. Note also that RET can be applied on the amplitude quantities, ( 12) and ( 13), but the difference 0and 0should be perfectly controlled. Indeed they create in the Fourier space an additional residual phase term which can be written as:

= 2 Δ (16)
For this phase term to be negligible, we must have 2 , that is:

2 Δ 2 ⇒ Δ with = 1 (17)
Since the maximum frequency is 1.5 THz, we obtain Δ 200 m. To summarize, this last condition should be satisfied when working with the amplitude coefficients rather than with the intensity coefficients.

Non dispersive analysis of the samples

In the Fourier space multiple reflections at each interface of a single layer classically allow to write the amplitude reflection and transmission coefficients of the layer as:

( ) = 0 + 0 0 1 ( 0 1 ) 2 ( +1) (18)
( ) = 0 1 ( 0 1 ) (2 +1) (19)
with and the Fresnel coefficients of interface ( ) between media (incident) and + 1, while , are the same coefficients when the illumination is reversed. The phase term is given by:

= , with = ( 2 1 -2 ) 1 2 and = 2 (20)
with the illumination wavelength in vacuum, 1 the complex index of the layer and its thickness. The spatial pulsation is given by : = 0 0 with 0 the illumination incidence angle. Note here that the substrate is air ( 2 = 0 = 1 ).

In formulae ( 18) and ( 19) if we assume negligible absorption (real indices ) and index dispersion, the frequency dispersion is only hold in the phase term , so that the inverse time Fourier transform is immediate at normal illumination ( = = / ), that is:

ˇ ( ) = 0 ( ) + 0 0 1 ( 0 1 ) [ -2( + 1) 1 ] (21)
ˇ ( )

= 0 1 ( 0 1 ) [ -(2 + 1) 1 ] (22)
with the Dirac function. These last formulae emphasize the time locations of the echoes at instants ( ) = 2( + 1) 1 / by reflection and ( ) = (2 + 1) 1 / by transmission. Taken into account the fact that the orders of magnitude of the Fresnel transmission ( ) and reflection ( ) coefficients are unity and a few percent, respectively, we also observe that the first two echoes by reflection have the same order of magnitude (they are proportional to 0 and 0 0 1 ), which is not the case of the echoes by transmission (proportional to 0 1 and 0 1 0 1 ). For that reason and when possible (depending on the receiver detectivity), one would measure the first echo by transmission, and the first two echoes by reflection.

A first evaluation of the sample optical thickness can be extracted from these time instants in reflection and transmission by measuring time delay peaks :

, = 2( + 1) 1 Reflection ( 23
) , = (2 + 1) 1 Transmission ( 24
)
For both modes, reflection and transmission, time delays between two peaks and + 1, is equal to:

Δ ( , +1) = 2 1 ( 25
)
Results are given in Table 1, 2 and 3 where optical thicknesses are extracted for each sample. Note that the agreement is high between values extracted by reflection and transmission. Moreover, the agreement is also high for values extracted from more than 2 echoes (case of silicon wafer sample). However, from these data, we only evaluated the optical thickness of each sample. In order to compare the results to the thickness values of the provider, we have to know the value of the refractive index, which is addressed in the next 2 sub-sections.

Dispersion estimation of refractive index

Though neglected until now, a few information can also be emphasized about the index dispersion. Indeed, minima of the reflection spectrum correspond to optical thicknesses related to a half wavelength in the form:

= 2 or = 2 (26)
with the location of the ℎ minimum of the reflection curve ( ) and the thickness of the sample. Actually in the frequency-domain and following the optical thicknesses given in Table 1, starts at 1 with the first minimum of the ( ). From this consideration we can define and evaluate the variation rate of the index Δ / = ( +1 -)/ as :

Δ = 2 1 +1 - 1 + 1 +1 (27)
Finally, after calculus, we obtain:

Δ = +1 + 1 -1 (28)
Results are depicted in Table 4, 5 and 6 and show for the three samples that the dispersion of the real index is negligible (lower than 3%).

Refraction index estimation

In a similar way, an estimation of the refractive index can be done with the amplitude of the maxima. Indeed the maxima occur for optical thicknesses in the form:

= (2 + 1) 4 ( 29
)
Provided that absorption can be neglected at low frequencies, means large wavelengths, the reflection maxima can be written as: That is: Notice that absorption at low frequencies or large wavelengths can be approximated by:

= 1 -2 1 + 2 ⇔ = ( 1 -2 1 + 2 ) 2 (30)
= 1 + 1 - (31)
1 ⇒ ≈| 1 + | 2 2 (32)
with = 2 and is the imaginary part of the permittivity. This relation shows why absorption is lower at low frequencies (under the assumption of negligible dispersion), which legitimates the use of (31) to estimate refraction index. Results are given in Table 7 and show realistic values of the real index and associated thickness. For example, if we consider the silicon wafer refractive index as 3.4, we find a thickness of the wafer sample ≈ 521 m. This value is in a very good agreement with manufacturing data since the thickness is announced to be 525 +/-20 m.

THz organization analysis of the sunflower leaf

RET

Here, we go further with a RET, a well-known procedure currently used in optical thin films. In regard to the previous sections, note also that absorption is now considered, though dispersion is still neglected.

In this section, the data are analyzed in the Fourier space while we could have chosen to do so in time-domain like some studies already did [START_REF] Peretti | THz-TDS Time-Trace Analysis for the Extraction of Material and Metamaterial Parameters[END_REF][START_REF] Krimi | Highly accurate thickness measurement of multi-layered automotive paints using terahertz technology[END_REF]. For each set of parameters, we define a merit function as the distance between the measurements ( , ) and the calculation ( , ) performed with these parameters, that is :

( , , ) = =1 {[ ( ) -( )] 2 + [ ( ) -( )] 2 } (33)
Notice here that the distance is not normalized. and are respectively the thicknesses and the complex indices of the stack layers, while is the number of frequencies.

As it is shown in Fig. 4, each layer is associated with 3 different parameters which are thickness , and the two parameters of the complex refractive index + . With layers, the number of parameters is 3 if index dispersion can be neglected.

Results on etalon samples

The algorithm was first tested on the 3 etalon samples we used in section 3. When compared to those of the preceeding section, we find the following relative differences:

Δ( ) ( ) ≈ 0.6% Silicon Wafer (37) Δ( ) ( ) ≈ 5.0% Microscope slide (38) Δ( ) ( ) ≈ 6.1% Cover slip (39)
These results are rather successfull and the slight differences can be attributed to absorption which was not considered in the previous section. The corresponding signals by reflection in the time and frequency domains are plotted in Figs. 567, and show a high agreement between calculation and measurement since the ending value of the normalized merit function (residual error) are lowest than 10% in each case. Hence these first results validate the inverse engineering procedure.

Menzel

Preliminary information on the sunflower leaf

Conversely to living tissues, investigations on vegetal physiology has just started few decades ago but is nowadays a domain largely studied with numerous different aspects. Many works on this topic have been published [START_REF] Hopkins | Physiologie végétale[END_REF][START_REF] Buckley | Stomatal responses to humidity in isolated epidermes[END_REF]. For example, photosynthesis mechanisms are clearly identified and understood by biological community as well as stoma functions in temperature and hydration regulation in plants [22]. However, other mechanisms like production of Abscisic acid in leaves or adaptive behaviour in front of different kinds of stress are only partially assimilated. The internal structure of plant leaves is classically described in books [START_REF] Hopkins | Physiologie végétale[END_REF] but few information can be found about the typical dimensions of their different layers. Thickness, roughness and hairiness of leaves could be a trait that would reveal the way the varieties protect themselves against stresses.

In opposition with the previous cases of single-layer reference samples, the leaf problem is more complex due to the multilayer structure. We will see that eight layers are required to approach the leaf structure. Given this complexity, one cannot deal anymore with multiple reflections at each interface, so that the time echoes and periodic oscillations (when existing in the time and frequency domains) cannot be attributed to a specific interface. Actually the data must be analyzed with the admittance formalism which takes account of the global structure of the stack. Furthermore the assumption of negligible dispersion cannot be justified a priori, so that the number of parameters may be excessive. For that reason preliminary information on the sunflower leaf is crucial before the inverse algorithm is applied on the THz data.

Procedures for scanning electron microscopy (SEM) are quite hard and delicate. For observations of samples have to be water free and covered by a thin conductive layer before being placed onto the sample holder of the system. However plant leaves are living tissues that contain almost 90 percent of water with various internal structures. This means that a simple drying would destroy the leaf organization and lead to useless measurements. Hence, instead of drying the leaf, it is necessary to fix the cell structure and substitute the water by an organic product such as osmium tetroxide OsO4. These two steps are very delicate and take several hours because of the fragile character of the leaves and their complex structures. Eventually a nanometer-thin gold layer is sprayed over the sample to allow electric interaction so that secondary electrons can be produced and detected [START_REF] Pathan | Sample preparation for scanning electron microscopy of plant surfaces-horses for courses[END_REF]. The resulting sample is then ready for the SEM analysis and the results are given in Fig. 8.

As expected, the leaf structure reveals strong heterogeneities but the structure can still be glimpsed. In fact, the protocol has roughly maintained the organization of the leaf together with part of the cells shape. The high-resolution picture emphasizes a specific elongated shape of the palisade mesophyll cells, and the multilayer organization. According to this SEM investigation, and though the protocol has largely reduced the typical dimensions and destroyed part of the inner structure, the sunflower leaf can roughly be described as a 4 layer stack organization, that is: the two upper and lower epidermis layers and the largest mesophyll ones which are called palisade mesophyll and spongy mesophyll. Both of them contain water and chlorophyll pigment but have not the same function. Palisade mesophyll presents a well-organized elongated cells which collect light, use it to photosynthesis reaction and guide it toward lower layers. Spongy mesophyll is composed by cells with different shapes and sizes completely desorganised and unsticked at each other. This typicall organization in spongy mesophyll permits gazes to move. Actually, thickness of palisade layer is about 2 times thinner than spongy one. Though the conclusions drawn from Fig. 8 cannot serve as a reliable basis for the leaf investigation, they give a first idea of the leaf heterogeneities responsible for the failure of optical techniques (due to complex speckle patterns). Also, they emphasize the symmetry properties within the multilayer structure of the sunflower sample, a symmetry that we use for RET in section 4.5. Eventually at this step we keep in mind that the thicknesses have been reduced by the preparation process and that more layers are expected to be revealed, as announced by the litterature (see Fig. 9 involving 6 layers). Note also that other techniques (environmental SEM or CryoSEM) could be tried to allow less damage on the samples.

Protocol for the sunflower leaf

Different sunflower plants (dicotyledonous plants) similar in size were watered regularly and used for the THz experiments. The plants were grown from seeds in green houses room under log-day conditions in single pots. Plants were enlightened with horticultural lamp during 8 hours per day (Fig. 10). Two leaves were cut from the petiole before being placed in the THz system, due to the difficulty to place the entire plant close the THz device. In order to minimize plant water stress, we wrapped up the petiole with a wet cotton. Then, to be measured, the samples are placed on a metallic sample holder with a hole at the centre. The leaf is delicately stretched and fixed with little magnet to hold the interested zone flat and stable in time. It is necessary to ensure that the measured area does not contain a large rib but only tertiary. These manipulations last about 5 minutes before measurements. The duration of acquisition can varies from 10 to 15 minutes for each acquisition depending on the number of scan we choose.

Reverse engineering on the sunflower leaf

The reflexion THz data on the leaf sample are given in Fig. 11, both in the time and frequency domains.

One key problem is the unicity of the solution, and the time-consuming algorithm. This is why we simplify the problem by considering the leaf structure as a symmetric organization (except mesophyll layers), with an a priori knowledge that was summarized in Fig. 12. Therefore the solutions will be explored in the vicinity of these values. Due to the high number of parameters, we will extract the best solutions (lowest merit functions) for different layer numbers ranging from 1 to 8. We have to face another difficulty in the case of leaves, due to the stationarity of their properties. Indeed when passing from the reflection to the transmission mode, the sample must be removed and replaced and we cannot guarantee that the same surface is illuminated with the required accuracy. In a near future a specific mechanical system will be available to solve this difficulty, but for the moment we are limited to reflection measurements, which reduces the data for RET. However the absence of transmission data can be partially compensated if we work with the amplitude data by reflection (rather than the modulus as done until now) and work with the real and imaginary parts of the reflection spectrum. This is possible because the geometry of the system is not modified between the reference measurement and the reflection measurement, so that the difference | 0 -| is much lower than 200 m, which satisfies relation [START_REF] Jepsen | Terahertz spectroscopy and imaging -Modern techniques and applications[END_REF]. To summarize, in the case of leaves we will use the RET simultaneously on the real and imaginary part of the reflected spectrum. = 2 Figure 13 is given to show that the sunflower leaf cannot be reduced to a double layer structure. Indeed the best solution is far from the measured data, both in the time and frequency domains. Hence we have to deal with the complexity of the leaf structure and go further with more layers. A good approximation is then reached with 8 layers. With this 8-layer solution (see Fig. 14) the two curves (theory and experiment) are quasi-superimposed, and the layer thicknesses of the two mesophyll layers are in agreement with what was expected.

This solution is quite interesting since it allows to recognize the characteristic cross-section of the leaf that was discussed in section 4. However their index value is low and this makes us think of the trichomes of the leaf (its hair) whose characteristic dimension is known to vary in the range (50 -200 m); such trichrome layers are seen here as homogenized layers and this would explain their low index value.

It is important to stress on the fact that the lowest value (18%) of the normalized merit function is reached for layer numbers greater than 6 ( < 6). Odd numbers were eliminated because of the symmetry conditions that we forced. Hence the 8 layers solution gives the minimum layer = 8 number which minimizes the merit function. Higher even numbers (10 or 12) also give a merit function around 18%, but the layer thicknesses are not realistic since they are all quasi-similar (around 60 micrometers including the cuticle layers) and do not provide a credible representation of the dicotyledonous layer stack.

Conclusion

Due to the leaf dimensions (around 500 m), a THz investigation is well adapted to recover the opto-geometrical multilayer organization of a leaf. Futhermore, these large THz wavelengths reduce all scattering phenomena so that the leaf behaves like a homogeneous multilayer. This balance allowed us to apply the well known (optical thin films) RET to analyze the organization of a sunflower leaf. Notice however that the reflected and transmitted leaf signals are here measured in the time-domain regime, while they are mostly measured in the frequency-domain in the optical thin film community. These time-domain signals provide the key advantage of an amplitude RET analysis, rather than an intensity one. They were recorded at a ps scale, and their Fourier spectra analyzed in the frequency range f < 1.5 THz.

In a first step our RET techniques were tested and validated with success on 3 etalon single layer inorganic samples. Then they were extended to the analysis of multilayer leaves. The leaf problem is more complex due to the high number of parameters (number of layers, thicknesses and complex indices for each layer, index dispersion). In order to reduce this complexity, the RET thickness solutions were explored in the vicinity of premilinary results issued from SEM microcopy (leaf cross-section). Absorption was considered but dispersion was neglected, thanks to the absence of resonances.

First results clearly show that the leaf organization cannot be reduced to one or 2 layers. A highly realistic solution is reached with 8 layers and allows to recognize the leaf cross-section, including the presence of two mesophyll layers, two epidermis layers, two cuticle layers and two trichrome layers, all of them showing the expected thickness magnitude orders. Notice that the trichrome layers ("leaf hair") lie at the top and the bottom of the leaf and exhibit the higher thicknesses and lower indices, so that they take account of a homogeneization process.

In a near future all results will be improved taking into account, when using the RET technique, the amplitude of both the reflected and transmitted signals on the leaves. This will first allow to approach a unicity in the solution, and then consider complex index dispersion. A specific motorized mechanical system with high positioning accuracy will be required for this application.

All results emphasize the fact that the multilayer leaf structure can be revealed with high accuracy thanks to THz data, and this may open the door to a complementary classification of leaves. In addition, the same investigations can be applied to different kinds of leaf stress, including the water stress.

Fig. 1 .

 1 Fig. 1. Input signals delivered in the absence of sample, (a) for reflection and (b) transmission modes. Signal in reflection has been acquired thanks to a gold mirror.

Fig. 2 .

 2 Fig. 2. Frequency-domain projection of incident beam (a) in the reflection mode ˆ ( ) and (b) transmission mode ˆ ( ).

  3 reference materials are one Corning microscope slide (provided by SIGMA-ALDRICH ref : CLS294775X50) of thickness in the range 900 to 1100 m, one Menzel Cover Slip (provided by Thermo Fisher Scientific ref : 24X30-1W-YO) of thickness around 150 m, and a high-resistivity silicon wafer double slide polished with a thickness around 525 m (provided by NEYCO). The time-domain signals and their Fourier spectra are given in reflection and transmission modes for these 3 samples in Fig.

Fig. 3 .

 3 Fig. 3. Time signals and Fourier spectra of the 3 etalon samples measured in reflection and transmission modes, (a-d) for Corning Microscope Slide and (e-h) for Menzel Cover Slip sample. For High-resistivity Silicon Wafer, (i) is time-domain signal on reflection mode and (j) its transfer function.

Fig. 4 .

 4 Fig. 4. Modelisation of a -layer stack. is a gathering matrix of the 3 parameters.

  1. Minimization has been done over a range of frequency from 0.1 to 1.5 THz. The parameters (real optical thickness and complex index) which give the lowest merit function for each sample are the following:

Fig. 5 .Fig. 6 .Fig. 7 .

 567 Fig. 5. Numerical fit in frequency-domain on reflection (a) and transmission (b) for < 1.5 THz on the cover slip sample. Reconstruction of the time-domain signal in reflection (c). The measured data are plotted in blue, while the calculation is in red.

Fig. 8 .

 8 Fig. 8. Scanning electron micrograph of the sunflower leaf cross-section showing a 4 layers stack organization.

Fig. 9 .

 9 Fig. 9. Expected natural structure of leaf cross-section.

Fig. 10 .

 10 Fig. 10. Growing chamber of sunflower plants equipped with horticultural lamp, ventilation system and tide tank.

Fig. 11 .

 11 Fig. 11. (a) squared modulus of frequency-domain projection and (b) the reflection signal acquired in time-domain.

Fig. 12 .

 12 Fig. 12. Symmetry assumptions of the layers for the RET algorythm implemented with Matlab.

Fig. 13 .

 13 Fig. 13. (a) numerical fit of the Fourier amplitude signal plotted in the complex plane for a 2-layer stack. (b) is a time-domain reconstruction by inverse Fourier transform.

3 .

 3 Figure 13 is given to show that the sunflower leaf cannot be reduced to a double layer structure. Indeed the best solution is far from the measured data, both in the time and frequency domains. Hence we have to deal with the complexity of the leaf structure and go further with more layers. A good approximation is then reached with 8 layers. With this 8-layer solution (see Fig.14) the two curves (theory and experiment) are quasi-superimposed, and the layer thicknesses of the two mesophyll layers are in agreement with what was expected.This solution is quite interesting since it allows to recognize the characteristic cross-section of the leaf that was discussed in section 4.3. The two mesophyll layers in the middle of the stack exhibit the expected magnitude orders, with thicknesses equal to 4 = 108.1 m and 5 = 176.3 m. Similar results are obtained for the two epidermis layers ( 3 = 15.3 m and 6 = 30 m) and the two cuticle layers ( 2 = 3.5 m , 7 = 1 m). Eventually the first and the last layer of the leaf are unexpected, since they exhibit high thicknesses ( 1 = 76.4 m, 8 = 134.7 m).However their index value is low and this makes us think of the trichomes of the leaf (its hair) whose characteristic dimension is known to vary in the range (50 -200 m); such trichrome layers are seen here as homogenized layers and this would explain their low index value.It is important to stress on the fact that the lowest value (18%) of the normalized merit function is reached for layer numbers greater than 6 ( < 6). Odd numbers were eliminated because of the symmetry conditions that we forced. Hence the 8 layers solution gives the minimum layer

Fig. 14 .

 14 Fig. 14. (a) numerical fit of the Fourier amplitude signal plotted in the complex plane for a 8-layer stack. (b) is a time-domain reconstruction by inverse Fourier transform. Table 9. Estimated parameters, , and of the eight-layer stack model. The color of the layer 4 and 5 are different since they represent the two mesophyll layers which are not considered symmetrical.

Table 1 .

 1 Time delays of the echoes measured by reflection and transmission for Menzel Cover Slip sample with the associated optical thicknesses.

			Menzel
	Time delay (ps)	0	1	Δ	= 1 ( )
	R	-18.51 -16.09 2.42	363
	T	-9.11	-6.71	2.4	360
	Average			2.41	362

Table 2 .

 2 Time delays of the echoes measured by reflection and transmission for Corning Microscope Slide sample with the associated optical thicknesses.

			Corning
	Time delay (ps)	0	1	Δ	= 1 ( )
	R	-18.47 -1.43 17.04	2556
	T	-4.41 12.63 17.04	2556
	Average			17.04	2556

Table 3 .

 3 Time delays of the echoes measured by reflection and transmission for High resistivity silicon Wafer sample with the associated optical thicknesses.

		Wafer
	Time delay (ps)	0	1	2	3
	R	-18.21	-6.53	5.09 17.03
	Time delay (ps)	<Δ >	1 ( )
	R	11.81	1772

Table 4 .

 4 Dispersion of the real indices for the Menzel Cover Slip etalon, measured between successive wavelengths of minimum reflection.

			Menzel
	(	)	0.44 0.87 1.32
	Δ / (%)	1.15 1.14

Table 5 .

 5 Dispersion of the real indices for the Corning Microscope Slide etalon, measured between successive wavelengths of minimum reflection.

			Corning
	(	)	0.06 0.12 0.18 0.24 0.30 0.36
	Δ / (%)	0.00 0.00 0.00 0.00	0	2.44
	(	)	0.41 0.48 0.53 0.60
	Δ / (%)	2.38 1.89 1.85

Table 6 .

 6 Dispersion of the real indices for the High resistivity silicon Wafer etalon, measured between successive wavelengths of minimum reflection.

				Wafer
	(	)	0.08 0.17 0.26 0.34 0.43 0.51
	Δ / (%)	0	1.92 1.96 1.16 1.18 0.83
	(	)	0.60 0.68 0.77 0.85 0.94
	Δ / (%)	0.84 0.65 0.65 0.53

Table 7 .

 7 Refraction index and thickness values extracted at frequencies of maximum reflection.

	Sample	1 max (THz) R (%)		( )
	Menzel	0.21	48	2.34	155	0.36.10 -3
	Corning	0.09	48	2.34	1092	2.56.10 -3
	Wafer	0.13	71	3.42	518	1.77.10 -3

Acknowledgments

This work was conducted within the OPTIPAG Project supported by the grant ANR-16-CE04-0010 from the French National Research Agency and by the Montpellier University of Excellence I-Site MUSE (PRIME@MUSE), the Région Occitanie, and by the "Gepeto Terahertz platform".

Funding information

French National Research Agency (ANR-16-CE04-0010)

Disclosures

The authors declare no conflicts of interest.