A Savoyant

M Rollo

M Texier

R E Adam

S Bernardini

O Pilone

O Margeat

O Nur

M Willander

S Bertaina

Light-induced High-spin State in ZnO Nanoparticles

Keywords: numbers: 32.30.Dx, 32.80.Fb, 75.75.-c, 76.30.-v

The effects of white-light irradiation on ∼ 20 nm diameter ZnO nanoparticles are investigated by means of electron paramagnetic resonance, near liquid-nitrogen and liquid-helium temperatures. Under dark conditions, usual core-and surface-defects are detected, respectively at g = 1.960 and g = 2.003. Under white-light illumination, the core-defect signal intensity is strongly increased, which is to be correlated to the light-induced conductivity's augmentation. Beside, a four-lines structure appears, with the same gravity center as that of the surface defects. Simulations and intensity power-dependence measurements show that this four-line-structure is very likely to arise from a localized high spin S = 2, induced by light irradiation, and subjected to a weak axial anisotropy. At 85 K, this high-spin state can last several hours after the light-irradiation removal, probably due to highly spin-forbidden recombination process. The possible excited resonant complexes at the origin of this signal are discussed. Other light-induced S = 1/2-like centers are detected as well, which depends on the nanoparticles growth conditions.

I. INTRODUCTION

Zinc Oxide (ZnO) is an extensively studied metal oxide semiconductor, due to its remarkable physical properties 1 .

It is a n-type wurtzite semiconductor with a direct bandgap of about 3.37 eV at room temperature and relatively large exciton binding energy of 60 meV. Naturally present in abundance, non toxic and of low cost, ZnO material is widely used in electronic and optoelectronic due to its high electron mobility 2 , photoelectric response 3 , excellent chemical and thermal stability 4,5 , piezoelectric applications 6 , gas sensors 7-9 , light emitting diode 10,11 , and so on. From the physical point of view, one of the most striking characteristic of ZnO is its strong interaction with light, resulting in photo-induced effets, such as photo-conductivity 12,13 , photocatalysis [14][15][16][17] and photo-luminescence 18,19 . All of these effects are mainly driven by excitonic and point-defects related properties, as well as coupling of both [20][21][22][23][24] . In nanosized ZnO systems, these photo-induced effects can be affected by the important surface-to-volume ratio, and/or by the presence of specific point-defects acting for example as new recombination or charge donor/acceptor centers. The study of these sizeinduced defects is thus a matter of great interest for tailoring efficient nanostructured ZnO-based devices.

Electron paramagnetic resonance (EPR) is a tool of choice to study point defects. This technique has been used for more than fifty years for characterizing intrinsic and extrinsic param-agnetic defects in bulk ZnO [25][26][27][28][29][30][31] . Regarding ZnO nanostructures, EPR studies have shown that two major kinds of intrinsic paramagnetic defects are almost always present within, resonating at g = 1.960 and g = 2.003, and arising respectively from core-and shell-defects [32][33][34] . However, the exact nature of these defects is still under debate. A mean of gaining insight into the EPR lines origin is to introduce an additional excitation during the experiment, such as an external electric field or an optical light beam. Indeed, the idea of applying a light irradiation during an EPR experiment goes back to the very early stage of its application on semiconductor materials, where light is likely to ionize point defects into intraband-gap states 35,36 . Regarding ZnO materials, most of these photo-EPR experiments have been performed on bulk single crystal [37][38][39] , showing in particular the apparition of some angle-dependent lines, attributed to excited oxygen and zinc vacancies 40 . For ZnO nanostructures, the main result of such an irradiation is to greatly increase the core-defect (g = 1.960) intensity and, to a lesser extent, that of the shell-defect (g = 2.003) 41 . However, such photo-EPR studies have not been exhaustively performed onto ZnO nanoparticles (NPs), in particular regarding their size, growth method and morphology, all characteristics known to greatly influence the population of defects.

In this work, we report on photo-EPR experiment performed on two types of ZnO NPs, near liquid-helium and liquid-nitrogen temperatures. The main new feature observed here is a four-line structure whose gravity center is that of the shell-defect (SD), g = 2.003, reversibly induced by white light irradiation. Microwave (MW) power-dependence study near liquid and helium temperature reveals sub-structures of this four-line group, as well as other photo-induced S = 1/2-like lines, not apparently linked together. This fourlines group appears to be highly saturated around 10 K and still partially saturated at 85 K. Simulations of the unsaturated spectrum strongly support a S = 2 effective-spin origin, subjected to a weak axial anisotropy. This effective high-spin state is likely to arise from particular excitonic complex, involving one or more point defects. In this paper, after explaining the experimental details of materials and method (Sec. II), we present the raw EPR results under dark (Sec. III A) and illuminated (Sec. III B) conditions, followed by spectra simulation based on a S = 2 spin state (Sec. III C). Finally, a qualitative discussion about the origin of this high-spin state is given in the last section (Sec. IV).

II. EXPERIMENTAL

The reported EPR measurements were performed on two kinds of hydrothermally grown ZnO NPs, NP1 and NP2, with slightly different preparation methods, but important morphological difference. NP1 were prepared as published elsewhere [START_REF] Diallo | Superlattices Microstruct[END_REF] . In a typical experiment, a solution of 0.4415 g of potassium hydroxide KOH (Aldrich Chemicals 99.99%) in 23 mL of anhydrous methanol (Alfa Aesar) was added dropwise to a solution of 0.8202 g of zinc acetate (Aldrich Chemicals 99.99%) in 42 mL of anhydrous methanol and 0.25 mL of distilled water at 60 • C under magnetic stirring. After 2h30 min at this temperature, the mixture was centrifuged (7800 rpm, 30 min) and the supernatant discarded. The product was washed with methanol and centrifuged, then finally dried under vacuum to obtain the ZnO nanoparticles as powder. For NP2 fabrication [START_REF] Adam | [END_REF] , Zinc acetate dehydrate (0.1 M) and sodium hydroxide (0.2 M) solutions were prepared separately and magnetically stirred overnight at room temperature. Then, the two solutions were mixed into one beaker and magnetically stirred at 750 rpm for 2 hours at temperature of 60 • C. The precipitated white milky solution was separated by centrifugation at 4500 rpm for 2 minutes. Afterward, the precipitation product was washed several times with deionized water and once time with acetone. Finally, the ZnO NPs in the form of powder was obtained by drying the precipitation product into a laboratory oven at 75 • C for 6 hours. Transmission electron spectroscopy (TEM) characterizations of both samples (see Fig. 1) show that NP1 are almost spherical in shape with typical diameter of 15 nm, while NP2 have spindle-like shape with 400 nm and 100 nm in long and short dimensions, respectively. This points to the great versatility for size and shape of ZnO NPs, and to its great sensitivity to growth parameters. Both EPR measurements have been performed on a conventional Bruker EMX continuous wave spectrometer operating at X-band (ν = 9.402 GHz) using a standard TE 102 mode cavity with optical access. The low (11 K) and intermediate (85 K) temperature conditions were obtained by cooling down the samples with liquid helium and nitrogen within a cryostat inserted in the EPR cavity. The microwave (MW) incident power was controlled by an attenuator, allowing power value from 0 to 200 mW (P 0 ). Improvement of the signal-to-noise ratio has been achieved by scanning EPR spectra from 4 to 81 times. White-light irradiation was performed by focusing the photon beam from halogen lamp toward the sample, via the cavity hole. The white-light source has been calibrated so that an optical power of approximatively 80 mW is reached on the sample. Additional experiment with violet laser mentioned in Sec. IV, was performed with a Radius 405 (nm) laser of 20 mW.

III. RESULTS

A. Light-off EPR spectra

In figure 2, we present the results of EPR experiments recorded on the NP1 sample at temperatures T = 9 K and T = 85 K, under dark conditions (light off). Two lines are observed, at g = 1.960 and g = 2.003, which have been reported many times in EPR measurement of ZnO nanostructures [32][33][34] . They are commonly identified as arising from core-and shell-defect, respectively. The g = 1.96 signal arises from most of the ZnO nanostructures 44,45 , regardless of growth method and morphology, but also from bigger microsized crystal as ceramics and bulk single crystal 25,37 . On the contrary, the g = 2.003 signal only appears in NPs with the smallest diameter (< 50 nm), when surface defects are assumed to play an important role. This is the foundation of the core-shell model (CSM) proposed by Erdem et al., which describe a very stable ZnO NP's core, coated by a shell subjected to interactions with the surroundings 32,33 . The NP2 sample spectrum (not shown) displays identical core-defect (CD) signal, while SD signal is not detected at none of the studied temperatures. This is consistent with the NPs bigger size in NP2, for which SD signal is not expected to occur. From the power-dependence of NP1 EPR lines intensity (insets (a) and (b) of fig. 2), qualitative information on the spin's dynamics is obtained. If not saturated, EPR intensities vary as the squareroot of incident the MW power, which results in straight dashed lines in the insets of figure 2. All departures from these lines indicate a (partially or totally) saturated spin resonance. It can be seen that the SD signal saturates at weaker power than the CD one, so that we can draw the rough picture that SD has spin-lattice relaxation time significantly greater than that of the CD. This may be understood within the coreshell model of NPs, in which the CD spin must be more strongly coupled/linked to the lattice than the SD spin.

B. Light-on EPR spectra

When the white-light source is turned on and focalized on samples, two main drastic changes are observed in the EPR spectra. First, the CD signal's intensity is greatly enhanced with time constant of several minutes at RT, depending on the power light intensity. The intensity reaches a saturation value which also depends on the light-power intensity. This first point is correlated to the observed photoconductivity in ZnO NPs 30 . When white-light is turned off, the enhanced CD signal lasts for several hours at 85 K, which is certainly related to the (problematic) persistent photoconductivity observed in the ZnO material 12,13,21 . The g = 1.960 signal thus certainly arises from a donor center, located within the NPs core. The second remarkable change in the EPR spectra is the apparition of a four-line structure (FLS), whose gravity center is that of the SD signal. This is the photo-induced EPR signal studied and discussed in this paper.

Figure 3 shows the result of white-light illumination during an EPR experiment on NP1 sample: an intense four line structure (FLS) appears, centered at the SD g-factor (2.003), and whose line are denoted L 1 to L 4 . Remarkably, these new lines have equal spacing (∆ = 2.25 mT) and equal peak-to-peak line width (L pp = 0.4 mT), pointing to a single origin, possibly related to the SD. Such a relation is not obvious because the SD g-factor is so close to its free-electron value (2.0023), where every isotropic 1/2 spin resonate. We note that the observed lines shape are clearly not those of a simple EPR line, due to unresolved lines in between and/or to saturation phenomena.

In order to get insight into the relaxation dynamics of these EPR transitions, the temperature is slightly increased (from 9 to 11 K) and the MW power is varied from 0.007 to 1.584 mW. As shown in figure 4 The spectra dependence on MW power (Fig. 4) shows that the three additional-lines intensity (L 5 to L 7 ) are almost power-independent, while the FLS intensity notably decreases with power. This signs a basic difference between the dynamics of these lines and the FLS, the formers having longer relaxation time than the latter. L 5 and L 7 displays no evidence for a relation with the FLS, whereas L 6 (g = 2.003) is likely to be the original SD defect itself, which is also fully saturated at 11 K. Indeed, in other kind of samples in which the FLS is observed under illumination (see NP2 Fig. 10 and Sec. IV), L 5 is also present with different g-valOmerue, while L 6 and FLS have the same position. It may be assumed that the photo-induced FLS is related to the SD, while being although of different nature, and that L 5 and L 7 arise from versatile excited defects, depending on growth method and conditions. The satellite lines of the FLS appear only at highest powers, and each of them is located at a distance value ∼ ∆/4 from the nearest FLS line. Interestingly, L 2 and L 3 appear to be accompanied by two of these satellite lines, while L 1 and L 4 by only one. These satellites are likely to arise from multiphotons absorption, forbidden transitions, or low symmetry fields. The increase of their intensity is mainly responsible for the line-shape deformation at high power.

Regarding the FLS itself, its intensity varies notably with MW power variation, while preserving a constant ∆ spacing and gravity center, and it becomes more and more asymmetric as the MW-power decreases. However, it can be seen on Fig. 5 that intensity of lines L 1 and L 4 do not vary anymore from P = 0.405 mW to P = 1.584 mW, contrarily to the center lines L 2 and L 3 . This indicates a slight difference in the dynamics of these transitions, which will be quantified in the following and discussed in Sec. IV. The FLS line widths are not affected by the MW power variation.

We now turn to the examination at higher temperature (85 K) of the FLS in order to observe fully unsaturated resonance, for it is the only case which can be treated by simulations. EPR spectra recorded on NP1 sample under illumination at T = 85 K are shown in Fig. 5, for different MW powers. The widths of the L 5 and L 6 lines are weaker than at T = 11 K, certainly due to partial desaturation, which allows for an almost complete resolution of them. On the contrary, the L 7 line is still not resolved. Line widths of FLS are reduced from 0.4 mT to 0.3 mT. The satellite lines are relatively much weaker, still only appearing for the highest MW powers. The interesting feature of the spectra plotted in figure 5 is the observation of three completly different EPR-line dynamics. First, L 1 to L 4 lines do not reach a total saturation regime, since their intensity always increases with increasing MW power. Beside, the L 5 line reaches a total saturation regime for MW power of 0.79 mW, while L 6 and L 7 lines appears to be always fully saturated, at any MW power. This is a srong indication that all these line groups arise from differents photo-excited paramagnetic centers. A key physical information to be extracted from these spectra is obtained from their intensity MW power dependence, which gives insight into the dynamics of each EPR line, as well as into the spin-resonance regime (saturated or unsaturated). Peak-to-peak height for lines L 1 to L 4 (proportional to absorbed intensities for constant peak-to-peak line width) have been plotted as functions of MW power in Fig. 6, at 11 and 85 K. It is to recall that an EPR line is said to be not saturated if by doubling √ P , the intensity is doubled. From this simple criterium, we see that at T = 85 K, the FLS is not saturated only for incident MW power as weak as P = 0.05 mW. At T = 11 K, whatever the weakness of the incident MW power is, the FLS structure is always saturated: the required power for a total unsaturated regime of the FLS is so low that the lines cannot be detected anymore in our experimental conditions.

More quantitative information about EPR dynamics is obtained by the fitting of the saturation curves with an appropriate function including physical parameters. The simplest modeling of the saturation phenomena is derived from the isotropic Bloch equations, including longitudinal and transverse relaxation times T 1 and T 2 . It results in the following power dependence of the EPR intensity [START_REF] Weil | Electron Paramagnetic Resonance[END_REF] :

I(P ) = K. √ P 1 + A.P (1)
The parameter A is proportional to the product of the relaxation times T 1 (spin-lattice) and T 2 (spin-spin). It determines the maximum of the saturation curve. The parameter K contains the square of the transition probability, the Boltzmann factor, and the T 2 time. The fitting of EPR intensities vs √ P is displayed in Fig. 6, which show an overall good agrement. At T = 85 K, parameter A is found to be almost identical for the four considered lines. However, at T = 11 K, saturation curves clearly show that L 1 and L 4 reach maximum at MW power substantially weaker than L 2 and L 3 . The fitting allows to determine that A 1,4 /A 2,3 ≃ 3. By assuming that T 2 is almost the same for all these lines, it is finally deduced that the spin-lattice relaxation time T 1 is notably longer for extreme lines of the FLS (T 1,4 1 /T 2,3 1 ≃ 3). This fact will be of importance in the following, when searching for the appropriate spin Hamiltonian to fit the FLS. Finally, we note that after light-irradiation removal the FLS remains detectable for several hour at 85 K, which will be discussed in Sec. IV.

C. Simulations

In order to gain insight into the origin of the photo-induced FLS, spectra simulations have been performed, based on assumed spin Hamiltonian and spin state. It is found that two radically different spin systems can give account for the observed FLS. This latter can arise either (1), from superhyperfine coupling of one electronic 1/2 spin with 3 adjacent protons, or (2), from a S = 2 spin under weak axial anisotropy. In the following, we present arguments which make the first scenario very unlikely to exist. Then, we show that the second one can give well account for the experimental spectra.

In fact, a complex made of one 1/2 spin (hole or electron) surrounded by 3 protons can exist in ZnO [START_REF] Van De Walle | [END_REF]48 . These protons can be provided by the ubiquitous unintentional hydrogen-doping, and the 1/2 spin can be carried by a point defect. However, in such a case, signals arising from simpler complexes (involving only 1 or 2 protons) must then occur with higher probabilities, which are indeed not detected. In addition, due to the identical ∆ spacing between lines, the superhyperfine coupling constant would have to be exactly the same for all protons, this constant (A = 66 MHz) being moreover very intense for such a distant interaction (A = 1.4 MHz for donor spin coupling with adjacent hydrogen in ZnO 49 ). It is also to be noted that the FLS has been observed in other kind of hydrothermally grown samples with exactly the same ∆ spacing: this indicates a stable and reproducible property, poorly compatible with random distribution of unintentional hydrogen impurities. Finally, the most important argument against this first scenario, and in favor of the second one, is provided by the saturation curves of the FLS (Fig. 6). These show that at T = 11 K the extreme lines L 1 and L 4 have not the same dynamics than that of the center lines, L 2 and L 3 . This has been quantitatively demonstrated by the fits of Fig. 6, leading to the ratio of spin-lattice relaxation times T 1,4 1 /T 2,3 1 ≃ 3. In the case of hyperfine couplings, however, all of the four lines would arise from | -1/2 ↔ | + 1/2 electronic-state transitions, and must then have the same spin-lattice coupling, thus saturating at the same MW power.

On the contrary, if the FLS arise from a S = 2 spin under weak axial anisotropy, this difference in dynamics is expected because, as shown in Fig. 7-b, the extreme lines L 1 and L 4 arise from transitions | ± 2 ↔ | ± 1 , while the center lines L 2 and L 3 do from transitions | ± 1 ↔ |0 . Different electronic spin states having different coupling with lattice's phonons, the above-mentioned transition must have different dynamics, consistently with our observations. In the following, based on the above-mentioned arguments, we focus on the S = 2 scenario and present the simulated spectra based on the spin Hamiltonian written in Figure 7-b. This latter contains the Zeeman interaction with the static magnetic field B, and an axial anisotropy term characterized by the D parameter (D < gµ B B). For θ = 90 • orientation, the spacing between lines simply equals D, in magnetic field units. The simulated powder spectrum is strongly asymmetric, and is dominated by the most occurring orientation, that is θ = 90 • . Then, the experimental ∆ spacing is a direct measure of the D parameter.

Regarding simulations, it must be kept in mind that only unsaturated resonance of one-photon transitions can be accounted for by the current EPR software as EASYSPIN 50 used here. Rigorously, the spectra to be simulated must be restricted to those recorded in unsaturated conditions, which is only the case at T = 85 K and P = 0.05 mW (see Fig. 6). This experimental spectrum displays the most asymmetric FLS, that well corresponds to the powder spectra of S = 2 with axial anisotropy (Fig. 7-a), and no satellite lines around the FLS.

In order to reproduce this unsaturated experimental spectrum, the powder spectra of S = 2 with axial anisotropy (Fig. 7-a) is used, involving only one energetic parameter (D).

Also the three lines L 5 to L 7 evoked before (Sec. III B) are added, plus another one (denoted L 8 ) which appears necessary to obtain an overall good agreement. The resulting simulation, presented in figure 8, appears very satisfactory. The axial anisotropy parameter is found to be |D| = 66 MHz (= 22 × 10 -4 cm -1 ), which corresponds to a very weak axial anisotropy. Its sign, which depends on the in-plane or out-of-plane character cannot be determined in our experimental conditions (too high temperature). The number of 1/2 spins used in this simulation (4) may appear as too important for drawing conclusion about the FLS. However, these additional 1/2 spins are clearly not directly related to the FLS 59 , so that their introduction in the simulation does not prevent to conclude about the actual presence of a S = 2 local electronic spin.

The unsaturated spectrum being well reproduced by simulation, the saturated case is now examined in order to understand the difference between both. Apart of the satellite lines apparition briefly discussed in Sec. III B (not taken into account in simulation), the main characteristic of this experimental saturated spectrum is the almost symmetric shape of the FLS, which is not compatible with a simple anisotropic S = 2 powder (see Fig. 7-a). The explanation can be that powder-spectrum simulation does not explicitly includes anisotropic line widths 60 , while it is likely to actually arise in ZnO saturated EPR. Surprisingly, it can be noticed that the saturated EPR spectrum tends to look alike the B ⊥ z simulation (Fig. 7-a, θ = 90 • ). This can be understood if the line width along the other direction (θ = 0 • ) becomes more and more important as the MW power is increased, so that its apparent height decrease. As a consequence, the powder spectrum is disturbed and tends to the oriented θ = 90 • direction as the saturation is reached. Figure 9 presents the comparison between the saturated experimental spectrum at 85 K and the simulation. This latter include the S = 2 powder spectrum, the four additional lines L 5 to L 8 , but also a contribution from the θ = 90 • orientated S = 2 simulation. We may consider this last contribution as an effective way of taking into account a strong anisotropy of the EPR line width, arising in saturated cases.

In the following, it is proven that the photoinduced FLS is reproducible in other kind of dried hydrothermally-grown sample, with different size and morphology (NP2 sample). In figure 10, the EPR spectrum recorded on NP2 under white-light illumination is shown, together with its corresponding simulated EPR spectrum. First, the photo-induced FLS is also observed, with exactly the same ∆ spacing (2.25 mT) and almost the same gravity center (g = 2.004 vs g = 2.003 in NP1). This points to the great stability of the paramagnetic center responsible for this EPR signal. Interestingly, the L 7 and L 8 lines observed in NP1 are not detected in NP2, while L 5 is now positioned at a different g-value. This support the assumption that although all these lines are photo-induced by white light, they arise from independent paramagnetic centers and can exist independently from each other. The L 6 line is still present at the gravity center of the FLS, but with greater relative intensity. Indeed, in this NP2 sample, the FLS is much weaker than in NP2 and high MW power is required in order to obtain acceptable signal-to-noise ratio. The presented experimental spectrum of figure 10 is thus recorded in saturated conditions (P = 6.32 mW). Regarding the simulation, and consistently with the above-presented results, such a saturated case needs the contribution of the θ = 90 • orientated S = 2 spin, for giving account of the symmetric shape of the FLS. The resulting simulation can be viewed as satisfactory.

It is then concluded that there actually exist a photo-excited S = 2 localized magnetic state in ZnO NPs, together but independently of other S = 1/2 photo-induced paramagnetic centers. In the following, the possible origins for this high-spin state are discussed.

IV. DISCUSSION

Regarding the existance of the S = 2 spin , it must be noticed that such an even high spin state is quite rarely observed in matter, where strong local anisotropy often prevent the related EPR transitions to occur (∆m S = ±1 required). The present observation of S = 2 spin with very weak axial anisotropy (|D| = 66 MHz) thus excludes the possibility for a transition metal ion origin, such as the 5 D ground multiplet of Fe 2+ ion, which would be expected to have larger axial anisotropy parameter D. For comparison, D = -710 MHz in ZnO:Mn 2+ (S = 5/2) 51 , and D = +82 743 MHz in ZnO:Co 2+ (S = 3/2) 52 .

Focusing on an intrinsic-defect origin for the observed photo-induced FLS, we recall that the assumed S = 2 spin-state must be viewed as effective, possibly made of orbital and spin magnetic momentum.

Without orbital contribution, the existence of a S = 2 spin requires the local coupling of at least four 1/2 spins (bear by electrons or holes). Such a situation can exist if two neutral vacancies defects in triplet state (S = 1) are ferromagnetically coupled. Indeed, it has been theoretically demonstrated that neutral zinc vacancy (two localized holes) has triplet ground state, while neutral oxygen vacancy (two localized electrons) has triplet excited state, lying at about 2-3 eV above its singlet ground state 53 . It is thus conceivable that, if two of these vacancies lay closed together enough, a photon absorption may causes this complex to transit to an S = 2 spin state. However, such a spin-pair would be expected to have different coupling constant, according to the relative position of each spin within the ZnO matrix, leading to a distribution in the axial anisotropy parameter D.

On the contrary, the observed S = 2 spin appears to arise from a unique local center, with stable and identical properties among different samples. A plausible explanation for such an intrinsic high-spin center created by light irradiation in ZnO can be related to excitons, which are known to be possibly bounded to defect 23,[54][55][56] . At least three different coupling schemes involving excitons can be thought about. (i) a tripletstate (S = 1) exciton can be ferromagnetically bounded to a triplet-state defect (neutral V Zn or V 0 discussed above), resulting in the observed quintet state. (ii) Two triplet excitons can be coupled together, forming a so-called high-spin state bi-exciton 57,[START_REF] Weiss L R | Dynamics are completly different and, moreover[END_REF] , potentially localized at the surface of NPs. (iii) In cubic crystal of local symmetry T d , single exciton can exist in a quintet state (S ≡ J = 2 effective spin), if the hole in the J = 3/2 (Γ 8 ) valence band and the electron in the J = 1/2 (Γ 6 ) conduction band are coupled by ferromagnetic exchange interaction. In wurtzite crystal (local symmetry C 6v ), the J = 3/2 valence band is moreover split into ±3/2 (Γ 9 , heavy hole) and ±1/2 (Γ 7 , light hole) band 56 . Again, in case of ferromagnetic exchange coupling between hole and electron , this situation can also give rise to an effective S ≡ J = 2 spin, but with weak axial anisotropy, as observed here. This last scenario is probably the most relevant, because it involved only one simple excitation. In any cases, even that of a FLS extrinsic origin, high local electric fields at the ZnO NPs surface certainly play a major role in the formation and localization of such excitonic high spin.

Among the experimental facts which can confirm or disconfirm the proposed scenarios, there may be the apparition of the satellite lines at high MW power, which are not taken into account by our simulations. These are likely to arise from multiphotons absorption, but their position are difficult to explain, in particular the single satellite of the extreme lines L 1 and L 4 , moreover outside of the FLS. Forbidden transitions or low symmetry Hamiltonian components mays also have to be considered to explain the satellite lines. In parallel, the study of the FLS signal decay after light removal gives a key information for its understanding as well as for its application. So far, it has been observed that the FLS has a very long life time, losing about 10 % of its intensity value two hours after the light removal, at 85 K. The only mean of speed up this decay is to heat the sample up to room temperature. This unusual very long life-time may be interpreted as the result of the high-spin state forbidden recombination, in a sort of "magnetic phosphorescence". This point is particularly interesting for potential application.

We now point to the intriguing fact that such an interesting photo-induced EPR signal has not been reported before to our knowledge. In particular, in their work onto EPR of ZnO NPs, Erdem al. have illuminated their samples with red (671 nm), green (532 nm) and blue (445 nm) laser lights, without seeing any photo-induced FLS 41 . The reason is certainly that these excitation energies are simply too low to give rise to the FLS. Indeed, we have performed the same photo-EPR experiment with violet laser light (405 nm), and observed the FLS just as with white light illumination. From these preliminary experiments, it can be assumed that this minimal energy lay between 2.79 eV (445 nm) and 3.06 eV (405 nm), which would be compatible with a highly bound exciton. A complete spectral study (up to UV energies) of this photo-induced FLS would be an important source of information for determining its origin. For example, if the FLS (L 1 to L 4 ) and the additional lines (L 5 to L 8 ) are really independent, their minimal excitation energy must certainly be different. Then, illuminating the sample with the lowest of these energies excitation would permit to separate the lines. Nevertheless, excitation energy is not the only crucial parameter for observing the photo-induced FLS. In fact, in some other samples of hydrothermally grown ZnO NPs, the photo-induced FLS is not observed. Morphology, solvent, drying and storage conditions, must have to be considered in a systematic study in order to find the exact con-ditions for favoring the FLS presence. However, we have observed it in several different samples, grown in two different labs, with exactly the same ∆ = 2.25 mT spacing between lines. We thus believe it to be a reproducible intrinsic property of ZnO NPs. In any cases, the control and the manipulation of such a high-spin state carried by a single NP open the way to promising possibility in terms of information storage or processing. Since a particular maximal wavelength must be able to switch on the high-spin state, it is very likely that some other wavelength would be able to switch it off. Thus, a light-writable and -erasable magnetic memory would be realized, readable by EPR. In such a context, the observed long life-time of the photo-induced FLS would be an advantage to preserve information with minimal refreshing. The presented results could path the way for reversibly encoding magnetic information within a stable single crystalline NP.

V. CONCLUSION

In this paper, we have reported on a new EPR signal, the FLS, reversibly induced by white-light irradiation in ZnO NPs. Power-dependence, temperature studies and spectra simulations show that this signal is very likely to arise from a localized high-spin state S = 2 with weak axial anisotropy. At T = 85 K, this spin state can last several hours after the light-irradiation removal, certainly due to highly spin-forbidden recombination process. The intrinsic character of this light-induced EPR signal is suggested, most probably related to boundexciton physics. Other photo-induced EPR lines without structure and with different spin-lattice dynamics are detected, not apparently linked to the FLS.

VI. ACKNOWLEDGEMENT

The authors acknowledge Prof. Bruno Guigliarelli for useful discussion and Clément Sébastio for help in the optical installation.

Figure 1 :

 1 Figure 1: (color online) Transmission electron spectroscopy images of NP1 and NP2 samples, showing differences in size and morphology.

Figure 2 :

 2 Figure 2: (color online) Dark EPR spectra of ZnO nanoparticles recorded at T = 9 and 85 K on NP1 sample, with different MW powers, showing the shelldefect (SD) and core-defect (CD) lines. Insets (a) and (b) show the saturation curve of CD at 85 K and those of CD and SD at 9 K, respectively.

Figure 3 :Figure 4 :

 34 Figure 3: (color online) Dark (up, black) and illuminated (down, red) EPR spectra of ZnO NPs, recorded at 9 K on NP1 sample. The dark spectrum is magnified by factor 3. ∆ = 2.25 mT

Figure 5 :

 5 Figure 5: (color online) EPR spectra under illumination recorded on NP1 sample at 85 K (colored, bottom) for different MW powers, and at 9 K (black, top) for P = 0.21 mW. Magnetic-field modulation is 2 G. Satellite line are indicated by black marks.

Figure 6 :

 6 Figure 6: (color online) Experimental saturation curves of the FLS (scatters) recorded on NP1 sample at T = 85 K (a) and T = 11 K (b), and fit from Eq. (1) (solid lines).

Figure 7 :

 7 Figure 7: (color online) (a) Simulated EPR spectra for S = 2 spin with axial anisotropy D for powder, B z (θ = 0 • ) and B ⊥ z (θ = 90 • ). (b) Corresponding Hamiltonian and levels scheme.

Figure 8 :

 8 Figure 8: (color online) Simulated (red) and experimental (black) EPR spectra recorded on NP1 at 85 K in unsaturated conditions (P = 0.05 mW) under illumination. Various contributions to the simulation are display (blue), all divided by a factor 3.

Figure 9 :

 9 Figure 9: (color online) Simulated (red) and experimental (black) EPR spectra recorded on NP1 at 85 K in saturated conditions (P = 3.18 mW) under illumination. Various contributions to the simulation are display (blue), all divided by a factor 3.

Figure 10 :

 10 Figure 10: (color online) Experimental and simulated spectra recorded at T = 85 K and P = 6.32 mW on the S2 sample. All simulation's contributions, with their respective weight, are divided by 3.

1

  Özgür U, Hofstetter D and Morkoç H 2010 Proc. IEEE, 98, 1255 2 Wang Z L 2008 ACS Nano 2 1987-1992 3 Kundu S and Nithiyanantham U 2014 Eng. Chem. Res. 53 13667 4 Roth P and Williums D F 1981 J. Appl. Phys. 52 6685 5 Suchea M, Christoulakis S, Moschovis K, Katsarakis N and Kiriakidis G 2006 Thin Solid Films 515 551 6 Ambia M G, Islam M N and Obaidul Hakim M 1992 Journal of Materials Science 27 19 5169 7 Zhu L and Zeng W 2017 Sensors and Actuators A:


➢ Euclidian distances [START_REF] Babel | Phonetic and social selectivity in speech accommodation[END_REF][START_REF] Ruch | Vowel convergence and divergence between two Swiss German dialects[END_REF] ➢ Discriminant analysis [START_REF] Delvaux | The influence of ambient speech on adult speech productions through unintentional imitation[END_REF] ➢ Moving averages [START_REF] Tamminga | The dynamics of variation in individuals[END_REF] ➢ GAMMs [START_REF] Sonderegger | The medium-term dynamics of accents on reality television[END_REF][START_REF] Solanki | Brains in dialogue: Investigating accommodation in live conversational speech for both speech and EEG data[END_REF] ➢ Euclidian distances [START_REF] Babel | Phonetic and social selectivity in speech accommodation[END_REF][START_REF] Ruch | Vowel convergence and divergence between two Swiss German dialects[END_REF] ➢ Discriminant analysis [START_REF] Delvaux | The influence of ambient speech on adult speech productions through unintentional imitation[END_REF] ➢ Moving averages (

  Adding other types of conversations, more interesting, more relevant, more promising, more harder ➢ Include additional factors and information for further quantitative and qualitative analysis

	➢ Corpus used: all peer-to-peer conversations in the corpus, ➢ Basic model run using the mgcv package in R (v3.6.1): ➢ Basic model run using the mgcv package in R (v3.6.1): 1) Speech rate (deviation in syll/sec) 2) Vowel duration (in ms, log) 3) F1 and F2 (one model per formant per vowel per pair; ➢ Take a step back and look at the whole dataset together ➢ Extracting the whole corpus with ISCAN
	i.e. 19 pairs or 38 speakers bam(dependent variable ~ fixed or random effects bam(dependent variable ~ fixed or random effects formant values have been Lobanov-normalized)
	+ s(time, by = speaker) + s(time, by = speaker) ➢ bam(speech_rate_dev ~ ➢ bam(phone_logdurms ~ phone_label + follseg + nsyl + ➢ Two options here: ➢
	90 s(phone_begin, by = speaker_name), data = pair) 00 speech_rate_dev + logfreq + ➢ bam(F2.normr ~ Preceding.POA + Following.POA + ➢ Treat pairs as a random factor to compare microvariation
	M ➢ Plotting using itsadug 4 F + 4 M ➢ Plotting using itsadug s(phone_begin, by = speaker_name), data = pair) 6 F + 4 M Target.duration.ms + trajectories by social factors (age, gender);
	Y example for speech rate deviation, pair 00Mf1 4 F + 4 M ➢ A lot of wigglyness within speakers s(phone_begin, by = speaker_name), data = pairvowel) ➢ or actually base the analysis of convergence/divergence 6 F + 6 M plot_smooth plot_diff ➢ There is variation within speakers but not as spiky as in patterns on pairs themselves?
	➢ All pairs have periods of significant difference between speech rate ➢ There is some variation within speakers but very smooth
	➢ Speech rate deviation, vowel duration for all speakers
	➢ Residuals linear plotting monophthongs, vowel quality for FLEECE, BOOT, COT ➢ Most pairs have periods of significant difference between ➢ Models were only run for 3 pairs (=18 models!):
	(Alexander 2018, Chevalier speakers (13 pairs out of 19) no consistence in significant difference periods in formant
	2018) values within pairs

[START_REF] Tamminga | The dynamics of variation in individuals[END_REF]

➢ GAMMs

[START_REF] Sonderegger | The medium-term dynamics of accents on reality television[END_REF][START_REF] Solanki | Brains in dialogue: Investigating accommodation in live conversational speech for both speech and EEG data[END_REF]


M Ackermann

M Ajello

A Allafort

P Schady

L Baldini

J Ballet

G Barbiellini

D Bastieri

R Bellazzini

R D Blandford

E D Bloom

A W Borgland

E Bottacini

A Bouvier

J Bregeon

M Brigida

P Bruel

R Buehler

S Buson

G A Caliandro

R A Cameron

P A Caraveo

E Cavazzuti

C Cecchi

E Charles

R C G Chaves

A Chekhtman

C C Cheung

J Chiang

G Chiaro

S Ciprini

R Claus

J Cohen- Tanugi

J Conrad

S Cutini

F D'ammando

F De Palma

C D Dermer

S W Digel

E Do Couto E Silva

A Domínguez

P S Drell

A Drlica-Wagner

C Favuzzi

S J Fegan

W B Focke

A Franckowiak

Y Fukazawa

S Funk

P Fusco

F Gargano

D Gasparrini

N Gehrels

S Germani

N Giglietto

F Giordano

M Giroletti

T Glanzman

G Godfrey

I A Grenier

J E Grove

S Guiriec

M Gustafsson

D Hadasch

M Hayashida

E Hays

M. S Jackson

T Jogler

J Kataoka

J Knödlseder

M Kuss

J Lande

S Larsson

L Latronico

F Longo

F Loparco

M N Lovellette

P Lubrano

M N Mazziotta

J E Mcenery

J Mehault

P F Michelson

T Mizuno

C Monte

M E Monzani

A Morselli

I V Moskalenko

S Murgia

A Tramacere

E Nuss

J Greiner

M Ohno

T Ohsugi

N Omodei

M Orienti

E Orlando

J F Ormes

D Paneque

J S Perkins

M Pesce-Rollins

F Piron

G Pivato

T A Porter

S Rainò

R Rando

M Razzano

S Razzaque

A Reimer

O Reimer

L C Reyes

S Ritz

A Rau

C Romoli

M Roth

M Sánchez-Conde

D A Sanchez

J D Scargle

C Sgrò

E J Siskind

G Spandre

P Spinelli

Łukasz Stawarz

D J Suson

H Takahashi

T Tanaka

J G Thayer

D J Thompson

L Tibaldo

M Tinivella

D F Torres

G Tosti

E Troja

T L Usher

J Vandenbroucke

V Vasileiou

G Vianello

V Vitale

A P Waite

B L Winer

K S Wood

M Wood

The Imprint of The Extragalactic Background Light in the Gamma-Ray Spectra of Blazars

University, Columbus, OH 43210, USA † majello@slac.stanford.edu, *buehler@stanford.edu, ‡ anita.reimer@uibk.ac.at

The light emitted by stars and accreting compact objects through the history of the Universe is encoded in the intensity of the extragalactic background light (EBL). Knowledge of the EBL is important to understand the nature of star formation and galaxy evolution, but direct measurements of the EBL are limited by Galactic and other foreground emissions. Here we report an absorption feature seen in the combined spectra of a sample of gamma-ray blazars out to a redshift of z∼1.6. This feature is caused by attenuation of gamma rays by the EBL at optical to UV frequencies, and allowed us to measure the EBL flux density in this frequency band.

The bulk of the intergalactic gas in the Universe must have been reionized between the epoch of cosmic recombination, when the Universe was only 300,000 years old (z∼1100), and 1 billion years later (z∼6) as indicated observationally by the spectra of distant quasi-stellar objects [START_REF] Fan | [END_REF]. However, the sources, modes and nature of this cosmic reionization are largely unknown because most of this redshift range has yet to be explored. Photoionization by UV radiation, produced by the first stars and galaxies of the Universe, represents the primary suspect for the ionizing process (2,3). Direct detection of the UV radiation fields is thus of fundamental importance, but at present extremely difficult (3).

An indirect but powerful means of probing the diffuse radiation fields is through γ-γ absorption of high-energy gamma rays (4)(5)(6). In this process, a gamma-ray photon of energy E γ and an EBL photon of energy E EBL annihilate and create an electron-positron pair. This process occurs for head-on collisions when (e.g.) E γ × E EBL ≥ 2(m e c 2 ) 2 , where m e c 2 is the rest mass energy of the electron. This introduces an attenuation in the spectra of gamma-ray sources above a critical gamma-ray energy of E crit (z) ≈ 170(1 + z) -2.38 GeV (7,8).

The detection of the gamma-ray horizon (i.e. the point beyond which the emission of gamma-ray sources is strongly attenuated) is one of the primary scientific drivers of the Fermi Gamma-ray Space Telescope [START_REF] Hartmann | The First GLAST Symposium[END_REF](10)[START_REF] Kashlinsky | The First GLAST Symposium[END_REF]. Several attempts have been made in the past but none detected the long-sought EBL attenuation (12)(13)(14). So far, limits on the EBL density have been inferred from the absence of absorption features in the spectra of individual blazars (13,15), distant galaxies with bright gamma-ray emission powered by matter accreting onto central, massive black holes. While this feature is indeed difficult to constrain for a single source, we show that it is detected collectively in the gamma-ray spectra of a sample of blazars as a cut-off that changes amplitude and energy with redshift. We searched for an attenuation of the spectra of blazars in the 1-500 GeV band using the first 46 months of observations of the Large Area Telescope (LAT) on board the Fermi satellite. At these energies gamma rays are absorbed by EBL photons in the optical to UV range. Thanks to the large energy and redshift coverage, Fermi-LAT measures the intrinsic (i.e. unabsorbed) spectrum up to ∼100 GeV for any blazar at z<0.2, and up to ∼15 GeV for any redshift.

The LAT has detected >1000 blazars to date (16). We restricted our search to a subset of 150 blazars of the BL Lacertae (BL Lac) type that are significantly detected above 3 GeV, because of the expected lack of intrinsic absorption (17). The sample covers a redshift range 0.03-1.6 [START_REF]Agenzia Spaziale Italiana[END_REF]19). The critical energy is therefore always ≥25 GeV, which means that the spectrum measured below this energy is unabsorbed and a true representation of the intrinsic spectrum of the source. We thus determined the intrinsic source spectrum relying on data between 1 GeV and the critical energy E crit and extrapolated it to higher energies. By combining all the spectra we were able to determine, the average deviation, above the critical energy, of the measured spectra from the intrinsic ones, which ultimately provides a measurement of the optical depth τ γγ .

The analysis was performed using the Fermi Science Tools [START_REF]The software to perform analysis of Fermi data is[END_REF]. We determined the spectral parameters of each blazar by maximizing the likelihood of a given source model. The model comprised the Galactic and isotropic diffuse components and all sources in the second Fermi LAT catalog [START_REF] Nolan | [END_REF] within a region of interest (ROI) of 15 • radius. We modeled the spectra of the sources in our sample as parabolic in the logarithmic space of energy and flux (see Eq. 2 in [START_REF] Nolan | [END_REF] for a definition). Their spectra were modified by a term e -τγγ (E,z) that describes the absorption of gamma-ray photons on the EBL. In the above we defined

τ γγ (E, z) = b • τ model γγ (E, z),
where the τ model γγ (E, z) is the optical depth predicted by EBL models (7,[START_REF] Kneiske | [END_REF](23)(24)(25) and b is a scaling variable, left free in the likelihood maximization. In particular, this allowed us to assess the likelihood of two important scenarios: i) there is no EBL attenuation (b=0), ii) the model prediction is correct (b=1).

We combined the data from all the ROIs in a global fit that determined the common parameter b for a given EBL model (see Table S1). All those models with a minimal EBL density based on (or compatible with) resolved galaxy counts (2,7,(24)(25)(26)(27) were found to be acceptable descriptions of the Fermi data (i.e are consistent with b=1 within ≈25 %, see also Figure 1) yielding a significance of the absorption feature of up to ∼6 σ. Models that predict a larger intensity of the EBL particularly in the UV [START_REF] Kneiske | [END_REF]23) would produce a stronger-than-observed attenuation feature and are therefore incompatible with the Fermi observations. Our measurement points to a minimal level of the optical-UV EBL up to redshift z≈1.6 which combined with the upper limits (15,28,29) derived at lower redshift (using observations of blazars at TeV energies) on the near-infrared EBL highlights the conclusion that most of the EBL intensity can be explained by the measured galaxy emission.

Our measurement relies on the accuracy of the extrapolation of the intrinsic spectra of the sources above the critical energy (30). This in turn depends on a precise description of the gamma-ray spectra by our source parametrization. To verify that this is the case and to exclude the possibility that the detected absorption feature is intrinsic to the gamma-ray sources (17), we performed the analysis in 3 independent redshift intervals (z<0.2, 0.2≤z<0.5, and 0.5≤z<1.6). The deviations from the intrinsic spectra in the three redshift intervals are displayed in Figure 2. In the local Universe (z<0.2), EBL absorption is negligible in most of the Fermi-LAT energy band (E crit ≥120 GeV). The lowest redshift interval therefore reveals directly the intrinsic spectra of the sources and shows that our spectral parametrization is accurate [START_REF]Agenzia Spaziale Italiana[END_REF]. The absorption feature is clearly visible above the critical energy in the higher redshift bins. Its amplitude and modulation in energy evolve with redshift as expected for EBL absorption. In principle, the observed attenuation could be due to a spectral cut-off that is intrinsic to the gamma-ray sources. The absence of a cut-off in the spectra of sources with z<0.2 would require that the properties of BL Lacs change with redshift or luminosity. It remains an issue of debate whether such evolution exists (31)(32)(33)(34). However, in case it were present, the intrinsic cut-off would be expected to evolve differently with redshift than we observe. To illustrate this effect, we fitted the blazar sample assuming that all the sources have an exponential cut-off at an energy E 0 . From source to source the observed cut-off energy changes because of the source redshift and because we assumed that blazars as a population are distributed in a sequence such as that proposed in (31)(32)(33)(34). E 0 was fitted to the data globally like b above. As apparent from Figure 2, it appears difficult to reconcile the observed feature with an intrinsic characteristic of the blazars' spectra. We therefore associate the spectral feature to the EBL absorption.

At energies ≤100 GeV, gamma rays observed at Earth and coming from redshift ≥1 interact mostly with UV photons of ≥5 electron volts. An UV background in excess of the light emitted by resolved galaxies can be produced locally by AGN or at higher redshift (z≈7-15) by low-metallicity massive stars (35). By comparing the results from the best-fit EBL models, we measured the UV component of the EBL to have an intensity of 3(±1) nW m -2 sr -1 at z≈1. A contribution to the UV background from AGN as large as the one predicted by (36) (i.e., ≈ 10 nW m -2 sr -1 ) and used in the EBL model of ( 22) is thus excluded by our analysis at high confidence. However, the recent prediction (37) of the UV background from AGN (≈ 2 nW m -2 sr -1 ) is in agreement with the Fermi measurement. Direct measurements of the extragalactic UV background are hampered by the strong dust-scattered Galactic radiation (38). The agreement between the intensity of the UV background as measured with Fermi and that due to galaxies individually resolved by the Hubble Space Telescope (39) (3±1 nW m -2 sr -1 versus 2.9-3.9 nW m -2 sr -1 , respectively) shows that the room for any residual diffuse UV emission is small. This conclusion is reinforced by the good agreement of the Fermi measurement and the estimate of the average UV background, at z≥1.7, of 2.2-4.0 nW m -2 sr -1 using the proximity effect in quasar spectra (40).

Zero-metallicity population-III stars or low-metallicity population-II stars are thought to be the first stars to form in the Universe and formally marked the end of the dark ages when, with their UV light, these objects started ionizing the intergalactic medium (41). These stars, whose mass might have exceeded one hundred times the mass of our Sun, are also believed to be responsible for creating the first metals and dispersing them in the intergalactic medium (42)(43)(44). A very large contribution of population-III stars to the near-infrared EBL had already been excluded by (15). Our measurement constrains, according to (45,46), the redshift of maximum formation of low-metallicity stars to be at z≥10 and its peak co-moving star-formation rate to be lower than 0.5 M ⊙ Mpc -3 yr -1 . This upper limit is already of the same order of the peak star-formation rate of 0.2-0.6 M ⊙ Mpc -3 yr -1 proposed by (47) and suggests that the peak star-formation rate might be much lower as proposed by (48).

Energy [GeV]

≈ z

Figure 1 Measurement, at the 68 % and 95 % confidence levels (including systematic uncertainties added in quadrature), of the opacity τ γγ from the best fits to the Fermi data compared to predictions of EBL models. The plot shows the measurement at z≈1 which is the average redshift of the most constraining redshift interval (i.e. 0.5≤z<1.6). The Fermi-LAT measurement was derived combining the limits on the best-fit EBL models. The downward arrow represents the 95 % upper limit on the opacity at z=1.05 derived in (13). For clarity this figure shows only a selection of the models we tested while the full list is reported in Table S1. The EBL models of (49), which are not defined for E≥250/(1 + z) GeV and thus could not be used, are reported here for completeness. Figure 2 Absorption feature present in the spectra of BL Lacertae objects as a function of increasing redshift (data points, from top to bottom). The dashed curves show the attenuation expected for the sample of sources by averaging, in each redshift and energy bin, the opacities of the sample (the model of (7) was used) and multiplying this average by the best-fit scaling parameter b obtained independently in each redshift interval. The vertical line shows the critical energy E crit below which ≤5 % of the source photons are absorbed by the EBL. The thin solid curve represents the best-fit model assuming that all the sources have an intrinsic exponential cut-off and that blazars follow the blazar sequence model of (32,33).

Supplements Source Selection

The second LAT AGN Catalog ( 16) contains many flat-spectrum radio quasars (FSRQs) and BL Lacertae (BL Lac) objects whose spectra extend significantly above 10 GeV. The classification relies on the conventional definition of BL Lac objects outlined in (50)(51)(52) in which the equivalent width of the strongest optical emission line is <5 Å and the optical spectrum shows a Ca II H/K break ratio C<0.4. Sources are then classified also according to the position of the peak of the synchrotron component as low-synchrotron-peaked (LSP, ν peak <10 14 Hz), intermediatesynchrotron-peaked (ISP, 10 14 <ν peak <10 15 Hz), and high-synchrotron-peaked (HSP, ν peak >10 15 Hz).

FSRQs generally have soft spectra (photon index greater than 2.3) at GeV energies, making it more difficult to detect absorption features at energies greater than 10 GeV. Additionally their spectra might suffer from intrinsic absorption of gamma rays by the photons of the broad line region or of the accretion disk (17). This makes them non-ideal candidates to constrain the EBL. We therefore focused on the BL Lac objects. For computational reasons, only 50 sources can be analyzed in one combined likelihood fit. In order to perform our analysis in three redshift intervals we selected the 150 BL Lacs that show the largest detection significance in the 3-10 GeV energy band (the significance is always larger than 3.5 σ) with no selection on spectral shape. Most of the sources (100 out of 150) were detected with a significance ≥3 σ above 10 GeV already after two years of observations [START_REF] Nolan | [END_REF]. In the second LAT AGN Catalog catalog (16) only ∼190 BL Lac have a redshift measurement, so our sample is a representative set of Fermi-detected BL Lacs. Figure S1 shows the redshift distribution of the BL Lac objects employed in this analysis.

Data Analysis

Each source in our sample is analyzed using 46 months 1 of Fermi observations using version v9r27 of the Science Tools 2 . The data were filtered, removing time periods in which the instrument was not in sky-survey mode, and removing photons whose zenith angle is larger than 100 • to limit contamination from the Earth limb emission. We consider only photons collected within 15 • of the source position with 1≤E≤ 500 GeV. We employ the P7SOURCE V6 3 instrumental response function (IRF) and perform a binned likelihood analysis. The Galactic and isotropic diffuse emissions are modeled using respectively the gal 2yearp7v6 v0.fits and iso p7v6source.txt templates.

We rely on the Composite Likelihood tool of the Fermi software to perform the likelihood maximization. This allows us to fit simultaneously the data from different ROIs with the aim of constraining the scaling parameter b that is applied to the opacity curves predicted by each 1 From August 4th 2008 to June 1st 2012 2 http://fermi.gsfc.nasa.gov/ssc/data/analysis/software/ 3 Our analysis is robust against change of the dataset and IRF. one of the EBL models we tested. The total number of free parameters across all ROIs in one redshift bin is approximately 1000. It is computationally unfeasible to fit all of these parameters simultaneously. We therefore proceeded in three steps. First we fit each ROI individually. The fit is performed on the entire energy band (1-500 GeV). The parameters (i.e. flux and photon index) of all the sources within 4 • of the target source, along with the parameters of the diffuse components, are left free to vary. More distant sources have parameters frozen at the values measured in the second Fermi LAT catalog [START_REF] Nolan | [END_REF], unless the inspection of the residual map showed that a given source underwent strong variability. In that case, the normalizations of those sources were left free in the fit. The spectra of the sources of interest are modeled using a LogParabola model [START_REF] Nolan | [END_REF]. We then proceed to re-optimize the parameters of the source of interest in each ROI up to that energy for which the EBL absorption becomes no longer negligible (see the definition of E crit in the main text). In the third step we fix the curvature of the LogParabola model to the value obtained in the previous fit. Finally, we use the Composite Likelihood to constrain the gamma-ray opacity. The Composite Likelihood allows the user to tie any parameters between any ROIs. In our case the only tied parameter is the renormalization factor b of the opacity of a given EBL model. The significance of our finding can be evaluated using the Test Statistic (TS) as T S = 2(logL (b) -logL (b = 0)), where L is the likelihood function and the b = 0 case (i.e. no EBL absorption) represents the null hypothesis. Since the null hypothesis is a special case of the hypothesis we test, we expect the T S to be distributed as a χ 2 with one degree of freedom (see also next sections). This allows one to transform the TS into the corresponding number of standard deviations of a Gaussian distribution as n σ = √ T S.

Redshift

Modeling the Intrinsic Blazar Spectrum

One of the basic assumptions of this analysis is that on average the spectrum of a BL Lac blazar can be adequately modeled using a LogParabola in the 1-500 GeV band. Figure S2 shows the residuals of the best fits to all the sources in the z<0.2 interval. It is apparent that the LogParabola provides a good representation of the spectra of blazars in our sample.

In an additional test we artificially decreased the critical energy E crit for all the sources in the z<0.2 interval from the typical ≥120 GeV (for z<0.2) to ∼40 GeV which is representative of the z≥0.5 case. Even in this case the results are unchanged, showing that the properties of the intrinsic spectrum can be determined over a more restricted energy range and can be safely extrapolated above the critical energy.

Blazars are however known to have complex spectra that cannot always be modeled using simple empirical functions. To further test our basic assumption, we fit the spectra of bright blazars with coverage in the GeV and TeV band. We rely on the published spectra of Mrk 421 (53), Mrk 501 (54) and RBS 0413 (55). We find that in the energy range relevant for the present analysis the LogParabola model provides an adequate fit to all these spectra. A further confirmation comes from our validation studies using simulations of broad-band blazar spectra (see next sections). We thus conclude that in a small energy range like the one adopted here (1-500 GeV), the LogParabola is an adequate representation of the intrinsic blazar spectrum.

Results

Table S1 reports the results of the joint-likelihood fit for all the predictions of the EBL models we tested. In many cases, a given EBL model (because of tunable parameters or uncertainties connected to e.g, galaxy evolution, star-formation rate etc.) comprises several different predictions for the opacity. In this case, we use these prediction at face value and test them against our observations. All the models with the exclusion of ( 22) 4 rely on a standard concordance cosmology (H 0 =70 km s -1 Mpc -1 , Ω M =1-Ω Λ =0.3). For each model we report the significance of the absorption feature detected once we allow the opacity to be renormalized by b. The significances have been obtained as a square root of the sum of the TS in the 3 redshift intervals while the b parameter as reported in Table S1 represents the weighted average5 of the 3 b values determined in the 3 redshift bins. For example, for the model of ( 7) we measure: b=1.18 +0.94 -0.81 (TS≈4 for z<0.2), b=0.82 +0. 41 -0.30 (TS≈7 for 0.2≤z<0.5), and b=1.29 +0.43 -0.36 (TS≈25 for 0.5≤z<1.6). The weighted average yields (as reported in Table S1) b = 1.02 ± 0.23 and TS≈36.

We also quantify the level of compatibility of the prediction of a given model with the Fermi data, by comparing the likelihood of the original model (b = 1) to the one of the best fit scenario (b free to vary). For b values significantly different than 1 a given EBL model predicts an attenuation larger than observed. Table S1 shows that the intensity of the optical-UV background in 4 of the 15 models we tested is not compatible with the Fermi data at ≥ 3 σ.

Figure S3 shows that the significance of the absorption feature is detected collectively6 among the BL Lacs considered and is not attributable to just a few sources. Indeed, it is apparent that the TS increases linearly with the number of sources as expected in a background-limited scenario.

As shown already in Figure 1, our analysis probes absorption signatures that are a factor ∼4 times weaker than those probed in (13) which relied on only one year of data. Even a conservative 95 % upper limit derived from this analysis (b U L = 1.5 referred to the model of (7), i.e. τ γγ <0.6 at z≈1 and E=77 GeV) is a factor ∼2 below the one reported by (13).

Our analysis relies on the assumption that BL Lac spectra (and HSP in particular) do not change dramatically across the z=0.2 redshift barrier. In the next section we address possible sources of systematic effects. a Models tested are implemented in the Fermi Science Tools. As an example the recent model of (49) which is not defined for E≥250/(1 + z) GeV could not be used, but its predictions are for completeness reported in Figure 1.

b Reference number in the 'References and Notes' section of the main text. c Significance, in units of σ, of the attenuation in the spectra of blazars when a given EBL model is scaled by the factor b. In this case b=0 (i.e. no EBL absorption) constitutes the null hypothesis.

d This column lists the maximum likelihood values and 1 σ confidence ranges for the opacity scaling factor.

e Here the b=1 case (i.e. EBL absorption as predicted by a given EBL model) constitutes the null hypothesis. This column shows the compatibility (expressed in units of σ) of the predictions of EBL models with the Fermi observations. Large values mean less likely to be compatible. 7) produced in the joint-likelihood fit (to the 0.5≤z<1.6 interval) while adding one source at a time. The sources have been sorted in redshift (from lowest to highest). The dashed line shows the best-fit linear increase of the TS with the number of sources. The inset shows the best-fit value of the renormalization parameter b applied to the opacity predicted by (7) (see text for details).

Validation of the Analysis using Simulations

We validated our analysis using simulated datasets starting from physically motivated spectral energy distributions (SEDs) of BL Lacs. We rely on SEDs that reproduce well the range of peak curvatures, peak frequencies (for both the synchrotron and high-energy component), and gamma-ray photon indices observed for the LAT-detected BL Lacs. These SEDs have been produced using the same numerical code presented in (57) and take into account all the important effects that contribute to determine the intrinsic curvature of the gamma-ray spectrum of the Fermi blazars. Of particular importance are those effects that depend on the curvature of the electron distribution, as well as those related to the Thomson (TH) to Klein-Nishina (KN) transition, in the inverse Compton process (IC). The latter effect is crucial since it is well known (57,58) that the TH/KN transition in the IC cross section can result in a steepening of the high-energy branch of the IC SED, compared to the low-energy one. Our simulations fully take into account this effect, allowing us to study the potential bias in the EBL estimate. We performed ∼500 simulations randomly selecting 50 SEDs from our template library. To each SED we assign a redshift (in the 0.5≤z<1.6 range) and a 0.1-100 GeV flux randomly extracted from those of BL Lacs in the second LAT AGN Catalog (16). We then proceeded to simulate the data expected from those sources based on Fermi's instrument response function and pointing history. Each one of those realizations was processed with the same analysis chain used on actual data.

Figure S4 presents the distribution of TS (left panel) and best-fit b parameters (right panel) for the simulated blazar population under the following two scenarios. First, no EBL attenuation was included, in order to check whether our analysis would pick a signal consistent with EBL attenuation when this was not present (i.e. a false positive). The dashed lines in the panels show that our analysis yields very small TS values and b values compatible with zero, as expected for a robust analysis.

In the second scenario, the SED models were attenuated using the opacities of the model of (7). The results summarized in Figure S4 show that TS values ≥20 and b values compatible with 1 (within the statistical uncertainty) are correctly retrieved. Moreover, the simulations show that our analysis is free of any major systematic uncertainty since the average values of b retrieved for both scenarios are compatible with the expected ones (0 and 1 respectively) within 1 %. Figure S4 Results from joint-likelihood fits to the simulated datasets. The left panel shows the TS of the detection of the attenuation produced by the EBL for the case in which no EBL was applied to the SEDs (dashed line) and the case for which an attenuation consistent with the one of (7) was used (solid line). The right panel shows the best-fit b parameter for the same two scenarios.

Systematic Uncertainties of the Instrument Response Function

In order to gauge the systematic uncertainties of this analysis we use the IRF bracketing method as described in (59). By deriving two different sets of IRFs and repeating the entire analysis we find that the systematic uncertainty on the opacity τ γγ is of the order ∼7 %. Moreover, the entire analysis has been repeated using the P7CLEAN V6 photon selection and IRF. The results presented in the above sections are fully confirmed and the systematic uncertainty on the opacity τ γγ due to changing event selection and IRF is of the order 3 %. Thus, we consider the total systematic uncertainty on τ γγ to be ∼10 %.

Source Intrinsic Effects

The blazar evolution with luminosity as described by the blazar sequence (31)(32)(33)(34) and the change of blazar type with redshift might in principle produce an intrinsic cut-off that changes with redshift. An intrinsic spectral cut-off might be detected if the part of the blazar spectrum sampled by Fermi corresponds to the tail of the electron distribution. In our case, this might happen only for the LSP sources and it has already been shown (see previous section) that their contribution to the total signal is negligible. On the other hand most of the signal is dominated by HSPs that are detected by Fermi right below or at the peak of the inverse Compton component, thus excluding the possibility that their emission is produced in the tail of the electron distribution.

However, in order to study the compatibility of the detected signal with an intrinsic origin, we assumed that all blazars have an exponential cut-off at a source-frame energy E 0 . We also assumed that the blazar spectrum moves to lower frequency for increasing source luminosity as dictated by the so-called blazar sequence (31)(32)(33)(34). This might represent the case in which the maximum electron energy depends on the jet power or the luminosity. Figure 2 shows that this model (E 0 is fitted to the data) provides a poor representation of the Fermi data. Even assuming that the maximum electron energy depends on a variable power of the luminosity does not improve the results. The result is even worse if the blazar sequence is not invoked. We thus conclude that the impact of source intrinsic effects on our analysis is likely minor.

Reprocessed Emission

The electron-positron pairs generated in the interaction of gamma rays and EBL photons can initiate a cascade by subsequent Compton scattering of the photons of the Cosmic Microwave Background. Typically the cascades originated by TeV gamma rays are reprocessed in the GeV energy range. In the case of a weak intergalactic magnetic field (IGMF), the pairs are not deflected out of the line of sight (60) and the reprocessed emission can be a substantial fraction of (or even dominate) the total source signal in the 100 GeV band. The intensity of the reprocessed emission depends primarily on the EBL density, on the intensity of the IGMF, on its coherent length and on the position of the peak of the high-energy component of a blazar (61). For typical coherent lengths of ∼1 Mpc and SED peaks located at or below 10 TeV, lower limits of B 10 -15 G have been inferred (61,62). For these intensities of the IGMF the reprocessed emission is expected to be subdominant with respect to the intrinsic component (63). While it is true that these estimates are based on an implied activity of the source of a few million years (64), the majority of the Fermi BL Lacs is expected to have the peak of the high-energy component located at ≤1 TeV (65) greatly reducing the amount of reprocessed emission in the GeV band (61). Similar considerations (i.e. suppression of the reprocessed emission by the IGMF) hold also for the case in which blazars are also sources of ultrahigh energy cosmic rays (66)(67)(68).

62 .

 62 Consorzio Interuniversitario per la Fisica Spaziale (CIFS), I-10133 Torino, Italy 63. Dipartimento di Fisica, Università di Roma "Tor Vergata", I-00133 Roma, Italy 64. Department of Physics, Center for Cosmology and Astro-Particle Physics, The Ohio State

  LAT best fit --2 sigma Franceschini et al. 2008 Finke et al. 2010 --model C Stecker et al. 2012 --High Opacity Stecker et al. 2012 --Low Opacity Kneiske et al. 2004 --highUV Kneiske et al. 2004 --best fit Kneiske & Dole 2010 Dominguez et al. 2011 Gilmore et al. 2012 --fiducial Abdo et al. 2010

Figure

  Figure S1Redshift distribution of the BL Lac objects used in this analysis.

Figure

  Figure S2Average of the residuals with respect to the best fit models for all the blazars in the z<0.2 redshift bin.

Figure

  Figure S3Increase in the TS value of the (renormalized) EBL model of(7) produced in the joint-likelihood fit (to the 0.5≤z<1.6 interval) while adding one source at a time. The sources have been sorted in redshift (from lowest to highest). The dashed line shows the best-fit linear increase of the TS with the number of sources. The inset shows the best-fit value of the renormalization parameter b applied to the opacity predicted by(7) (see text for details).

Table S1 .

 S1 Joint-likelihood results for different EBL models.

	Model a	Ref. b Significance of b=0 Rejection c	b d	Significance of b=1 Rejection e
	Stecker et al. (2006) -fast evolution (23)	4.6	0.10±0.02	17.1
	Stecker et al. (2006) -baseline	(23)	4.6	0.12±0.03	15.1
	Kneiske et al. (2004) -high UV	(22)	5.1	0.37±0.08	5.9
	Kneiske et al. (2004) -best fit	(22)	5.8	0.53±0.12	3.2
	Gilmore et al. (2012) -fiducial	(27)	5.6	0.67±0.14	1.9
	Primack et al. (2005)	(56)	5.5	0.77±0.15	1.2
	Dominguez et al. (2011)	(25)	5.9	1.02±0.23	1.1
	Finke et al. (2010) -model C	(24)	5.8	0.86±0.23	1.0
	Franceschini et al. (2008)	(7)	5.9	1.02±0.23	0.9
	Gilmore et al. (2012) -fixed	(27)	5.8	1.02±0.22	0.7
	Kneiske & Dole (2010)	(26)	5.7	0.90±0.19	0.6
	Gilmore et al. (2009) -fiducial	(2)	5.8	0.99±0.22	0.6

The model of[START_REF] Kneiske | [END_REF] adopts a value of H 0 =65 km s -1 Mpc -1 . We checked that neglecting the change of H 0 for the model of (22) introduces an error on the opacity of <7 % which has negligible influence on our analysis..

The weights are 1/σ 2 b where σ b is the uncertainty on the b parameters.

In FigureS3the sources have been sorted in redshift from low to high redshift.

The sample used in the previous sections contains for z<0.2 35 HSPs, 10 ISPs and 5 LSPs. The 0.2≤z<0.5 sample comprises 27 HSPs, 18 ISPs, 5 LSPs, while the z≥0.5 sample contains 10 HSPs, 19 ISPs, and 21 LSPs.

In the z<0.2 bin there are respectively 41 HSPs, 10 ISPs and 5 LSPs; while in the z≥0.2 there are 50 HSPs, 30 ISPs, and 23 LSPs.

Marco Ajello acknowledges generous support from the Fermi guest investigator program (proposals ID 31117 and 51258) and the Swift and the GROND teams for observing ∼hundred Fermi blazars in an effort to constrain their redshifts. The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique / Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d' Études Spatiales in France.

by a grant from the

Additional Tests

A number of additional tests have been performed to test the stability of the results presented in the previous sections. These are documented below.

Sample Selection

The selection criteria used in the previous sections were chosen a priori to yield a uniform source population that could be used to probe in a sensitive and least-biased way the EBLinduced absorption feature. Several tests have been performed to check this assumption.

The BL Lac population is known to evolve in redshift with HSP objects detected predominantly at low redshift and LSPs detected predominantly at larger redshifts with ISPs bridging the gap 7 . This could artificially induce a decline in the combined (produced by all subtypes) BL Lac spectrum at an energy that decreases with increasing redshift, and hence could mimic the effect of EBL-caused absorption. To test whether such bias affects our results we run the analysis separately for the three source classes (HSPs, ISPs and LSPs). For each of them we define two redshift bins: z<0.2 and z≥0.2. In the lowest redshift bin, due to the low statistics 8 the joint likelihood yields no significant detection of EBL softening. For the HSPs in the low redshift bin, the joint likelihood yields b=1.46 +0.89 -0.78 and a TS=4.1. At higher redshift, the results are: b =1.04 +0. 32 -0.28 (TS=24), b =2.13 +1.11 -0.98 (TS=9), and b =0.43 +0.67 -0.71 (TS=1), for HSPs, ISPs and LSPs respectively. The weighted average of the above values yields b =1.04 +0. 27 -0.24 and a total TS of ∼38. These results are in very good agreement with the ones presented in the previous sections and show that high-redshift LSPs are not responsible for the observed spectral softening, and HSPs are mainly responsible for it. Moreover, since the signal is dominated by only one source population (namely HSPs) the probability that selection effects from changing the sample composition with redshift strongly affect our result is rather low.

We also checked if the highest redshift HSPs are driving the result. To this extent we isolated (from the above sample) the 6 HSPs with the highest redshift. The joint likelihood yields a TS=6 and b =1.99 +1.83 -1.10 indicating that the detection of the EBL cut-off is not caused by a few high redshift sources.

We also tested the population of BL Lacs that did not pass our selection criteria. This comprises 32 objects detected with a rather low significance in the 3-10 GeV over a redshift range 0.03-1.9 (average of 0.4). The joint likelihood yields b=1.93 +2.73 -1.39 and TS=3.5 compatible with our result. Thus, including low significance sources marginally improves the results of our analysis.

Binning in Redshift

We bin the sources in redshift to cope with the maximum number of free parameters (i.e. 100) that can be optimized in a single maximum likelihood fit. Our analysis is however independent of the redshift binning used. This is already evident from the tests reported above using separately HSPs, ISPs and LSPs and a different redshift binning (with respect to the one used in the main analysis). However, we performed an additional test and divided the most constraining redshift bin (i.e. 0.5≤z<1.6) into halves (at z≈0.75) containing 25 sources each. When fitting the Fermi data using the model of ( 7), we obtain b=1.71 +1.12 -0.82 (TS=10) and b=1.17 and TS=25.1 obtained for the full 0.5≤z<1.6 redshift bin. Our analysis is thus robust against choices of the redshift bins.

Source Variability

Blazars are inherently variable objects with variability in flux of up to a factor 10 or more. Throughout this work we use spectra of blazars accumulated over 46 months of observations and blazar variability might in principle bias our result. In the Fermi sample, the most variable BL Lac objects are those belonging to the LSP and ISP class and it has been shown already that the entire analysis can be confirmed without using those two classes. To test possible biases deriving from variability we selected the 30 most variable, according to (16,[START_REF] Nolan | [END_REF], and most significant (in the 3-10 GeV band) BL Lac objects. Repeating the entire analysis using this sample, we obtain b = 1.20 +0. 62 -0.52 and a TS=9. These results are in good agreement with the ones presented in the previous section and show that blazar variability does not affect our results.

Critical Energy

Our measurement relies on the extrapolation of the intrinsic (i.e. unabsorbed) spectrum above the critical energy E crit . Aggressive choices of this energy (i.e. large values) might bias the result towards a low level of opacity. We show that this is not the case for our measurement and that our analysis is robust against conservative choices of E crit . In order to demonstrate it, we imposed E crit =10 GeV for any source in the most constraining redshift interval: i.e. 0.5≤z<1.6. For most EBL models, with the exception of those presented in (23), the opacity is negligible at 10 GeV up to z≈1.6. When fitting the Fermi data using the model of (7), our analysis yields a value of b=1.24 +0. 44 -0.37 and TS≈22 for the constant E crit =10 GeV case. This can be compared to b=1.29 +0. 43 -0.36 and TS≈25.1 for the variable E crit case presented in the main text. We thus conclude that our analysis is robust against any conservative choice of E crit .