Astronomy & Astrophysics manuscript no. MAGIC_1ES1011_EBL_Feb2016_toArX c©ESO 2016 February 19, 2016 MAGIC observations of the February 2014 flare of 1ES 1011+496 and ensuing constraint of the EBL density M. L. Ahnen1, S. Ansoldi2, L. A. Antonelli3, P. Antoranz4, A. Babic5, B. Banerjee6, P. Bangale7, ∗, U. Barres de Almeida7,26, J. A. Barrio8, J. Becerra González9,27, W. Bednarek10, E. Bernardini11,28, B. Biasuzzi2, A. Biland1, O. Blanch12, S. Bonnefoy8, G. Bonnoli3, F. Borracci7, T. Bretz13,29, E. Carmona14, A. Carosi3, A. Chatterjee6, R. Clavero9, P. Colin7, E. Colombo9, J. L. Contreras8, J. Cortina12, S. Covino3, P. Da Vela4, F. Dazzi7, A. De Angelis15, B. De Lotto2, E. de Oña Wilhelmi16, C. Delgado Mendez14, F. Di Pierro3, D. Dominis Prester5, D. Dorner13, M. Doro7,30, S. Einecke17, D. Eisenacher Glawion13, D. Elsaesser13, A. Fernández-Barral12, D. Fidalgo8, M. V. Fonseca8, L. Font18, K. Frantzen17, C. Fruck7, D. Galindo19, R. J. García López9, M. Garczarczyk11, D. Garrido Terrats18, M. Gaug18, P. Giammaria3, N. Godinović5, A. González Muñoz12, 35, ∗, D. Guberman12, A. Hahn7, Y. Hanabata20, M. Hayashida20, J. Herrera9, J. Hose7, D. Hrupec5, G. Hughes1, W. Idec10, K. Kodani20, Y. Konno20, H. Kubo20, J. Kushida20, A. La Barbera3, D. Lelas5, E. Lindfors21, S. Lombardi3, F. Longo2, M. López8, R. López-Coto12, A. López-Oramas12,31, E. Lorenz7, P. Majumdar6, M. Makariev22, K. Mallot11, G. Maneva22, M. Manganaro9, K. Mannheim13, L. Maraschi3, B. Marcote19, M. Mariotti15, M. Martínez12, D. Mazin7, U. Menzel7, J. M. Miranda4, R. Mirzoyan7, A. Moralejo12, ∗, E. Moretti7, D. Nakajima20, V. Neustroev21, A. Niedzwiecki10, M. Nievas Rosillo8, K. Nilsson21,32, K. Nishijima20, K. Noda7, R. Orito20, A. Overkemping17, S. Paiano15, J. Palacio12, M. Palatiello2, D. Paneque7, R. Paoletti4, J. M. Paredes19, X. Paredes-Fortuny19, M. Persic2,33, J. Poutanen21, P. G. Prada Moroni23, E. Prandini1,34, I. Puljak5, W. Rhode17, M. Ribó19, J. Rico12, J. Rodriguez Garcia7, T. Saito20, K. Satalecka8, C. Schultz15, T. Schweizer7, S. N. Shore23, A. Sillanpää21, J. Sitarek10, I. Snidaric5, D. Sobczynska10, A. Stamerra3, T. Steinbring13, M. Strzys7, L. Takalo21, H. Takami20, F. Tavecchio3, P. Temnikov22, T. Terzić5, D. Tescaro9, M. Teshima7, 20, J. Thaele17, D. F. Torres24, T. Toyama7, A. Treves25, V. Verguilov22, I. Vovk7, J. E. Ward12, M. Will9, M. H. Wu16, and R. Zanin19 (Affiliations can be found after the references) Received ; accepted ABSTRACT Context. In February-March 2014, the MAGIC telescopes observed the high-frequency peaked BL Lac 1ES 1011+496 (z=0.212) in flaring state at very-high energy (VHE, E>100GeV). The flux reached a level more than 10 times higher than any previously recorded flaring state of the source. Aims. Description of the characteristics of the flare presenting the light curve and the spectral parameters of the night-wise spectra and the average spectrum of the whole period. From these data we aim at detecting the imprint of the Extragalactic Background Light (EBL) in the VHE spectrum of the source, in order to constrain its intensity in the optical band. Methods. We analyzed the gamma-ray data from the MAGIC telescopes using the standard MAGIC software for the production of the light curve and the spectra. For the constraining of the EBL we implement the method developed by the H.E.S.S. collaboration in which the intrinsic energy spectrum of the source is modeled with a simple function (≤ 4 parameters), and the EBL-induced optical depth is calculated using a template EBL model. The likelihood of the observed spectrum is then maximized, including a normalization factor for the EBL opacity among the free parameters. Results. The collected data allowed us to describe the flux changes night by night and also to produce differential energy spectra for all nights of the observed period. The estimated intrinsic spectra of all the nights could be fitted by power-law functions. Evaluating the changes in the fit parameters we conclude that the spectral shape for most of the nights were compatible, regardless of the flux level, which enabled us to produce an average spectrum from which the EBL imprint could be constrained. The likelihood ratio test shows that the model with an EBL density 1.07 (-0.20,+0.24)stat+sys, relative to the one in the tested EBL template (Domínguez et al. 2011), is preferred at the 4.6 σ level to the no-EBL hypothesis, with the assumption that the intrinsic source spectrum can be modeled as a log-parabola. This would translate into a constraint of the EBL density in the wavelength range [0.24 µm,4.25 µm], with a peak value at 1.4 µm of λFλ = 12.27+2.75 −2.29 nW m−2 sr−1, including systematics. Key words. gamma rays – cosmic background radiation – BL Lacertae objects Use \titlerunning to supply a shorter title and/or \authorrunning to supply a shorter list of authors. Article number, page 1 of 8 ar X iv :1 60 2. 05 23 9v 2 [ as tr o- ph .H E ] 1 8 Fe b 20 16 A&A proofs: manuscript no. MAGIC_1ES1011_EBL_Feb2016_toArX 1. Introduction 1ES 1011+496 (RA: 10h 15m 04.1s, DEC: +49◦ 26m 01s) is an active galactic nucleus (AGN) classified as a high-frequency peaked BL Lac (HBL), located at redshift =0.212 (Albert et al. 2007). HBLs have spectral energy distributions (SED) charac- terized by two peaks, one located in the UV to soft X-ray band and the second located in the GeV to TeV range, which makes it possible to detect them in very-high-energy (VHE, E>100GeV) γ rays. 1ES 1011+496 was discovered at VHE by the MAGIC Collaboration in 2007 following an optical high state reported by the Tuorla Blazar Monitoring Programme (Albert et al. 2007). The observation of a bright source at intermediate redshift, like 1ES 1011+496, provides a good opportunity to constrain the impact of the Extragalactic Background Light (EBL) on the propagation of γ rays over cosmological distances. The EBL is the diffuse radiation that comes from the contributions of all the light emitted by stars in the UV-optical and near infrared (NIR) bands. It also contains the infrared (IR) radiation emitted by dust after absorbing the starlight, plus a small contribution from AGNs (Hauser & Dwek 2001). VHE γ rays from extragalac- tic sources interact with the EBL in the optical and NIR bands, producing electron-positron pairs, which causes an attenuation of the VHE photon flux measured at Earth (Gould & Schréder 1967). Measuring directly the EBL is a challenging task due to the intense foreground light from interplanetary dust. For the opti- cal band strict lower limits to the EBL have been derived from galaxy counts (Madau & Pozzetti 2000; Fazio et al. 2004; Dole et al. 2006). At NIR wavelengths, one way to access the EBL is by large-scale anisotropy measurements (e.g. Cooray et al. 2004; Fernandez et al. 2010; Zemcov et al. 2014). Making rea- sonable assumptions on the intrinsic VHE spectra of extragalac- tic sources (e.g. the limit in the hardness of the photon spec- tra of 1.5, coming from theoretical limits in the acceleration mechanisms), upper limits to the EBL density can be derived (e.g. Stecker & de Jager 1996; Aharonian et al. 2006; Mazin & Raue 2007). More recently, extrapolations of data from the Fermi Large Area Telescope have been used to set constraints to the intrinsic VHE spectra of distant sources, which, in combina- tion with Imaging Atmospheric Cherenkov Telescopes (IACT) observations of the same objects, have also provided upper lim- its to the EBL density (Georganopoulos et al. 2010; Orr et al. 2011; Meyer et al. 2012). The Fermi collaboration employed a different technique to constrain the EBL density using a likelihood ratio test on LAT data from a number of extragalactic sources (Ackermann et al. 2012). SEDs from 150 BL Lacs in the redshift range 0.03 - 1.6 were modeled as log parabolae in the optically-thin regime (E < 25 GeV), then extrapolated to higher energies and compared with the actually observed photon fluxes. A likelihood ratio test was used to determine the best-fit scaling factor for the optical depth τ(E, z) according to a given EBL model, hence provid- ing constraints of the EBL density relative to the model predic- tion. Several EBL models were tested using this technique (e.g. Stecker et al. 2006; Finke et al. 2010), including the most widely and recently used by IACTs by Franceschini et al. (2008) and Domínguez et al. (2011). They obtained a measurement of the UV component of the EBL of 3 ± 1 nW m−2 sr−1 at z ≈ 1. The H.E.S.S. collaboration used a similar likelihood ratio test to constrain the EBL taking advantage of their observations of distant sources at VHE (Abramowski et al. 2013). The EBL absorption at VHE is expected to leave an imprint in the ob- served spectra, coming from a distinctive feature (an inflection point in the log flux–log E representation) between ∼ 100 GeV and ∼5-10 TeV, a region observable by IACTs. This feature is due to a peak in the optical region of the EBL flux density, which is powered mainly by starlight. The H.E.S.S collaboration mod- eled the intrinsic spectra of several AGNs using simple func- tions (up to 4 parameters), then applied a flux suppression factor exp(−α × τ(E, z)), where τ is the optical depth according to a given EBL model and α a scaling factor. A scan over α was per- formed to achieve the best fit to the observed VHE spectra. The no-EBL hypothesis, α = 0, was excluded at the 8.8 σ level, and the EBL flux density was constrained in the wavelength range between 0.30 µm and 17 µm (optical to NIR) with a peak value of 15 ± 2stat ± 3sys nW m−2 sr−1 at 1.4 µm. In Domínguez et al. (2013), data from 1ES 1011+496 was used as part of a data set from several AGNs to measure the cosmic γ-ray horizon (CGRH). The CGRH is defined as the energy at which the optical depth of the photon-photon pair production becomes unity as function of energy. Using multi- wavelenght (MWL) data, Domínguez et al. modeled the SED of each source, including 1ES 1011+496, doing a extrapolation to the VHE band. Then they made a comparison with the observed VHE data. In the case of 1ES 1011+496, they modeled the SED using the optical data from 2007 (Albert et al. 2007) and X-ray data (from the X-Ray Timing Explorer) taken in 2008 May, and compared it with the VHE data taken in 2007 by MAGIC. Their prediction was below the observed VHE data, which led to no optical-depth information. The prediction may have failed due to the lack of simultaneity in the data. A similar approach was presented by Mankuzhiyil et al. (2010), where they modeled the SED of PKS 2155-304 making a prediction for the VHE band and compared it to the observed data to give attenuation limits. After the discovery of 1ES 1011+496 in 2007 (Albert et al. 2007), two more multi-wavelength campaigns have been organ- ised by MAGIC: the first one between 2008 March and May (Ah- nen et al. 2016a) and a second one divided in two periods, from 2011 March to April and 2012 from January to May (Ahnen et al. 2016b). In all previous observations (including the discov- ery) the source did not show evidence of flux variability within the observed periods and the observed spectra could be fitted with simple power-law functions, with photon indices ranging between 3.2 ± 0.4stat and 4.0 ± 0.5stat, and integral fluxes, above 200 GeV, between (0.8±0.1stat)×10−10 and (1.6±0.3stat)×10−11 photons cm−2s−1. In this paper we present the analysis of the extraordinary flare of 1ES 1011+496 in 2014 February-March observed by the MAGIC telescopes, and apply a technique based on Abramowski et al. (2013) for constraining the EBL. The observations and the data reduction are described in Sect. 2, the results in Sect. 3, the procedure for constraining the EBL in Sect. 4, the inclusion of the systematic uncertainty is shown in Sect. 5, and the results of the EBL constraint are discussed in Sect. 6. 2. Observations & Analysis MAGIC is a stereoscopic system of two 17 m diameter Imag- ing Atmospheric Cherenkov Telescopes (IACT) situated at the Roque de los Muchachos, on the Canary island of La Palma (28.75◦N, 17.86◦W) at a height of 2200 m above sea level. Since the end of 2009, it has been working in stereoscopic mode with a trigger threshold of ∼50 GeV. During 2011 and 2012, MAGIC underwent a series of upgrades which results in a sensitivity of (0.66 ± 0.03)% of the Crab nebula flux above 220 GeV in 50 hours at low zenith angles (Aleksić et al. 2015a,b). Article number, page 2 of 8 M. L. Ahnen et al.: MAGIC observations of the February 2014 flare of 1ES 1011+496 and ensuing constraint of the EBL density Time [MJD] 56700 56710 56720 ] fo r E > 2 0 0 G e V ­1 s ­2 F lu x [ c m 0 0.05 0.1 0.15 0.2 0.25 ­9 10× Fig. 1. 1ES 1011+496 light curve between February 6th and March 7th 2014 above an energy threshold of 200 GeV with a night-wise binning. The blue dashed line indicates the mean integral for the MAGIC ob- servations of 2007 and the red dotted line the MWL campaign between 2011 and 2012. On February 5th 2014, VERITAS (Weekes et al. 2002) is- sued an alert for the flaring state of 1ES 1011+496. MAGIC per- formed target of opportunity (ToO) observations for 17 nights during February-March 2014 in the zenith range of 20◦−56◦. After the quality cuts, 11.8 hrs of good data were used for fur- ther analysis. The data were taken in the so-called wobble-mode where the pointing direction alternates between four sky posi- tions at 0.4◦ away from the source (Fomin et al. 1994). The four wobble positions are used in order to decrease the systematic un- certainties in the background estimation. The data were analyzed using the standard routines in the MAGIC software package for stereoscopic analysis, MARS (Zanin et al. 2013). 3. Results After background suppression cuts, 6132 gamma-like excess events above an energy of 60 GeV were detected within 0.14◦ of the direction of 1ES 1011+496. Three control regions with the same γ-ray acceptance as the ON-source region were used to estimate the residual background recorded together with the signal. The source was detected with a significance of ∼ 75 σ, calculated according to Li & Ma (1983, eq. 17). Fig. 1 shows the night by night γ-ray light curve for ener- gies E > 200 GeV between February 6th and March 7th 2014. The emission in this period had a high night-to-night variability, reaching a maximum of (2.3 ± 0.1) × 10−10 cm−2s−1, ∼14 times the mean integral flux measured by MAGIC in 2007 and 2008 for 1ES 1011+496 (Albert et al. 2007; Reinthal et al. 2012) and ∼29 times the mean integral flux from the observation in 2011- 2012 (Ahnen et al. 2016b). For most of the nights the exposure time was ∼40 minutes, except for two nights (February 8th and 9th) in which the observations were extended to ∼2 hours. No significant intra-night variability was observed. The gap between observations seen in Fig. 1 was due to the strong moonlight pe- riod. The average observed spectral energy distribution (SED) is shown in Fig. 2. The estimated intrinsic spectrum, assum- ing the EBL model by Domínguez et al. (2011), can be fit- ted with a simple power-law function (PWL) with probability 0.35 (χ2/d.o.f.= 13.2/12) and photon index Γ = 2.0 ± 0.1 and normalization factor at 250 GeV f0 = (5.4 ± 0.1) × 10−11 Energy (GeV) 210 3 10 E 2 d ϕ /d E [T e V c m − 2 s − 1 ] ­1210 ­1110 ­1010 Fig. 2. Spectral energy distribution (SED) of 1ES 1011+496 for the 17 nights of observations between February 6th and March 7th 2014. The black dots are the observed data and the blue triangles are the data after EBL de-absorption. The red line indicates the fit to a broken power-law with transition region function of the observed SED whereas the blue line indicates the fit to a power-law function of the de-absorbed SED. cm−2s−1TeV−1. The observed spectrum is clearly curved. Several functions were tried to parametrize it: power-law with an expo- nential cut-off (EPWL), log-parabola (LP), log-parabola with ex- ponential cut-off (ELP), power-law with a sub/super-exponential cutoff (SEPWL) and a smoothly-broken power-law (SBPWL). Of these, only the SBPWL, dF dE = f0 ( E E0 )−Γ1 [ 1 + ( E Eb )g] Γ2−Γ1 g (1) achieves an acceptable fit (P = 0.17, χ2/d.o.f.= 12.8/9), though with a sharp change of photon index by ∆Γ = 1.35 within less than a factor 2 in energy. For the SBPWL, the normal- ization factor at E0 = 250 GeV is f0 = (4.2 ± 0.2) × 10−11 cm−2s−1TeV−1, the first index is Γ1 = 0.35 ± 0.01, the second index Γ2 = 1.7 ± 0.1, the energy break Eb = 298 ± 21 GeV and the parameter g = 12.6 ± 1.5. Among the other, smoother functions, the next-best fit is provided by the LP (shown in Fig. 2), with P = 1.7 × 10−3 (χ2/d.o.f.= 29.8/11). The photon index for the LP is Γ = 2.8 ± 0.1, curvature index β = 1.0 ± 0.1 and normalization factor at E0 = 250 GeV f0 = (3.6 ± 0.1) × 10−11 cm−2s−1TeV−1. This non-trivial shape of the observed spectrum, and its simplification when the expected effect of the EBL is cor- rected, strongly suggests this observation has high potential for setting EBL constraints. The night-wise estimated intrinsic spectra could all be fitted with power-laws, and the evolution of the resulting photon in- dices is shown in Fig. 3. In the latter part of the observed period, the activity of the source was lower, resulting in larger uncer- tainties for the fits. There is no evidence for significant spectral variability in the period covered by MAGIC observations, de- spite the large variations in the absolute flux. 4. EBL constraint We follow the procedure described in Abramowski et al. (2013) for the likelihood ratio test. The absorption of the EBL is de- scribed as e−ατ(E,z) where τ(E, z) is the optical depth predicted Article number, page 3 of 8 A&A proofs: manuscript no. MAGIC_1ES1011_EBL_Feb2016_toArX Time [MJD] 56695 56700 56705 56710 56715 56720 56725 P o w e r­ la w p h o to n i n d e x 0.5 1 1.5 2 2.5 3 3.5 Fig. 3. Evolution of the photon index from power-law fits to the de- absorbed night-wise spectra of 1ES 1011+496 between February 6th and March 7th 2014. The error bars are the parameter uncertainties from the fits. The red line represents the fit to a constant value, for which the probability is 10%. by the model, which depends on the energy E of the γ-rays and the redshift z of the source. With the optical depth scaled by a factor α, the observed spectrum is formed as:( dF dE ) obs = ( dF dE ) int × exp(−α × τ(E, z)) (2) where (dF/dE)int is the intrinsic spectrum of the source. The emission of HBLs, like 1ES 1011+496, is often well described by basic synchrotron self-Compton (SSC) models (e.g. Tavec- chio et al. 1998). A population of electrons is accelerated to ultrarelativistic energies with a resulting power-law spectrum with index Γe of about 2. The high energy electrons are cooled faster than the low energy ones, resulting in a steeper Γe. These electrons produce synchrotron radiation with a photon index Γ = Γe+1 2 = 1.5. In the Thomson regime the energy spectrum in- dex of the inverse-Compton scattered photons is approximately the same as the synchrotron energy spectrum, whereas in the Klein-Nishima regime, the resulting photon index is even larger. These arguments put serious constraints to the photon index of the energy spectrum of VHE photons. Additionally, in most of the SSC models, the emission is assumed to be originated in a single compact region, which results in a smooth spectral energy distribution with two concave peaks. The shape of the individ- ual peaks could be modified in a multizone model, where the emission is a superposition of several one-zone emission regions. However the general two-peak structure is conserved. For the modeling of the intrinsic source spectrum we have used the same functions as in Mazin & Raue (2007) and Abramowski et al. (2013) which were also used to fit the ob- served spectrum: PWL, LP, EPWL, ELP and SEPWL. We have added the additional constraint that the shapes cannot be convex, i.e. the hardness of the spectrum cannot increase with energy, as this is not expected in emission models, nor has it been ob- served in any BL Lac in the optically-thin regime. In particular, the un-absorbed part of BL Lac spectra measured by Fermi-LAT are well fitted by log-parabolas (Ackermann et al. 2012). The PWL and the LP are functions that are linear in their parameters in the log flux–log E representation (hence well- behaved in the fitting process), and both can model pretty well the de-absorbed spectrum found in Sect. 3. The EPWL, ELP and SEPWL have additional (non-linear) parameters that are physi- cally motivated, e.g. to account for possible internal absorption at the source. Note that these functions (except the PWL) can also mimic the overall spectral curvature induced by the EBL over a wide range of redshifts, but will be unable to fit the in- flection point (in the optical depth vs. log E curvature) that state- of-the-art EBL models predict around 1 TeV. We therefore ex- pect an improvement of the fit quality as we approach the true value of the scaling factor α, hence providing a constraint of the actual EBL density. The chosen spectral functions, however, do not exhaust all possible concave shapes. Therefore the EBL constraints we will obtain are valid under the assumption that the true intrinsic spectrum can be well described (whitin the uncer- tainties of the recorded fluxes) by one of those functions. As we saw in section 3, the 5-parameter SBPWL (not included among the possible spectral models) provides an acceptable fit to the ob- served spectrum; if considered a plausible model for the intrinsic spectrum, it would severely weaken the lower EBL density con- straint. On the contrary, the upper constraint (arguably the most interesting one VHE observations can contribute) from this work would be unaffected, as we will see below. To search for the imprint of the EBL on the observed spec- trum, a scan over α was computed, varying the value from 0 to 2.5. In each step of the scan, the model for the intrinsic spec- trum was modified using the EBL model by Domínguez et al. (2011), with the scaled optical depth using the expression (2) and then was passed through the response of the MAGIC telescopes (accounting for the effective area of the system, energy recon- struction, observation time). Then the Poissonian likelihood of the actual observation (the post-cuts number of recorded events vs Eest, in both the ON and OFF regions) was computed, af- ter maximizing it in a parameter space which includes, besides the intrinsic spectral parameters, the Poisson parameters of the background in each bin of Eest 1. The maximum likelihood was thus obtained for each value of alpha. This likelihood shows a maximum at a value α = α0, indicating the EBL opacity scal- ing which achieves a best fit to the data. A likelihood ratio test was then performed to compare the no-EBL hypothesis (α = 0) with the best-fit EBL hypothesis (α = α0). The test statistics TS = 2 log(L(α = α0/L(α = 0)), according to Wilks theorem, asymptotically follows a χ2 distribution with one degree of free- dom (since the two hypotheses differ by just one free parameter, α). Despite changing the flux level, the EBL determination method should work properly as long as the average intrinsic spectrum in the observation period can be described with one of the tested parameterizations. Assuming that is the case for the different states of the source, it will also hold for the average spectrum if the spectral shape is stable through the flare. A sim- ple way to check the stability of the spectral shape is fitting the points on Fig. 3 to a constant value. The χ2/d.o.f. of this fit is 23.5/16 and the probability is 10%, so there is no clear signature of spectral variability —beyond a weak hint of harder spectra in the second half of the observation period. A varying spectral shape would in any case need quite some fine tuning to repro- duce, in the average spectrum, a feature like the one expected to be induced by the EBL. 1 Note that in the Poissonian likelihood approach we have included the point at E ∼ 55 GeV which was shown just as an upper limit in Fig. 2, since it has an excess of just around 1 standard deviation above the background. The fit performed with the Poissonian likelihood approach has therefore one more degree of freedom than the χ2 fits reported in section 3, and the 55 GeV point is included in Fig. 7 Article number, page 4 of 8 M. L. Ahnen et al.: MAGIC observations of the February 2014 flare of 1ES 1011+496 and ensuing constraint of the EBL density αOpacity normalization 0 0.5 1 1.5 2 p ro b a b ili ty 2 χ 0 0.1 0.2 0.3 0.4 0.5 LP EPWL ELP PWL SEPWL Fig. 4. χ2 probability distributions for the average spectrum of the Feb- March flare of 1ES 1011+496 for the 5 models tested. PWL in blue (dashed line), LP in red (solid line), EPWL in green (dash-dot line), ELP in pink (dotted line) and SEPWL in light blue (long dash-dot line). αOpacity normalization 0 0.5 1 1.5 2 2 χ F it 0 20 40 60 80 100 120 140 160 180 200 LP EPWL ELP PWL SEPWL Fig. 5. Fit χ2 distributions for the average spectrum of the Feb-March flare of 1ES 1011+496 for the 5 models tested. PWL in blue (dashed line), LP in red (solid line), EPWL in green (dot-dash line) and ELP in pink (dotted line) and SEPWL in light blue (long dash-dot line). The LP red line is overlapping ELP pink line. Notice how all curves converge after reaching the minimum. Table 1. χ2 probabilities (P) for the cases of α = 0.0 and α = 1.07 Function P(α = 0.0) P(α = 1.07) LP 6.0 × 10−4 0.38 PWL 7.0 × 10−34 0.46 EPWL 4.5 × 10−9 0.38 ELP 3.2 × 10−4 0.30 SEPWL 6.2 × 10−5 0.30 Fig. 4 shows the χ2 probabilities for the five tested mod- els, also listed in Table 1 for the no-EBL case (α = 0.0) and the best-fit α = 1.07. The model that gives the highest proba- bility in the scanned range of α is the PWL. Following the ap- proach in Abramowski et al. (2013) would lead us to choose the PWL as model for the intrinsic spectrum, as the next models in complexity (LP and EPWL) are not preferred at the 2 σ level. However, choosing a PWL as the preferred model is rather ques- tionable, since would not allow any intrinsic spectral curvature, αOpacity normalization 0 0.5 1 1.5 2 2.5 T e s t S ta ti s ti c ( T S ) 0 5 10 15 20 Fig. 6. Test statistics distribution for the data sample for the 2014 Feb- March flare of 1ES 1011+496. The vertical lines mark the maximum and the uncertainty corresponding to 1 σ. meaning that all curvature in the observed spectrum will be at- tributed to the EBL absorption. If this procedure is applied to a large number of spectra, as in Biteau & Williams (2015), in- dividual < 2 σ hints of intrinsic (concave) curvature might be overlooked and accumulate to produce a bias in the EBL estima- tion. In our case, the assumption of power-law intrinsic spectrum for 1ES 10111+496 would lead the likelihood ratio test to prefer the best-fit α value to the no-EBL hypothesis by as much as 13 σ. We prefer to adopt a more conservative approach, choosing the next-best function, the LP. Note however that at the best-fit α, all the tested functions become simple power-laws, hence the fit probabilities at the peak depend only on the number of free parameters. At the best-fit α = 1.07, the parameters of the PWL are: photon index Γ = 1.9 ± 0.1 and normalization factor at 250 GeV f0 = (9.2±0.2)×10−10 cm−2s−1TeV−1. The other functions have the same values for these parameters. Going deeper in the behaviour of the fits for the five models, it can be seen in the Fig. 5 that after reaching the minimum, the χ2 are identical for all models. This happens because of the con- cavity restriction imposed to the functions. After reaching the point where the EBL de-absorption results in a straight power- law intrinsic spectrum, all three functions converge, and the de- absorbed spectra becomes more and more convex as α increases (so no concave function can fit it any better than a simple power- law). The shape of the spectrum observed by MAGIC is thus very convenient for setting upper bounds to the EBL density, un- der the adopted assumption that convex spectra are “unphysical”. Given the arguments in the previous lines, we take the LP as our model for the intrinsic spectrum. For the data sample from the 2014 February-March flare of 1ES 1011+496, the test statis- tics has a maximum of TS = 21.5 at α0 = 1.07+0.09 −0.13 (Fig. 6). This means that the EBL optical depth from the model of Domínguez et al. (2011) scaled by the normalization factor α0 is preferred over the null EBL hypothesis with a significance of 4.6 σ. Us- ing the EBL model of Franceschini et al. (2008) as template (as in Abramowski et al. (2013)) the test statistic using the LP as model for the intrinsic spectrum has a maximum of TS=20.6 at α0 = 1.14+0.09 −0.14, which is compatible with the result using Domínguez et al. (2011) within statistical uncertainty. We again remark that allowing for other concave spectral shapes, like the SBPWL, would severely affect our lower EBL bound. This would also be the case for earlier published EBL Article number, page 5 of 8 A&A proofs: manuscript no. MAGIC_1ES1011_EBL_Feb2016_toArX lower constraints based on gamma-ray data —especially those in which the PWL is among the allowed models for the intrin- sic spectrum. For the observations discussed in the present pa- per, the SBPWL would achieve an acceptable fit even in the no- EBL assumption. This and earlier constraints on the EBL den- sity through its imprint on γ-ray spectra hence rely on somewhat tentative assumptions on the intrinsic spectra —but assumptions which, as far as we know, are not falsified by any observational data available on BL Lacs. On the other hand, the upper bound we have obtained is robust, since it is driven by the fact that convex spectral shapes (completely unexpected for BL Lacs at VHE) would be needed to fit our observations if EBL densities above the best-fit value are assumed. 5. Systematic uncertainty The MAGIC telescopes has a systematic uncertainty in the abso- lute energy scale of 15% (Aleksić et al. 2015b). The main source of this uncertainty is the imprecise knowledge of the atmospheric transmission. In order to assess how this uncertainty affects the EBL constraint, the calibration constants used to convert the pixel-wise digitized signals into photoelectrons were multiplied by a scaling factor (the same for both telescopes) spanning the range -15% to +15% in steps of 5%. This procedure is similar as the one presented by Aleksić et al. (2015b). For each of the scaling factors the data were processed in an identical manner through the full analysis chain, starting from the image cleaning, and using in all cases the standard MAGIC MC for this observation period. In this way we try to asses the effect of a potential miscalibration between the data and the MC simulation. For all scaled data samples, χ2 profiles for α between 0 and 2.5 were computed. From the 1 σ uncertainty ranges in α ob- tained for the different shifts in the light scale, we determine the largest departures from our best-fit value α0, arriving to the final result α0 = 1.07 (-0.20,+0.24)stat+sys. 6. Discussion The relation of the γ-ray of energy Eγ from the source (measured in the observed frame) and the EBL wavelength at the peak of the cross section for the photon-photon interaction is given by: λEBL(µm) = 1.187 × Eγ(TeV) × (1 + z)2 (3) where z is equal or less than the redshift of the source. The en- ergy range used for our calculations was between 0.06 and 3.5 TeV. However, the constraining of the EBL following the method from Abramowski et al. (2013) is based in the fact that after de- absorbing the EBL effect, with the right normalization, the fea- ture between ∼ 100 GeV and ∼ 5-10 TeV is suppressed. In Fig. 7 we show a comparison between two cases where the residu- als were computed (ratio between the observed events and the expected events from the model). The plot on the left shows the residuals for the null EBL hypothesis α = 0, while the right pad shows the same plot for the case of the best fit EBL scaling α = 1.07. The differences start to show after 200 GeV, a region where the EBL introduces a feature (an inflection point) that can- not be fitted by the log-parabola. This is the feature that drives the TS value on which the EBL constraint is based. We therefore calculate the EBL wavelength range for which our conclusion is valid from the VHE range between 0.2 and 3.5 TeV. Energy [GeV] 210 3 10 o b s e rv e d /e x p e c te d 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 = 0.0α Energy [GeV] 210 3 10 = 1.07α Fig. 7. Ratio between the observed events and the expected events from the model of the intrinsic spectrum for two normalization values of the EBL optical depth, α = 0 to the left and α = 1.07 to the right, which corresponds to the normalization where the maximum TS was found. In both plots the line corresponding to a ratio=1 is shown. The energy range has to take into account the redshift de- pendency in Eq. (3) since the interaction of the γ-ray and the EBL photons can happen in any point between the Earth and the source. The range is between [(1 + z)2Emin, Emax], correspond- ing to a wavelength range of the EBL where the interaction with the γ-ray can take place along the entire path between the source and the Earth. In Fig. 8 we show the contours from the statistical + systematic uncertainty of the EBL flux density, derived scal- ing up the EBL template model by Domínguez et al. (2011) at redshift z = 0 . The wavelength coverage is in the so-called cos- mic optical background (COB) part of the EBL, where we found the peak flux density λFλ = 12.27+2.75 −2.29 nW m−2 sr−1 at 1.4 µm, systematics included. 7. Conclusions We have reported the observation of the extraordinary outburst from 1ES 1011+496 observed by MAGIC from February 6th to March 7th 2014 where the flux reached a level ∼ 14 times the observed flux at the time of the discovery of the source in 2007. The spectrum of 1ES 1011+496 during this flare displays little intrinsic curvature over > 1 order of magnitude in energy, which makes this an ideal observation for constraining the EBL. Al- though the source showed a high flux variability during the ob- served period, no significant change of the spectral shape was observed during the flare, which allowed us to use the aver- age observed spectrum in the search for an imprint on it of the EBL-induced absorption of γ rays. Such EBL imprint can be seen in the fit residuals of the best-fit achieved under the no- EBL assumption (Fig. 7, left). In the approach that we followed along this work, for the description of the intrinsic spectrum at VHE we restricted ourselves to smooth concave functions which have shown in the past to provide good fits to BL Lac spec- tra whenever the expected EBL absorption was negligible. Un- der this assumption, the best-fit EBL density at λ = 1.4µm is Fλ = 12.27+2.75 −2.29 nW m−2 sr−1, which ranks among the strongest EBL density constraints obtained from VHE data of a single source. Acknowledgements. We would like to thank the Instituto de Astrofísica de Ca- narias for the excellent working conditions at the Observatorio del Roque de los Muchachos in La Palma. The financial support of the German BMBF and Article number, page 6 of 8 M. L. Ahnen et al.: MAGIC observations of the February 2014 flare of 1ES 1011+496 and ensuing constraint of the EBL density m]µ [λ 1 10 210 3 10 ] ­1 s r ­2 [ n W m λ F λ 1 10 210 Dominguez et al. 2011 This work Abramowski et al. 2013 (HESS) Ackermann et al. 2012 (Fermi) Bernstein et al. 2002, 2005 Berta et al. 2010 (Herschel/PEP) thermin et al. 2010 (SPITZER)eB thermin et al. LL 2010 (SPITZER)eB Brown et al. 2000 (HST/STIS) sy et al. 2001 (DIRBE/2MASS)eCambr Dole et al. LL 2004 (SPITZER) Dole et al. 2006 (SPITZER) Dube 1979/Leinert 1998 UL Dwek & Arendt 1998 (DIRBE) Dwek & Arendt UL 1998 (DIRBE) Edelstein et al. 2000 (Shuttle/UVX) Elbaz et al. 2002 (ISO) Fazio et al. 2005 (SPITZER) Finkbeiner et al. 2000 (DIRBE) Frayer et al. 2006 (SPITZER) Gardner et al. 2000 (STIS) Gorjian et al. 2000 (DIRBE) Hauser et al. 1998 (DIRBE/FIRAS) Hauser et al. UL 1998 (DIRBE/FIRAS) Kashlinksy et al. 1996 Keenan et al. 2010 Lagache et al. 2000 (DIRBE) Lagache et al. UL 2000 (DIRBE) Levenson et al. 2007 (DIRBE/2MASS) Levenson & Wright 2008 (SPITZER) Madau & Pozzetti 2000 (HST) Martin et al. 1991 (Shuttle/UVX) Mattila et al. UL 1990 Mattila et al. 2011 Matilla et al. UL 2011 Matsumoto et al. 2015 (IRTS) Matsuoka et al. 2011 (Pioneer 10/11) Matsuura et al. 2010 (AKARI) Metcalfe et al. 2003 (ISO) Papovich et al. 2004 (SPITZER) Penin et al. 2012 Schlegel et al. 2000 Thompson et al. 2007 (NICMOS) Toller 1983/Leinert UL 1998 Voyer et al. 2011 Wright & Reese 2001 (DIRBE) Xu et al. 2005 (GALEX) Zemcov et al. 2014 Fig. 8. Extragalactic background light intensity versus wavelength at z = 0. The solid black line is the EBL template model (Domínguez et al. 2011) that we used for our calculations. The azure shaded area spans the wavelength range for which our constraint is valid, scaled from the EBL template. The width of the shaded area includes the statistical and systematic uncertainties. The orange area is the EBL constraint by Abramowski et al. (2013) and the green shaded area is the constraint by Ackermann et al. (2012). As a comparison we include the direct measurements by Dwek & Arendt (1998), Hauser et al. (1998), Finkbeiner et al. (2000), Lagache et al. (2000), Gorjian et al. (2000), Cambrésy et al. (2001), Wright (2001), Krick & Bernstein (2005), Matsumoto et al. (2005), Matsuoka et al. (2011), and Matsumoto et al. (2015). Also galaxy-count data is included, from Brown et al. (2000), Gardner et al. (2000), Madau & Pozzetti (2000), Elbaz et al. (2002), Metcalfe et al. (2003), Dole et al. (2004), Papovich et al. (2004), Fazio (2005), Xu et al. (2005), Frayer et al. (2006), Levenson et al. (2007), Thompson et al. (2007), Levenson & Wright (2008), Béthermin et al. (2010), Berta et al. (2010), Keenan et al. (2010) and Voyer et al. (2011). The upper limits shown are from Dube et al. (1979), Martin & Rouleau (1991), Mattila & et al. (1991), Kashlinsky et al. (1996) and Edelstein et al. (2000). Also shown are large-scale anisotropy measurements from Pénin et al. (2012) and Zemcov et al. (2014). MPG, the Italian INFN and INAF, the Swiss National Fund SNF, the ERDF under the Spanish MINECO, and the Japanese JSPS and MEXT is grate- fully acknowledged. This work was also supported by the Centro de Exce- lencia Severo Ochoa SEV-2012-0234, CPAN CSD2007-00042, and MultiDark CSD2009-00064 projects of the Spanish Consolider-Ingenio 2010 programme, by grant 268740 of the Academy of Finland, by the Croatian Science Foundation (HrZZ) Project 09/176 and the University of Rijeka Project 13.12.1.3.02, by the DFG Collaborative Research Centers SFB823/C4 and SFB876/C3, and by the Polish MNiSzW grant 745/N-HESS-MAGIC/2010/0. We thank the anonymous referee for a thorough review and a very constructive list of remarks that helped improving the quality and clarity of this manuscript. References Abramowski, A., Acero, F., & et al. 2013, A&A, 550, A4 Ackermann, M., Ajello, M., Allafort, A., Schady, P., & et al. 2012, Science, 338, 1190 Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., & et al. 2006, Nature, 440, 1018 Ahnen, M. L., Ansoldi, S., Antonelli, L. A., & et al. 2016a, submitted to MNRAS Ahnen, M. L., Ansoldi, S., Antonelli, L. A., & et al. 2016b, accepted to A&A Albert, J., Aliu, E., Anderhub, H., Antoranz, P., & et al. 2007, ApJ, 667, L21 Aleksić, J., Ansoldi, S., Antonelli, L. A., et al. 2015a, APh, in press Aleksić, J., Ansoldi, S., Antonelli, L. A., et al. 2015b, APh, in press Berta, S., Magnelli, B., Lutz, D., et al. 2010, A&A, 518, L30 Béthermin, M., Dole, H., Beelen, A., & Aussel, H. 2010, A&A, 512, A78 Biteau, J. & Williams, D. A. 2015, ApJ, 812, 60 Brown, T. M., Kimble, R. A., Ferguson, H. C., et al. 2000, AJ, 120, 1153 Cambrésy, L., Reach, W. T., Beichman, C. A., & Jarrett, T. H. 2001, ApJ, 555, 563 Cooray, A., Bock, J. J., Keatin, B., Lange, A. E., & Matsumoto, T. 2004, ApJ, 606, 611 Dole, H., Lagache, G., Puget, J.-L., et al. 2006, A&A, 451, 417 Dole, H., Rieke, G. H., Lagache, G., et al. 2004, ApJS, 154, 93 Domínguez, A., Finke, J. D., Prada, F., et al. 2013, ApJ, 770, 77 Domínguez, A., Primack, J. R., Rosario, D. J., et al. 2011, MNRAS, 410, 2556 Dube, R. R., Wickes, W. C., & Wilkinson, D. T. 1979, ApJ, 232, 333 Dwek, E. & Arendt, R. G. 1998, ApJ, 508, L9 Edelstein, J., Bowyer, S., & Lampton, M. 2000, ApJ, 539, 187 Elbaz, D., Flores, H., Chanial, P., et al. 2002, A&A, 381, L1 Article number, page 7 of 8 A&A proofs: manuscript no. MAGIC_1ES1011_EBL_Feb2016_toArX Fazio, G. G. 2005, Neutrinos and Explosive Events in the Universe, 47 Fazio, G. G., Ashby, M. L. N., Barmby, P., & et al. 2004, ApJS, 154, 39 Fernandez, E. R., Komatsu, E., Iliev, I. T., & Shapiro, P. R. 2010, ApJ, 710, 1089 Finkbeiner, D. P., Davis, M., & Schlegel, D. J. 2000, ApJ, 544, 81 Finke, J. D., Razzaque, S., & Dermer, C. D. 2010, ApJ, 712, 238 Fomin, V. P., Stepanian, A. A., Lamb, R. C., et al. 1994, Astroparticle Physics, 2, 137 Franceschini, A., Rodighiero, G., & Vaccari, M. 2008, A&A, 837 Frayer, D. T., Huynh, M. T., Chary, R., et al. 2006, ApJ, 647, L9 Gardner, J. P., Brown, T. M., & Ferguson, H. C. 2000, ApJ, 542, L79 Georganopoulos, M., Finke, J. D., & Reyes, L. C. 2010, ApJ, 714, L157 Gorjian, V., Wright, E. L., & Chary, R. R. 2000, ApJ, 536, 550 Gould, R. J. & Schréder, G. P. 1967, Physical Review, 155, 1408 Hauser, M. G., Arendt, R. G., Kelsall, T., et al. 1998, ApJ, 508, 25 Hauser, M. G. & Dwek, E. 2001, ARA&A, 39, 249 Kashlinsky, A., Mather, J. C., Odenwald, S., & Hauser, M. G. 1996, ApJ, 470, 681 Keenan, R. C., Barger, A. J., Cowie, L. L., & Wang, W.-H. 2010, ApJ, 723, 40 Krick, J. E. & Bernstein, R. A. 2005, American Astronomical Society Meeting Abstracts, 37, 177.14 Lagache, G., Haffner, L. M., Reynolds, R. J., & Tufte, S. L. 2000, A&A, 354, 247 Levenson, L. R. & Wright, E. L. 2008, ApJ, 683, 585 Levenson, L. R., Wright, E. L., & Johnson, B. D. 2007, ApJ, 666, 34 Li, T. & Ma, Y. 1983, ApJ, 317 Madau, P. & Pozzetti, L. 2000, MNRAS, 312, L9 Mankuzhiyil, N., Persic, M., & Tavecchio, F. 2010, ApJ, 715, L16 Martin, P. G. & Rouleau, F. 1991, Extreme Ultraviolet Astronomy, 341 Matsumoto, T., Kim, M. G., Pyo, J., & Tsumura, K. 2015, ArXiv e-prints Matsumoto, T., Matsuura, S., Murakami, H., et al. 2005, ApJ, 626, 31 Matsuoka, Y., Ienaka, N., Kawara, K., & Oyabu, S. 2011, ApJ, 736, 119 Mattila, K. & et al. 1991, The Early Observable Universe from Diffuse Back- grounds, 133 Mazin, D. & Raue, M. 2007, A&A, 471, 439 Metcalfe, L., Kneib, J.-P., McBreen, B., et al. 2003, A&A, 407, 791 Meyer, M., Raue, M., Mazin, D., & Horns, D. 2012, A&A, 542, A59 Orr, M. R., Krennrich, F., & Dwek, E. 2011, ApJ, 733, 77 Papovich, C., Dole, H., Egami, E., et al. 2004, ApJS, 154, 70 Pénin, A., Lagache, G., Noriega-Crespo, A., et al. 2012, A&A, 543, A123 Reinthal, R., Rügamer, S., Lindfors, E. J., et al. 2012, Journal of Physics Con- ference Series, 355, 012017 Stecker, F. W. & de Jager, O. C. 1996, Space Science Reviews, 75, 401 Stecker, F. W., Malkan, M. A., & Scully, S. T. 2006, ApJ, 774 Tavecchio, F., Maraschi, L., & Ghisellini, G. 1998, ApJ, 509, 608 Thompson, R. I., Eisenstein, D., Fan, X., Rieke, M., & Kennicutt, R. C. 2007, ApJ, 657, 669 Voyer, E. N., Gardner, J. P., Teplitz, H. I., Siana, B. D., & de Mello, D. F. 2011, ApJ, 736, 80 Weekes, T. C., Badran, H., Biller, S. D., et al. 2002, Astroparticle Physics, 17, 221 Wright, E. L. 2001, ApJ, 553, 538 Xu, C. K., Donas, J., Arnouts, S., et al. 2005, ApJ, 619, L11 Zanin, R., Carmona, E., Sitarek, J., & et al. 2013, Proc of the 33rd ICRC, Rio de Janeiro, 2 Zemcov, M., Smidt, J., Arai, T., et al. 2014, Science, 346, 732 1 ETH Zurich, CH-8093 Zurich, Switzerland 2 Università di Udine, and INFN Trieste, I-33100 Udine, Italy 3 INAF National Institute for Astrophysics, I-00136 Rome, Italy 4 Università di Siena, and INFN Pisa, I-53100 Siena, Italy 5 Croatian MAGIC Consortium, Rudjer Boskovic Institute, University of Rijeka, University of Split and University of Zagreb, Croatia 6 Saha Institute of Nuclear Physics, 1\AF Bidhannagar, Salt Lake, Sector-1, Kolkata 700064, India 7 Max-Planck-Institut für Physik, D-80805 München, Germany 8 Universidad Complutense, E-28040 Madrid, Spain 9 Inst. de Astrofísica de Canarias, E-38200 La Laguna, Tenerife, Spain; Universidad de La Laguna, Dpto. Astrofísica, E-38206 La Laguna, Tenerife, Spain 10 University of Łódź, PL-90236 Lodz, Poland 11 Deutsches Elektronen-Synchrotron (DESY), D-15738 Zeuthen, Germany 12 IFAE, Campus UAB, E-08193 Bellaterra, Spain 13 Universität Würzburg, D-97074 Würzburg, Germany 14 Centro de Investigaciones Energéticas, Medioambientales y Tec- nológicas, E-28040 Madrid, Spain 15 Università di Padova and INFN, I-35131 Padova, Italy 16 Institute for Space Sciences (CSIC\IEEC), E-08193 Barcelona, Spain 17 Technische Universität Dortmund, D-44221 Dortmund, Germany 18 Unitat de Física de les Radiacions, Departament de Física, and CERES-IEEC, Universitat Autònoma de Barcelona, E-08193 Bel- laterra, Spain 19 Universitat de Barcelona, ICC, IEEC-UB, E-08028 Barcelona, Spain 20 Japanese MAGIC Consortium, ICRR, The University of Tokyo, De- partment of Physics and Hakubi Center, Kyoto University, Tokai University, The University of Tokushima, KEK, Japan 21 Finnish MAGIC Consortium, Tuorla Observatory, University of Turku and Department of Physics, University of Oulu, Finland 22 Inst. for Nucl. Research and Nucl. Energy, BG-1784 Sofia, Bulgaria 23 Università di Pisa, and INFN Pisa, I-56126 Pisa, Italy 24 ICREA and Institute for Space Sciences (CSIC\IEEC), E-08193 Barcelona, Spain 25 Università dell’Insubria and INFN Milano Bicocca, Como, I-22100 Como, Italy 26 now at Centro Brasileiro de Pesquisas Físicas (CBPF\MCTI), R. Dr. Xavier Sigaud, 150 - Urca, Rio de Janeiro - RJ, 22290-180, Brazil 27 now at NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA and Department of Physics and Department of Astronomy, University of Maryland, College Park, MD 20742, USA 28 Humboldt University of Berlin, Istitut für Physik Newtonstr. 15, 12489 Berlin Germany 29 now at Ecole polytechnique fédérale de Lausanne (EPFL), Lau- sanne, Switzerland 30 also at INFN, Padova, Italy 31 now at Laboratoire AIM, Service d’Astrophysique, DSM\IRFU, CEA\Saclay FR-91191 Gif-sur-Yvette Cedex, France 32 now at Finnish Centre for Astronomy with ESO (FINCA), Turku, Finland 33 also at INAF-Trieste 34 also at ISDC - Science Data Center for Astrophysics, 1290, Versoix (Geneva) 35 now at Instituto de Física, Universidad Nacional Autónoma de Méx- ico, Apartado Postal 20-364, México D. F. 01000, México 36 * Corresponding authors: A. González Muñoz (adiv.gonzalez@fisica.unam.mx), A. Moralejo (moralejo@ifae.es) and P. Bangale (priya@mpp.mpg.de) Article number, page 8 of 8 1 Introduction 2 Observations & Analysis 3 Results 4 EBL constraint 5 Systematic uncertainty 6 Discussion 7 Conclusions