Agnès Lewden
email: agnes.lewden@univ-brest.fr

Chelsea Ward

Aude Noiret

Sandra Avril

Lucie Abolivier

Caroline Gérard

Tracey L Hammer

Émilie Raymond

Jean-Patrice Robin

Vincent A Viblanc

Pierre Bize

Antoine Stier
email: antoine.stier@iphc.cnrs.fr

Surface temperatures are influenced by handling stress independently of glucocorticoid levels in wild king penguins

Keywords: Thermal imaging, Stress, Corticosterone, Heterothermy, Bird, Beak, Eye Ambient air temperature: Ta Glucocorticoid hormones HPA Maximum eye region surface temperature

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Measuring physiological stress is of central interest in wild animals to understand how they cope with environmental change [START_REF] Ellenberg | Elevated hormonal stress response and reduced reproductive output in Yellow-eyed penguins exposed to unregulated tourism[END_REF][START_REF] Romero | Physiological stress in ecology: lessons from biomedical research[END_REF]. Traditionally, ecologists have evaluated individual stress by assessing the activation of the Sympathetic-Adrenal-Medullary (SAM) or Hypothalamic-Pituitary-Adrenal (HPA) axes, either at baseline levels or in response to acute disturbances in the environment, by metrics such as increased heart rate [START_REF] Cabanac | Temperature and heart rate as stress indicators of handled common eider[END_REF][START_REF] Viblanc | Coping with social stress: heart rate responses to agonistic interactions in king penguins[END_REF], ventilation rate [START_REF] Carere | Shy and bold great tits (Parus major): body temperature and breath rate in response to handling stress[END_REF] or glucocorticoid hormone (GC) secretion [START_REF] Sapolsky | How Do Glucocorticoids Influence Stress Responses? Integrating Permissive, Suppressive, Stimulatory, and Preparative Actions[END_REF]. Stress exposure is also known to trigger core body and surface temperature changes [START_REF] Oka | Chapter 35 Stress-induced hyperthermia and hypothermia[END_REF]. As the latter can be measured using thermal imaging [START_REF] Mccafferty | Applications of thermal imaging in avian science[END_REF], it led to a recent burgeoning of thermal imaging studies to measure the stress response of captive and wild endotherms [START_REF] Mccafferty | Advances in Thermal Imaging[END_REF]. Both the SAM and HPA can potentially influence changes in body surface temperature, but the importance of these two pathways in mediating the stress response(s) measured through thermal imaging remains largely unknown. This information is however essential for the proper interpretation of the physiological stress response measured by thermal imaging. SAM responds within seconds to a stressor by releasing catecholamines hormones (adrenaline and noradrenaline) into the blood stream, which induces an immediate increase in heart rate (tachycardia) and vasoconstriction of peripheral blood vessels [START_REF] Wingfield | Tempests, Poxes, Predators, and People[END_REF]. This redistribution of peripheral blood flow to essential internal organs generally leads to a response called stress-induced hyperthermia along with a decrease of surface temperature [START_REF] Cabanac | Temperature and heart rate as stress indicators of handled common eider[END_REF][START_REF] Herborn | Skin temperature reveals the intensity of acute stress[END_REF][START_REF] Oka | Chapter 35 Stress-induced hyperthermia and hypothermia[END_REF]. The SAM response is followed by the slower (within a few minutes) activation of HPA that releases GC into the blood stream [START_REF] Wingfield | Tempests, Poxes, Predators, and People[END_REF]. GC are metabolic hormones that helps responding to a stressor by mobilizing energy resources and by triggering an increase in behavioural activity and metabolic rate [START_REF] Wingfield | Tempests, Poxes, Predators, and People[END_REF]. Circulating GC levels are usually considered at baseline (i.e. without an inducing stressor) or at stress-induced levels, which are often interpreted as mirroring chronic stress/metabolic demand and acute stress response respectively. Yet, the actual interpretation of circulating GC levels is likely more complex [START_REF] Macdougall-Shackleton | Glucocorticoids and "Stress" Are Not Synonymous[END_REF]Lemmonier et al. 2022). GC-mediated changes in metabolism may cause an increase in internal heat production that, in turn, may lead to changes in body surface temperature [START_REF] Oka | Chapter 35 Stress-induced hyperthermia and hypothermia[END_REF]. Additionally, previous evidence suggest that glucocorticoids might be necessary for catecholamine-induced thermogenic effects such as shivering, free fatty acid mobilization or vasoconstriction [START_REF] Deavers | The function of glucocorticoids in thermogenesis[END_REF].

Therefore, both the SAM and HPA axes could trigger, independently from each other, or together, a change in body core and/or surface temperatures [START_REF] Ouyang | Infrared thermography is an effective, noninvasive measure of HPA activation[END_REF]. The immediate activation of SAM could lead to a decrease of peripheral temperature independently of CORT release. On the other hand, the short to long-term elevation of baseline CORT might lead to increased metabolic rate, heat production and, ultimately, higher peripheral temperature (to facilitate heat dissipation and maintain homeothermy) independently of SAM. To date, the use of thermal imaging as a non-invasive tool to measure stress in unmanipulated animals, as well as the roles of SAM and HPA axes in mediating changes in surface temperature, remains debated. For instance, data in captive and wild birds show that acute stress exposure is often leading to a decrease in body surface temperature (reviewed in Table 1, but note some discrepancies between species and/or body parts), as predicted if these changes are driven by the SAM axis. Evidence for an impact of chronic stressors on body surface temperature (Herborn et al., 2018;Robertson et al., 2020aRobertson et al., , 2020b;;[START_REF] Winder | Body surface temperature responses to food restriction in wild and captive great tits[END_REF] or for a correlation between circulating GC levels and body surface temperature [START_REF] Giloh | Skin surface temperature of broiler chickens is correlated to body core temperature and is indicative of their thermoregulatory status[END_REF][START_REF] Jerem | Eye region surface temperature reflects both energy reserves and circulating glucocorticoids in a wild bird[END_REF][START_REF] Weimer | Broiler stress responses to light intensity, flooring type, and leg weakness as assessed by heterophil-to-lymphocyte ratios, serum corticosterone, infrared thermography, and latency to lie[END_REF], however, remains mixed. Therefore, the question of whether acute or chronic elevation of GC causally induces changes in surface temperature remains largely untested (but see [START_REF] Ouyang | Infrared thermography is an effective, noninvasive measure of HPA activation[END_REF]. Since both HPA and SAM axes provide strong and very reactive responses to stress, teasing apart their contribution to changes in peripheral temperature is challenging.

In this study, we took up the challenge of distinguishing the role of HPA and SAM on changes in body surface temperatures using an experimental approach in adult king penguins (Aptenodytes patagonicus). The king penguin is a heterothermic bird whose body temperature largely varies between core and peripheral tissues (Lewden et al., 2017a(Lewden et al., , 2017b) ) . Its thick body plumage provides a very efficient insulation, while counter-current vascular heat exchangers in the appendages allow to decrease peripheral temperature loss [START_REF] Thomas | Penguin heat-retention structures evolved in a greenhouse Earth[END_REF]. These adaptations allow birds to conserve heat when foraging in cold Antarctic waters [START_REF] Handrich | Hypothermia in foraging king penguins[END_REF], and to rapidly recover to normothermia when exiting the water (Lewden et al., 2020). Thermal imaging has been frequently used in the last years to investigate acute (see an overview in Table 1) and chronic [START_REF] Herborn | Surface temperature elevated by chronic and intermittent stress[END_REF]Robertson et al., 2020aRobertson et al., , 2020b;;[START_REF] Tabh | Stress-induced changes in body surface temperature are repeatable, but do not differ between urban and rural birds[END_REF][START_REF] Winder | Body surface temperature responses to food restriction in wild and captive great tits[END_REF] stress-response in birds. Because of their dense and watertight plumage, the use of thermal imaging in penguins is restricted to measuring bare parts of the body, such as the eye (periorbital region) and beak (Lewden et al., 2020). These regions are highly vascularised and provide valuable information on body surface temperatures.

Breeding king penguins were subcutaneously implanted with either a corticosterone (CORT) implant or a placebo implant. The use of CORT implants in king penguin (see methods and Fig. 1), as in many other bird species (see [START_REF] Torres-Medina | Corticosterone implants make stress hyporesponsive birds[END_REF], leads to higher baseline GC (or CORT) levels, while inhibiting the acute CORT responses in the days following implantation. Indeed, CORT-implanted birds no longer exhibit high blood CORT release in response to acute stress induced by a handling stress protocol. Hence, the use of CORT implants in conjunction with handling stress allows for the separation of the contributions of the HPA and SAM axes to acute changes in body surface temperatures. Here, we therefore measured the effects of our CORT treatment, together with a handling stress protocol, on surface temperatures in the king penguin. If changes in surface temperatures are mainly driven by activating the HPA axis, we expect to find higher initial body surface temperatures in CORT-implanted individuals (with elevated baseline CORT) prior to capture and handling. However, because CORT implants shut-down the release of CORT in response to handling, we expect to find no change in body surface temperature in response to handling stress in CORT individuals if such changes are mediated by the HPA only. On the contrary, if changes in surface temperatures are mainly driven by the activity of the SAM axis, we expect our CORT treatment to have no additional effect on body surface temperatures (both before or after capture and handling). We also measured egg surface temperature on the same day as the body surface temperatures of implanted parents and the surface temperatures of chicks at day 20 after hatching, before they were fully thermally emancipated. Based on experimental data in Adélie penguin showing that elevated baseline CORT led to lower parental care (Thierry et al., 2013a), we hypothesized that increasing CORT levels may lower the surface temperatures of eggs and/or chicks during incubation and brooding respectively, as parents may pay less attention in keeping the egg/chick covered and warm under their brood pouch.

Material and Methods

Study site and species

The study was conducted in 2018-2019 in a king penguin colony of ca 22,000 breeding birds at 'La Baie du Marin' on Possession Island [START_REF] Barbraud | Population trends of penguins in the French Southern Territories[END_REF], Crozet Archipelago, in the Southern Indian Ocean (46°26' S, 51°52' E). We followed 49 breeding pairs from courtship (early November) until the onset of the Austral winter (April). All adults were identified by a hair dye mark on the breast feathers. During the breeding season, male and female alternate between periods on-land caring for their single-egg or chick and periods foraging at sea [START_REF] Weimerskirch | The breeding biology and population dynamics of king penguins Aptenodytes patagonica on the Crozet Islands[END_REF]. The male takes care of the first incubation shift while the female forages at sea. The female returns ~15 days later to relieve her partner, and the parents continue to alternate shifts throughout the incubation (~53 days) and early chick brooding.

Chicks become thermally emancipated from the parents around one month of age, allowing both parents to go back at sea simultaneously to forage [START_REF] Weimerskirch | The breeding biology and population dynamics of king penguins Aptenodytes patagonica on the Crozet Islands[END_REF].

Use of sub-cutaneous CORT and SHAM implants

We randomly divided the 49 breeding pairs into two groups of either female-treated pairs (N = 23) or male-treated pairs (N = 26). In the female-treated pairs, we randomly implanted females with either a CORT (NG-111, 50mg corticosterone, 90-days release; N = 12 females) or a Placebo (NC-111, same vehicle but no CORT; N = 11 females) implant. Subcutaneous implants were purchased from Innovative Research of America (Sarasota, FL, USA).

Females were implanted during their first incubation shift, 3 days after returning from sea.

Similarly, in the male-treated pairs, we randomly implanted 13 males with a CORT and 13 males with a Placebo implant, 3 days after returning from sea. In males, the first incubation shift is immediately after egg laying, and males have already been fasting for ~15 days during courtship. Therefore, we chose to implant the males during their second incubation shift so that they would have comparable fasting times to the females, and therefore both sexes would be in a comparable physiological state during measurements. Male and female partners in, respectively, the female-treated and male-treated pairs were not implanted. Hereafter, they are referred to as "Partner-CORT" or "Partner-Placebo" regardless of their sex. Partner penguins were only considered in the chick analysis (see below).

On the day of the implantation, adults were captured while incubating their egg and restrained by one experimenter with a hood covering the head to keep them calm. Subcutaneous CORT or placebo implants were inserted in the upper part of the back by a veterinary surgeon (SA) under local anaesthesia (ca. 0.5mg/kg xylocaine and 0.0001mg/kg adrenalin, Aspen Pharma). Immediately after implantation, the small incision was closed using 3 sterile surgical staples, and a prophylactic dose of antibiotics (cephalexin ca. 50mg/kg, Rilexine®, Virbac) was injected to prevent any risk of infection. The wound was checked 3 days later, when taking the thermal images without any sign of infection on individuals, nor implant rejection, and staples were removed 6 days after implantation.

To verify the efficacy of our CORT treatment, blood samples were taken 6 days after implantation in less than 4 min (baseline CORT) and after 30 min of standardized handling (stress-induced, see [START_REF] Stier | Oxidative stress and mitochondrial responses to stress exposure suggest that king penguins are naturally equipped to resist stress[END_REF][START_REF] Viblanc | An integrative appraisal of the hormonal and metabolic changes induced by acute stress using king penguins as a model[END_REF]. Plasma total CORT levels were determined by immunoassay according to guidelines provided by the manufacturer (Corticosterone EIA Kit, Arbor Assay, USA). Intra-plate coefficient of variation based on duplicates was 8.17 ± 0.90%, and inter-plate coefficient of variation based on one repeated sample was 8.08%. Implants were successful in raising circulating baseline CORT levels and in inhibiting stress-induced CORT release (Fig. 1).

Fig. 1: Baseline and stress-induced plasma corticosterone levels of breeding king penguins (N = 20 individuals per group) implanted either with a Placebo, or a CORT subcutaneous implant (50 mg, 90 days release) during incubation. Birds were blood sampled 6 days after implantation. CORT levels were influenced by the interaction between handling stress (baseline: < 4 min vs. stress-induced: 30 min of handling stress) and CORT treatment (LMM on log-transformed plasma CORT: F1,37.2 = 62.4, p < 0.001). Letters indicate significant differences according to post-hoc tests with Tukey adjustments (p < 0.030).

Thermal image collection and analysis

Thermal pictures of implanted adults (Fig. 2A) and their egg (Fig. 2B) were taken 3 days after implantation and thermal pictures of chicks (Fig. 2C) were taken at day 20 after hatching (40 to 60 days after CORT or Placebo implant of the focal parent, i.e. within the 90-days release given by the CORT pellet manufacturer). Thermal pictures were taken using FLIR E8 thermal imaging camera, Lewden et al., 2020). When approaching the adults and/or chicks, we started a stopwatch when the focal individual showed the first signs of alarm behaviour to measure the duration of our disturbance at the time the different pictures were taken. In adults, we is known to increase by ca. 100% above baseline levels in king penguins [START_REF] Viblanc | An integrative appraisal of the hormonal and metabolic changes induced by acute stress using king penguins as a model[END_REF].

Just after capturing the adult, we removed the single-egg from the brood pouch and immediately took a thermal picture at a distance of ca. 0.5 m. We then measured egg length and breadth using callipers and estimated egg volume (cm 3 ) following [START_REF] Narushin | Egg geometry calculation using the measurements of length and breadth[END_REF]. When chicks were 20 days old, we took a first set of pictures (1 to 2 pictures per individual) at a distance of ca. 0.5 m immediately after they were removed from the brood pouch (i.e. before handling: 44 ± 5 sec after the first sign of alarm behavior from the parent) and a second set of pictures (1-2 pictures per individual) was taken after handling, 5.4 ± 0.3 min later. Thermal images of the chicks were only taken in non-rainy days, which reduced the sample size from 38 chicks that survived until day 20, to 16 chicks with thermal images available. Chick sex was unknown and body mass was not significantly influenced by the treatment (CORT: 1.14 ± 0.11 vs. Placebo 1.32 ± 0.09, N = 8/group, t = -1.28, p = 0.22). Chicks were kept at ambient temperature (14.9 ± 0.9°C) between the first and second set of pictures, a temperature at which they are known to be able to maintain their internal body temperature constant [START_REF] Duchamp | Ontogeny of thermoregulatory mechanisms in king penguin chicks (Aptenodytes patagonicus)[END_REF]. In adults and chicks, we defined the first set of measures taken < 1 minute as 'before handling' surface temperatures and the second set of measures taken after > 5 minutes of handling stress as 'after handling' surface temperatures. 2021;[START_REF] Tabh | Infrared thermography as a technique to measure physiological stress in birds: Body region and image angle matter[END_REF]. For each image we set the emissivity at 0.98 and controlled in the analyses for daily variation in air temperature (Ta) and relative humidity using daily climatic measures from a permanent weather station located 2 km inland from the colony (https://rp5.ru/Archives_m%C3%A9t%C3%A9o_sur_la_base_Alfred-Faure). We extracted measures of surface temperatures on maximum eye temperature (Teye) and minimum beak temperature (Tbeak) in adults and chicks. We used maximum Teye as recommended by [START_REF] Jerem | Thermal Imaging to Study Stress Non-invasively in Unrestrained Birds[END_REF][START_REF] Jerem | Eye region surface temperature dynamics during acute stress relate to baseline glucocorticoids independently of environmental conditions[END_REF] and the minimum Tbeak to gain insight on the maximum state of vasoconstriction for this body surface. Finally, we estimated egg surface temperature (Tegg) as the mean of average length and width surface temperature.

Ethical note

All the procedures were approved by the French Ethical Committee (APAFIS#16465-2018080111195526 v4) and the Terres Australes et Antarctiques Françaises (Arrêté TAAF A-2018-118).

Statistical analysis

We ran separate analyses for adults, eggs and chicks. In adults, we investigated variations in Teye and Tbeak using two separated linear mixed models (LMMs) where treatment (CORT vs. Placebo), handling stress (before vs. after), Sex and Ta were specified as fixed effects.

Relative humidity was initially tested but removed from final models since it was never significant. Two-ways interactions were also initially included but removed from final model if p > 0.10 using a backward stepwise procedure. The p-values (just before removal in the backward stepwise procedure) for the focal interaction between treatment and handling stress are reported in tables. Bird identity was included as a random intercept to control for repeated measures. We investigated variation in Tegg by entering parental CORT treatment (CORT vs. Placebo), the sex of the parent, Ta and egg volume as fixed effects. We investigated variation chick Teye and Tbeak by entering parental treatment (CORT vs. Placebo), the category of the brooding parent at the time of measurement (implanted individual vs. non-implanted partner), handling stress (before vs. after), and Ta as fixed effects, and chick identity as a random intercept. Interactions were treated as described above for adults. All statistical analyses were conducted using R Studio (version 3.6.2) and lme4, emmeans and pbkrtest packages. Results are reported as least-square means ± SE.

Results

Adult surface temperatures

B

Treatment (CORT vs. Placebo), either alone or in interaction with handling stress, had no significant effect on eye (Teye) or beak (Tbeak) surface temperature in incubating adults (Table 2, Fig. 3A and3B). Handling stress significantly affected Teye and Tbeak, in opposite directions (Table 2): over an 8-min handling stress, Teye increased by 0.8 ± 0.3°C, on average (Fig. 3A), while Tbeak decreased by -1.5 ± 0.4°C, on average (Fig. 3B). The increase in Teye appeared to be sex-specific (marginally significant interaction Handling * Sex: p = 0.060), with only females showing a significant increase in Teye in response to handling stress (Fig. 3C; males: t = 0.63, p = 0.53; females: t = 3.08, p = 0.003). As expected, Teye and Tbeak were significantly and positively associated with Ta (Table 2).

Egg surface temperature

Neither parental treatment (F1,45 = 1.38, p = 0.25; CORT: 34.5 ± 0.2 vs. Placebo: 34.9 ± 0.3°C; N = 46), sex (F1,45 = 0.47, p = 0.49), ambient temperature (F1,45 = 0.45, p = 0.50), nor egg volume (F1,45 = 0.01, p = 0.93) significantly affected Tegg. 3 for details). Least-square means ± SE from final statistical models (Table 3) are presented. Different letters indicate significant differences. The interaction between parental treatment (CORT vs. Placebo) and the category of the parent brooding the chick (implanted vs. non-implanted partner) at the time of measurement significantly affected chick Teye (F1,11 = 6.37, p = 0.028, Table 3) and chick Tbeak (F1,11 = 5.49, p = 0.039, Table 3). Chicks brooded by the partner of CORT-implanted individuals had higher overall Teye and Tbeak than chicks brooded by CORT parents, or by Placebo and associated partners (Fig 4A and4B). While we found no significant interaction between parental treatment and Handling (Table 3), both chick Teye (Fig 4C ) and chick Tbeak (Fig 4D)

Chick surface temperatures

significantly decreased after handling (by -0.9 ± 0.3°C and -3.3 ± 0.6°C respectively, Table 3).

As expected, chick Teye and Tbeak were significantly positively associated with Ta (Table 3).

Discussion

By experimentally manipulating CORT levels in adult breeding king penguins, we show that glucocorticoids are unlikely to play a major role in determining changes in surface temperatures (both before or in response to handling stress) in this large seabird species.

Indeed, we found no significant difference in baseline surface temperature between CORT and Placebo implanted individuals. However, whatever the parental CORT treatment, our handling stress protocol led to an increase in Teye and a decrease in Tbeak in response to capture and handling. The contrasting results between parental treatment and handling stress suggest that changes in surface temperature are probably driven primarily by the activity of the SAM axis in king penguins (but see [START_REF] Ouyang | Infrared thermography is an effective, noninvasive measure of HPA activation[END_REF], for a potential causal involvement of CORT in house sparrows Passer domesticus). Contrary to our prediction, parental treatment had no significant effect on egg temperature. Yet, we found an unexpected pattern on chick surface temperatures: chicks brooded by the non-implanted partner of CORT individuals had higher surface temperatures (both Teye and Tbeak) than chicks brooded by Placebo (or associated nonimplanted partner) or CORT individuals. Finally, in the thermally non-emancipated chicks, both

Teye and Tbeak decreased with handling, but irrespectively of the parental treatment.

CORT and adult surface temperatures

Previous correlative studies have considered baseline and stress-induced CORT levels to investigate the relationship between GC signalling and changes in body surface temperature in birds (blue tit, Cyanistes caerulescens, [START_REF] Jerem | Eye region surface temperature reflects both energy reserves and circulating glucocorticoids in a wild bird[END_REF]; domestic chicken, Gallus gallus domesticus, [START_REF] Giloh | Skin surface temperature of broiler chickens is correlated to body core temperature and is indicative of their thermoregulatory status[END_REF]; [START_REF] Weimer | Broiler stress responses to light intensity, flooring type, and leg weakness as assessed by heterophil-to-lymphocyte ratios, serum corticosterone, infrared thermography, and latency to lie[END_REF]). To go deeper in the involved mechanisms, we experimentally manipulated CORT in king penguin to help to causally distinguish the contribution of the HPA and SAM axes to changes in body surface temperature. Since we did not detect any significant effect of CORT manipulation on baseline or stress-induced surface temperatures, our results provide little support for a causal involvement of the HPA in mediating stress-related changes in surface temperature.

Our results contrast with a recent study showing that an acute stimulation of the HPA axis with ACTH induced a decrease in Teye within 5 min in captive-held house sparrows [START_REF] Ouyang | Infrared thermography is an effective, noninvasive measure of HPA activation[END_REF]. While we used an experimental approach and the lack of significant effect of baseline CORT on body surface temperatures appears to be robust, we cannot fully exclude that our handling stress protocol for thermal imaging was too short to detect a potential role of stress-induced CORT in influencing changes in surface temperatures, since in 8 minutes of handling, CORT levels are expected to double in king penguin, while the maximum CORT levels are only attained after ca. 80 minutes of handling in our study species [START_REF] Viblanc | An integrative appraisal of the hormonal and metabolic changes induced by acute stress using king penguins as a model[END_REF].

Overall, our results suggest that changes in surface temperature are probably more likely mediated by the SAM rather than the HPA in adult king penguins. Future studies testing for a direct acute stimulation of the HPA axis (using ACTH) or SAM axis (using catecholamines) in the king penguin and other bird species could be useful to gain greater insight on the contribution of the HPA and SAM axes in mediating changes in body surface temperature [START_REF] Ouyang | Infrared thermography is an effective, noninvasive measure of HPA activation[END_REF].

Stress-induced changes in surface temperatures

In breeding adults, we observed a mild (sex-specific) significant increase of Teye and a significant decrease of Tbeak in response to acute capture-handling stress. Differences in temperature change between body regions in response to acute stress are not rare (see Table 1). For instance, Moe et al. ( 2017) measured a decrease of Teye contrasting with an increase of Tfeet in chicken, whereas in captive budgerigars (Melopsittacus undulatus) Teye increased transiently after exposure to a stressor while Tfeet did not change [START_REF] Ikkatai | Eye surface temperature detects stress response in budgerigars (Melopsittacus undulatus)[END_REF]; Table 1). Anatomical and functional reasons can likely explain the difference of stressinduced changes in surface temperature between body regions in penguins, with for instance a potential interest in maintaining or increasing eye blood flow to maintain or enhance visual acuity when acutely threatened by an environmental challenge.

Teye has also been shown to positively correlate with internal (cloacal) temperature in various avian species (budgerigars: [START_REF] Ikkatai | Eye surface temperature detects stress response in budgerigars (Melopsittacus undulatus)[END_REF]; chicken: Cândido et al.

(2020)). Consequently, it is possible that the increase in Teye we observe in adults could reflect an increase in internal body temperature [START_REF] Oka | Chapter 35 Stress-induced hyperthermia and hypothermia[END_REF]. Consistently, in closely related Emperor penguins (Aptenodytes forsteri), acute stress leads to an approximate 1.5°C increase in internal (stomach) temperature [START_REF] Regel | Effect of human disturbance on body temperature and energy expenditure in penguins[END_REF]. Such an increase in Teye in response to handling stress has been reported two times in avian species to the best of our knowledge [START_REF] Ikkatai | Eye surface temperature detects stress response in budgerigars (Melopsittacus undulatus)[END_REF][START_REF] Jakubas | Eye Region Surface Temperature and Corticosterone Response to Acute Stress in a High-Arctic Seabird, the Little Auk[END_REF], while most previous studies have shown either a significant decrease (6 studies) or no significant change (2 studies; Table 1). Those contrasted findings might be explained by at least three factors: 1. variation in body size, since smaller individuals/species are expected to favour internal heat conservation by reducing more markedly peripheral blood flow and surface temperature compared to large ones (thermoprotective hypothesis; Robertson et al., 2020a); 2. differences in the thermal environment, since it has been shown that during acute stress response, heat conservation is favoured below the thermoneutral zone while heat dissipation is favoured above its upper limit (Robertson et al., 2020a); 3. the various delays at which Teye was measured in response to acute stress, since for instance the increase in Teye found by [START_REF] Ikkatai | Eye surface temperature detects stress response in budgerigars (Melopsittacus undulatus)[END_REF] was short-lived (visible 5 min after the stressor, but not later on). Although the interaction between handling stress and sex was only marginally significant (p = 0.060), it appeared that only females exhibited a noticeable increase in Teye in response to handling stress in our study.

There are known sex-differences in stress physiology (Handa and McGivern, 2017), but sex differences in surface temperature changes induced by acute stress have been rarely investigated (but see Robertson et al. (2020a), for a result opposite to ours between sexes).

We previously observed no sex effect in HPA responsiveness between males and females king penguins [START_REF] Viblanc | Mutually honest? Physiological "qualities' signalled by colour ornaments in monomorphic king penguins[END_REF], but as mentioned above SAM is the likely driver of acute changes in surface temperatures. Consequently, this suggests that sexes might differ in the stress-sensitivity of their SAM, although this would need to be confirmed by direct measurements of the SAM activity.

The decrease we observe in Tbeak in response to acute stress in adults appears like the typical peripheral vasoconstriction response previously reported for instance by [START_REF] Tabh | Infrared thermography as a technique to measure physiological stress in birds: Body region and image angle matter[END_REF]. Such peripheral vasoconstriction enables the redistribution of the peripheral blood circulation to internal organs and tissues, such as the brain or muscles, favouring their oxygenation and nutrition to sustain the fight-or-flight response. King penguins breed in a highly dense and aggressive colonial environment (up to 100 aggressive interactions per hour; [START_REF] Côté | Aggressiveness in king penguins in relation to reproductive status and territory location[END_REF]. Aggressive interactions, frequently leading to injuries, are known to result in increased heart rates [START_REF] Viblanc | Coping with social stress: heart rate responses to agonistic interactions in king penguins[END_REF]. A decrease of Tbeak during acute stress response likely reflects peripheral vasoconstriction, that is likely widespread to other peripheral body parts, which according to the haemoprotective hypothesis could help to reduce blood loss in case of injury (Robertson et al., 2020a).

The changes observed in surface temperatures in chicks differ from those found in adults. Indeed, we observed a significant decrease in chick's Teye whereas Teye increased in adult females and stayed stable in adult males. While Tbeak significantly decreased as observed in adults, the decrease was even more pronounced in chicks (-3.3 ± 0.6°C vs. -1.5 ± 0.4°C, in chicks and adults, respectively). This is likely explained by greater peripheral thermal losses in chicks suddenly exposed to an ambient temperature drop (ca. 35°C under the brood pouch of the adults vs. 15°C when taken out for measurement and handling), despite their ability to maintain their internal temperature at this stage and ambient temperature [START_REF] Duchamp | Ontogeny of thermoregulatory mechanisms in king penguin chicks (Aptenodytes patagonicus)[END_REF]. The more pronounced drop in surface temperature we observed in chicks during acute handling could reflect their greater need to conserve heat at a relatively mild Ta compared to adults.

Impact of parental CORT on egg and chick surface temperature

Contrary to our prediction, we observed no significant impact of parental CORT treatment on egg temperature, suggesting little alteration of incubation behaviour by increased CORT levels, contrary to a previous study in Adelie penguin reporting a reduction of 1.3 ± 0.2°C in Tegg (Thierry et al., 2013a). It is possible that incubating two eggs in the colder environment of Adélie penguins is more challenging (and thus more likely to be impacted by high CORT levels) than incubating a single egg under the milder climate experienced by king penguins. Alternatively, it is also possible that king penguins are more resilient to increased CORT levels than Adelie penguins, which is supported by the good reproductive success of king penguins implanted with CORT (Stier et al. unpublished), contrary to what has been observed in Adelie penguins (Thierry et al., 2013b). Measurements of incubation quality using a dummy egg with temperature and rotation sensors (Thierry et al., 2013a) may provide more accurate data on this question, but it is important to note that the Tegg measured in this study (Placebo: 34.9 ± 0.3°C) was close to the Tegg measured in the same penguin colony using internal sensors in dummy eggs (35.7 ± 0.4°C; [START_REF] Groscolas | Refeeding signal in fasting-incubating king penguins: changes in behavior and egg temperature[END_REF]. Our results therefore suggest that further studies may benefit from the use of minimally invasive thermal imaging to measure incubation quality, for instance in the context of parental behaviour and climate change [START_REF] Cook | Parenting in a warming world: thermoregulatory responses to heat stress in an endangered seabird[END_REF].

The measure of chick surface temperatures revealed some surprising patterns. While we expected CORT chicks to have lower surface temperatures due to poorer parental care (less efficient brooding), we found the opposite result, but only when chicks were brooded by the non-implanted partner of the CORT-parent. This result should be considered with caution considering the limited sample size available for chick surface temperatures (N = 16). One possibility is that partners from CORT-implanted individuals somehow perceived the 'stress' levels of their partner (as shown between siblings; [START_REF] Noguera | Family-transmitted stress in a wild bird[END_REF] and somehow compensated for parental care through more efficient brooding. This hypothesis requires further study to determine whether the concept of 'family-transmitted stress' [START_REF] Noguera | Family-transmitted stress in a wild bird[END_REF] applies to the king penguin. In this species, both parents must rely on each other over more than one year to successfully fledge their single chick, which itself is entirely dependent on its parents for its food supply throughout its growth. Hence, the king penguin could be promising system to investigate the 'family-transmitted stress' hypothesis.

Conclusion

Evaluating stress levels and responses of free-living animals in a non-invasive manner using thermal imaging is a burgeoning field of research [START_REF] Jerem | Thermal Imaging to Study Stress Non-invasively in Unrestrained Birds[END_REF][START_REF] Tabh | Stress-induced changes in body surface temperature are repeatable, but do not differ between urban and rural birds[END_REF].

However, little is still known about the underlying pathways of stress physiology that influence the changes in surface temperature. Our experimental study in king penguin points toward a likely preponderant role of the SAM. Additionally, by summarizing the current evidence for acute changes in surface temperatures during stress exposure (Table 1), we highlight the relative complexity and inconsistency between the effects observed by different studies.

Further experimental studies related to SAM and HPA involvement are thus required to clarify what stress component(s) we are measuring through non-invasive thermal imaging.

Declarations

Author contribution

Data availability

The datasets used in this manuscript are available on FigShare through a private link for review purposes (https://figshare.com/s/73b39dff978b0dfd2545), and public access with the following doi (10.6084/m9.figshare.20134181) will be enabled pending acceptance of the manuscript.

Acknowledgments

We are grateful to the French Polar Institute (IPEV) and the Terres Australes et Antarctiques

Françaises for providing financial and logistical support for this study through the polar program #119 (ECONERGY). We wish to thank the Zone Atelier Antarctique et Terres Australes AS was supported by a 'Turku Collegium for Science and Medicine' Fellowship and a Marie Sklodowska-Curie Postdoctoral Fellowships (#894963).

Competing interests

Antoine Stier is an Associate Editor of J Thermal Biol but took no part in the peer review and decision-making processes for this research article.

  heat and food stress are not included). Area corresponds to the body region being measured, response to the direction of the temperature response (=: no significant change, ↘: significant decrease and ↗: significant increase, Δ T°C to the amplitude of the temperature response, stressor the nature and duration of the experimental stressor. Some studies measured temperature responses at different time points, and when the effects of stress exposure change through time, the responses are shown with consecutive symbols (e.g. =, ↗: initial absence of change followed later on by an increase in surface temperature).

  set of pictures (1 to 2 pictures per individual) before handling, on average (mean ± SE) 38 ± 4 sec after the first sign of alarm behaviour at a distance of ca. 2 meters from the focal individuals. Adults were then captured-restrained, wound healing was checked, and we took a second set of pictures (1 to 2 pictures per individual) just after terminating handling stress (mean ± SE handling stress duration was 8.3 ± 0.3 min; a timeframe within which CORT

Fig. 2 :

 2 Fig. 2: Infrared images of (A) breeding adult, (B) incubated egg and (C) 20 days-old chick

Fig. 3 :

 3 Fig. 3: Adult king penguin surface temperature responses to an experimental manipulation of glucocorticoid levels (CORT vs. Placebo subcutaneous implants) and handling stress: (A) Maximum eye temperature response to handling stress, (B) Minimum beak temperature response to handling stress, and (C) Sex-specific response of eye temperature to handling stress. In panels A and B, CORT (blue) vs. Placebo (green) individual responses (raw data) are presented. Since CORT treatment had no significant effect either alone or in interaction with handling stress (see Table2for details), the overall mean ± SE is shown in black. Least-square means ± SE from final statistical models (Table2) are presented. Different letters indicate significant differences in panels A and B, and ** in panel C represent the significant increase of Teye observed only in females.

Fig. 4 :

 4 Fig. 4: King penguin chick surface temperature in response to parental glucocorticoid manipulation and handling: (A) Maximum eye temperature according to parental treatment and category of the parent brooding the chick (i.e. partners are not implanted), (B) Minimum beak temperature according to parental treatment and category of the parent brooding the chick, (C) Maximum eye temperature response to handling, and (D) Minimum eye temperature response to handling. In panels C and D, CORT (blue) vs. Placebo (green) chick's individual responses (raw data) are presented, but only the global mean ± SE are presented since there was no significant interaction between parental CORT treatment and handling (see Table3for details). Least-square means ± SE from final statistical models (Table3) are presented. Different letters indicate significant differences.

  Study design: AS, VAV, PB and AL. Funding acquisition: JPR, PB, VAV, AS. Data collection in the field: AS, SA, TH, LA, CG, and JPR. Data collection in the lab: ER and AS. Data analysis: CW, AL, PB, AS. Writing original draft: AL, PB and AS. Writing review and editing: VAV, JPR, AN, SA, TH, LA, ER and CG.

(

  ZATA) from the CNRS for financial support, and the members of the Alfred Faure field station for their help and support in the field. AL was supported by ISblue project, Interdisciplinary graduate school for the blue planet (ANR-17-EURE-0015) and co-funded by a grant from the French government under the program "Investissements d'Avenir" embedded in France 2030.

Table 1 : Summary table of avian studies testing the impact of acute stress exposure on changes in body surface temperatures measured with thermal imaging

 1

(studies related to

  Table 2 for details), the overall mean ± SE is shown in black. Least-square means ± SE from final statistical models (Table 2) are presented. Different letters indicate significant differences in panels A and B, and ** in panel C represent the significant increase of Teye observed only in females.

				A
	Maximum eye temperature	(T eye °C)	25 30 35 40	a Mean ± SE	CORT	b Placebo
				20
					Before handling		After handling
	Minimum beak temperature	(T beak °C)	10 15 20 25 30 35	Mean ± SE a		CORT	Placebo b
				5
					Before handling		After handling
				C
	Maximum eye temperature	(T eye °C)	28 29 30 31 32	Females Males Handling * Sex: p = 0.060 ns **
				27
					Before handling Before handling	After handling After handling

Table 2 : Summary of linear mixed models

 2

	Random effects:		Variance
	Bird ID	Intercept	5.04
	Residual		3.21
	Fixed effects:	Estimate Std. Error	t	df	p
	Intercept	24.94	1.71	14.54	47.9	<0.001
	Treatment (CORT)	-0.10	0.72	-0.15	45.8	0.89
	Handling (after)	1.45	0.47	3.04	89.3	0.003
	Sex (Male)	0.76	0.78	0.97	63.7	0.34
	Ta	0.29	0.13	2.33	45.0	0.024
	Handling*Sex	-1.19	0.63	-1.90	89.7	0.060
	(Treatment x Handling)					(0.44)
	B. Tbeak (n = 138 observations; N = 49 individuals)
	Random effects:		Variance
	Bird ID	Intercept	15.91
	Residual		6.40
	Fixed effects:	Estimate Std. Error	t	df	p
	Intercept	13.84	2.92	4.74	45.0	<0.001
	Treatment (CORT)	-1.13	1.24	-0.91	44.5	0.37
	Handling (after)	-1.47	0.45	-3.29	84.6	0.001
	Sex (Male)	0.90	1.23	0.73	44.6	0.47
	Ta	0.54	0.22	2.45	43.8	0.018
	(Treatment x Handling)					(0.61)

(LMMs) investigating the determinants of (A) maximum eye temperature (Teye) and (B) minimum beak temperature (Tbeak) in adult king penguins. A. Teye (n = 138 observations; N = 49 individuals)

Table 3 : Summary of linear mixed models (LMMs) investigating the determinants of (A) maximum eye

 3

	A. Teye (n = 63 observations; N = 16 chicks)
	Random effects:		Variance
	Bird ID	Intercept	0.81
	Residual		1.18
	Fixed effects:	Estimate Std. Error	t	df	p
	Intercept	29.47	1.39	21.26	11.1	<0.001
	Treatment (CORT)	-0.88	0.68	-1.29	11.3	0.22
	Parent (implanted)	-2.84	0.82	3.48	11.1	0.005
	Handling (after)	-0.85	0.28	-3.08	47.2	0.003
	Ta	0.21	0.09	2.41	10.9	0.035
	Treatment*Parent	2.78	1.10	2.53	11.0	0.028
	(Treatment x Handling)					(0.79)
	B. Tbeak (n = 63 observations; N = 16 chicks)
	Random effects:		Variance
	Bird ID	Intercept	5.53
	Residual		4.89
	Fixed effects:	Estimate Std. Error	t	df	p
	Intercept	13.55	3.42	3.96	11.1	0.002
	Treatment (CORT)	1.15	1.68	0.69	11.3	0.51
	Parent (implanted)	-6.82	2.02	-3.38	11.1	0.006
	Handling (after)	-3.26	0.56	-5.77	46.9	<0.001
	Ta	0.57	0.21	2.65	11.0	0.023
	Treatment*Parent	6.37	2.72	2.34	11.0	0.039
	(Treatment x Handling)					(0.94)

temperature (Teye) and (B) minimum beak temperature (Tbeak) in king penguin chicks. Only one parent was implanted with CORT or Placebo implant, thus we considered both parental treatment (CORT vs. Placebo), the category of the parent brooding the chick at the time of measurement (implanted parent vs. non-implanted partner) and their interaction in our analyses.


R Moreau

T Calligaro

L Pichon

B Moignard

S Hermon

I Reiche

Development of a metaspectral scanner coupling XRF, UV-Vis-NIR Photoluminescence and Vis-NIR-SWIR Reflectance Imaging Spectroscopies

Keywords: cultural heritage, instrumental development, Luminescence Imaging Spectroscopy (LIS), Reflectance Imaging Spectroscopy (RIS), scanning x-ray fluorescence (MA-XRF)

Scanning X-Ray Fluorescence (MA-XRF), Reflectance Imaging Spectroscopy (RIS) and Luminescence Imaging Spectroscopy (LIS) are scientific imaging methods that are extensively used in the investigation and mapping of painted works of art materials in Cultural Heritage Science. Nevertheless, no instrument combining these three highly complementary methods for a simultaneous acquisition exists. In this paper, we present a scanning system combining MA-XRF, RIS and LIS, generating datacubes sharing the same spatial resolution and spatially aligned. The nature and processing of the obtained data is described, in particular their exploration and investigation through the coupled exploitation of datacubes. The operation of this instrument, which is the first to combine XRF with RIS and LIS mapping in an extended domain (400-2500 nm and 200-1000 nm respectively), is exemplified on an anonymous test painting, studied previously. A focus on a restored area of this painting allows to illustrates the complementarity of the three aforementioned methods, for the visualization and characterization of pigments, varnish and binders.

Introduction

The development of scanning X-Ray Fluorescence (MA-XRF) during the past decades has triggered a revolution in the scientific investigation of paintworks and many instruments were developed by different institutions [START_REF] Alfeld | Optimization of Mobile Scanning Macro-XRF Systems for the in Situ Investigation of Historical Paintings[END_REF][START_REF] Hocquet | Elemental 2D Imaging of Paintings with a Mobile EDXRF System[END_REF][START_REF] Alfeld | A Mobile Instrument for in Situ Scanning Macro-XRF Investigation of Historical Paintings[END_REF][START_REF] Romano | Real-Time Elemental Imaging of Large Dimension Paintings with a Novel Mobile Macro X-Ray Fluorescence (MA-XRF) Scanning Technique[END_REF][START_REF] Alberti | A Fast and Reconfigurable Macro X-Ray Fluorescence Scanner for in-Situ Investigations of Polychrome Surfaces[END_REF][START_REF] Ravaud | Development of a Versatile XRF Scanner for the Elemental Imaging of Paintworks[END_REF][START_REF] Pouyet | Development of a Highly Mobile and Versatile Large MA-XRF Scanner for in Situ Analyses of Painted Work of Arts[END_REF] . It has allowed the imaging of pigment distribution over large areas by means of associated elemental signatures [START_REF] Alfeld | Strategies for Processing Mega-Pixel X-Ray Fluorescence Hyperspectral Data: A Case Study on a Version of Caravaggio's Painting Supper at Emmaus[END_REF] , providing better insights in artist's painting materials, evidencing restoration areas and materials, and often revealing underlying and hidden compositions. Nevertheless, despite the short time spend per pixel for XRF data acquisition, MA-XRF scanning operation remains time consuming, requiring minutes to hours and sometimes days. Based on these observations, we developed a new MA-XRF system from a previous versatile scanning system designed by our group [START_REF] Ravaud | Development of a Versatile XRF Scanner for the Elemental Imaging of Paintworks[END_REF][START_REF] Moreau | Application of a MA-XRF/RIS/PL Scanner to Paintwork Studies[END_REF] to acquire during XRF scanning Reflectance Imaging Spectroscopy (RIS) in the Vis-NIR-SWIR range (400-2500 nm) and Luminescence Imaging Spectroscopy (LIS) in the UV-Vis-NIR range (200-1000 nm). The RIS and LIS spectroscopic imaging techniques are highly complementary to MA-XRF, and owing to their quick and non-destructive character, appears suitable for simultaneous implementation during MA-XRF acquisition [START_REF] Ricciardi | Use of Visible and Infrared Reflectance and Luminescence Imaging Spectroscopy to Study Illuminated Manuscripts: Pigment Identification and Visualization of Underdrawings[END_REF][START_REF] Delaney | Use of Imaging Spectroscopy, Fiber Optic Reflectance Spectroscopy, and X-Ray Fluorescence to Map and Identify Pigments in Illuminated Manuscripts[END_REF][START_REF] Delaney | Visible and Infrared Imaging Spectroscopy of Paintings and Improved Reflectography[END_REF][START_REF] De Viguerie | XRF and Reflectance Hyperspectral Imaging on a 15th Century Illuminated Manuscript: Combining Imaging and Quantitative Analysis to Understand the Artist's Technique[END_REF][START_REF] Radpour | Identification and Mapping of Ancient Pigments in a Roman Egyptian Funerary Portrait by Application of Reflectance and Luminescence Imaging Spectroscopy[END_REF] . The coupled exploitation of the signatures can help understanding the diverse materials generally encountered in cultural heritage paintworks, from pigments to binders, through varnishes and restoration materials, each technique complementing and corroborating each other. For instance, while the XRF elemental analysis generally allows identifying many pigments based on the elemental signature, reflectance in the Vis-NIR provides information on electronic bands [START_REF] Elias | The Colour of Ochres Explained by Their Composition[END_REF][START_REF] Cloutis | Identification of Historic Artists' Pigments Using Spectral Reflectance and X-Ray Diffraction Properties I. Iron Oxide and Oxy-Hydroxide-Rich Pigments[END_REF] , photoluminescence on crystal defects [START_REF] Gonzalez | Revealing the Origin and History of Lead-White Pigments by Their Photoluminescence Properties[END_REF] while SWIR spectral features provides additional molecular information [START_REF] Vagnini | FT-NIR Spectroscopy for Non-Invasive Identification of Natural Polymers and Resins in Easel Paintings[END_REF] (e.g. overtones of hydroxyl or carbonate bands) out of reach of the previous techniques. The datacubes generated with this instrument are naturally aligned in term of position and spatial resolution, facilitating coupled data processing, each analyzed point having its own XRF, RIS and LIS signatures. In addition, the orthoscopic imaging obtained (without optical deformation) also avoids the complex image corrections required to merge images from hyperspectral cameras. The instrument was designed to be carried on the field and to scan curved and non-vertical surfaces, and offers an attractive and cost-effective instrument in the field of heritage applications. The following sections describe the specification of the instruments, data processing steps and its application to a case study where the complementarity of the three presented methods is demonstrated.

Experimental

The present MA-XRF instrument is major evolution of a previous MA-XRF system whose details and operation are described elsewhere [START_REF] Ravaud | Development of a Versatile XRF Scanner for the Elemental Imaging of Paintworks[END_REF][START_REF] Moreau | Application of a MA-XRF/RIS/PL Scanner to Paintwork Studies[END_REF] . Using stepped-motor linear stages in a X/Y/Z stack (LTM80, OWIS GmbH) mounted on a sturdy geared head (Manfrotto 400), the instrument can map curved and non-vertical surfaces. It can be transported on the field, the analytical head weighting 20 kg and the entire instrument with all accessories 60 kg (Fig. 1).

The analytical head (Fig. 2) can scan an area extending up to 300 mm (horizontal) × 300 mm (vertical) × 120 mm (depth), following a raster scanning trajectory at constant speed. It is kept in focus during scanning using the Z movement and a telemetric laser pointer. The XRF section is built upon a microfocus X-ray generator (50 kV, 0.7 mA, Rh tube, iMOXS HELMUT-FISCHER) producing a sub-millimeter X-ray beam. The beam shape is defined either using a changeable collimator ( 0.2-0.8 mm) to produce a parallel beam for analysis at long working distance with broad excitation energy range (Rh L to 50 keV) or a using a polycapillary mini-lens for higher fluxes at shorter distances but narrower energy excitation range (4 keV-Rh K). X-rays are collected using a 70-mm² SDD X-ray detector (X123 FastSDD AMPTEK) with USB connection. XRF spectra are acquired on 512 channels, in order to maximize statistics per channel and to minimize download time, and 1-µs shaping time allowing counts rates of 200 k counts/sec with acceptable pile-ups and short dead-times. During the dwell time of the XRF spectrum acquisition, four triggerable light sources (UV-C, UV-B, red and white lasers) are sequentially illuminating the probed spot and optical spectra in the 200-1000 nm range collected using a UV-Vis-NIR spectrometer. In parallel with these acquisitions, a XRF spectrum is acquired and the probed spot is continuously shined by an IR light source and the resulting reflectance spectrum in the 1000-2500 nm range is acquired using a dedicated SWIR spectrometer.

As shown in Figs. 2 and3, UV illuminations are obtained using two LED sources: UV-C (250 nm, 11 nm FWHM, 1 mW optical power), UV-B (365 nm, 10 nm FWHM, 15 mW optical power). The UV-C LED integrating a ball lens directly illuminates the sample at close distance (15 mm) to compensate for its limited optical power. The UV-B illumination is transported through an UVsolarization resistant fiber and focused using a quartz lens located at 45° in the vertical plane. Both UV illuminations are filtered out using a UV band pass filter to block spurious emission of the UV-LEDs in the 400-700 nm range. The telemetric laser pointer is used to illuminate the analyzed spot in the red (655 nm). A supercontinuum laser (360 -1000 nm, 10 mW optical power, SAMBA-UV, LEUKOS) provides the white illumination for RIS. The UV-C, UV-B, red laser and white emissions are collected using a single quartz lens placed at 45° in the horizontal plane to avoid specular reflection and transported through fiber optics to a CCD spectrometer (200-1000 nm, 1044 channels, 100 µm slit, 3 nm resolution Ocean Insight QEPRO). A stabilized tungsten light source providing IR illumination and filtered out to preserve the low intensity luminescence in the NIR range is sent through a fluoride fiber and optics to the sample surface. The IR emitted signals are collected using a fluoride lens and brought through fluoride fiber optics to a SWIR spectrometer using an InGaAs sensor (1000-2500 nm, 256 channels, 50 µm, 13 nm resolution spectrometer, Avaspec NIR-256 2.5 HSC evo, AVANTES).

The synchronization of the illumination alternation is driven by a homemade electronic sequencer that count the number of acquisitions carried out by the QEPRO spectrometer, ensuring a perfect synchronization between the active illumination and its datacube. As a MS Visual Basic software running on a portable computer through USB connections, i.e. without additional electronics, insure the synchronization of the stage X & Y with the six spectra acquisition and the storage in their respective datacube, the reading of the distance and the control the Z stage to keep all the system in focus.

Such optimized optical components and geometry allows collecting acceptable UV-Vis-NIR and SWIR spectra despite the short collection time (e.g. 20 ms and 80 ms respectively). Due to the low intensity of the signals, PL spectra are better acquired in the dark or attenuated light. The scanner delivers in a single pass six datacubes containing all spectra at each pixel: XRF, photoluminescence in the UV-Vis-NIR under excitation at at 250 nm, 365 nm, and 655 nm and reflectance in the Vis-NIR and SWIR ranges.

Data processing

Data formats

The hyperspectral datacubes (XRF, UV-PL at 255/365 nm, NIR-PL at 655 nm and RIS) are stored in EDF format (ESRF Data Format) during scanning since it was the native format for data visualization and investigation using the PyMCA software package [START_REF] Solé | A Multiplatform Code for the Analysis of Energy-Dispersive X-Ray Fluorescence Spectra[END_REF] . A global consensus of the cultural heritage MA-XRF and hyperspectral communities towards a common data format has not yet emerged, but the structured and "self-describing" format Hierarchical Data Format version 5 (HDF5) [START_REF] Collette | Python and HDF5[END_REF] is arising as a very convenient option. For the instrument presented in this article, conversion from EDF format to HDF5 files is achieved after acquisition. Conversion towards other data formats (ENVI dataset, SPECTRONON dataset, etc…) is also possible.

Preprocessing operation

The usual preprocessing operation related to the acquisition of such dataset are realized through homemade Python routines, and available upon request to the corresponding author. Data stitching of multiple scans and correction for small misalignments of the entire datacube are realized through an interactive GUI, while removal of spectral artefacts such as spurious emission lines in the supercontinuum laser and normalization operation of the RIS datacube are carried out with in-line command. Correction for intrinsic luminescence 21 (Kubelka-Munk) and wavelength to energy Jacobian transformation for the emission spectra [START_REF] Mooney | Get the Basics Right: Jacobian Conversion of Wavelength and Energy Scales for Quantitative Analysis of Emission Spectra[END_REF] can also be integrated into the preprocessing workflow of the datacube, again with home-written Python scripts.

Processing software

The MAXRF dataset is first visualized and explored through the versatile and broadly acknowledged PyMCA software initially designed for the processing of XRF spectra. The Region Of Interest (ROI) tool allows displaying the MAXRF acquired maps, for a quick display of the different element distributions over the map by spectral ROI selection, or to show the mean spectrum of a selection of point in the map. The fast linear fitting of the entire dataset enables the quick extraction and display of the different elemental maps. In addition, a basic display and coupled analysis of the other dataset (LIS and RIS) can be achieved thanks to PyMCA ROI Imaging Tool by loading stacks of data of the same shape. When several stacks of data are loaded, a master/slave interaction allows the user to select a spatial region across one map and this region will be selected on the other stacks of data. The mean spectral signature of the region over all the loaded dataset can be displayed, allowing association and visualization of spectral features from the different spectroscopic methods. The dataset acquired with the scanner presented here are particularly well suited for this operation as each dataset feature the same size, spatial resolution and are aligned, i.e. to each pixel correspond a set of XRF, FORS and PL signatures.

To process the RIS and LIS datasets according to usual trends of spectral unmixing in the field of cultural heritage it is possible to use different software, such as ENVI (Harris) [START_REF]Image Processing & Analysis Software | Geospatial Image Analysis Software | ENVI®[END_REF] , SPECTRONON (RESONON) [START_REF]Hyperspectral Scanning System | Imaging Software | Resonon[END_REF] , or DataHandlerP [START_REF] Datahandlerp | [END_REF] , the first one being proprietary and the latter two open source. General workflows of spectral unmixing involves dimensional reduction of the dataset, endmembers selection and abundance estimation of the selected or found endmembers [START_REF] Delaney | Use of Imaging Spectroscopy, Fiber Optic Reflectance Spectroscopy, and X-Ray Fluorescence to Map and Identify Pigments in Illuminated Manuscripts[END_REF][START_REF] Radpour | Identification and Mapping of Ancient Pigments in a Roman Egyptian Funerary Portrait by Application of Reflectance and Luminescence Imaging Spectroscopy[END_REF][START_REF] Delaney | Macroscale Multimodal Imaging Reveals Ancient Painting Production Technology and the Vogue in Greco-Roman Egypt[END_REF] . The different endmembers and their distributions can then be compared with spectral database to identify materials.

Case study

The scanning operation is exemplified on an anonymous test painting dated from the 18 th century, which has already been intensively studied in the frame of the development of emerging nuclear methods for forensics in Arts [START_REF] Calligaro | Emerging Nuclear Methods for Historical Painting Authentication: AMS-14C Dating, MeV-SIMS and O-PTIR Imaging, Global IBA, Differential-PIXE and Full-Field PIXE Mapping[END_REF] (Fig. 4). The painting is an untitled oil on canvas composition measuring 32 cm x 23 cm, depicting a rural landscape, with a cliff and trees on the left side and few houses around a lake with two small characters on the bottom. Most of the composition is devoted to the depiction of sky at sunset, with various shades of blue, pink, and yellow. Curators from the National Center for Research and Restauration in French Museums (C2RMF) suggested that such a style could correspond to a French school composition of the 17 th -18 th century. The instrument described here was used to scan the painting with a 400-µm resolution and 80 milliseconds of dwell time for a total scanning time of 8 hours, with the X-Ray generator operating at 40 kV, 0.5 mA. In Fig. 5 are presented elemental maps obtained after the MA-XRF scanning operation. The elemental maps of Pb, Hg, Fe, As, and Sb shows that original composition is mainly painted with pigments such as lead white (PbCO3)2•Pb(OH)2, vermilion (HgS), iron oxide (FeO•OH and Fe2O3), orpiment (As2S3) and Naples yellow (Pb2Sb2O7) [START_REF]Artists' Pigments: A Handbook of Their History and Characteristics[END_REF] . A focus is made on a specific area (red insert on Fig. 4), which went under an undated restoration operation, covering what seems to be a church steeple or a tower. The restoration operation appears to have been carried out with a mix of various pigments typical of the 19 th and the 20 th century: cadmium yellow (CdS) and cadmium sulfo-selenide red (CdS + CdSe), titanium white (TiO2), and zinc white (ZnO) pigments [START_REF]Artists' Pigments: A Handbook of Their History and Characteristics[END_REF] . Various elemental maps from the selected area are shown in Fig. 6, Pb-L, Pb-M, Hg-L, Cd-L, Fe, Zn, Ti, Se, all shown in linear scale. The Pb and Hg maps show the abundancy of these elements in the original composition, and Cd, Zn, Ti and Se maps clearly display the composition of the restauration area.

All spectra, in Vis-NIR and SWIR ranges, were calibrated using a Spectralon® SRT-MS-050 reference target. Reflectance signatures from different areas are presented in Fig. 7. Inspection of elemental map tends to indicate that the yellow tint of the sky is made out of a mix of Orpiment [START_REF]Artists' Pigments: A Handbook of Their History and Characteristics[END_REF] due to the presence of As, Naples yellow [START_REF]Artists' Pigments: A Handbook of Their History and Characteristics[END_REF] due to the presence of Sb, although reflectance spectrum of this area does not seems to correspond to such a mix. This spectral signature could correspond to the varnish and painting material combination. The restored area display a very specific signature with two sharp absorption front, with inflexion point at 500 nm and 580 nm respectively, with a straight plateau in-between. This shape can be the result of a mix of yellow pigment (Cadmium yellow, Chromium Yellow) and very little red pigment such as Cadmium red, which would be coherent with the distribution maps of such elements seen before (Fig. 6). A comparison between the reflectance spectrum of the restored area and a simulated spectrum of a Chromium Yellow and Cadmium Red mix, derived from Kubelka-Munk theory [START_REF] Kubelka | An Article on Optics of Paint Layers[END_REF][30][START_REF] Liang | Advances in Multispectral and Hyperspectral Imaging for Archaeology and Art Conservation[END_REF] is shown in Fig. 7. b. The actual reflectance spectrum and the simulated one display very similar shape, with a slight shift of 7 nm, corroborating the possible pigment identification from XRF data only.

Exploration of the UV-365nm LIS dataset reveal different emission characteristic of various materials. ROI spectral distribution of certain signals observed are shown in Fig. 8 with their associated spectral signatures. Large band centered around 500 nm is characteristic of aged organic material emission such as binders and varnishes [START_REF] De La Rie | Fluorescence of Paint and Varnish Layers (Part II)[END_REF][START_REF] Thoury | Nondestructive Varnish Identification by Ultraviolet Fluorescence Spectroscopy[END_REF] , whereas on the restauration area, with presumably more recent varnish or no varnish at all and different binder, this large band is not present. The restored area display a very specific luminescence signature that can be associated to emissions of modern materials, presumably Cadmium and Zinc based pigments [START_REF] Cesaratto | Analysis of Cadmium-Based Pigments with Time-Resolved Photoluminescence[END_REF][START_REF] Bellei | Photoluminescence Properties of Zinc White: An Insight into Its Emission Mechanisms through the Study of Historical Artist Materials[END_REF] . Looking closely at Fig. 8.b and c, the shape of the restauration area observed is different compared to the MA-XRF results. Indeed, the generated maps are display of emission signals from varnishes fluorescence but also from aforementioned pigments combined with varnish absorption, this also being an example of how challenging interpretation of LIS can be.

The RIS-SWIR dataset provide information regarding the nature of binders and varnishes used [START_REF] Vagnini | FT-NIR Spectroscopy for Non-Invasive Identification of Natural Polymers and Resins in Easel Paintings[END_REF][START_REF] Dooley | Mapping of Egg Yolk and Animal Skin Glue Paint Binders in Early Renaissance Paintings Using near Infrared Reflectance Imaging Spectroscopy[END_REF][START_REF] Dooley | Standoff Chemical Imaging Finds Evidence for Jackson Pollock's Selective Use of Alkyd and Oil Binding Media in a Famous 'Drip' Painting[END_REF][START_REF] De Viguerie | Mapping Pigments and Binders in 15th Century Gothic Works of Art Using a Combination of Visible and near Infrared Hyperspectral Imaging[END_REF] . The SWIR datacube was processed using ENVI software. The global image in Mininum Noise Fraction in the hourglass wizard clearly show distinct MNF components in the painting, especially for the restored part, that are barely distinguished in XRF or PL (Fig. 9.a and 9.b). The SWIR spectrum of the original composition and of the restored area are presented in Fig 10 .a. The distribution of specific band are presented in Fig. 10.c. and Fig. 10. e. To obtain these maps, the SWIR dataset was loaded into ROI Imaging Tool widget of PyMCA, selecting spectral Region of Interest, with Background subtraction PyMCA feature, by drawing a straight line between the spectral ROI limits, and using inverted greyscale, so that stronger absorption bands would appear brighter than a smaller absorption bands. One can observe a signal at 1450 nm (Fig. 10.b), attributed to hydroxyl (-OH) stretching from hydrocerussite [START_REF] Delaney | Mapping the Pigment Distribution of Vermeer's Girl with a Pearl Earring[END_REF] , which makes sense as the associated distribution map of this signal (Fig. 10.c) follows the one of Pb distribution (Fig. 6). We can also observe a doublet at 1707 nm and 1735 nm (Fig. 10.d), that is present only for the area of the original composition. As the restored area display a more recent varnish layer or no varnish at all, we can assume that this signal can be attributed to the varnish, which is coherent with results from litterature, where this spectral signature can be found in NIR spectra of different varnishes [START_REF] Vagnini | FT-NIR Spectroscopy for Non-Invasive Identification of Natural Polymers and Resins in Easel Paintings[END_REF][START_REF] Dooley | Standoff Chemical Imaging Finds Evidence for Jackson Pollock's Selective Use of Alkyd and Oil Binding Media in a Famous 'Drip' Painting[END_REF] .

As briefly discussed previously, one can easily see that the instrument presented here generated large datasets from techniques that differ in interpretation, meaning that the exploration of the results requires broad knowledge, not to say expertise, in very different domains. That interpretation is user dependent, meaning that different user could come up with different results based on their experience and knowledge. Moreover, this operation of interpretation can be very time consuming, consideration the very large amount of data generated.

Thus, new approaches can be necessary to explore the dataset and identify the materials probed: data fusion can be considered as the most promising approach. From the work of Llinas et al. [START_REF] Llinas | An Introduction to Multi-Sensor Data Fusion[END_REF] , Ramos et al. [START_REF] Ramos | Micro-Raman and X-Ray Fluorescence Spectroscopy Data Fusion for the Classification of Ochre Pigments[END_REF] identify two major strategies suitable to be used with such datasets: data level fusion, also called low-level fusion, and feature level fusion, also called mid-level fusion. In low-level fusion, the raw data are fused, and the fused dataset are considered as a single metaspectra datacubes. In mid-level fusion, the fusion is realized on the most relevant features of each dataset, introducing a previous step of variable selection. They applied both procedures to XRF and Raman spectral dataset of different ochre pigments in order to classify them, and showed that the fused dataset allowed a more accurate classification of these pigments compared to using the XRF or Raman dataset only. Alfeld et al. [START_REF] Alfeld | Joint Data Treatment for Vis-NIR Reflectance Imaging Spectroscopy and XRF Imaging Acquired in the Theban Necropolis in Egypt by Data Fusion and t-SNE[END_REF] fused MAXRF and RIS datasets from two different instruments, an operation that needed a computational spatial alignment of the datasets, and performed t-SNE clustering on this joined dataset. The authors observed that the t-SNE classification of the fused dataset allowed a better visualization of the results, as the generated clusters highlights the presence of sub-layers surfaces.

As the instrument presented here generates spatially aligned dataset with the same spatial resolution, analytical workflow based on data fusion can be considered as of great interest of development. Indeed, it could allow to visualize and explore the relationship existing between pigments, binder and varnishes for each pixels of the generated maps. Moreover, a relative notion of stratigraphy could be visualized as MA-XRF and RIS-SWIR probe the matter with a relative depth whereas RIS-Vis-NIR and LIS only reach the first opaque surface encountered. Data fusion processing could thus reveal new information regarding the layered organization of the studied area.

Conclusion and perspectives

To the authors knowledge, the instrument described here constitutes the first instrument dedicated to the study of cultural heritage that delivers six datasets (MAXRF, RIS Vis-NIR, RIS SWIR, UV-LIS at 255 nm, UV-LIS at 365 nm and NIR-LIS at 655nm) in a single pass. These datasets are naturally aligned and with the same resolution, constituting a single datacube that can be process as one, or separately. A wide range of materials generally encounter in cultural heritage fields can thereby be probed in a non-invasive and non-destructive way. Application of this instrument was exemplified on an anonymous painting and the exploration of the generated datasets allowed to visualize and characterize the materials of a restored area. Among the perspectives related to the use of this instrument, the development of routines algorithms and automated process workflow for the rapid identification of materials coupling the various signatures obtained after a scanning operation is currently in progress.

Figures

Figure 1 :

 1 Figure 1: Multi-sensor scanner in operation on the field on a mural painting.

Figure 2 :Figure 3 :

 23 Figure 2: Close up of the analytical head showing the implementation of RS and PL system  X-Ray collimator X-Ray detector SWIR illumination SWIR collection UV 250-nm White Laser UV-Vis-NIR collection UV 365-nm Telemetric red laser Camera

Figure 4 :

 4 Figure 4: Test painting, untitled and anonymous composition on canvas. The red dotted squared corresponds to the scanned area.

Figure 5 :

 5 Figure 5: MAXRF elemental maps: a) Pb-L, b) Fe, c) Hg-L, d) As-Kb, e) Sb-L, f) Cr, g) Ti, h) Zn, i) Cd-L, all shown in linear scale. See text for comments.

Figure 6 :

 6 Figure 6: MAXRF elemental maps of the scanned area: a) Visible Image of the area, b) Pb-L, c) Pb-M, d) Hg-L, e) Cd-L, f) Fe, g) Zn, h) Ti, i) Se. See text for comments.

Figure 7 :

 7 Figure 7: a) Reflectance signature Vis-NIR of the original composition (red), and restored area (black) b) Reflectance signature of the restored area and simulated spectrum derived from Kubelka-Munk theory, with a mix ratio 0.98/0.02 of chromium yellow and cadmium red c) First derivative of the reflectance spectrum of the restored area and the simulated spectrum.

Figure 8 :

 8 Figure 8: a) Emission spectra of the original composition (red) and restored area (black) b) distribution map of the large band of luminescence observed centered at 500 nm, attributed to organic materials such as binder and varnish c) Distribution map of the large band of luminescence observed between 750 nm and 900 nm, attributed to Cadmium and Zink based pigments.

Figure 9 :Figure 10 :

 910 Figure 9: a) RGB display on the RIS SWIR map of the first three MNF components contribution b) Plot of the 8 first MNF components of the RIS SWIR datacube.


Crocherie Axel

Boulenc Pierre

Jérôme Vaillant

Hirigoyen Flavien

Hérault Didier

Tavernier Clément

From photons to electrons: a complete 3D simulation flow for CMOS image sensor

Keywords: Image sensors, Simulation, FDTD methods, Optical propagation, Semiconductor devices modeling, TCAD

I. INTRODUCTION

S pixel size shrinks in CMOS image sensors markets, the crosstalk effects become more and more critical. They can be divided in two components: the optical crosstalk due to the light propagation inside the pixel itself (including color filters transmission curves), and the electrical crosstalk related to the minority carrier diffusion and collection efficiency inside the substrate. Since these crosstalks are strongly interdependent, coupled Optical/TCAD simulations are mandatory to accurately reproduce this interaction and thus real pixel performances.

The next section deals with the methodology we used: on one hand for the Optical simulation tool, and on the other hand for the Electrical simulation tool. The third section presents the whole simulation flow calibration and the comparison between obtained results and experimental measurements.

II. METHODOLOGY

A. Optical part 1) FDTD

For small pixels that are used nowadays in CMOS image sensors, diffraction effects can substantially affect light propagation and photon collection. Thus, ray-tracing description is not accurate anymore [START_REF] Hirigoyen | FDTD-based optical simulations methodology for CMOS image sensor pixels architecture and process optimization[END_REF] and we must use a more fundamental description to simulate these optical effects. We chose to adopt an electromagnetic simulation tool based on Finite Difference Time Domain (FDTD) [START_REF] Yee | Numerical Solution of Initial Boundary Value Problems Involving Maxwell's Equations in Isotropic Media[END_REF][START_REF] Taflove | Computational Electrodynamics : the finite-difference time-domain method[END_REF], available from Lumerical Solutions [4], to describe light propagation and photon collection inside the pixels while correctly simulating diffraction effects.

FDTD is a fully vectorial method that gives both timedomain and frequency domain information by exploiting Fourier Transforms, i.e. when a broadband pulse is used as the source, the response of the system over a wide range of wavelengths could be obtained in a single simulation. Thus, this method is well adapted to our topic as we will be interested in a frequency averaging around the wavelength of interest (450nm, 532nm, and 633nm). We will discuss later about this subject.

2) Pixel modeling

Lumerical software allows us to import the GDS layouts created by designers to generate a 3-Dimensional representation of the pixels. Besides, specific shape of microlenses could also be imported in the software from the AFM data. This allows a complete and real modeling of the pixels as we could see on Fig. 1 below for a 1.75µm pixel. In the simulation, we are interested in the response of the blue, the green, and the red pixels of the structure. As the layout is periodic for the four pixels of the Bayer pattern, we will only simulate these pixels with periodic boundary conditions at the four lateral sides of the structure. Besides, absorbing boundary conditions are used at the top and the bottom of the structure to avoid unwanted reflections.

Materials used, like oxides, nitrides, silicon, color filters, are based on STMicroelectronics experimental database. Refractive indexes of these materials are parameterized in Lumerical tool to rigorously simulate their dispersion on the whole visible spectrum of the simulation [4,[START_REF] Crocherie | Three-dimensional broadband FDTD optical simulations of CMOS image sensor[END_REF]. Finally, metal for the interconnections are considered as perfect conductors (in this way, any wave reaching the metal line is reflected with the same incident angle and without any attenuation).

3) Accuracy of simulation: mesh step

The mesh step affects the accuracy of the simulation: a smaller mesh step will give more accurate results but requires more resources. One defines a parameter N λ that represents the number of mesh cells per wavelength as shown by:

n x N . ∆ = λ λ (1)
with ∆x the mesh step along one direction of the grid, λ the smallest wavelength of the simulation source, and n the highest Axel Crocherie, Pierre Boulenc, Jérôme Vaillant, Flavien Hirigoyen, Didier Hérault, Clément Tavernier From photons to electrons: a complete 3D simulation flow for CMOS image sensor STMicroelectronics, 850 rue Jean Monnet, 38920 Crolles, France A refractive index of the simulation structure. The mesh step has been fixed to N λ =10 here with the adaptive grid, a good tradeoff between mesh accuracy and resources requirement for 3D simulations [START_REF] Crocherie | Three-dimensional broadband FDTD optical simulations of CMOS image sensor[END_REF].

4) Source modeling

The source development is a key part of this methodology. We need a source that reproduces a product-like illumination, keeping reasonable computation time and memory. Regarding scaling problems between the objective-lens (several millimeters) and the pixels (several micrometers), we chose to adopt a local approach by simulating a group of pixels receiving the same uniform illumination from the exit pupil of the objective-lens (see Fig. 2). We have demonstrated in a previous paper [START_REF] Vaillant | Uniform illumination and rigorous electromagnetic simulation applied to CMOS image sensor[END_REF] that this kind of diffuse-like source could be represented by the incoherent sum of angularly uniformly distributed plane waves with incidence angles limited by the fnumber of the objective-lens. Interference inside the sensor back-end-of-line (BEOL) causes oscillations on the transmission spectrum. In order to smooth them but also to take into account the dispersion that exists in process for layers' thicknesses, optical simulations are made on a narrowband spectrum (±10nm with a 1nm sample step) around the wavelengths of interest. Then results are incoherently averaged.

Finally, we have to set the power of the source. As we make here an electro-optical coupling tool, we must simulate a realistic number of electrons in the silicon, i.e. a realistic number of photons at the input of the optical simulator. We consider here a pixel of 1.75µm pitch with a full well capacity of 9000 electrons [START_REF] Cohen | Fully Optimized Cu based process with dedicated cavity etch for 1.75µm and 1.45µm pixel pitch CMOS Image Sensors[END_REF]. At mid-saturation, with a quantum efficiency around 45%, this leads to roughly N photons =10,000 photons at the source input. Then, we could calculate the source power I λ (in W.m -2 ) using the following equation:

pixel photons A t h N I . . int υ λ = (2)
with hυ the photon energy, t int the integration time (equal to 66.7ms, i.e. 15 frames/second), and A pixel the area of the pixel.

5) Optical Generation

The information we want to extract from the optical simulations is the Optical Generation, i.e. the photogenerating rate of electron/hole pairs in the silicon in m -3 .s -1 , given by:

( ) υ σ υ h E h P G opt 2 . . Re r r r = ∇ - = (3)
with P the Poynting vector in W.m -2 . The next step is to inject G opt into the electrical tool to simulate the behavior of carriers in the substrate. As we want to understand the diffusion, recombination, and collection of all generated electrons, the depth in the silicon is parameterized to ensure that all the photons are absorbed (3µm in blue, 6µm in green and 10µm in red).

B. TCAD simulations 1) Principle

We use the Synopsys Sentaurus TCAD [START_REF]Synopsys TCAD tools[END_REF] suite to firstly perform process simulations that deal with modeling the frontend-of-line (FEOL) steps in the sensor manufacturing. Then, device simulations are the electro-optical behavior of the system in operating conditions. Both simulations are required to allow the most advanced description of complex physical, optical and electrical phenomenon encountered in CMOS Image Sensors.

2) 3D Bayer modeling

The FEOL part of the 3D Bayer has been built thanks to the following methodology. We firstly perform layout based process simulations of 2D domains that are perpendicular to the pixels transfer gates: these different cuts consider topology that might significantly impact doping distribution inside the pixels such as the presence / absence of Shallow-Trench-Isolation (STI). Then, the obtained end-of-process 2D doping distributions are extruded thanks to Sentaurus Structure Editor [START_REF]Synopsys TCAD tools[END_REF] along the third dimension to fill the pixel regions. Additional extrusions of 1D end-of-process specific implants are incorporated to fill the 3D doping distribution outside the domain covered by the previous 2D cuts. This strategy based on 2D process modeling ensures an adequate level of accuracy while keeping the simulation time compatible with large design of experiments.

Once the full 3D doping distribution is generated, the structure is meshed prior to device simulations. Our meshing strategy is to concentrate the smallest elements (20nm) in the photodiode / sensing-node regions to catch the strong doping gradients. Coarser (80nm) isotropic mesh is then used for photodiode surroundings (down to 1.1µm under the surface). Another coarser mesh (160nm) is applied between 1.1µm and 3µm for 450nm, 6µm for 532nm and 8µm for 633nm simulations. Finally, the bottom of the structure has the biggest elements (900nm).

Since the previously simulated Optical Generation data field is mapped on a Lumerical mesh, we need to interpolate it on a Sentaurus grid format. This mandatory step is performed on an equivalent tensor grid thanks to a routine with a maximum interpolation error that is smaller than 0.1%. Finally, the Optical Generation data is once again interpolated on the doping dependent mesh which will be responsible for another interpolation error.

A compromise must be found between high accuracy on Optical Generation gradients, doping distribution gradients and the final number of elements to allow reasonable CPU time. Both strongest Optical Generation and doping gradients are found in the same regions -i.e. close to the Silicon surface -therefore, a doping dependent mesh provides a satisfactory description. Our final structure has 340,000 elements (Fig. 3) and the resulting Optical Generation interpolation relative errors are respectively 6%, 3% and 2% at 450nm, 532nm and 633nm.

Fig. 3 3D Optical Generation in the Bayer interpolated on device simulation mesh (from left to right: λ=450nm, λ=532nm and λ=633nm). Half of the pixels has been blanked in the x and y directions to show light focalization inside the pixels.

3) TCAD device simulation setup

Sensor device simulations have been performed thanks to Sentaurus Device tool [START_REF]Synopsys TCAD tools[END_REF] which resolves the Poisson equation and both electrons and holes continuity equations in the structure.

The Image Sensor is a periodic array of Bayer structures, therefore, the use of periodic boundary conditions in both the x and y directions is mandatory to accurately consider electrical crosstalk i.e. electrons that may flow from one pixel to all its neighbors.

As in light capture operational conditions, the device simulation experimental protocol consists of the three following steps: pixels Reset, Illumination, Readout. The Reset is performed by initializing electron Fermi level to an arbitrary high value in the pixels in order to empty all the photodiodes. Then, thermal generation of electron-hole pairs (Shockley-Read-Hall model with Scharfetter doping dependence [START_REF] Roulston | Modeling and Measurement of Minority-Carrier Lifetime versus Doping in Diffused Layers of n+-p Silicon Diodes[END_REF][START_REF] Fossum | Computer-Aided Numerical Analysis of Silicon Solar Cells[END_REF][START_REF] Fossum | A Physical Model for the Dependence of Carrier Lifetime on Doping Density in Nondegenerate Silicon[END_REF][START_REF] Fossum | Carrier Recombination and Lifetime in Highly Doped Silicon[END_REF]) fills the photodiodes by dark current. Because of the sudden character of this artificial Reset operation, we let the system go back to a linear regime by waiting for 1ms. During that time, only dark current fills the photodiodes. Light is then turned on with the Lumerical Optical Generation map in the 3D Bayer and integration is performed for t int =66.7ms. Finally, the 3D electron density is integrated in each pixel 1ms after light is stopped to obtain the number of collected electrons. Another simulation without light is performed to obtain the number of dark current electrons to be subtracted in the final pixel quantum efficiency (QE) calculation ( 4 where λ is the wavelength, n is the electron density, t int is the integration time and Ф is the incident photon flux.

III. COUPLED OPTICAL/TCAD SIMULATIONS RESULTS

1) Measurements and calibration

For the calibration, we choose to limit the angular distribution of the source to limit the computation time and memory requirement. So we used plane wave for simulation and nearly collimated (f-number=200), narrow band beam for 1.75µm pixels QE characterization. The experimental setup is based on a stabilized halogen light source. The intensity is adjusted by inserting neutral density filters. The bands of interest (450nm±5nm, 532nm±5nm, 633nm±5nm) are obtained with interference filters. The illumination is recorded in real-time by a calibrated photodiode. Sensor under test provides images under standard operating conditions. Standard pre-processing (image average, dark subtraction) is performed in order to minimize the measurement noise.

Process simulations calibration has been performed on the basis of STMicrolectronics SIMS profiles and TCAD devices simulations database. It shows a good agreement with a large number of standard electrical characteristics on various devices. Furthermore, we add a constant doping profile in a 1µm range from the surface to account for 3D effects including Well-Proximity-Effects that are not taken into account in our 2D simulations. Such effects reduce 3D photodiode extension, thus photodiode efficiency.

The coupled Optical/TCAD simulations results are hence used to compute the QE(λ) accounting for both optical and electrical crosstalks. Following results aim at demonstrating that our coupled simulation strategy is calibrated and reproduces various BEOL and FEOL trials.

2) Bayer without color filters Fig. 4 shows the quantum efficiency of a Bayer without color filters. Moreover, this Bayer corresponds to a specific FEOL design where voluntary different doping distributions have been realized inside each pixel. Although at 450nm all pixels behave in a similar manner, this is no longer the case at 532nm and 633nm. Indeed, Green and Red pixels catch more electrons than the Blue one because of its specific isolation scheme. The small difference between Green and Red pixels is explained by a smaller Red photodiode extension. Note that at each wavelength, a scaling factor is applied on simulated Blue, Green and Red QE(λ). It is calculated as the ratio between measured Bayer average QE(λ) and simulated Bayer raw average QE(λ). This factor reflects imperfections that are not considered in optical simulations (like interface roughness) and is consequently wavelength dependent.

3) Bayer with color filters

Adding color filters in the BEOL strongly affects the overall Bayer response because of color filters transmission (Fig. 5). FEOL process #2 has identical doping distribution in Green and Red pixels. The strong response of Red pixel at 532nm and Green pixels at 633nm indicate the presence of crosstalk between neighboring pixels. In order to evidence that our simulation tool well describes the evolution of electrical crosstalk with respect to various pixels doping distributions, we simulated two additional FEOL process changes that can be compared with measurements. FEOL process #3 is a slight modification of the doping concentration in depth for all pixels. Its response (Fig. 6) shows it is possible to improve Green QE while reducing Red pixel QE at 532nm by acting on electron diffusion in depth. This behavior is confirmed by the reduction of Green QE at 633nm which is well reproduced by our simulations.

Accordingly, our simulations also well predict that FEOL process #4 (Fig. 7), which adds a cut-off in depth, reduces Red QE at 633nm in particular. As expected, this also strongly reduces both 532nm and 633nm electrical crosstalks.

IV. CONCLUSION

In very small pixels, crosstalk becomes one of the main performance limiting factors, with strong 3D and neighboring effects. Moreover crosstalk arises from two distinct phenomena (optical and electrical) that are in strong interaction. Therefore, an approach that couples dedicated simulations tools is mandatory. Our strategy is thus to link Optical simulations to TCAD process and device modeling.

The present study shows the ability to describe both BEOL and FEOL process changes in agreement with sensor characterization data. Although some adjustments are required, like global optical stack transmission scaling and 3D corrections to 2D doping distributions, electro-optical performances of four different sensors can be reproduced within a few percents error.

This methodology, validated on 1.75um pixels, can easily be implemented on more complex Image Sensors technologies.

Fig. 1

 1 Fig. 1 3D modeling of pixels with the optical simulation tool

Fig. 2

 2 Fig. 2 Light shape in the case of a uniform pixel illumination provided by an objective-lens

Fig. 4

 4 Fig.4Pixels QE(λ) without color filters in the BEOL and with FEOL process #1. Scaling factors are 0.84 (λ=450nm), 0.89 (λ=532nm) and 0.88 (λ=633nm). Normalization is performed vs. Green QE at 532nm.

Fig. 7

 7 Fig.[START_REF] Cohen | Fully Optimized Cu based process with dedicated cavity etch for 1.75µm and 1.45µm pixel pitch CMOS Image Sensors[END_REF] Pixels QE(λ) with color filters in the BEOL and with FEOL process #4. Scaling factors are 0.82 (λ=450nm), 0.86 (λ=532nm) and 0.69 (λ=633nm). Normalization is performed vs. Fig.5Green QE at 532nm.