© ESO, 2017

1. Introduction

η Carinae is the most luminous massive binary system of our galaxy and the first one to have been detected at very high energies without hosting a compact object¹. It is composed by one of the most massive stars known (η Car A) with an initial mass estimated above M_A ≳ 90 M_⊙ (Hillier et al. 2001) and of a companion (η Car B) believed to be an O supergiant or a WR star.

η Carinae has been studied across the whole electromagnetic spectrum from the radio (Murphy et al. 2007) to the TeV (H.E.S.S. Collaboration et al. 2012), passing through the infrared (Whitelock et al. 1994, 2004), optical (Damineli 1996; Damineli et al. 2000, 2008), X-rays (Corcoran et al. 2001; Okazaki et al. 2008; Leyder et al. 2008, 2010), and γ-rays (Tavani et al. 2009; Abdo et al. 2010; Farnier et al. 2011). The presence of the companion star so far has been only indirectly inferred from the effect of the wind-wind collision (in particular by the X-ray emission from a multi keV plasma) and from the variable ultraviolet emission photoionizing nearby circumstellar clouds (Iping et al. 2005; Mehner et al. 2010).

η Car A is accelerating a very dense wind with a mass-loss rate of ~8.5 × 10^-4M_⊙ yr^-1 and a terminal wind velocity of ~420 km s^-1 (Groh et al. 2012). Its companion probably emits a fast low-density wind at 10^-5M_⊙ yr^-1 reaching a velocity of 3000 km s^-1 (Pittard & Corcoran 2002; Verner et al. 2005; Parkin et al. 2009).

The regular modulation detected in the X-ray light curves suggests that the two stars are located in a very eccentric orbit (Corcoran et al. 2001; Okazaki et al. 2008). The estimated orbital period at the epoch of the Great Eruption that happened between 1837–1856 was ~5.1 yr and since then has increased up to the current ~5.54 yr (Whitelock et al. 2004; Corcoran 2005; Damineli et al. 2008) owing to the huge quantities of mass and energy dissipated during the past century. During the Great Eruption, η Carinae experienced a huge outburst ejecting an impressive quantity of mass estimated as 10–40 M_⊙ (Gomez et al. 2010) at an average speed of ~650 km s^-1 (Smith et al. 2003), giving rise to the formation of the Homunculus Nebula and becoming one of the brightest stars of the sky. The energy released in such a catastrophic event (10⁴⁹⁻⁵⁰ erg) was comparable with a significant fraction of the energy emitted by a supernova explosion.

Given the high eccentricity of the orbit, the relative separation of the two stars varies by a factor ~20, reaching its minimum at periastron, when the two objects pass within a few AU of each other; the radius of the primary star is estimated as 0.5 AU. In these extreme conditions their supersonic winds interact forming a colliding wind region of hot shocked gas where charged particles can be accelerated via diffusive shock acceleration up to high energies (Eichler & Usov 1993; Dougherty et al. 2003; Reimer et al. 2006). As these particles encounter conditions that vary with the orbital phase of the binary system, one can expect a similar dependency in the γ-ray emission.

The hard X-ray emission detected by INTEGRAL (Leyder et al. 2008) and Suzaku (Okazaki et al. 2008), with an average luminosity (4–7) × 10³³ erg s^-1, suggested the presence of relativistic particles in the system. The following year AGILE detected a variable source compatible with the position of η Carinae (Tavani et al. 2009). Other γ-ray analyses followed, which reported a luminosity of 1.6 × 10³⁵ erg s^-1 (Abdo et al. 2010; Farnier et al. 2011; Reitberger et al. 2012), and suggested the presence of hard component in the spectrum around periastron, which subsequently disappeared around apastron. Such a component has been explained through π⁰ decay of accelerated hadrons interacting with the dense stellar wind (Farnier et al. 2011), or interpreted as a consequence of γ-ray absorption against an ad hoc distribution of soft X-ray photons (Reitberger et al. 2012).

An alternative acceleration scenario suggested by Ohm et al. (2010), which associated particle acceleration to the blast wave of the 1843 Great Eruption and foresaw a constant flux emission, was ruled out by the variability detected in the Fermi-LAT light curves.

At even higher energies, the observation by the H.E.S.S. Collaboration et al. (2012) did not lead to any significant detection, raising only an upper limit at energies ≳500 GeV. This in turn would imply a sudden drop in the γ-ray flux, which could be related to a cut-off in the accelerated particle distribution or to severe γ−γ absorption.

Our study starts with a new analysis of the Fermi-LAT data, including 25% more data than published previously and the latest version of the pipeline processing. We then present the results of a simulation of particle acceleration in the colliding wind, based on detailed hydrodynamic simulations of Parkin et al. (2011), and compare these with the observations. These comparisons are very successful and lead to a number of conclusions reinforcing our previous interpretations (Leyder et al. 2008; Farnier et al. 2011).

2. Fermi-LAT data analysis

Launched on 2008 June 11, the Large Area Telescope (LAT) on-board the Fermi Gamma-ray Space Telescope is the most sensitive γ-ray telescope to date, covering an energy range from 20 MeV to 300 GeV (Atwood et al. 2009). The LAT is characterized by large field of view (2.4 sr at 1 GeV) and collecting area (~6500 cm² at 1 GeV), a low deadtime (<100 μs per event), a high time resolution (<10 μs), and an energy dependent point-spread function (PSF) improving from ~5° (68% containment) at 100 MeV to ~0.1° at 40 GeV. The LAT consists of a charged particle tracker, a calorimeter, and an anti-coincidence system.

The electron-positron pair conversion tracker is made of 36 layers of silicon strip detectors to track charged particles, interleaved with 16 layers of tungsten foil to facilitate the conversion of γ-rays to pairs, these layers are comprised of 12 thin layers in the front section followed by 4 thick layers in the back, of 0.03 and 0.18 radiation length, respectively. The calorimeter, located at the bottom of the instrument, has ~8.5 total radiation lengths of caesium iodide to measure the total event energy. Given the intense background of charged particles from cosmic rays and trapped radiation at the orbit of the Fermi satellite, the instruments are protected by a segmented anti-coincidence detector used to reject charged-particle background events. More information about the LAT is provided in Atwood et al. (2009), the LAT in-flight calibration is described in Abdo et al. (2009), Ackermann et al. (2012a), and Ackermann et al. (2012b).

We performed the analysis of the recent Fermi-LAT data since the beginning of the regular survey-mode observations on 2008 August 4, until 2015 July 1, more precisely mission elapsed time (MET): 239557417 to 457485024 s. We used only source class data (i.e. reconstructed events with high probability of being photons) that have been reprocessed with the PASS8² pipeline and subsequently analysed using the Fermi Science Tool v10r0p5 package³. With respect to the previous PASS7, the new PASS8 introduced improvements in the reconstruction of the event direction, energy measurement, event selection, ghost handling, and track/anti-coincidence detector matching information. All these yield to an enhancement in the entire performance of the telescope. These enhancements include higher acceptance, larger field of view, smaller PSF, better energy resolution, deeper differential sensitivity, smaller systematic uncertainties, and the introduction of four subgroups of PSF and energy dispersion to improve even further the containment radius and energy resolution at the expense of the statistics. We selected data with evclass = 128 corresponding to the source class and evtype = 3 indicating that both categories of photons converted in the front and in the back part of the LAT were selected. We followed the analysis recommendations of the Fermi-LAT team for the time and event selection, rejecting photons with apparent zenith angle greater than 90° in order to minimize the background due to the atmospheric γ-rays originated from the Earth’s limb, which lies at a zenith angle of ~113°. We did not perform a zenith angle cut based on the region of interest (ROI) with gtselect, but did correct the exposure with the zmax = 90 option in the gtltcube tool.

With its orbital period of 96.5 min, its maximal rocking angle of 60° and its precession period of 53.4 days, Fermi-LAT spends more than 80% of its operative time in survey mode providing a uniform coverage of the sky every two orbits (~3 h). As a consequence, our target of interest is not always visible, so knowing the exact position and orientation of the satellite from the spacecraft files at each time, we ran gtmktime to obtain all the good time intervals in which to perform our analysis; we also excluded the readout dead-time (~9%) and the time intervals corresponding to the south Atlantic anomaly (~13%) when data taking is suspended. Finally data are filtered to accept only those events flagged with (DATA_QUAL==1) && (LAT_CONFIG==1), which exclude bad quality events and instrument configuration not recommended for scientific analysis, respectively.

Given the relative distance between the apparent position of the Sun during the year and the nominal position of η Carinae, we can clearly neglect the γ-ray contribution from our local star. The same is also valid for the Moon.

The instrument response function (IRF) of the Fermi-LAT, i.e. the description of the instrument performance provided for data analysis, strongly depends on the energy (Ackermann et al. 2012a); thus to better exploit its performance we made a separate analysis for photons above and below 10 GeV. This threshold was chosen because of the spectral shape of η Carinae, as explained in the third paragraph of Sect. 2.2.

2.1. High-energy emission: 10–300 GeV

For the high-energy analysis, we took into account only photons arriving within 3° from the nominal position of η Carinae (RA = 161.264775, Dec = –59.684431), as the total (front+back) 95% PSF containment angle above 10 GeV is smaller than 1° on-axis. We created a sky model using the 3FGL Fermi-LAT four-year catalogue (Acero et al. 2015) including all sources up to 1° outside of the ROI. We used the same source spectral models as indicated in the catalogue, leaving the normalization parameter free to vary for all those sources within 2° from η Carinae and with an average test statistic⁴ (TS) reported from the catalogue greater than 100 (corresponding to a detection significance ~10σ) or presenting a variability index higher than 72.44, which indicates a 99% confidence probability that the source is variable on a monthly timescale (Acero et al. 2015). The only exception to the model was made for the source 3FGL J1043.6-5930 (hereafter J1043), which does not satisfy the condition TS> 100 but is only 0.27° away from η Carinae (for a 68% containment radius of 0.23° above 10 GeV, Acero et al. 2015). This source has been modelled in the 3FGL catalogue as a power law (PL) with index Γ = 2.07 ± 0.11 and we decided to leave its normalization parameter free as well. For completeness we also checked if there was any source from the FAVA (Ackermann et al. 2013) weekly flare list⁵ present in our model. As the closest source lies more than 6° away from the binary system, we did not take any particular precaution. Finally, with the only exception of η Carinae, no other source is mentioned in the second Fermi-LAT catalogue of high-energy sources (Ackermann et al. 2016).

η Carinae lies on the tangential projection of the Carina-Sagittarius Arm of the Milky Way, less than 1° away from the Galactic plane. Consequently a correct description of the diffuse Galactic γ-ray background plays a key role in the analysis. To reproduce this emission we used the gll_iem_v06.fits⁶ (Acero et al. 2016) Galactic background model. Unfortunately, for several reasons the representation of the diffuse emission is a serious concern in the Carina region, becoming more and more delicate towards lower energies (see Sect. 2.2). As a matter of fact a significant number of sources from the catalogue (Acero et al. 2015) are flagged with “c” in that region, indicating that they are considered to be potentially confused with the Galactic diffuse emission. In order to obtain a more realistic representation of the Galactic diffuse emission in such a small region we decided to let its normalization free. The isotropic background component has been represented using the model iso_P8R2_SOURCE_V6_v06.txt⁶.

In the high-energy range we described the emission coming from η Carinae (its signal-to-noise is very limited) with a simple PL model dN/d$\mathit{E}\mathrm{=}{\mathit{N}}_{{\mathit{E}}_{\min }\mathrm{-}{\mathit{E}}_{\max }}\mathrm{(}\mathrm{\Gamma }\mathrm{+}\mathrm{1}\mathrm{)}{\mathit{E}}^{\mathrm{\Gamma }}\mathit{/}\mathrm{(}{\mathit{E}}_{\max }^{\mathrm{\Gamma }\mathrm{+}\mathrm{1}}\mathrm{-}{\mathit{E}}_{\min }^{\mathrm{\Gamma }\mathrm{+}\mathrm{1}}\mathrm{)}$, where N_Emin−Emax indicates the integrated photon flux over the fitting energy range [ E_min,E_max ] and Γ represents the PL index.

Exploiting this model we performed a maximum-likelihood analysis (Mattox et al. 1996) on the previously described sample of photons, using the P8R2_SOURCE_V6 IRF and the gtlike tool. Checking the results of the first fit, we found that the significance of some of the sources, left free in our model, were below 5σ. This does not contradict the catalogue results because even though we are increasing the statistic covering a longer period with our data sample, we previously took the average significance over the whole energy band in our choice of the model, while here we are performing the analysis only above 10 GeV. Thus we proceeded to freeze those sources parameters to the catalogue values. We did the same for all the “c” flagged sources, as they might be just artefacts of a bad representation of the diffuse Galactic model at lower energies. Finally we let free the normalization parameter of the Galactic diffuse emission and isotropic background, and the fit resulted in an integrated flux for η Carinae of F_{10−300 GeV} = (5.06 ± 0.52) × 10^-10 ph cm^-2 s^-1 and index Γ = −2.28 ± 0.15 with a TS value of 366 (significance ≳19σ). Those results represent an average over the whole time period of nearly seven years and are in good agreement with that reported in Reitberger et al. (2015). The slightly lower flux value could be explained by the lower than average flux in the data of the previous year (Corcoran et al. 2015a).

The resulting normalization factor of the Galactic diffuse emission differs significantly from unity. This is a known problem in the Carina region. Running again the fit, forcing the normalization of the Galactic diffuse emission component to unity, we obtain a 22% higher integrated flux for η Carinae, while the Γ index remains unchanged. On the contrary, the isotropic background does not present specific issues on the Galactic plane, so we decided to freeze it to unity and let vary just the Galactic component normalization. Running again the likelihood, the fit yields a Galactic normalization of 1.95 ± 0.06, a slightly lower integrated flux for η Carinae F_10−300GeV = (5.00 ± 0.51) × 10^-10 ph cm^-2 s^-1 and an index Γ = −2.27 ± 0.15, for a significance of nearly 19σ. Finally, as we do not expect any temporal variation of the Galactic diffuse emission we kept its normalization frozen to the seven-year average for the subsequent fit.

We then refined our analysis by dividing the photon sample into more time intervals. The duration of the bins was chosen to obtain a clear detection (TS ≳ 16) of η Carinae in all time intervals and to obtain a sequence repeating from an orbit to the next. We split the light curve into two periods, one for periastron and one for apastron, and subsequently we applied a static binning in each of them. The final eight time bins are reported in Table 1.

We attributed to each photon a phase calculated from periastron times using JD = 2 450 799.792 + N_p·(2024 ± 2) (Corcoran 2005), where N_p counts the successive periastrons. A more recent analysis (Damineli et al. 2008) suggests a shorter orbital period of 2022.7 ± 1.3 days, but such a small variation does not have any impact on our results. We then ran another likelihood analysis for each time bin, but given the shorter exposure, we detected a very low significance for the source J1043 in the 4th, 5th, and 8th bin and, consequently, we decided to fix its contribution in those bins to the average value reported in the 3FGL catalogue.

The light curve of the high-energy flux of η Carinae obtained from the likelihood analysis is reported in Fig. 1. After the first periastron passage of 2009 the flux of η Carinae decreased slightly towards apastron. The flux did however not increase again toward the periastron of 2014. The last two bins of this light curve can be directly compared with the first two, having the same exposures and orbital phases. In order to search for any faster variability we reduced the temporal size of the bins, which confirmed the absence of any excess during the second periastron.

Fig. 1 Seven-year high-energy flux light curve of η Carinae obtained from the binned analysis, using the binning as reported in Table 1 (blue points) and with a smaller binning (red points). Error bars are 1σ and superposed upper limits are 95%. For comparison we plotted an arbitrarily rescaled X-ray light curve (grey line).

When performing the fit in each bin, letting free the spectral index Γ, we can observe that the resulting values are constant within the uncertainties. When performing the same analysis fixing Γ to its average value, we found that the resulting fluxes were changing by a few percent at most, i.e. much less than the statistical uncertainties. Therefore the light curve presented in Fig. 1 does not depend significantly on the exact spectral shape assumed.

Given the relatively low event statistics at high energy, as a countercheck we also performed an unbinned analysis (Mattox et al. 1996) on the same sample of photons. The results that we obtained are F_{10−300 GeV} = (4.74 ± 0.49) × 10^-10 ph cm^-2 s^-1 and index Γ = −2.29 ± 0.15. Both the flux light curve and the spectral index trend are in very good agreement with the results of the binned analysis. For comparison we also ran an analysis with a smaller bin as reported in Fig. 1. These results show a bigger error but are consistent with the previous analysis within the uncertainties, even if on a few occasions the smaller statistics did not yield a firm detection of the source and only reached a 95% upper limit.

We gave special attention in our analysis to J1043. In the Fermi 3FGL catalogue a relatively hard spectrum is reported for that source. At low energy its flux is more than one order of magnitude lower than that of η Carinae, and it reaches nearly one-third of that flux above 10 GeV. Given its small distance (≲2r₆₈ PSF), we analysed the impact that a wrong representation of this source could have on the flux of η Carinae. To set an upper limit on the possible systematic error introduced, we considered the following two cases. We first ran a likelihood analysis keeping the J1043 parameters frozen to its seven-year average value, which gave us a light curve for η Carinae with flux values reduced by 5%. Then as a counter check, we completely removed J1043 from the model, obtaining a biased light curve for η Carinae with values up to 11% higher. All those results are in agreement with the TS maps we obtained for each single bin (see Fig. 2).

Fig. 2 High-energy TS maps for each of the time intervals, corrected for the small differences of exposure times between the 8 phase bins such that the images illustrate the source variability. Each image has the same width (0.77°). Rows a) and b) show TS maps obtained from the binned analysis, respectively including and excluding J1043 from the model. Rows c) and d) represent the same as a) and b) but from the unbinned analysis. The linear colour map spans TS from 0 to 100. The green and purple crosses are the positions of η Carinae and J1043.

The average seven-year TS maps perfectly match the PSF of the instrument, while in the shorter time bins the emission from η Carinae is often broadened. So far in our analysis we have always kept fixed the location of η Carinae, using its nominal coordinates. Exploiting the unbinned analysis and the gtfindsrc tool, we left the spatial coordinates as a free parameter in the fit and looked for the best coordinates to maximize the likelihood. Such an analysis has been performed in each single bin. The results are shown in Fig. 3 and reported in Table 2. The eight error circles of η Carinae are represented for each time bin, while the black circle indicates the result obtained running gtfindsrc on the whole seven-year data sample. The radius of each circle represents a 1σ error. The nominal position of η Carinae, shown by the black cross, is very well in agreement with the results of the Fermi analysis, being within the error circle more than 78% of the time.

Fig. 3 Error circles of 1σ for η Carinae derived from the analysis of each time bin. Labels identify the time bin as in Table 1; black circle refers to the whole 7 years of data. The black cross shows the nominal position of η Carinae. As for Fig. 2, the image is in Galactic coordinates with north up and longitude increasing towards the left.

2.2. Low-energy emission: 0.3–10 GeV

Extending the analysis to lower energies, the PSF becomes broader and the effective area, acceptance, and energy resolution worsen, making the analysis more challenging. In particular it requires us to enlarge the ROI and consequently to increase drastically the number of sources, the number of free parameters, the uncertainties (also related to the non-perfect Galactic diffuse emission model), and the computation requirements. For these reasons we chose a lower bound for our analysis of 300 MeV, as a compromise. This choice also keeps the flux systematic uncertainty below 5%⁷.

As the photon statistic is much better at lower energies, we used only a binned analysis, similar to that described in Sect. 2.1 but with the following adaptations. As the total (front+back) 95% containment angle above 300 MeV is now smaller than 7° (on-axis), we took into account all photons arriving within 14° from the nominal position of η Carinae and created a sky model including all sources in the ROI enlarged by 7°. We let the normalization parameters free to vary for all sources within 7° from η Carinae with an average TS> 100 or presenting a variability index higher than 72.44 (see Sect. 2.1 for explanation). As the flux of J1043 is much smaller than that of η Carinae at low energy, we kept frozen all its parameters to the values given in the 3FGL catalogue. Even though the extended pulsar wind nebula HESS J1303-631 lies more than 16° away from η Carinae, we included an appropriate extended⁸ representation in the model. We checked again in the FAVA (Ackermann et al. 2013) weekly flare list and found seven events within our ROI and verified that all were effectively associated with sources whose normalization were left free to vary.

To obtain a rough idea of the spectral shape of η Carinae, we started to split our analysis into separate energy bins, where the spectrum could be approximated locally as a simple PL. At the same time we kept a sufficient statistic in order to detect our source with at least 5σ. We defined five logarithmically equal energy bins, from 300 MeV up to 100 GeV, and performed a separate likelihood analysis in each of them for the different orbital phase intervals. In the most energetic bin [30–100 GeV] η Carinae did not always reach the required TS, and in these cases we merged the 4th and 5th energy bins together in the analysis. We computed the integral energy flux for each bin, converted them to luminosities (assuming a distance of 2.3 kpc), and plotted them in Fig. 7 (later described in Sect. 3) for two orbital phase bins (0.92–1.05 and 0.39–0.59), corresponding to periastron and apastron, respectively. In these plots, the centre of each point is not a logarithmic average, but is computed making a weighted average using the energy dependent function resulting from each fit as a weight. The result of this analysis indicates that the low-energy spectrum of η Carinae features some curvature that could be represented locally, for example by a cutoff power law or a broken power law and that an excess could be observed in some spectra above 10 GeV.

We therefore performed the complete analysis, in band 0.3–10 GeV, assuming a power law spectrum with an exponential cutoff (PLEC) dN/dE = N_{1 GeV}(E/ 1 GeV)^Γe^{− E/E_c} for η Carinae, where N_{1 GeV} is the normalization in units of ph cm^-2 s^-1 MeV^-1, Γ is the power law index, and E_c is the cutoff energy.

To perform the analysis in the orbital phase bins we always fixed the Galactic and extragalactic diffuse emission normalization to their seven-year averaged values (N_Gal = 0.962 ± 0.002 and N_Exgal = 1.20 ± 0.03). Leaving these two normalizations free affects the flux of η Carinae by less than 5%. The difference of normalization obtained between the high and low-energy analysis for the diffuse emission is related to the very different sizes of the respective ROI.

Fig. 4 Integrated flux light curve of η Carinae assuming different spectral models: PLEC with three (green) or two (blue) free parameters and BPL with three (red) or two (magenta) free parameters. Error bars are 1σ.

The values of the photon index and energy cutoff are locally correlated and not very meaningful. Figure 4 compares the integrated flux obtained with the PLEC model with all parameters free or fixed to the average of Γ_PLEC = −2.12 and with a broken power law (BPL) model with the low-energy spectral index left free or fixed to the average of Γ_BPL = −2.14. The low-energy γ-ray flux light curve does not change significantly and is therefore a good measure of the emission of η Carinae. Figure 5 shows the combined result of the binned high-energy analysis obtained in Sect. 2.1 together with the results derived here for the low-energy band (assuming a distance of 2.3 kpc).

As the data above 10 GeV of the second periastron fails to reproduce the high flux levels detected during the first periastron, the light curve of the low-energy band is now of prime importance to confirm the orbital variability of the γ-ray emission. The probability of obtaining an orbital flux variation by chance, as reported in Fig. 4, is lower than 5 × 10^-9 (5.9σ).

Finally, as the low-energy γ-ray variability appears similar for the two periastrons (see Fig. 4) and the spectra are well compatible within 2σ (see Fig. 7), we used the good statistics available to perform a merged analysis of these two periods on even shorter time bins. We merged the data covering the period from phase [0.92–1.19] and [1.92–2.19] and split the data according to phase bin intervals of 2% (i.e. ~40 days). We performed a single likelihood analysis in each time bin, representing the emission of η Carinae with a simple PL. η Carinae was detected above 12σ in every bin, increasing the detection significance up to 50% compared to a single periastron analysis. We also performed the analysis by shifting the central phase point of each bin by half of its width to obtain a better sampling of the variability. The results of both analyses are shown in Fig. 6. We also performed another similar analysis increasing the lower energy threshold to 600 MeV and 1 GeV to reduce the energy range, trying to exploit the better PSF at higher energy and increasing at the same time the robustness of the model approximations. The results showed exactly the same variability trend and similar likelihoods were obtained. We added in Fig. 6 the results obtained previously for two broad bins adjacent to the periastron period.

3. Comparison with simulations

Parkin et al. (2011) presented three-dimensional hydrodynamical simulations of η Carinae including radiative driving of the stellar winds (Castor et al. 1975), optically thin radiative cooling (Kaastra & Mewe 2000), gravity, and orbital motion. The main aim of these simulations was to reproduce the X-ray emission by analysing the emissivity and self-obscuration of the stellar wind. The simulations reproduced the observed X-ray spectra and light curves reasonably well, excepting the post-periastron extended X-ray minimum, where flux was overestimated and the wind collision disruption was inhibited. Additional gas cooling, for example by particle acceleration and inverse-Compton processes, could increase the cooling and disruption of the central wind collision zone.

Fig. 5 Simulated and observed X-ray and γ-ray light curves of η Carinae. The black and purple lines and bins show the predicted inverse-Compton and neutral pion decay light curves. The green and red points show the observed Fermi-LAT light curves at low (0.3–10 GeV) and high (10–300 GeV) energies. The dim grey light curves show the observed (continuous) and predicted (dash, without obscuration) thermal X-ray light curves. Error bars are 1σ.

Parkin et al. (2011) provided us with the results of their simulations, i.e. temperature, density, and three-dimensional velocities in the cells of the adaptive mesh for various orbital phases. To estimate the non-thermal emission we first calculated the maximum energies that could be reached by electrons and hadrons (as in Farnier et al. 2011) cell by cell assuming a di-polar magnetic field at the surface of the main star, perpendicular to the orbital plane (reality is probably more complex with the two stars contributing). The magnetic field is the only additional parameter, which can be tuned. We calculated shock velocities and mechanical power in every cell, including those outside the shock region. As expected, most of the shock power is released on both sides of the wind collision zone and in the cells downstream the wind-collision region (Reimer et al. 2006). The increasing shock area compensates for the loss of the released energy density up to a relatively large distance from the centre of mass, explaining why the X-ray luminosity at apastron is about a third of the peak emission at periastron.

The energy available in electrons and hadrons were then summed in the ranges 0.3 < E_e < 10 GeV and E_p > 20 GeV, respectively, to match the spectral bands observed by Fermi-LAT. The local cell physical properties can be used to easily estimate pion production as long as the Larmor radius is similar to the cell size. The minimum size of the cells in the simulation is ~10¹¹ cm, which is larger than the proton Larmor radius for Lorentz gamma factor up to 10⁵. Only one-third of the power accelerating protons is available to produce γ-rays through the neutral pion channel. Electron cooling and pion decay occur instantaneously when compared to other timescales.

To consider the possible effects of photon-photon opacity we calculated the X-ray thermal emission in each cell and evaluated the optical depth along different lines of sight. As the current orientation of the binary system, with respect to the Earth, still presents some uncertainties (Madura & Groh 2012), we used several possible directions; this provided optical depth τ that varied between ~10^-6 at apastron and ~10^-2 at periastron. This excludes explaining the 1–100 GeV spectral shape by the effects of photon-photon absorption (Reitberger et al. 2012).

Fig. 6 Merged Fermi LAT analysis (0.3–10 GeV) of the two periastrons for narrow time bins. The two broad bins and the black curve are the same as in Fig. 5.

Thermal emission increases towards periastron. The mechanical luminosity available in the shock also increases towards periastron and almost doubles in the phase range ≈1.05–1.15. The latter peak corresponds to a bubble with reverse wind conditions developing because of the orbital motion, effectively doubling the shock front area during about a tenth of the orbit (see Fig. 9 of Parkin et al. 2011). The density of this bubble is low so its thermal emission (∝density²) does not contribute significantly to the X-ray light curve. The mechanical luminosity shows a local minimum between phases 1.0 and 1.05 when the central part of the wind collision zone is disrupted.

Electron cooling, through inverse-Compton scattering, is very efficient and such γ-rays are expected to peak just before periastron. A secondary inverse-Compton peak could be expected above phase 1.05, although its spectral shape could be very different as the UV seed thermal photons are of lower density when compared to the location of the primary shock close to the centre of the system. In our simplified model we assumed that the spectral shape of the seed photons is the same in all cells of the simulation (r^-2 dependency is taken into account), and that these soft photons are sufficient to cool down all the relativistic electrons. The relative importance of the second peak, however, depends on the magnetic field geometry; radiation transfer, which is neglected in our model; obscuration; and details of the hydrodynamics, which do not represent the soft X-ray observations very well in this phase range. These details are not well constrained by the available observations and we did not try to refine them.

The situation is different for hadrons. Unless the magnetic field is very strong (>kG) hadronic interactions mostly take place close to the centre and a single peak of neutral pion decay is expected before periastron.

Figure 5 shows the X and γ-ray light curves predicted by the simulations for a magnetic field of 500 G and assuming that 1.5% and 2.4% of the mechanical energy is used to accelerate electrons and protons, respectively. To ease the comparison between observations and simulations, the results of the latter were binned in the same way as the observed data.

The thermal X-ray emission matches the observations pretty well (by construction, Parkin et al. 2011). For the simulated curve in Fig. 5, we use the orbit-IA model by Parkin, which uses instantaneous acceleration and not radiative driving of the winds. In addition, the thermal X-ray light curve does not take self-obscuration into account and therefore does not match the observations around periastron. The predicted γ-ray emission induced by the hadrons and electrons are also at the right level, although significant differences exist between simulations and observations.

Fig. 7 Electrons and photons luminosity spectra at periastron (phase 0.92–1.06; left) and at apastron (phase 0.39–0.59; right). The top panels show the spectra (arbitrary units) of the electrons accelerated in the wind of the primary (blue) and secondary (green) stars and their sum (red). The lower panels show the inverse-Compton emission of both components and the total emission, under the highly simplified assumption that the inverse-Compton parameters (geometry and soft photon spectra) are the same in all cells. The black and grey points are the broadband fluxes derived from Fermi data for the first and the second orbital cycle, respectively. The simulation results were averaged over the orbital phase range corresponding to the periastron observation, as the electron spectra vary quickly during that interval.

Both the predicted inverse-Compton emission and the observed (0.3–10 GeV) LAT light curve show a broad peak extending on both sides of periastron, as expected from the evolving shock geometry. The amplitude of the variability in the simulation depends on the number/size of those cells where particles can be accelerated up to relevant energies, which in turn depends on the magnetic field. Probing the range suggested by Walder et al. (2012), a surface magnetic field that is larger than 400 G provides a good match to the observations, while lower fields produce variations that are too large. In this work we did not considered any magnetic field amplification at the shock, which in turn could obviously scale down the surface magnetic field required to get equivalent results. Assuming a field of 500 G for the rest of the discussion, the predicted flux at phase 1.1 is two times larger than observed. This discrepancy largely comes from the energy released in the inverted wind bubble after periastron. The ratio of the emission generated in the shocks on both sides of the wind collision zone is relatively constant along the orbit except at phase 1.1, where much more power is generated in the shock occurring in the wind of the secondary star. The inverted bubble might either be unstable in reality or might produce a significantly different inverse-Compton spectrum.

Relativistic electrons immersed in such a high magnetic field produce a synchrotron radiation at low energy. The ratio of the energy that electrons lose via synchrotron and inverse-Compton processes is equal to the ratio of the magnetic field energy density over the photon field energy density, i.e. P_synch/P_IC = U_B/U_rad = B²4πR²c/ 8πL ≈ 7.8 × 10^-31·B [ G ] ²·R [ cm ] ². If we know the inverse-Compton spectrum, we can estimate the synchrotron peak luminosity, which around apastron results to be several orders of magnitude (~10⁶) fainter than the inverse-Compton peak. The synchrotron emission peak should reach its maximum in the optical band only very close to periastron, as it is only two orders of magnitude fainter than the inverse-Compton peak. Those limits are in agreement with the estimated radio upper limit (Duncan & White 2003).

Since the low-energy spectra during both periastrons are sufficiently in agreement (see Fig. 7 described later), we analysed simultaneously the Fermi LAT low-energy data derived from the two periastrons, binned in shorter time intervals (Fig. 6). These data show a peak at periastron, a minimum at phase 1.02, and a second broad peak at phase 1.1. This is very similar to the prediction of the simulation for the inverse-Compton luminosity. The only notable exception is that the observed second broad peak is slightly shifted towards earlier phases and has a lower luminosity when compared to the simulation. The similarities between the observations and the simulation for the γ-ray peak and minimum with consistent duration and amplitude are very encouraging. The phase difference could be related to the eccentricity (ϵ = 0.9) assumed in the simulation, which is not well constrained observationally (Damineli et al. 2000; Corcoran et al. 2001), and this has an important effect on the inner shock geometry.

Figure 7 shows that the distribution of γ_e, weighted by the emissivity, is relatively smooth and that the expected photon distribution is very smooth. The difference in the electron spectral shape on both sides of the wind collision zone cannot explain the two components γ-ray emission as suggested by Bednarek & Pabich (2011), who assumed a simplified geometry. We obtain a good match between the observed low-energy γ-ray spectrum and the predictions of the simulations at periastron, even though some discrepancy can be observed at apastron where an excess is observed between 2 and 10 GeV.

The inverse-Compton emission peaks slightly below 1 GeV and does not extend beyond 10 GeV at a level that is consistent with the observations during the first periastron in contrast with the conclusions from Ohm et al. (2015), which attribute the full Fermi LAT detection to hadronic emission. Their simulations predict a smaller variation between periastron and apastron, a longer flare around periastron, and a deeper minimum when compared to the observed data. Such discrepancies might be due to the simplified geometry assumed by the authors and by the artificially reduced particle acceleration at periastron. Inverse-Compton emission and neutral pion decay (Farnier et al. 2011) remains therefore a very good candidate to explain the Fermi observations. The fraction of the shock mechanical luminosity accelerating electrons appears to be slightly smaller than the fraction that accelerates protons. These results differ from the efficiencies derived from simulations of particle acceleration in supernova remnants (Park et al. 2015), but those simulations involve low magnetic field, radiation energy, and particle densities, i.e. very different physical conditions than found in η Carinae. An instrument sensitive in the 1–100 MeV band would be able to discriminate between our model and the one proposed by Ohm et al. (2015).

The simulated pion induced γ-ray light curve and its variability amplitude show a single peak of emission centred at periastron, which is in good agreement with the Fermi LAT observations of the first periastron. The results of the observations of the second periastron are different in that they have a weaker emission. It has been suggested that the change of the X-ray emission after that periastron, where a significant decrease can be observed in Fig. 5, (see also Corcoran et al. 2015b), was the signature of a change of the wind geometry possibly because of cooling instabilities. A stronger disruption or clumpier wind after the second periastron could perhaps induce a decrease of the average wind density and explain that fewer hadronic interactions and fewer thermal emission took place without affecting inverse-Compton emission much.

Fig. 8 Protons luminosity spectra (arbitrary units) at periastron (red; phase 1.0), apastron (blue; phase 0.5), and accelerated on average along the orbit (black).

Figure 8 shows the proton spectra obtained from the simulation at apastron, periastron, and averaged over the orbit. Protons could be accelerated up to 10¹⁵ eV around periastron and reach 10¹⁴ eV on average. The choice of a lower magnetic field reduces those energies at apastron to ~6 × 10¹² eV and ~2 × 10¹² eV and at periastron to ~5.6 × 10¹⁴ eV and ~1.9 × 10¹⁴ eV for 300 G and 100 G, respectively. η Carinae can therefore probably accelerate particles close to the knee of the cosmic-ray spectrum. The spectra and the maximum particle energy depend of course on several assumptions, in particular the magnetic field. The highest energy γ-rays are photo-absorbed and orbital modulation could be expected in the TeV domain. The duration of the periastron bin [0.92–1.06] corresponds to more than 260 days and is longer than the interaction timescale of the protons responsible for the flux variability.

γ-ray observations can probe the magnetic field and shock acceleration in detail, however the quality of the current data above 1 GeV does not yet provide enough information to test hydrodynamical models including detailed radiation transfer (inverse-Compton, pion emission, and photo-absorption). The interplay between disruption and obscuration does not yet account for the X-ray minimum and orbit to orbit variability. More sensitive γ-ray observations will provide a wealth of information and allow us to test the conditions and physics of the shocks at a high level of details, making of η Carinae a perfect laboratory to study particle acceleration in wind collisions.

4. Conclusions

We have used the hydrodynamic simulation of Parkin et al. (2011), which was developed to reproduce the thermal soft X-ray light curve of η Carinae, and have estimated leptonic and hadronic Fermi acceleration, inverse-Compton emission, and neutral pion decay cell by cell assuming a di-polar magnetic field at the surface of the primary star. The results of the simulation were compared with the light curves and spectra observed by Fermi LAT between mid-2008 and mid-2015. We increased the data sample by ~25% with respect to previous analyses and exploited the much better performance of the new PASS8 Fermi-LAT pipeline and of the updated instrument responses. We performed a low-energy and high-energy analysis, from 300 MeV to 10 GeV and from 10 GeV up to 300 GeV, respectively, using the binned (Cash 1979) and unbinned analyses (Mattox et al. 1996). We used different time bins and also performed a low-energy merged analysis combining data with the same orbital phases, when possible, to increase the signal-to-noise ratio. We looked for high and low-energy flux variability of η Carinae and analysed its spectral variations at different orbital phases.

We found a good match between the accuracy of the simulation, even if simplified, and the signal-to-noise of the observations. The comparison between simulation and observations led to several results:

• 1.

  The centroid of the γ-ray source observed by Fermi LAT is compatible with the position of η Carinae within less than 1 arcmin. The low-energy (0.3–10 GeV) γ-ray light curve is modulated along the orbit and shows a very similar and highly significant modulation (5.9σ) during the periastrons of 2009 and 2014, indicating that it is driven by the orbital motion of the system.

• 2.

  Around periastron the low-energy (0.3–10 GeV) γ-ray flux varied by nearly a factor 2 in less than 40 days. A significant fraction of the γ-rays are therefore emitted by a source smaller than the Homunculus Nebula in contrast with the hypothesis of Ohm et al. (2010).

• 3.

  The maximum over the minimum flux ratio observed at low energy (0.3–10 GeV) is 1.53 considering broad phase bin and 1.92 considering bins of 40 days. This matches the results of the simulations assuming that the magnetic field at the surface of the primary star is larger than ~400 G. Smaller values of the magnetic field shorten the volume where electrons could be accelerated to sufficient energies, increase the expected variability amplitude beyond the observed one, and decrease the expected γ-ray luminosity.

• 4.

  A surface magnetic field larger than ~1 kG would produce a secondary peak of emission after periastron that is stronger than the periastron peak, which is not observed. A large part of the secondary peak, observed in the data, is linked with a bubble with reversed wind conditions created after periastron and lasting for about a tenth of the orbit. We note that γ-ray observations together with improved simulations should allow us to constrain the magnetic field in the system even more accurately.

• 5.

  The primary maximum observed just before periastron perfectly matches the prediction of the simulation (amplitude, phase, and duration). The secondary peak occurs slightly earlier and with a lower amplitude than predicted. We assume that these discrepancies come from an inaccurate eccentricity and from the extremely simplified treatment of inverse-Compton scattering. The γ-ray observations should allow us to constrain the eccentricity of the orbit of η Carinae more accurately than possible with current optical observations.

• 6.

  The amplitude and pattern of the low-energy (0.3–10 GeV) γ-ray variability correspond in general very well with the predictions. The luminosity of the pion decay depends on the density and a larger variability is expected. The low-energy γ-rays are therefore very likely to be emitted by inverse-Compton emission in contrast with the claims of Ohm et al. (2015).

• 7.

  The match between the electron distribution predicted by the simulation and the observed cutoff energy, as well as the negligible photon-photon opacity due to the hot shocked gas in the wind collision as computed along different lines of sight, are strong arguments against the scenario suggested by Reitberger et al. (2012).

• 8.

  The γ-ray spectrum observed at apastron shows a discrepancy with the predictions assuming a simplified inverse-Compton treatment. This very likely indicates that the seed soft photon spectrum is not identical everywhere, as currently assumed by the simulations. Spectral variability therefore provides additional constraints on the shock geometry that can be used by more accurate simulations.

• 9.

  The high-energy (>10 GeV) γ-ray component is poorly constrained by the observations. It was well detected during the periastron of 2009, but only weakly detected during the periastron of 2014 and at apastron. The amplitude of variability and the level of the emission however match the expectations for pion decay; inverse-Compton emission is ruled out at such energies, in contrast with the claims of Bednarek & Pabich (2011). Both the high-energy γ-rays and the thermal X-ray emission were weaker during the second periastron, while the inverse-Compton emission was not affected much. This indicates that something peculiar happened in the densest region of the wind collision zone in 2014. Observation of the next periastrons with Fermi-LAT and by the Cherenkov Telescope Array (CTA) are required to probe the high-energy component and the wind and shock geometry further through γ-ray pair conversion.

• 10.

  With the constraints derived on the magnetic field, the simulations predict that η Carinae should be a Pevatron, as this object is able to accelerate protons nearly up to ~10¹⁵ eV. Assuming that for each photon originated via hadronic processes we also have the production of one neutrino, we derive a neutrino flux above 10 TeV that might reach 10^-9 GeV s^-1 cm^-2 on average, which is of the order of the IceCube neutrino sensitivity for several years of observations (Aartsen et al. 2017). Stacking some months of periastron data over many orbits should in principle allow the detection of one PeV neutrino, well above the atmospheric background.

• 11.

  The lack of statistics at high energy does not allow us to constrain any physical information about the hadronic spectrum. But it is evident that a pure leptonic scenario is not able to reproduce the high-energy spectrum observed during the first periastron. A strong γ-ray variability is expected above 100 GeV. Depending on the assumed soft energy photon distribution and the consequent γ-γ absorption at very high energy, η Carinae could be detected by the CTA southern array (including four large size telescopes) at more than 10σ in spectral bins of ΔE/E = 20% for exposures of 50 h, which would be sufficient to measure separately the variability along the orbit of the high-energy component and of photo absorption (Acharya et al. 2013). η Carinae could yield to 10⁴⁸⁻⁴⁹ erg of cosmic-ray acceleration, which is a number close to the expectation for an average supernova remnant (Becker Tjus et al. 2016).

η Carinae is a wonderful laboratory to study particle acceleration in wind collisions. We have demonstrated that the data from Fermi match the simulation expectations, confirming that Fermi acceleration takes place and providing a new tool to diagnose magnetic fields, shock processes, and a complex geometry. Hadronic acceleration is likely but the ultimate proof requires further observations. The evolution of the geometry along the orbit of η Carinae provides a wealth of constraints that future observations and simulations will profit from.

Acknowledgments

We thank Ross Parkin for making the hydrodynamic simulation results available to us and Etienne Lyard for helping us interpret these files.

References

All Tables

All Figures

	Fig. 1 Seven-year high-energy flux light curve of η Carinae obtained from the binned analysis, using the binning as reported in Table 1 (blue points) and with a smaller binning (red points). Error bars are 1σ and superposed upper limits are 95%. For comparison we plotted an arbitrarily rescaled X-ray light curve (grey line).
In the text

Fig. 2 High-energy TS maps for each of the time intervals, corrected for the small differences of exposure times between the 8 phase bins such that the images illustrate the source variability. Each image has the same width (0.77°). Rows a) and b) show TS maps obtained from the binned analysis, respectively including and excluding J1043 from the model. Rows c) and d) represent the same as a) and b) but from the unbinned analysis. The linear colour map spans TS from 0 to 100. The green and purple crosses are the positions of η Carinae and J1043.

In the text

	Fig. 3 Error circles of 1σ for η Carinae derived from the analysis of each time bin. Labels identify the time bin as in Table 1; black circle refers to the whole 7 years of data. The black cross shows the nominal position of η Carinae. As for Fig. 2, the image is in Galactic coordinates with north up and longitude increasing towards the left.
In the text

	Fig. 4 Integrated flux light curve of η Carinae assuming different spectral models: PLEC with three (green) or two (blue) free parameters and BPL with three (red) or two (magenta) free parameters. Error bars are 1σ.
In the text

Fig. 5 Simulated and observed X-ray and γ-ray light curves of η Carinae. The black and purple lines and bins show the predicted inverse-Compton and neutral pion decay light curves. The green and red points show the observed Fermi-LAT light curves at low (0.3–10 GeV) and high (10–300 GeV) energies. The dim grey light curves show the observed (continuous) and predicted (dash, without obscuration) thermal X-ray light curves. Error bars are 1σ.

In the text

	Fig. 6 Merged Fermi LAT analysis (0.3–10 GeV) of the two periastrons for narrow time bins. The two broad bins and the black curve are the same as in Fig. 5.
In the text

Fig. 7 Electrons and photons luminosity spectra at periastron (phase 0.92–1.06; left) and at apastron (phase 0.39–0.59; right). The top panels show the spectra (arbitrary units) of the electrons accelerated in the wind of the primary (blue) and secondary (green) stars and their sum (red). The lower panels show the inverse-Compton emission of both components and the total emission, under the highly simplified assumption that the inverse-Compton parameters (geometry and soft photon spectra) are the same in all cells. The black and grey points are the broadband fluxes derived from Fermi data for the first and the second orbital cycle, respectively. The simulation results were averaged over the orbital phase range corresponding to the periastron observation, as the electron spectra vary quickly during that interval.

In the text

	Fig. 8 Protons luminosity spectra (arbitrary units) at periastron (red; phase 1.0), apastron (blue; phase 0.5), and accelerated on average along the orbit (black).
In the text