Norhane Benkahla Mediatron

Sup'com

Hajer Tounsi
email: hajer.tounsi@supcom.tn

Ye-Qiong Song
email: ye-qiong.song@loria.fr

Mounir Frikha
email: m.frikha@supcom.tn

Enhanced ADR for LoRaWAN networks with mobility

Keywords: IoT, LoRaWAN, ADR, Spreading factor, Positioning, mobility

LoRa is becoming an attractive low cost and low power WAN solution for many real-world IoT applications. LoRa has been designed for static end-devices to individually use the optimal configuration through an adaptive data rate mechanism (ADR), thanks to the possibility to choose a set of LoRa physical layer transmission parameters. However a large class of IoT applications (e.g. connected farm) also includes mobile nodes with specific mobility patterns. For those applications, the current ADR control algorithm may not be efficient when the radio channel attenuation rapidly changes because of the node mobility. This paper contributes to enhance the ADR mechanism by taking into account the position of the mobile devices and their trajectories in order to have a dynamic allocation. The Enhanced-ADR (E-ADR) minimizes the transmission time and energy consumption as well as packet loss for mobile devices. The testbed-based experiments show that E-ADR improves the quality of service (QoS) of the overall networks.

I. INTRODUCTION

The Internet of Things (IoT) community is currently focusing on the design of large-scale network infrastructures targeting the coverage of massive-scale city-wide scenarios using LPWAN technologies such as LoRaWAN. The LoRaWANs architecture is built on a star topology where multiple LoRa End-Devices (ED) are interconnected to one or many Gateways (GW). A message transmitted by an end-node can be received by multiple close and far-away gateways, which in turn forward the collected messages to a Network Server.

Communication between end-devices and gateways is spread out on different frequency channels and data rates. The selection of the data rate is a trade-off between communication range and message transmission time [START_REF] Sornin (semtech | LoRaWAN Specification v1.1[END_REF]. Communications with different data rates do not interfere with each other thanks to the choice of different spreading factors and frequency channels. The LoRa network server can manage the data rate and the transmission power of each device individually to maximize both the battery life and the overall network capacity by means of the Adaptive Data Rate (ADR) scheme which allows to choose an optimal combination (mode) among a range of parameters: Spreading factor (SF), Bandwidth (BW), Code Rate (CR), Transmission power (TP), Data Rate (DR). Accordingly, there is a trade-off between SF and communica-tion range. The higher the SF, the slower the transmission and the longer the communication range [START_REF] Mikhaylov | Analysis of capacity and scalability of the lora low power wide area network technology[END_REF].

The main challenge in this case is to determine the right configuration that allows a reliable communication with a low energy consumption.

However, the basic ADR scheme may not be efficient in the case of mobile nodes. Indeed the current rate adaptation is performed only after the reception of a set of frames (20 frames in default ADR) [START_REF] Sornin (semtech | LoRaWAN Specification v1.1[END_REF]. It does not always allow to take into account the degradation of the signal due to either the node mobility or the presence of moving obstacles. Our experiments revealed two issues: waiting for collecting 20 frames may be too long for quickly adapting to new situations and transmitting 20 frames may be too many, leading some times to increasing the Time on Air (ToA) beyond authorized 1% duty-cycle imposed by the regulation authorities [START_REF] Sornin (semtech | LoRaWAN Specification v1.1[END_REF].

In this paper, we propose E-ADR, an enhanced ADR scheme aiming at self-adapting to the node mobility through quickly reconfiguration of the mode. E-ADR is based on estimating the next position of a mobile device and its predefined trajectory. Thus the network server could define the best configuration mode for the device to reduce power consumption, channel occupancy time and packet loss rate. E-ADR has been implemented using Waspmote SX1272. It has been tested and compared to the basic ADR in a connected farm scenario.

The rest of the paper is organized as follows. In Section 2, we present the LoRa/LoRaWAN technology and the ADR strategy in LoRaWAN networks. Section 3 highlights related works. Section 4 describes E-ADR in detail. Our implementation and our experimentation scenario will be discussed in Section 5 by showing how the proposal is adapted to the mobile and fixed devices. We conclude in Section 6.

II. LORA/LORAWAN A. LoRa/LoRaWAN background

A LoRa network is based on two components, namely LoRa and LoRaWAN. LoRa is a proprietary physical layer developed by Semtech's Corporation. LoRaWAN (medium access control protocol) is described in an open specification developed by the LoRa Alliance [START_REF] Alliance | LoRaWAN What is it? A technical overview of LoRa and LoRaWAN[END_REF]. The LoRa physical layer uses Chirp-Spread Spectrum (CSS) modulation to enable long distance and low power communications. LoRa operates in the sub-GHz ISM band. Each LoRa transmission is characterized by five parameters: Spreading factor, transmission power, code rate, data rate and bandwidth [START_REF] Libelium Comunicaciones Distribuidas | Waspmote LoRa 868 MHz 915 MHz SX1272 networking guide[END_REF]. These settings affect communication range, data rate, robustness to interference or noise and the ability of a receiver to decode the signal. The configuration of spreading factors allows us to adjust the data rate and the accessible distance. To maximize both the battery life of the end-devices and overall network capacity, LoRa can manage the data rate and transmission power for each end-device individually by means of an adaptive data rate (ADR) scheme [START_REF] Alliance | LoRaWAN What is it? A technical overview of LoRa and LoRaWAN[END_REF]. A device interested in adapting its data rate activates the ADR flag in any up-link MAC frame header. When this is enabled, the network will be optimized to use the fastest possible data rate (to reduce ToA, energy consumption and increase the QoS). When the Network Server is unable to control the data rate of a device due to fast changes, the device's application layer should control it. This latter tries to minimize the aggregated ToA used given the network conditions [START_REF] Sornin (semtech | LoRaWAN Specification v1.1[END_REF].

B. Adaptive data rate scheme

The ADR mechanism in LoRa runs in two parts: on the LoRa node and on the network server. The server provides the most complex part to keep the nodes as simple as possible. The objective of the part running on the node is only to decrease the data rate for increasing the radio coverage if the up-link transmission does not reach the gateway (loss of connection). A counter defining the number of frames sent in up-link without being acknowledged is triggered. If this counter reaches a certain limit, then the node increases the SF (decreases the configuration mode) which could increase the probability of reaching the gateway. The part executed on the server makes it possible to change the transmission power and increase the data rate for the up-links. The server collects the SNRs (Signal-to-noise ratio) of the 20 frames received after activation of the ADR flag and estimates on this basis the new parameters for the future transmissions using equations ( 1) and (2) until the next ADR activation.

SN R margin = SN R max -SN R Req -10 (1)
N step = f loor(SN R margin /3) (2)
In [START_REF] Sornin (semtech | LoRaWAN Specification v1.1[END_REF], SN R max is the maximum SNR of the received frames, SN R Req is the corresponding SNR of the last received packet's mode and 10 refers to a margin constant [START_REF] Sornin (semtech | LoRaWAN Specification v1.1[END_REF]. In (2), N step is the adjustment step. In the case it is greater than 0 and the minimum SF (SF min ) is not yet reached, the server increases the data rate. If the SF is equal to SF min , the server decreases the Transmission Power. In the case where N step is less than 0 and the max power is not reached, the server increases the power otherwise it does not make change [START_REF]The Things Network[END_REF].

III. RELATED WORKS

Recent research on LoRa / LoRaWAN was focused on assessing LoRa performance in terms of capacity, life duration and coverage. These studies were carried out using real experiments in [START_REF] Petaejaejaervi | Performance of a low-power wide-area network based on lora technology: Doppler robustness, scalability, and coverage[END_REF], mathematical models in [START_REF] Bankov | Mathematical model of lorawan channel access[END_REF] or simulations in [START_REF] Abdelfadeel | Poster: A fair adaptive data rate algorithm for lorawan[END_REF]. Among works done for the performance evaluation, [START_REF] Georgiou | Low power wide area network analysis: Can lora scale?[END_REF] shows the exponential dependence of ToA on the SF. [START_REF] Bor | Do lora low-power widearea networks scale?[END_REF] confirms also that the throughput is limited either by the collision rate or by the duty-cycle limitation. As a solution, the authors in [START_REF] Bor | Lora transmission parameter selection[END_REF] implement a dynamic transmission scheme, i.e. ADR in LoRaWAN, and densify the infrastructure by adding additional gateways.

In addition, the authors in [START_REF] Bor | Lora transmission parameter selection[END_REF] were interested in the LoRa transmission parameters. In fact, LoRa device can be configured to use different spreading factors, bandwidth settings, coding rates and transmission powers, resulting in overall 6720 possible settings. The authors in [START_REF] Bor | Lora transmission parameter selection[END_REF] have developed a link probing mechanism with different settings to determine the optimal configuration for each device based on transmission energy. However, due to its high complexity, the trade-off to be found for this solution is the reduction of the number of probes to determine the optimal parameters.

In [START_REF] Slabicki | Adaptive configuration of lora networks for dense IoT deployments[END_REF], the authors considered a distribution of spreading factors based on the median of the SNR values received at the gateway instead of the maximal SNR. However, the case of obstacles and mobility is not considered and waiting for collecting 20 frames may be too long for quickly adapting to these cases.

In [START_REF] Cuomo | EXPLoRa: Extending the performance of LoRa by suitable spreading factor allocations[END_REF], the authors proposed two approaches to allocate SF, EXP-SF (Extending the performance using SF) and EXP-AT (Extending the performance using Air time). The goal is to achieve a high overall throughput. The principle of EXP-SF is to divide the nodes into 6 groups and assign each group a different configuration among the predefined SFs (SF7 to SF12). The first group corresponding to the nodes having the highest RSSIs are assigned the SF 7, the second group with the following RSSIs is assigned the SF8, and so on. The procedure is repeated until all subsets are served. The second approach EXP-AT is more dynamic than the first, where the allocation of the SF theoretically equalizes the transmission time of the nodes. Both approaches are interesting and have been integrated in the processing of the server in the case of the ADR scheme is activated. But in the case of EXP-SF, we can have use cases where the attributed SF is not consistent with the limitations provided in LoRaWAN which degrades performance.

In [START_REF] Kim | Adaptive data rate control in low power wide area networks for long range IoT services[END_REF], authors proposed a new ADR algorithm for LoRa networks at the nodes. The principle of this algorithm is to determine the level of congestion. It performs a learning method using the different transmission parameters. This learning method uses a logistic regression algorithm. In the case where congestion is predicted, the scheme adjusts latency instead of reducing throughput. The disadvantage is that the algorithm requires an active acknowledgment for every transmission. However, this mechanism would decrease the delivery ratio as down-link traffic has been demonstrated to have an impact on up-link throughput. In summary the above reviewed works aim at improving the performance by proposing either enhancements or new ADR, however, they have not considered the mobility case and we think that the degree of their reactivity compared to frequent mobility situations is very low.

IV. CONTRIBUTION

Our goal is to find the optimal configuration of parameters that will respond to the device's needs in terms of QoS even with node mobility. The choice of a configuration has a direct impact on the range of the transmission, the ToA and the energy consumption. Thus, we propose a configuration allocation model driven by the network server. Contrary to the basic ADR, the server can either increase or decrease the configuration mode. This model is suitable for both fixed nodes and the mobile ones that will be localized to adjust the configuration they will use. The main idea of our mechanism is to define a new configuration mode according to the predicted position of a device based on its previous ones and the trilateration method. [START_REF] Mardeni | Node Positioning in ZigBee Network Using Trilateration Method Based on the Received Signal Strength Indicator (RSSI)[END_REF]. Once the new position is defined, the server calculates the corresponding RSSI and looks for the best RSSI interval in which it can be located to determine the most suitable new mode.

A. Position estimation

There are several localization techniques. [START_REF] Kulaib | An overview of localization techniques for wireless sensor networks[END_REF] gives an overview of the localization techniques for wireless sensor networks. The authors have shown that the most accurate solution is the use of GPS. However, its disadvantage is the additional cost in terms of energy consumption. [START_REF] Zhu | APT: Accurate outdoor pedestrian tracking with smartphones[END_REF] has presented a GPS-free solution using an accelerometer and a compass. The disadvantage of this solution is that the error is cumulative over time. Other approaches that are low-cost and more energy efficient are based on the RSSI. [START_REF]Sigfox Geolocation[END_REF] used trilateration approach that consists to define the device's position based on the information collected from three gateways. Each gateway defines the distance that separates it from the device based on the RSSI of the received packet. The distance between a device and a gateway represents the radium of the circle having as center the gateway's position. The intersection of the three circles of the gateways represents the adjusted position of the device. [START_REF] Tsuji | Experiential evaluation of outdoor radio source localization using spatial information of array[END_REF] has used angle of arrival (AoA) technique where the location of a node is estimated by 2 base stations equipped with antenna array. The Time Difference Of Arrival (TDOA) technique [START_REF] Fargas | GPS-free geolocation using LoRa in low-power WANs[END_REF] has used three or more gateways with precise time references.

In our work, we choose the trilateration technique for its very low complexity. We assume that the network server knows the positions of the gateways. To deduce the distance separating the node from the gateway from the measured RSSI, we adopted a tendency curve presented in Fig. 1 that we have obtained experimentally. We draw the curve by moving our device from the gateway by a distance of 1 m for each test and measuring the RSSIs of 100 packets sent. An average RSSI is taken with a 95% confidence interval. We note that this calibration phase should be conducted every time we have a new deployment area. How to make it more general is part of our future work.

Avg displacement (X, Y ) = n-1 i=0 (P i+1 -P i ) n -1 (3)
The position of the device sending the (i+1) th packet is defined as P i+1 (X i+1 , Y i+1 ) is equal to (X i +X Avgdisplacement , Y i + Y Avgdisplacement ). For Each received packet, the server calculates the device's position and compares it to the covered area limits to eventually adjust the estimated position. In fact, our algorithm takes in account the change of the direction of the mobile device and adjusts the position with a distance calculated according to the mobility model. For example in the case of a Zigzag model (Fig. 2), the position is adjusted each time Y i exceeds Y min or Y max or X i exceeds Gap (Fig. 2).

Fig. 2: Zigzag mobility model example

According to the n th estimated position, the server will decide if the device should modify its configuration depending on the calculated change rate that will be defined in the next section.

B. Configuration change rate

We assume that the configuration change should be done at each n packets (1 < n < 20). Indeed, the choice of n very small (equal to 1) generates a significant extra cost in the configuration change, while n high ( for example 20) causes the loss of packets sent with inappropriate parameters. However, the choice of n depends on the devices mobility. After the estimation of the n new positions of a device, the server calculates a configuration change rate according to equation 4, where RSSI real (n) is the real RSSI measured for the n th received packet, B inf (n) corresponds to the lower bound of the RSSI limit range of the current mode and RSSI estim (2n) is the estimated RSSI of the 2n th packet.

R change = |RSSI real (n)-RSSIestim(2n)| |B inf (n)-RSSIestim(2n)| (4)
If the configuration change rate is equal to or greater than a fixed threshold α(α = 0.5 in our case), it means that the device is in a new zone and it is necessary to change its configuration mode. Otherwise, the mode does not change.

C. Choice of new configuration mode

To decide to which mode a device configuration should change, we calculate a transition rate R A→B from mode "a" to each other mode "b" according to equation ( 5). This rate defines the degree of belonging of the estimated RSSI of the 2 th packet in the interval [B inf , B sup ] of each configuration mode. B inf defines the lower bound of the old mode (a) and B sup is the upper bound of the target mode (b). Finally, the server will choose the mode with the biggest rate in the raw presented in equation [START_REF]The Things Network[END_REF]. The new configuration (mode) will be communicated to the device through the LinkADRReq command message.

R A→B = 1 /2×|Bsup(b)-B inf (a)| |Bsup(b)-RSSIestim(2n)|+|B inf (b)-RSSIestim(2n)| (5) R A→B = max(R a,1 , R a,2 , • • • , R a,9 , R a,10 ) (6)

V. EXPERIMENTATION RESULTS

To evaluate our proposed E-ADR, we have conducted several experiments using the Waspmote-SX1272 devices and gateways [START_REF] Libelium Comunicaciones Distribuidas | Waspmote LoRa 868 MHz 915 MHz SX1272 networking guide[END_REF], [START_REF] Libelium | Waspmote data frame programming guide[END_REF]. We have integrated a field in the LoRa Join Request frame where the device indicates to the server what mobility model will be used and the area to cover. We have tested several use cases (Cleaning robot, Drones for inspections, Monitoring robot, Feeding system and temperature sensor) in a smart farm environment with fixed and mobile devices.

We assume that we have 3 gateways (A, B, C) whose positions are respectively: G A [18.02 , 8.62], G B [-28.31 , 45.14], G C [-24.25 , 115.04]. We consider the case of 5 devices and evaluate the network performance in terms of ToA and energy consumption using equation ( 7) [START_REF] Peura | LoRaWAN optimization for a batery powered sensor network[END_REF], where I tx is the transmit current for the transmission time T tx , I rx is the receive current for the RX window duration T rx , I sleep is the sleep current for the time spent in sleep T sleep and U the voltage. We evaluate also the packet loss in the case of using E-ADR and basic ADR respectively.

E = (T tx × I tx + I w1w + T rx1w × I rx1w + I w2w + I rx2w + T sleep × I sleep ) × U (7)
We assume that the server assigns the devices the mode 1 in the join answer to use for the first packets. Also we choose to adapt the configuration every n = 3 packets, which gives good results according to our context. In our use case, most of the time the device needs to change its mode every n = 3 positions. If we check the allocation every one received packet, it increases the transmission delay and the energy consumption. If we increase the change frequency for a number higher than n = 3, we risk more packet losses. In the following, we present and analyze the experimental results for several uses cases.

A. Cleaning robot use case

The farmer needs to have the animal shelters cleaned and checked using a cleaning robot. we take the case of a robot that traverses an area of dimensions ([Y min , Y max ] = [24.6, 188.1],[X min , X max ] = [-20, 20]) and using the zigzag mobility model in Fig. 2. The robot sends 24 messages of 35 bytes each indicating the animal shelter state. Fig. 3 shows the comparison between the basic ADR and the E-ADR in terms of time on air and energy consumption.

Comparing the basic ADR (red line) to the proposed E-ADR (blue line), we can see the minimization of ToA in Fig. 3a for the cleaning robot. We have minimized 80.26% of the transmission time of the basic ADR. In the same way, the gain on energy consumption is increased about 60.23% (Fig. 3b). This is due to the reduced waiting packets (3 instead of 20) in E-ADR and its quicker configuration re-adaptation capacity.

B. Drone inspecting the parcels

Drones are set up to fly over the parcels using the "Square mobility" model in Fig. 4, covering ([Y min , Y max ] = [24.6, 188.1],[X min , X max ] = [-20, 20]) and sending 55 messages of 143 Bytes containing information about the deteriorated areas, the estimation of damage and the amount of missing product on these parcels. The positioning algorithm will help the farmer to know the positions of the drone in order to detect the damage's zones. Fig. 5 compares the basic ADR to E-ADR in terms of time on air and energy consumption.

Fig. 5 shows that E-ADR decreases the time on air about 81.04% and the energy consumption about 67.82% comparing to the basic ADR. Thus, it saves the time allowed by the duty cycle limitation for the device and the node will have more chance to send more data. Hence, we remark that E-ADR decreases the packet loss. In fact, our proposal avoids exceeding the duty cycle restriction (36s/cycle). In Fig. 5, we see that in the case of basic ADR, the drone can no longer send its data from the 8 th packet as it exhausts its allowed duty cycle. This represents a loss ratio of of 87% (we observed only 7 received packets over 55 sent packets). So, E-ADR in this case optimizes the QoS by choosing the right and the optimal configuration.

C. Vegetable and fruit maturity monitoring robot

We take a case of a robot that travels the parcels using the zigzag mobility model presented in Fig. 2 For this application, the basic ADR allocates a new configuration after each set of 20 received frames when the ADR is set. But this allocation will be done depending on the maximal RSSI received during these 20 frames. We remark that the new parameters allocated after the 20 th frame which correspond to mode 7 is not efficient to send the 34 th , 35 th and the 36 th frames (the robot will move away from the 3 gateways and exceeds the limit range of mode 7). The same problem happens after the new allocation done by the server after the 40 th frame which corresponds to mode 8. We remark that the 43 th , 44 th and the 45 th frames are not received because the new allocated parameters does not correspond to the robot's needs.

E-ADR however allows the robot to minimize its transmission time about 78.21% and its energy consumption about 45.55% comparing to the basic ADR and helps us to avoid the risk of packet loss that is about 13% in the case of basic ADR (39 packets received over 45 packets sent). D. Feeding system and temperature sensor An intelligent feeding system is also installed to send 15 messages of 50 Bytes containing the amount of fed animals and their identifiers. We also assume that a temperature sensor is installed for checking the temperature variation in the farm's house by sending 2 messages of 30 Bytes when the temperature becomes higher or lower than a targeting range. Fig. 7 compares the basic ADR to E-ADR in terms of energy consumption and time on air for the feeding system. We remark that E-ADR minimizes the ToA about 60.01% and the energy consumption about 50.7% (See Fig. 7). VI. CONCLUSION ADR may be used to reconfigure node parameters for adapting to node mobility in LoRa networks, but suffers from low adaptation speed and low performance. E-ADR that we proposed and implemented in Waspmote SX1272 solves those issues. Through experimental tests of different scenarios we have shown that E-ADR outperforms basic ADR in terms of time on air and energy consumption. Moreover, it has reduced or eliminated the packet losses since it minimizes the use of the allowed time limited by the duty cycle. So E-ADR allows LoRa networks to supporting mobility, with known mobility pattern. We are interested in a future work to compare our approach (E-ADR) to other proposals particularly EXP-SF and EXP-AT. In addition, we intend to work on a better accuracy of the device position estimation given its mobility model.

Fig. 1 :

 1 Fig. 1: LoRa range measurement vs. RSSI To estimate the next position of a mobile device, we calculate an average displacement based on the last n measured positions P i = (X i , Y i ) according to equation (3).

Fig. 4 :

 4 Fig. 3: Case of a mobile Cleaning Robot

  , covering ([Y min , Y max ] = [24.6, 188.1],[X min , X max ] = [-30, 30]) and sending 45 messages of 17 Bytes containing information about the maturity of the fruits and vegetables. The farmer saves the coordinates of the areas where the fruits and vegetables are not mature enough, to not send the employees to this zone. Fig. 6 compares the basic ADR to E-ADR in terms of time on air and energy consumption.

Fig. 5 :

 5 Fig. 5: Case of a mobile drone

  (a) Evaluation of Time on Air (ToA) (b) Evaluation of Energy Consumption

Fig. 7 :

 7 Fig. 7: Case of a feeding system

TABLE I

 I

	represents the
	different configuration modes used in Lora Waspmote devices
	[4] [5], where each mode n corresponds to a range of RSSI
	(Received Signal Strength) bounded by [B inf (n), B sup (n)]
	(except mode 10).

TABLE I :

 I Configuration modes[START_REF] Libelium Comunicaciones Distribuidas | Waspmote LoRa 868 MHz 915 MHz SX1272 networking guide[END_REF]

	mode	BW (kHz) CR	SF	RSSI (dB)	DR (Kbps)
	1	125	4/5	12	[-134, -131]	0,293
	2	250	4/5	12	[-131, -129]	0,585
	3	125	4/5	10	[-129, -128]	0,976
	4	500	4/5	12	[-128, -126]	1,718
	5	250	4/5	10	[-126, -125.5]	1,953
	6	500	4/5	11	[-125.5, -123]	2,148
	7	250	4/5	9	[-123, -120]	3,515
	8	500	4/5	9	[-120, -117]	7,031
	9	500	4/5	8	[-117, -114]	12,50
	10	500	4/5	7	[≥ -114]	21,875


Eugenia Mylona
email: eugenia.mylona@inserm.fr

Martin Ebert

Angela Kennedy

David Joseph

James Denham

Allison Steigler

Stéphane Supiot

Oscar Acosta

Renaud De Crevoisier

Stephane Supiot

Rectal and Urethro-Vesical Subregions for Toxicity Prediction After Prostate Cancer Radiation Therapy: Validation of Voxel-Based Models in an Independent Population

published or not. The documents may come   L'archive ouverte pluridisciplinaire

A c c e p t e d m a n u s c r i p t Introduction Gastrointestinal (GI) and genitourinary (GU) toxicities following prostate cancer radiotherapy are the most common radiation-induced complications, with 5-year grade ≥2 rates ranging from 6 to 38% and from 10 to 46%, respectively [START_REF] Wortel | Late Side Effects After Image Guided Intensity Modulated Radiation Therapy Compared to 3D-Conformal Radiation Therapy for Prostate Cancer: Results From 2 Prospective Cohorts[END_REF][START_REF] Jolnerovski | Intensity-modulated radiation therapy from 70Gy to 80Gy in prostate cancer: six-year outcomes and predictors of late toxicity[END_REF][START_REF] Zelefsky | Long term tolerance of high dose three-dimensional conformal radiotherapy in patients with localized prostate carcinoma[END_REF]. Traditionally, dose-volume histogram (DVH) parameters extracted from the whole organ-at-risk are used to analyze the dose-effect relationships.

Recently, novel approaches, commonly referred as voxel-based methods, enabled the analysis of the entire 3D dose distribution at fine spatial scales, providing evidence of variable intra-organ radiosensitivity [START_REF] Palma | A Voxel-Based Approach to Explore Local Dose Differences Associated With Radiation-Induced Lung Damage[END_REF][START_REF] Dréan | Identification of a rectal subregion highly predictive of rectal bleeding in prostate cancer IMRT[END_REF][START_REF] Mylona | Voxel-Based Analysis for Identification of Urethrovesical Subregions Predicting Urinary Toxicity After Prostate Cancer Radiation Therapy[END_REF][START_REF] Mcwilliam | Radiation dose to heart base linked with poorer survival in lung cancer patients[END_REF]. Such approaches rely on deformable image registration to align patient anatomies and map dose distributions to a common space of reference, enabling voxel-by-voxel dosimetric comparisons between patients with/without toxicity [START_REF] Acosta | Voxel-based population analysis for correlating local dose and rectal toxicity in prostate cancer radiotherapy[END_REF].

Two studies, aiming at identifying organ subregions predictive of GI and GU toxicities, applied voxel-based analyses on a dataset of 254 patients with localized prostate cancer, treated with IMRT/IGRT at total dose 78/80Gy, from two prospective trials [START_REF] De Crevoisier | Daily Versus Weekly Prostate Cancer Image Guided Radiation Therapy: Phase 3 Multicenter Randomized Trial[END_REF][START_REF] Catton | Randomized Trial of a Hypofractionated Radiation Regimen for the Treatment of Localized Prostate Cancer[END_REF]. A subregion located in the inferior-anterior hemi-anorectum predictive of late rectal bleeding [START_REF] Dréan | Identification of a rectal subregion highly predictive of rectal bleeding in prostate cancer IMRT[END_REF] and five urethro-vesical subregions predictive of grade ≥1 GU toxicities (6) were identified. However, the predictive value of these subregions was assessed in the same cohort used for model development (herein called "development" cohort). In order to translate these findings into clinical practice, the improvement in toxicity prediction when considering these subregions, compared to the whole organ, needs to be demonstrated in independent populations. In addition, the superiority of the subregions should be confirmed for more clinically relevant toxicity events (grade ≥2 ).

The objective of this study was, thus, to validate in a large independent population (herein called "validation" cohort), the predictive value of these rectal and urethro-vesical subregions for grade ≥2 GI and GU toxicity prediction.

MATERIALS AND METHODS

VALIDATION COHORT CHARACTERISTICS

The validation cohort comprises 450 patients from the prospective XXX(BLIND) trial, treated with 3DCRT at a total dose of 66, 70 or 74 Gy (2Gy/fraction) (11). This trial is registered with A1 (Supplementary material). The median follow-up was 72 months. Rectal and urinary toxicities were assessed using the LENT/SOMA and IPSS systems. Grade ≥2 toxicity rates were 23% for rectal bleeding and ranged from 2 to 16% for urinary toxicities.

The radiation treatment technique and organ dose constraints for the development and validation cohorts are described in Supplementary material (Appendix A). More details can also be found in previous publications [START_REF] Dréan | Identification of a rectal subregion highly predictive of rectal bleeding in prostate cancer IMRT[END_REF][START_REF] Mylona | Voxel-Based Analysis for Identification of Urethrovesical Subregions Predicting Urinary Toxicity After Prostate Cancer Radiation Therapy[END_REF](11)(12)(13)(14).

Figure 1 shows the workflow of the study.

DEFINITION OF RECTAL AND URINARY TOXICITY SUBREGIONS IN THE VALIDATION COHORT

Concerning rectal toxicity, a voxel-wise analysis of the dose distributions in the rectum allowed the identification of a generic, geometrically-defined subregion, predictive of 3-year rectal bleeding (grade ≥1 and grade ≥2) (Figure 1A) [START_REF] Dréan | Identification of a rectal subregion highly predictive of rectal bleeding in prostate cancer IMRT[END_REF]. In the validation cohort, dedicated software was used to automatically segment this subregion in each individual's anatomy (Figure 1B).

Concerning urinary toxicity, a voxel-wise analysis of the dose distribution in the bladder and the urethra allowed the identification of five subregions predictive of acute retention and incontinence and late retention, dysuria and hematuria (Figure A1) [START_REF] Mylona | Voxel-Based Analysis for Identification of Urethrovesical Subregions Predicting Urinary Toxicity After Prostate Cancer Radiation Therapy[END_REF]. In the validation cohort, the prostatic urethra was, first, segmented on planning CT images using a multi-atlas-based segmentation algorithm [START_REF] Acosta | Multi-atlas-based segmentation of prostatic urethra from planning CT imaging to quantify dose distribution in prostate cancer radiotherapy[END_REF]. Then, the five urethro-vesical subregions, defined in a templatepatient, were non-rigidly registered to each individual's anatomy. Registration accuracy was assessed using the Dice score (0.88±0.04 for the prostate and 0.92±0.06 for the bladder).

Figure 2A displays the six subregions identified in the development cohort.

STATISTICAL ANALYSIS

All the analyses were performed using the physical dose. DVHs were computed for the rectum, the bladder, the urethra and the six subregions. Wilcoxon test was performed bin-wise to compare patients with and without grade ≥2 toxicities. For each endpoint studied, patients with relevant baseline symptoms were excluded from the analysis (Supplementary material; Table A1). The discriminatory power of the DVHs was assessed in univariate analysis using logistic regression and the area under the ROC curve (AUC) for acute toxicities while Cox regression and the timedependent AUC (tAUC) (16) were used for late toxicities. Bootstrapping with 1000 replicates was performed to determine AUC, t-AUC, and 95% confidence intervals. Considering the subregions, significant differences were found for late rectal bleeding (V44-V74), acute urinary incontinence (V68-V72), late dysuria (V56-V68) and late retention (V14-V64).

Considering the whole organs, significant differences were found for late rectal bleeding (V63-V74), acute urinary incontinence (V67-V72) and late hematuria (V65-V70).

Table 1B shows the discriminatory power of the DVHs for the bladder, the rectum, the urethra and the six subregions. The DVHs of the whole rectum and the rectal subregion demonstrated equal discriminatory power (V68; tAUC=0.61[95%CI: 0.55-0.66]). The doses to three urethro-vesical subregions were predictive for grade ≥2 acute incontinence (V71; AUC=0.69 [95%CI:0.60-0.80]), late retention (V39; tAUC=0.74[95%CI:0.62-0.87]) and late dysuria (V65; tAUC=0.66[95%CI: 0.59-0.73]). The dose to the whole urethra was predictive only for grade ≥2 acute incontinence (V71; AUC=0.68 [95%CI:0.60-0.79]). The dose to the whole bladder was predictive for acute incontinence (V71; AUC=0.66 [95%CI:0.56-0.75]) and late dysuria (V68; tAUC=0.59 [95%CI:0.50-0.69]). Figure 4 shows the discriminatory ability of the DVHs (bin-wise) for the subregions and the whole organs.

DISCUSSION

Voxel-wise analysis is a novel approach allowing the identification of subregions in organsat-risk potentially predictive of toxicity. External validation of previously developed voxel-based spatial models [START_REF] Dréan | Identification of a rectal subregion highly predictive of rectal bleeding in prostate cancer IMRT[END_REF][START_REF] Mylona | Voxel-Based Analysis for Identification of Urethrovesical Subregions Predicting Urinary Toxicity After Prostate Cancer Radiation Therapy[END_REF] was performed, confirming that symptom-specific local dose-effect patterns exist in the urethra for acute urinary incontinence and in the bladder for late retention and dysuria.

By considering the dose to these subregions, instead of the whole bladder, the ability of the models to discriminate patients with and without grade ≥2 toxicities was improved.

In practice, sparing all three subregions in treatment planning, although promising for reducing multiple urinary side-effects, may not always be feasible. In this context, some structures can be considered as more important than others. Even if subregions have been validated for clinically A c c e p t e d m a n u s c r i p t relevant toxicity endpoints (grade ≥2), late symptoms, namely retention and dysuria, are expected to affect patients' quality of life more severely because of their prolonged duration. Late retention, in particular, has higher prevalence in both populations and the corresponding predictive model exhibits superior discrimination. Thus, sparing the posterior part of the bladder (above the trigone) corresponding to late retention, should be prioritized in treatment planning. Apart from the bladder, current findings also indicate a dose-effect relationship in the entire urethral volume for acute urinary incontinence. Even so, reducing the dose to the urethra requires more advanced treatment techniques, such as urethra-sparing stereotactic body radiation therapy (US-SBRT) (17).

On the other hand, two bladder subregions previously identified as predictive for acute retention and late hematuria, as well as the subregion in the rectum associated with rectal bleeding, were not confirmed in the validation cohort as more predictive than the whole organ. Differences between population settings (Tables A1), including the treatment technique, toxicity rates and geographic/temporal aspects, raise questions about a potential "cohort-effect" when performing population-based analyses (18). Noteworthy differences include: significant differences in bladder volume (Table A1), due to the different bladder filling protocols applied in the two populations; higher toxicity rates in the validation cohort ( see Table 1A). For instance, rates of grade ≥2 rectal bleeding were four times higher than in the development cohort, possibly due to the use of 3DCRT [START_REF] Wortel | Late Side Effects After Image Guided Intensity Modulated Radiation Therapy Compared to 3D-Conformal Radiation Therapy for Prostate Cancer: Results From 2 Prospective Cohorts[END_REF]. Overall, these differences may have implications on the generalizability of some subregions.

A limitation of this work is that the statistical DVH-based predictive models, obtained from the development cohort, could not be directly applied to the validation cohort as a part of the external validation process (19). This is due to the differences in the prescribed doses between the two populations. Moreover, the validation cohort was treated with 3D-CRT with prescribed doses ≤74 Gy, which is not representative of the current clinical practices, where higher and conformal doses are delivered to the tumor. Another weakness arises from the consideration of the planning dose distribution which may differ from the actual delivered dose due to intra-and inter-fraction anatomical variations (20, 21). Finally, clinical and biological parameters were not tested in a multivariate setting.

In conclusion, the present study validates three bladder and urethra subregions as more predictive of urinary toxicity, than the whole bladder. These findings highlight the need to go beyond the whole organ-based dosimetric constraints, towards optimized, patient-tailored suborgan constraints. Clinical trials, however, are needed to demonstrate that these optimized planning Red and green curves represent the average DVHs of the patient with/without toxicity symptoms (grade ≥ 2 toxicity), respectively. Continuous lines represent the subregions and dash lines the whole organ. The shadowed regions indicate dose bins with significant differences between patients with/without toxicity (Wilcoxon test; p ≤0.05). DVHs were computed for grade ≥2 toxicity with the exception of hematuria for which the analysis was performed only for grade ≥1, due to the low number of grade ≥2 events.

A c c e

  p t e d m a n u s c r i p t number XXX(BLIND) in ClinicalTrials.gov. Patient and treatment characteristics are summarized in Table

Figure 3

 3 Figure 3 displays organ and subregion DVHs for patients with and without toxicities.

Figure

  Figure 2B displays the three subregions confirmed as more predictive than the whole organs in the validation cohort.

A c c e

  p t e d m a n u s c r i p t strategies are able to significantly reduce urinary toxicity in the case of prostate cancer radiotherapy. A c c e p t e d m a n u s c r i p t A c c e p t e d m a n u s c r i p t

Figure 1 .

 1 Figure 1. Study workflow(A) Rectal and urethro-vesical subregions predictive of GI and GU toxicities, were identified in previous voxel-based studies[START_REF] Dréan | Identification of a rectal subregion highly predictive of rectal bleeding in prostate cancer IMRT[END_REF][START_REF] Mylona | Voxel-Based Analysis for Identification of Urethrovesical Subregions Predicting Urinary Toxicity After Prostate Cancer Radiation Therapy[END_REF]. (B) In the present study, these subregions were segmented in a validation cohort and their ability to discriminate patients with/without grade ≥2 toxicities was evaluated.

Figure 2 .

 2 Figure 2. Subregions identified as predictive of GI and GU toxicities in the development cohort (A) and subregions confirmed as predictive in the validation cohort (B)Three out of six subregions were confirmed as more predictive than the whole organ. All the subregions in the development cohort (A) were identified for grade ≥1 toxicity endpoints. All the subregions confirmed in the validation cohort (B) correspond to grade ≥2 toxicity endpoints.

Figure 3 .

 3 Figure 3. Dose-volume histograms of the subregions and the whole organs for patients with and without toxicities in the validation cohort

  A c c e p t e d m a n u s c r i p t

Figure 4 .

 4 Figure 4. Bin-wise grade ≥ 2 discriminatory power of the DVHs for the subregions and the whole organs Logistic regression was used to compute the discriminatory power of the DVH for acute toxicity (incontinence) and Cox regression for late toxicity (rectal bleeding, dysuria and hematuria). X-axis represents the dose level (x) corresponding to the volume (Vx) of the DVH. Continuous and dashed lines represent the p-values≤0.05 and p-values>0.05, respectively. Blue and orange lines represent the discriminatory power of the DVH in the subregions and the whole organs (bladder/rectum), respectively.

A c c e

  p t e d m a n u s c r i p t

  16. Chambless LE, Diao G. Estimation of time-dependent area under the ROC curve for longterm risk prediction. Stat. Med. 2006;25:3474-3486. 17. Sato M, Hirose K, Hatayama Y, et al. The Potential of the Urethra-Sparing Intensity Modulated Radiation Therapy for Localized Prostate Cancer for Preventing Genitourinary Toxicity: An Interim Analysis. Int. J. Radiat. Oncol. Biol. Phys. 2018;102:e140. 18. Mylona E, Cicchetti A, Rancati T, et al. Local dose analysis to predict acute and late urinary toxicities after prostate cancer radiotherapy: Assessment of cohort and method effects. Radiother. Oncol. 2020;147:40-49. 19. Moons Altman DG, Reitsma JB, et al. Transparent reporting of a multivariable prediction model for individual prognosis or diagnosis (TRIPOD): Explanation and elaboration. Ann. Intern. Med. 2015;162:W1-W73. 20. Nassef M, Simon A, Cazoulat G, et al. Quantification of dose uncertainties in cumulated dose estimation compared to planned dose in prostate IMRT. Radiother. Oncol. 2016;119:129-136. 21. Fuchs F, Habl G, Devečka M, et al. Interfraction variation and dosimetric changes during image-guided radiation therapy in prostate cancer patients. Radiat. Oncol. J. 2019;37:127.

Research This work was supported by the French Institut National du Cancer, a French government grant (through the CominLabs excellence laboratory and managed by the National Research Agency in the "Investing for the Future" program, under reference ANR-10-LABX-07-01) and the PHC FASIC mobility grant. The PROFIT study was funded by Programme Hospitalier de Recherche Clinique, French Ministry of Health. We also acknowledge funding from the Australian National Health and Medical Research Council (grants 300705, 455521, 1006447, 1077788), the Hunter Medical Research Institute, the Health Research Council (New Zealand), the University of Newcastle, the Calvary Mater Newcastle, Abbott Laboratories and Novartis Pharmaceuticals. We gratefully acknowledge the support of the Sir Charles Gairdner Hospital.

data are stored in an institutional repository and will be shared upon request to the corresponding author.

Tables & Figures

(A) Late toxicity rates and 95% confidence intervals were computed using the Kaplan-Meier method. (B) Logistic and Cox regression were used to compute the discriminatory power of DVHs for acute and late toxicities, respectively. Analysis was performed for grade ≥2 toxicities with the exception of hematuria for which the analysis was performed only for grade ≥1, due to the low number of grade≥2 events.


Morgane Gerardin
email: gerardin.morgane@gmail.com

Nicolas Holzschuch

Alain Ibanez

Bernard Schmitt

Pauline Martinetto

Correlation between micro-structural features and color of nanocrystallized powders of hematite

Keywords: hematite, pigment, reflectance, crystal structure

Pigments are quite complex materials whose appearance involves many optical phenomena. Here, we focused on hematite as it is a traditional pigment, whose origin of coloration has been well discussed in the literature. Pure nanocrystallized -FeFe 2 O 3 hematite powders have been synthesized using different synthesis routes. These powders have been characterized by X-Feray powder diffraction and scanning electronic microscopy. The color of the samples has been studied by visible-FeNIR spectrophotometry. We obtained hematite with both various grain morphologies and noticeably different shades going from orange-Fe red to purple. Colorimetric parameters in CIE L*a*b* color space and diffuse reflectance spectra properties were studied against the structural parameters. For small nanocrystals, hue is increasing with the grain size until a critical diameter of about 80 nm, where the trend is reversed. We believe a combination of multiple physical phenomena occurring at this scale may explain this trend.

Correlation between micro-structural features and color of nanocrystallized powders of hematite

EXPERIMENTS

Synthesis of Hematite powders: As the preparation of hematite has been widely studied in the literature, it exists a lot of different routes to control the nucleation and growth of -FeFe 2 O 3 crystallites in solution under standard or hydrothermal conditions. Our study has been done over thirteen samples, chosen to cover a wide range of morphologies and colors, and to mimic the color variations observed on archaeological rock paintings (Table 1).

Our samples numbered 1 to 6 have been respectively prepared following the methods 1 to 6 described by [START_REF] Schwertmann | Iron oxides in the laboratory: preparation and characterization[END_REF]. The samples 1 2 and 1 3 are obtained following the same route as 1 with adding an aging of one month after synthesis. The authors also propose a route to obtain purple hematite, but it resulted as goethite FeO(OH) which transformed to hematite after being annealed (here 16h at 270°C). It corresponds to the sample 'Annealed Goethite' (AG). This production process is similar to the one of some hematite-Febased pigments used by Palaeolithic artists for rock paintings (Pomiès et al. (1998); [START_REF] Gialanella | Artificial or natural origin of Hematite-based red pigments in archaeological contexts: the case of Riparo Dalmeri (Trento, Italy)[END_REF]). We generated peanuts-Felike hematite particles P 1 and P 2 following the routes proposed by Sugimoto et al. (1993). The sample HT has been prepared under hydrothermal conditions following the route of Sugimoto et al. (1993) for platelet particles. Finally, C is a commercial Iron(III) oxide powder sample purchased from Puratronic®.

Characterization: X-FeRay powder diffraction (XRPD) data were collected on the dried powders with a Bruker Endeavour D8 diffractometer operated at Cu Kα radiation (λ = 1.5404 Å). DIFFRAC.EVA's (Bruker) search/match module, which performs searches on the PDF4+ (2018) reference database, has been used for phase identification. All analyses of the X-Feray powder diffraction patterns were performed with the FP_Suite softwares [START_REF] Rodríguez-Fecarvajal | Recent advances in magnetic structure determination by neutron powder diffraction[END_REF]). Rietveld refinements [START_REF] Rietveld | A profile refinement method for nuclear and magnetic structures[END_REF]) were conducted to obtain the structural and micro-Fe structural parameters of each powder and the mass proportion of hematite and goethite in case of mixture. Morphology (size and shape) of the grains composing the dried samples were studied by scanning electron microscopy (SEM, Zeiss Ultra+ microscope, Néel Institute, Grenoble, France).

The diffuse reflectance of the hematite powders was measured using the spectro-Fegonio radiometer SHADOWS developed by [START_REF] Potin | SHADOWS: a spectro-Fegonio radiometer for bidirectional reflectance studies of dark meteorites and terrestrial analogs: design, calibrations, and performances on challenging surfaces[END_REF]. The measurements were performed on dried powder layer sample sufficiently thick so that the substrate does not contribute to the measured reflectance (>1mm), from 400nm to 1000nm with a step of 1 nm, with a normal incidence and an observation angle fixed at 30°. We used a Spectralon® as a white reference.

RESULTS AND DISCUSSION

XRPD: Hematite is identified in all of our samples (Figure 1). All samples are pure, except samples 5, 6 and HT, which contain a sufficiently small to be negligible amount of goethite (Figure 1 and Table 1). Crystallographic data and micro-Fe structural parameters have been obtained from the Rietveld powder structure refinement analysis: unit cell parameters, atomic positions/occupancies and average crystallite size. Only the latter parameter seems to significantly distinguish the different samples. It is calculated from the width of diffraction lines and corresponds to the size of domains over which diffraction is coherent. If the powder consists of single-Fecrystals, crystallite size and grain size (entities seen under the microscope) are identical. SEM: Except for 3 which is composed of two different size populations (a smaller one 3 s and and a larger one 3 l ), the monodispersity of all the samples is well controlled. Various shapes and sizes going from few nanometers to microns are observed (see Figure 2). HT is of big interest because this hexagonal platelet shape is quite similar to the morphology of some natural hematite-Febased pigment found in Lascaux's painting [START_REF] Chalmin | Les blasons de Lascaux[END_REF]). 1, 1 2 , 2, 5 and 6 all show a diamond shape of different but nanometric sizes. The observed grains also differ from one another from their surface appearance: some of them such as P 1 , P 2 and 1 3 look quite rough. It may also be the case for 4 but the powder is too fine to confirm it. The description of the sample C is quite complex because it can be done at two different scales: small spheres seem to have been annealed so that they melted and assembled to form micrometric grains. The average dimensions of the grains were deduced from SEM images and reduced to only one dimensional parameter per sample by computing the radius of an equivalent sphere. In order to better take into account the shape anisotropy of certain samples (AG and HT), this radius has been computed so that the ratio of the surface over the volume of the grain remains constant. Results are listed in Table 1. Grain size and average crystallite size are in good agreement for all samples: it shows that the powders consist of single-Fecrystals. Larger discrepancies are found for samples AG, C, P 1 and P 2 for which average crystallite size are very smaller than grain size. In this case, grains are therefore agglomerates of small single-Fecrystals, as the SEM image of sample P 2 seems to show (Figure 2f).

Reflectance: As all of our samples are quite pure hematite, they all have the same spectral signature: in the studied wavelength range, hematite has three major absorption bands at around 600nm, 700nm and 850nm whose origins are well described by Pailhé et al. (2008), and a fourth one near 430nm that is too saturated to be of such an interest. Depending on the sample under study, we either note a shift in wavelength, or a change of the intensity, the width and the slope of the bands (Figure 3). To achieve a colorimetric study and visualize the difference between the obtained colors, these spectra are converted into color coordinates in the CIE L*a*b* space.

Correlation between micro-structural features and color of nanocrystallized powders of hematite

The thirteen syntheses of hematite lead to crystals (or agglomerates of crystals) with different morphologies and/or sizes and to as many powders of different colors that are summarized in Table 1. The micro-Festructural parameters, spectra characteristics and color parameters have been studied alongside with each other to find a correlation between them. In every case, it is possible to distinguish different trends for the very small particles (<80 nm in diameter) and the bigger ones. The most striking results are obtained for the study of the hue against the grain sizes (see Figure 4). The hue, defined as the angular component of the color in the polar representation (see Figure 3), is first increasing for small particles, then the trend is reversed after reaching a critical diameter (~80nm). Numerous physical phenomena depending on the grain size that are at stake and becoming predominant over one another at this scale should be the cause of such a trend, absorption and scattering being the more important. Simulations of different scattering models are underway in order to fit the reflectance spectra, using numerous phase functions such as Rayleigh, Mie or Henyey-FeGreenstein. (Bohren and Huffman (1983), [START_REF] Henyey | Diffuse radiation in the galaxy[END_REF]). To faithfully describe the observed color of the powders, a precise measure of the hematite refractive index would be necessary. As far as we know, only one measurement of this index has been done by [START_REF] Querry | Optical constants[END_REF], and not any other ever since.

CONCLUSION

Thirteen nanocrystallized hematite powders showing various morphologies and/or sizes of crystals (or agglomerates of crystals) and a wide range of colors have been synthesized. By comparing the micro-Festructural features, the characteristics of measured reflectance spectra and the color parameters in the CIE L*a*b space, it has been observed that some parameters are highly correlated. The study of the hue versus the grain size of the powders shows that the trend following by these parameters is reversed when reaching a critical diameter of about 80nm. We believe that the influence of different scattering phenomena at stake is changing at this scale causing this change of trend.

Correlation between micro-structural features and color of nanocrystallized powders of hematite

Simulations of reflectance spectra with well described scattering models are needed to correctly understand the color of each sample according to the grain size and shape. An accurate measure of the refractive index would be useful to produce more reliable results. We will also aiming to model the bi-Fedirectional reflectance distribution function (BRDF), that would be of great interest for digital representation.

In the future we would like to extend this study to other inorganic pigments, like ochers, mixtures of hematite and clay.

Figure 1 :

 1 Figure 1: Rietveld refinement results (experimental pattern: red data points, calculated pattern: black full line, difference: blue full line) for samples 5. The vertical ticks indicate the Bragg positions for hematite and goethite respectively.

Figure 2 :

 2 Figure 2: SEM images of the different hematite powders. a) Crystals of samples 1, 12, 2, 5 and 6 show pretty similar morphology: same diamond shape and close sizes. Only 1 is shown here. b) Sample 3 presents two bi-pyramid shaped crystal populations different in size. c) Sample 4 is composed of very small spherical particles. d) Sample 13 consists of crystals with a rough surface e) Sample AG (flat needle shape whose thickness is much smaller than its width and its length). f) Samples P1 and P2 have the same shape and different sizes. They seem to consist of agglomerates of crystals, forming grains. Only P2 is shown here. g) Sample HT (large-sized hexagonal platelet shape). h) Sample C presents sticks composed of melted spheres of much smaller size.

Figure 3 :

 3 Figure 3: Reflectance spectra of the samples, calibrated by a Spectralon®, and position of the corresponding colors in the CIE L*a*b* diagram for L*=30.

Figure 4 :

 4 Figure 4: Hue against the equivalent sphere diameter deduced from SEM experiments. The hue is increasing for very small particles until a critical diameter (~80nm, dot line)) where the trend is reversed.

AIC

  Interim Meeting Natural Colours-FeDigital Colours | 20, 26-Fe28 November 2020 | Avignon, France | www.aic2020.org

2

  Diamond 20.0 ± 3.28 20.1 ± 3.00 (27.2, 26.2, 19.0) 1 2 Diamond 24.0 ± 3.22 25.4 ± 4.90 (25.0, 23.2, 15.9) 1 Diamond 25.0 ± 5.90 24.5 ± 4.

Table 1 :

 1 Features of the hematite crystals (or agglomerates of crystals) and colors of the powders. The samples have been organized by increasing crystallite size). Grain shapes and sizes are determined through SEM images. Average crystallite sizes and mass proportions of goethite are optimized through Rietveld refinements. Colors are calculated from the reflectance spectra in the CIE L*a*b* color space.

AIC Interim Meeting Natural Colours-FeDigital Colours | 20, 26-Fe28 November 2020 | Avignon, France | www.aic2020.org

ACKNOWLEDGEMENTS

The authors would like to thank S. Pairis, S. Giraud (Insitut Néel) for their help in the manipulation of the scanning electron microscope, and the chemical syntheses respectively, as well as C. Félix (Institut Néel) and O. Brissaud (IPAG) for their help in the reflectance measurements. This work was financially supported by the

French National Research Agency in the framework of the Investissements d