Free Access
Issue
A&A
Volume 558, October 2013
Article Number A47
Number of page(s) 22
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/201220963
Published online 02 October 2013

© ESO, 2013

1. Introduction

Blazars are a class of active galactic nuclei (AGN) characterised by extreme variability, a high degree of polarization, and a strong non-thermal continuum observed in the optical/UV spectrum (Stein et al. 1976; Moore & Stockman 1981). According to the unified AGN model, they are considered to be radio-loud AGN with the relativistic jet pointing in the direction of the observer (Urry & Padovani 1995). Their spectral energy distribution (SED) is thus dominated by the non-thermal emission from the jet, enhanced by relativistic effects. The blazar class is divided into the two subclasses of flat-spectrum radio quasars (FSRQs) and BL Lacertae (BL Lac) objects according to the strength of the non-thermal continuum emission relative to the thermal emission from the accretion disc that is partially reprocessed in the broad-line region (BLR). A blazar is then classified as a BL Lac object if the optical/UV spectrum is dominated by the continuum emission, while it is classified as an FSRQ if emission lines and/or the big blue bump are observed (the standard threshold between the two subclasses being an equivalent width of the emission lines equal to 5 Å, Angel & Stockman 1980).

The blazar SED is composed of two components, the first one peaking between infra-red and X-rays, the second one peaking in γ-rays (see e.g. Abdo et al. 2010). The position of the peak frequency of the first component is used to differentiate among subclasses of BL Lac objects: we can then identify a BL Lac object as high-frequency-peaked BL Lac (HBL) if the peak frequency is located in UV/X-rays, or as low-frequency-peaked BL Lac (LBL) if the peak frequency is located in infrared/visible light (see Padovani & Giommi 1995). The FSRQs always show a low-frequency peak. The position of the first peak frequency seems to be inversely correlated with the luminosity (the so-called blazar sequence, Fossati et al. 1998), even though there is not a general consensus on this point (see e.g. Nieppola et al. 2008; Padovani et al. 2012).

The origin of the first bump is generally ascribed to synchrotron emission from a non-thermal population of electrons and positrons in the emitting region, while the second bump is often ascribed to inverse Compton processes of the same leptons with their own synchrotron emission (synchrotron-self-Compton model, SSC, Konigl 1981) or with an external photon field (external-inverse-Compton, EIC, Dermer et al. 1992; Dermer & Schlickeiser 1993; Sikora et al. 1994) such as the BLR, the accretion disc or the external torus. While the SSC model successfully describes the SED of HBLs, an external Compton component is required for LBLs and FSRQs (see e.g. Ghisellini et al. 2011).

It is worth recalling that other kinds of models exist, in particular hadronic models in which the emitting region also contains relativistic protons, which are responsible (together with the secondary particles coming from p-γ interactions) for the high energy component of the SED (see e.g. Mannheim 1993; Mücke & Protheroe 2001).

The one-zone SSC modelling of HBLs is usually considered to be degenerate: the number of free parameters is higher than the number of observational constraints, and a best-fit solution cannot be given. The standard approach is thus to find one set of parameters that successfully describes the SED, without evaluating the uncertainties on these parameters. In this paper we present a new algorithm to constrain the parameter space of the model, using the information from the GeV and TeV γ-ray emission, as observed by Fermi-LAT (Atwood et al. 2009) and ground-based imaging atmospheric Cherenkov telescopes (IACTs), respectively.

Before describing our algorithm, we recall here the basis of the one-zone SSC model by Katarzyński et al. (2001), on which this study is based. The algorithm can easily be adapted to similar models. The emitting region is considered to be a spherical blob of plasma (with radius R), moving in the relativistic jet with a Doppler factor δ. The source is supposed to be filled with a tangled, homogeneous magnetic field B and a non-thermal population of leptons (e±). The particle energy distribution is described as a broken power-law function (a spectral break is expected in the presence of synchrotron radiation, see e.g. Inoue & Takahara 1996), defined by the two indexes α1 and α2; the three Lorentz factors γmin, γbreak, and γMax; and the normalization factor K. The free parameters of the model are then nine: three for the emitting region (R, δ, and B) and six for the particle population (γmin;break;Max, α1;2, and K). In general there are less than nine observables to constrain the parameters of the model, which remains thus degenerate. The basis of the constraints on the SSC model are described in the work by Tavecchio et al. (1998, hereafter T98. The idea is to determine a system of equations linking parameters and observables, and to solve it analytically. In their work, T98 used six observables: the frequency and luminosity of the synchrotron peak, the frequency and luminosity of the Compton peak, and the measured X-ray and γ-ray spectral slopes (ΓX, γ, linked directly to α1 and α2 through α1,2 = 2Γγ, X−1). The number of free parameters was reduced from nine to seven by considering a very low value of γmin and a very high value of γMax; additional constraints were added from information about the variability time-scale and from the requirement that the break Lorentz factor be consistent with the one expected from synchrotron cooling.

In this paper, we present a new algorithm to find the complete set of solutions for the interpretation of SEDs of HBLs with a one-zone SSC model. Our method is based on the approach by T98, but is fully numerical. In the next section we provide the details of our algorithm, showing in Sect. 3 a first application to the HBL 1RXS J101015.9 - 311909.

2. Description of the algorithm

Our algorithm can be seen as a numerical extension of the work done by T98. The basic idea is to define a set of equations linking SSC parameters and observables that, in our case, are obtained numerically. The algorithm is composed of three separate steps (summarized in Fig. 1):

thumbnail Fig. 1

Flow diagram for the new numerical algorithm for constraining the SSC model parameters.

Open with DEXTER

  • we simulate SEDs for a grid of parameter values (spanningthe region of interest of the SSC parameters p j ), and for each SED we compute the value of the corresponding observables (O i);

  • we then parametrize each observable O i as a function of the SSC parameters (1) where α j is the result of a fit;

  • we finally solve the set of equations for a given set of observables Oi, to obtain the solution of our model pj. If the observables Oi include an intrinsic uncertainty (σi), we iteratively solve the system for Oi ∈ [Oiσi,Oi + σi], thus producing a set of possible solutions.

We reduce the number of free parameters from nine to seven by fixing a reasonably low and high value for γmin (fixed to 100) and γMax (fixed to 5 × 106), respectively.1 The main difference with respect to the work by T98 lies in the choice of the observables. The synchrotron peak is relatively well constrained. Even when it is located in the observational gap between the ultra-violet and X-ray bands, its position and luminosity can be reasonably well estimated by extrapolating the low and high energy data. The position of the inverse Compton peak is, however, much more uncertain, and so we decided to use the actually observed GeV and TeV spectral slopes, and their flux normalizations, defined in the following as ΓGeV, ΓTeV, νFν;GeV, and νFν;TeV. The number of observables we use is thus seven: together with the four γ-ray observables, we consider the frequency and intensity of the synchrotron peak (νs and νFν;s), and the measured X-ray spectral slope (ΓX). We thus have seven free parameters and seven independent observables, and the problem can be solved.

It is important to underline that the four γ-ray observables are not degenerate: the shape of the inverse Compton component is not symmetrical with respect to its peak, given the fact that the TeV part is affected by the transition to the Klein-Nishina regime, the internal absorption by γ-γ pair production (which depends on the synchrotron emission, see e.g. Aharonian et al. 2008) and the absorption on the extra-galactic background light (EBL, Salamon & Stecker 1998). It is clear that a study of this kind is carried out more easily with a numerical approach (a purely theoretical derivation of the expression of the TeV spectral slope assuming all the effects described above would be non-trivial). Another improvement with respect to the work by T98 is that we do not use the simple approximation that the GeV slope is uniquely related to the electron spectral index below the break (when γmin is reasonably low): this approximation is true only if the inverse Compton component accounts for the GeV spectrum well before the peak; when the Compton peak is located at lower energies, the spectrum measured with Fermi-LAT is significantly softer.

In the SSC code used for this study (Katarzyński et al. 2001), a modification has been introduced concerning the definition of the normalization factor K of the electron distribution, which is originally defined as the number density at γ = 1. With this definition, the value of K depends strongly on the value of α1: in particular, assuming that the model correctly fits the synchrotron peak, but we want to modify α1, this imposes a correction of the value of K as well. For a simulation of a multitude of SEDs, this definition is not appropriate, imposing a huge range of values for K. We thus redefine K as the number density at γ = γbreak, (2)In this way, if the model correctly describes the synchrotron peak, one can modify the value of α1 without affecting the other parameters and, more importantly, the range of values of K′ to study becomes narrower. As another small modification, the EBL absorption is computed using the model by Franceschini et al. (2008).

Table 1

Summary of the constrained parameters for the SSC modelling of 1RXS J101015.9 - 311909.

In order to minimize the computing time, the code has been parallelized using OpenMP2. The parallelization is not at the level of the computation of the synchrotron and inverse Compton components, but at the higher level of the sampling of the parameter space. The distribution of the input parameters per thread is done dynamically, and the communication between threads takes place only at the end of each computation when the results are merged.

Iterating on the values of δ, B, K′, R, γbreak, and α1, SEDs have been simulated for a grid of parameter values (grid production in Fig. 1). The value of α2 is fixed, constrained by the measured X-ray spectral slope ΓX, and equal to 2ΓX − 1. For each modelled SED, we first identify the position of the synchrotron peak and its intensity. To compute the expected fluxes and spectral slopes detected in the GeV and TeV ranges, we simply fit the modelled SED in the range of the GeV and TeV detection with a power-law function, using as minimum and maximum energies the values adopted in the spectral fitting of the Fermi-LAT and IACT detection. This defines the two indexes ΓGeV;TeV, and the two fluxes (of the fitted power-law, not of the model) νFν;GeV;TeV at the respective decorrelation energies (which are measured quantities).

For the application that will be presented in the next section, we sampled five different values for each free parameter, resulting in the computation of 56 = 15   625 SEDs. As an indication, the computing time on our 16-core machine is roughly one hour. The result of this stage is a grid containing the corresponding set of observables for each set of parameters.

This grid is then fitted using the six observables as dependent variables, and the six free parameters as independent variables (grid fitting in Fig. 1). While for the synchrotron peak frequency and flux the relation between observables and parameters is simple (and consistent with the analytical expression), for the fluxes and slopes measured in the GeV and TeV energy bands, it becomes more complicated. In particular, the simple relation given in Eq. (1) no longer applies, and we need to consider a more complex relation of the form (3)More explicitly, we have to perform a fit considering all the possible polynomials of our parameters.

The problem is then to find a parametrization for a dependent variable, which is a function of six parameters, without knowing it a priori (nonparametric regression). This problem has been solved using the root software3, namely the TMultiDimFit class. This class provides, for each term of the fit function, the fitted coefficient and the weight of the term in the overall fit, providing a listing of the different terms in order of relative importance.

The computing time for this step is a few seconds for the synchrotron observables and for νFν;GeV;TeV, while for the GeV and TeV slopes it can be significantly longer, taking up to an hour each if we want to study all the possible polynomial functions (in order to be sure not to miss a high-order polynomial which might play an important role in the fit). The choice of the last polynomial considered in the fit is a free parameter of the algorithm: it is chosen by listing the fit terms according to their contribution to the χ2, and defining a threshold above which the contribution of further terms becomes negligible. The goodness of the fit has thus to be verified, as shown in Fig. A.1 for the application presented in the next section.

Once we have obtained the six equations relating the observables to the free parameters, the final stage of the algorithm is to solve this system of equations (system solving in Fig. 1). This task is carried out with the Mathematica software4, using the numerical FindInstance command, and the solutions have been reduced to the real domain. The computing time for this step depends strongly on the form of the equations in our system. In order to reduce it, we have fixed one of the parameters (α1 for the application presented in the next section), and searched the solutions of the system for different given values of this frozen parameter. For each given value, we solve a system of five equations plus two inequalities, corresponding to the minimum and the maximum of ΓGeV5. An additional inequality is added to the system, relating the size and the Doppler factor to the variability time-scale τvar, (4)where c is the speed of light and z is the redshift of the source.

To take into account the uncertainty on the five remaining observables, we iterate the solution of the system spanning the range Oi ∈ [Oi − σi,Oi + σi] to produce a set of solutions. We sampled three different values of the frozen parameter, and we span the five observables 12 times each, leading to 3 × 125 = 746   496 systems studied. The Mathematica code has been parallelized as well, and run on a computer grid, using up to 200 cores simultaneously. For the expressions used for the application presented in the next section, the computing time is roughly nine hours.

The three steps of the analysis (grid production, grid fitting, and system solving, summarized in the plot shown in Fig. 1) can then take several hours, and this imposes a compromise between the resolution adopted for the grid (which affects the computing time of steps one and two), the number of terms used in the fit (which affects the computing time of steps two and three), and the number of iterations done in the system solution (which affects the computing time of step three).

The sets of parameters found with our algorithm are then used to produce contour plots, showing the confidence regions as a function of two different parameters, as in Fig. 4. As a final step, we compute all the SSC models corresponding to the found solutions, and we evaluate for each of them the reduced chi-square () with respect to the observational data. This check allows us to estimate the true best-fit solution (the one with the minimum value), and to verify that the contour plots actually correspond to 1-σ contours (i.e. that all the solutions are characterised by with respect to the best-fit value).

thumbnail Fig. 2

SED of 1RXS J101015.9 - 311909 (Abramowski et al. 2012b, ; the H.E.S.S. spectrum is represented by the green bow-tie and the blue points, the Fermi-LAT spectrum by the orange bow-tie and the red empty circles; Swift-XRT data are shown by the pink crosses, Swift-UVOT data by the red stars, ATOM data by the blue open boxes, and archival data from the NED in grey). All the SSC models which describe the SED, as found with our algorithm, are plotted in grey, while the solid black curve represents the best-fit solution with . It is characterised by an extreme value of δ = 96.83, B = 0.015 G, R = 1.3 × 1016 cm, α1 = 2.0, K′ = 8.94 × 10-8 cm-3, and γbr = 5.31 × 104. The three different families of solutions, which can be distinguished in the range between 1011 and 1014 Hz, correspond to α1 = 1.6, 1.8, and 2.0, as discussed in Sect. 3. The infrared and visible data can be reproduced by taking into account the host-galaxy contribution.

Open with DEXTER

3. Application: 1RXS J101015.9 - 311909

In this section we present a first application of our numerical algorithm to the typical HBL 1RXS J101015.9 - 311909. The TeV emission from this source, located at a redshift of 0.143 (Piranomonte et al. 2007), has been detected with H.E.S.S. (Abramowski et al. 2012b). The multi-wavelength study presented in the H.E.S.S. discovery paper shows that there are two major systematic uncertainties in the SED of this object. The first one comes from a possible absorption effect in X-rays, affecting the evaluation of the position of the synchrotron peak; the second one comes from a problem in the determination of the Fermi spectrum at low energies due to a possible contamination by diffuse γ-ray emission from the Galactic foreground. In the following we consider that the synchrotron peak is located at an energy between the available UV and X-ray data (case B in Abramowski et al. 2012b), and the Fermi spectrum is evaluated using a 1 GeV low-energy threshold. This does not represent a limitation of the algorithm, which can be successfully applied to the other cases discussed by Abramowski et al. (2012b) as well.

The value of α2 is uniquely constrained by the X-ray spectrum measured with Swift-XRTX = 2.5 ± 0.1) and it has been fixed at 4.0. The grid of simulated SEDs has been produced sampling the other six free parameters in the ranges: δ ∈ [20,100], B ∈ [0.005,0.05] G, α1 ∈ [1.5,2.2], γbreak ∈ [3 × 104,1.3 × 105], K′ ∈ [2 × 10-9,2 × 10-6] cm-3, and R ∈ [4 × 1015,1017] cm. As laid out in the previous section, the grid has been produced sampling five different values for each parameter, logarithmically spaced (apart from α1 which has been sampled linearly) between the minimum and the maximum values. It is important to underline that solutions found outside of the sampled grid are not considered, because the system of equations is non-linear, and the fitted expression of the Oi is appropriate only for the sampled parameter space. This could be seen as a limit of our approach, in the sense that we preselect a fraction of the parameter space. However, the analytical constraints defined by T98 allow us to beforehand exclude regions of the parameter space that are not accessible in the one-zone SSC model. In particular we decided not to study values of δ higher than 100, as this would require extreme values of the bulk Lorentz factor of the emitting region.

The result of the parametrization of the six observables is presented in Tables A.1 (for νs and νFν;s), A.2 (for νFν;GeV,TeV), A.3 (for ΓGeV), and A.4 (for ΓTeV). As can be seen, while the expressions for νs and νfν;s are simple, and consistent with what is expected from analytical considerations, the expressions for the Fermi and H.E.S.S. observables are more complicated, and we are obliged to consider second-order terms in the fitting polynomial with the extreme case of ΓTeV, which is described satisfactorily only when using more than one hundred terms. In Fig. A.1 we show the relation between the values of each observable, and the reconstructed values (after fit): in a perfect fit, the values should follow exactly the linear relation represented by the thick blue line.

The system of equations was then solved for three different values of α1 (equal to 1.6, 1.8, and 2.0) and for twelve different values of each observable (apart from ΓGeV, which is defined through two inequalities) reaching from Oi − σi to Oi + σi. The range of values for α1 corresponds to expectations from standard acceleration scenarios. The values adopted for the observables are (in logarithm) νs ∈ [16.15,16.25], νfν;s ∈ [−10.74, −10.68], νfν;GeV ∈ [−11.98, −11.78], νfν;TeV ∈ [−12.36, −12.10], ΓGeV ∈ [−1.94, −1.48], and ΓTeV ∈ [−3.55, −2.61]. For the γ-ray observables, the considered error does include the systematic error on the observable (summed in quadrature to the statistic error). The system of equations includes an inequality for the variability time-scale (τvar = 24 h) plus a limit on the value of the Doppler factor (δ < 100, which corresponds to the limit of the sampled grid).The results of our algorithm defining the range of SSC parameters which correctly fits the data are shown in Table 1 and the corresponding histograms are shown in Fig. 3.

To compare our solutions with the observational data, in Fig. 2 we show the SED of 1RXS J101015.9 - 311909, together with all the modelled SEDs found with our algorithm. The of each modelled SED is estimated taking into account X-ray and γ-ray data only (i.e. excluding optical/UV and archival data), and the best-fit solution, shown in black in Fig. 2, has .

In Fig. 4 we show the contour plots corresponding to our solutions for four different pairs of parameters: B-δ, R-δ, γbr-B2R, and ue-uB (the particle and magnetic field energy densities). The different sets of solutions that can be seen in the contour plots are the result of the system of equations solved for discrete values of the observables in the range ± 1σ. The only contour which is statistically significant is the most extended one, that represents the 1-σ confidence region.

Our algorithm, unlike the standard modelling approach, permits us to provide general, quantitative constraints on the physical parameters of the emitting region. These are discussed in the following for the case of 1RXS J101015.9 - 311909.

3.1. Doppler factor

The first important result is that we can define a minimum Doppler factor that solves the system: δmin ≈ 31. In Fig. 4 (top-left plot) we show the contour plots of all our solutions in the B-δ plane, and we compare them to the analytical constraints computed following Tavecchio et al. (1998). It becomes clear that with our numerical approach we significantly narrow down the accessible parameter space.

The constraints on δ come mainly from the TeV slope and the variability time-scale; solutions with a lower Doppler factor exist, but would imply a variability time-scale higher than the one observed. In particular, comparing our result with the one presented by Abramowski et al. (2012b), the fact that they provide a solution with δ = 30 is due to a variability time-scale slightly higher than what has been assumed here (τvar ≈ 25 h) and to a rather soft TeV slope of the SSC model.

3.2. Size of the emitting region

The size of the emitting region is limited by the variability time scale (Eq. (4)), following the usual considerations in SSC modelling. The maximum allowed size is ≈1017 cm (at the edge of the values used to sample the grid). In Fig. 4 (top-right plot) we show the contour plots in the R-δ plane, and the limit corresponding to Eq. (4). As discussed above, solutions with a larger size of the emitting region and a lower Doppler factor exist, but violate the variability constraint. For the case of 1RXS J101015.9 - 311909, the variability constraint is relatively weak (no flaring behaviour observed), but for rapidly flaring sources this would permit us to constrain the parameter space even more.

thumbnail Fig. 3

Histograms showing the values of the SSC model parameters for the case of 1RXS J101015.9 - 311909. From top to bottom, and left to right, we show the distribution of the solutions for δ, B (in G), R (in cm), K′ (in cm-3), γbreak, and α1. We note that the value of α1 has been frozen and studied for three different cases (1.6, 1.8, and 2.0). The three colours represent the different solutions for the three values of α1 studied: violet-1.6, pink-1.8, and yellow-2.0.

Open with DEXTER

thumbnail Fig. 4

Contour plots of the solutions found for 1RXS J101015.9 - 311909. The contours are expressed in arbitrary units representing the number of solutions corresponding to a given pair of parameters. White represents zero. As discussed in Sect. 2, the most extended contour corresponds to a 1-σ contour with respect to the best-fit solution, as determined a posteriori. Top left: contours in the Bδ plane. The shadowed regions represent the exclusion regions defined following T98 (the region filled with diagonal lines is computed from Eq. (4), while the region filled with dots is computed from Eq. (11), considering an emitting region size equal to (1 + z)cτvar/δ and (1 + z)cτvar/30δ). We do not include the constraint on the consistency of γbreak with synchrotron cooling). Top right: contours in the Rδ plane. The filled region corresponds to a variability time-scale higher than 24 h. Bottom left: contours in the γbreakB2R plane. The black lines represent the expected value of γbreak in the presence of pure synchrotron cooling, assuming a value of βesc equal to (from left to right) 1/100, 1/50, and 1/6. Bottom right: contours in the ueuB plane. The black lines correspond to ue/uB values equal to (from bottom to top) 10, 50, and 100, respectively.

Open with DEXTER

3.3. Energy budget

An open question in blazar physics concerns the energy budget of the emitting region. For each solution we compute the magnetic energy density (uB = B2/8π in CGS units) and the particle energy density (ue = mec2dγ γN(γ)). In Fig. 4 (bottom-right plot) all the solutions we found are shown in the ue-uB plane, indicating that the system is significantly out of equipartition, with a ue/uB ratio comprised between 10 and 200. The evaluation of ue depends on our assumptions on the low-energy part of the particle population (which is not constrained by the data), and in particular on the value of γmin. A higher value of γmin would lower the value of ue, reducing the equipartition-ratio. This effect has not been tested here, because high values of γmin affect the evaluation of the GeV slope, and would require a dedicated grid of SED model curves. However, for values of α1 lower than 2.0, we do not expect that high values of γmin affect significantly this result, while it should be important when α1 ≥ 2.0.

The fact that the SSC modelling of HBLs requires that the emitting region be out of equipartition is consistent with the results obtained by several authors (see e.g. Abdo et al. (2011b) for the case of Mrk 421, Abdo et al. (2011a) for the case of Mrk 501, and Abramowski et al. (2012a) for the case of PKS 2155-304).

3.4. Particle energy distribution

The study of the solutions for different values of α1 reveals that the slope of the electron distribution is not as well constrained by the GeV data as it is generally assumed to be (and as in T98). The solutions for the different α1 values are quite similar. An additional constraint on the values of the index of the particle population can be provided by optical/UV measurements. In the case of 1RXS J101015.9 - 311909, however, the low-energy data cannot be used to define a unique α1 value because of the uncertainty on the correction by absorption and by the host-galaxy contamination (see Abramowski et al. 2012b, and Fig. 2). For this source we can use the optical and UV data only as upper limits, excluding slopes softer than 2.2. For other blazars, with a better estimate of the AGN emission in this part of the spectrum, a more precise value of α1 can be determined.

The values of α1 that we have considered point to one of the open questions of blazar physics: synchrotron cooling predicts a spectral break equal to one (i.e. α1 = 3.0, for α2 = 4.0, as constrained by the X-ray observations), which is clearly excluded by the data. This might indicate that the homogeneous one-zone model is too simple, and that more complex effects should be taken into account, such as non-linear inverse Compton cooling, energy-dependent acceleration and escape from the emitting region, non-homogeneity of both particles, and magnetic field. This justifies the choice of considering the two particle slopes α1 and α2 as independent parameters.

If the spectral break were a result of pure synchrotron cooling, the energy of the break in the particle distribution could be computed theoretically, and linked to the other free parameters. Assuming that particles are injected in the emitting region, escape from it at a speed βescc6, and lose their energy while emitting synchrotron radiation, the value of γbreak can be expressed as (see T98,Eq. (30)) (5)In Fig. 4 (bottom-right plot) we show the values of γbreak that we found, as a function of B2R, comparing them to the theoretical expectations. As can be seen, our solutions are consistent with a pure synchrotron cooling only if the escape speed is lower than c/6. This value can be increased by considering additional cooling terms, such as the energy loss by inverse Compton scattering. This point strengthens the fact that the observed spectral break is not consistent with pure synchrotron cooling, and that more complex effects have to be considered to explain the particle spectrum in the stationary state.

In this application we have fixed the values of γmin and γMax, which cannot be constrained in the particular case of 1RXS J101015.9 - 311909. For other sources with different observational constraints it would be possible to fix other parameters (such as α1, in the presence of precise measurements at lower frequencies) and to study γmin or γMax. In particular, this approach would be interesting for some other HBLs, for which a value of γmin higher than the one that has been assumed here is required in order to describe the SED with a SSC model (see e.g. Katarzyński et al. 2006; Kaufmann et al. 2011).

4. Conclusions

In this paper we have described a new algorithm to constrain the parameter space of the one-zone SSC model, that can be successfully applied to any HBL with simultaneous spectral measurements in the X-ray, GeV, and TeV energy-range. The algorithm follows the idea developed by T98 (i.e. the definition of a set of equations linking the free model parameters and the blazar observables), but is fully numerical, introduces as new observables the properties of the GeV and the TeV detection, and allows the derivation of a new set of equations.

The algorithm cannot be applied to the subclasses of LBLs and FSRQs because for these sources an external inverse Compton component is required, leading to a higher number of free parameters in the model.

As a first test, we have applied our algorithm to the typical HBL 1RXS J101015.9 - 311909 and derived the set of solutions that correctly describes its SED. It is important to underline that the range of solutions found here is much narrower than the one obtained with the analytical approach, showing that the information coming from the GeV–TeV detection permits the model to be better constrained. In this way, a true best-fit solution is provided.

Other numerical algorithms have been proposed in the literature to constrain the SSC model parameter space. Both Finke et al. (2008) and Mankuzhiyil et al. (2011) have proposed a χ2 minimization algorithm to fit the SEDs of BL Lac objects, showing as well that, thanks to the TeV detection, it is possible to obtain a best-fit solution. Neither of these two methods make use, however, of the recent Fermi data. On the other hand, their fits are strongly constrained by the low energy (optical and UV) photons, implicitly assuming that the γ-ray emitting region is responsible for the totality of this flux. Our method, based on the properties of the Fermi detection and the position of the synchrotron peak, allows us to relax this hypothesis, and to model more complex scenarios in which the low-energy flux is contaminated by emission from other, farther regions of the jet, or the host galaxy. It is worth noting as well that the method described by Finke et al. (2008) introduces the requirement that the energy budget of the emitting region is minimized, and pre-selects a set of solutions.

Several improvements to the algorithm we have presented here can be considered. First of all, the production of the grid of modelled SEDs, as was discussed here for the application to 1RXS J101015.9 - 311909, limits the parameter space: a larger grid can be produced, extending the region of the parameter space under investigation. To reduce the computing time we have also frozen the value of α1, solving the system for different, given values of this parameter; however, the system can, in principle, be solved leaving all the parameters free to vary, studying α1 as a free parameter as well.

In most of the current publications, the standard approach in blazar physics is still to consider the one-zone SSC model as degenerate (i.e. the number of free parameters is higher than the number of constraints), and a fit of the SED is usually not performed. It is thus common to explore only a few solutions to decide whether the data can be described in a satisfactory way, without studying the entire parameter space nor evaluating the errors on the parameters. We have shown here that, at least for GeV–TeV HBLs, it is possible to explore the parameter space more systematically and to significantly improve the constraints on the parameters of the model. Our algorithm can be applied to the continuously increasing sample of HBLs detected in γ-rays, and will help us improve our comprehension of this extreme class of blazars.


1

These values are consistent with the ones adopted by T98. For the particular case of 1RXS J101015.9 - 311909 presented in the next section, values of γmin lower than 100 would overestimate the radio measurements.

5

There is no constraint on the choice of the frozen parameter (α1 in this case), nor on the choice of the equation which is turned into two inequalities (ΓGeV in this case) to maintain the same number of equations and variables.

6

The particle escape term can be considered as well as an adiabatic term (i.e. associated to the expansion of the emitting region).

Acknowledgments

The work has been performed using the computation facilities of the Paris Observatory and the Harvard-Smithsonian Center for Astrophysics. The authors wish to thank Hélène Sol for fruitful discussions, as well as the anonymous referee for the comments and remarks which improved the present work.

References

Appendix A: Tables and figures containing the fit results

We present here the results of the fitting of the grid for the particular case of 1RXS J101015.9 - 311909. For each observable we provide the terms of the fit listed by increasing polynomial order of the parameters composing each term. The plot showing the goodness of each fit (comparison between the sampled and the reconstructed values) are given in Fig. A.1.

thumbnail Fig. A.1

Comparison between the sampled values of the six observables considered in this study and their reconstructed values, expressed as a function of the SSC-model parameters, for the case of 1RXS J101015.9 - 311909. In a perfect fit, the points would follow a linear relation (thin solid line). The six subplots are in the order (from top to bottom, and left to right), the synchrotron peak frequency (νs), the synchrotron peak flux (expressed as νFν;s), the Fermi flux (measured at the decorrelation energy and expressed as νFν;GeV), the H.E.S.S. flux (measured at the decorrelation energy and expressed as νFν;TeV), the measured Fermi photon index (ΓGeV), and the measured H.E.S.S. photon index (ΓTeV).

Open with DEXTER

Table A.1

Coefficients of the fit performed to obtain an expression of νs (left) and νFν;s (right) for the case of 1RXS J101015.9 - 311909.

Table A.2

Coefficients of the fit performed to obtain an expression of νFν;GeV (left) and νFν;TeV (right) for the case of 1RXS J101015.9 - 311909.

Table A.3

Coefficients of the fit performed to obtain an expression of ΓGeV for the case of 1RXS J101015.9 - 311909.

Table A.4

Coefficients of the fit performed to obtain an expression of ΓTeV for the case of 1RXS J101015.9 - 311909.

Appendix B: Mathematica code used in this work

We present here the details of the Mathematica code used to determine and solve the system of equations to constrain the SSC-model parameter-space. The code takes as input a grid containing the information of the simulated SSC model and the corresponding values of the observables. For the specific case of 1RXS J101015.9 - 311909 we fixed the values of γmin = 100; γMax = 5 × 106 and α2 = 4.0, and the grid contains only the information on the remaining free parameters. As an example, here is a part of the resulting grid used for the case of 1RXS J101015.9 - 311909 (all the values are expressed as log ; the grid file-name, as read by the Mathematica code is constraints_1RXSJ1010_log.dat).

α 1 γ br K B δ R ν syn νFν;syn νFν;GeV νFν;TeV ΓGeV ΓTeV
0.176091 4.5 -8.69897 -2.30103 1.30103 15.6021 14.6667 -18.3045 -22.1208 -23.4220 -2.15305 -3.79803
0.176091 4.5 -8.69897 -2.30103 1.47577 15.6021 14.7778 -17.6049 -21.4126 -22.5768 -2.08191 -3.75843
0.176091 4.5 -8.69897 -2.30103 1.65052 15.6021 15.0000 -16.9060 -20.7162 -21.7401 -2.01315 -3.71274
0.176091 4.5 -8.69897 -2.30103 1.82526 15.6021 15.1111 -16.2074 -20.0308 -20.9150 -1.94735 -3.65919
0.176091 4.5 -8.69897 -2.30103 2.00000 15.6021 15.3333 -15.5078 -19.3558 -20.1017 -1.88443 -3.59492
0.176091 4.5 -8.69897 -2.05105 1.30103 15.6021 14.8889 -17.8037 -21.6272 -22.8920 -2.08529 -3.84209
0.176091 4.5 -8.69897 -2.05105 1.47577 15.6021 15.0000 -17.1059 -20.9311 -22.0402 -2.00853 -3.79903
0.176091 4.5 -8.69897 -2.05105 1.65052 15.6021 15.2222 -16.4057 -20.2470 -21.1980 -1.93690 -3.74820
0.176091 4.5 -8.69897 -2.05105 1.82526 15.6021 15.4444 -15.7077 -19.5739 -20.3691 -1.87098 -3.68696
0.176091 4.5 -8.69897 -2.05105 2.00000 15.6021 15.5556 -15.0079 -18.9109 -19.5539 -1.81029 -3.61265
0.176091 4.5 -8.69897 -1.80104 1.30103 15.6021 15.1111 -17.3037 -21.1522 -22.3739 -2.01313 -3.89216

• Reading the grid containing model parameters and observables

• Parametrizing the observables as a function of the model parameters

For each observable, the best-fit formula has been estimated using the TMultiDimFit class in root (see http://root.cern.ch/root/html/TMultiDimFit.html for details).

(1) νsyn

(2) Fν;syn

(3) Fν;GeV

(4) Fν;TeV

(5) ΓGeV

(6) ΓTeV

• Solving the system of equations

(1) Solutions for α1=2.0

(2) Solutions for α1=1.8

(3) Solutions for α1=1.6

• Histogram plotting

• Exporting the solutions as text file

All Tables

Table 1

Summary of the constrained parameters for the SSC modelling of 1RXS J101015.9 - 311909.

Table A.1

Coefficients of the fit performed to obtain an expression of νs (left) and νFν;s (right) for the case of 1RXS J101015.9 - 311909.

Table A.2

Coefficients of the fit performed to obtain an expression of νFν;GeV (left) and νFν;TeV (right) for the case of 1RXS J101015.9 - 311909.

Table A.3

Coefficients of the fit performed to obtain an expression of ΓGeV for the case of 1RXS J101015.9 - 311909.

Table A.4

Coefficients of the fit performed to obtain an expression of ΓTeV for the case of 1RXS J101015.9 - 311909.

All Figures

thumbnail Fig. 1

Flow diagram for the new numerical algorithm for constraining the SSC model parameters.

Open with DEXTER
In the text
thumbnail Fig. 2

SED of 1RXS J101015.9 - 311909 (Abramowski et al. 2012b, ; the H.E.S.S. spectrum is represented by the green bow-tie and the blue points, the Fermi-LAT spectrum by the orange bow-tie and the red empty circles; Swift-XRT data are shown by the pink crosses, Swift-UVOT data by the red stars, ATOM data by the blue open boxes, and archival data from the NED in grey). All the SSC models which describe the SED, as found with our algorithm, are plotted in grey, while the solid black curve represents the best-fit solution with . It is characterised by an extreme value of δ = 96.83, B = 0.015 G, R = 1.3 × 1016 cm, α1 = 2.0, K′ = 8.94 × 10-8 cm-3, and γbr = 5.31 × 104. The three different families of solutions, which can be distinguished in the range between 1011 and 1014 Hz, correspond to α1 = 1.6, 1.8, and 2.0, as discussed in Sect. 3. The infrared and visible data can be reproduced by taking into account the host-galaxy contribution.

Open with DEXTER
In the text
thumbnail Fig. 3

Histograms showing the values of the SSC model parameters for the case of 1RXS J101015.9 - 311909. From top to bottom, and left to right, we show the distribution of the solutions for δ, B (in G), R (in cm), K′ (in cm-3), γbreak, and α1. We note that the value of α1 has been frozen and studied for three different cases (1.6, 1.8, and 2.0). The three colours represent the different solutions for the three values of α1 studied: violet-1.6, pink-1.8, and yellow-2.0.

Open with DEXTER
In the text
thumbnail Fig. 4

Contour plots of the solutions found for 1RXS J101015.9 - 311909. The contours are expressed in arbitrary units representing the number of solutions corresponding to a given pair of parameters. White represents zero. As discussed in Sect. 2, the most extended contour corresponds to a 1-σ contour with respect to the best-fit solution, as determined a posteriori. Top left: contours in the Bδ plane. The shadowed regions represent the exclusion regions defined following T98 (the region filled with diagonal lines is computed from Eq. (4), while the region filled with dots is computed from Eq. (11), considering an emitting region size equal to (1 + z)cτvar/δ and (1 + z)cτvar/30δ). We do not include the constraint on the consistency of γbreak with synchrotron cooling). Top right: contours in the Rδ plane. The filled region corresponds to a variability time-scale higher than 24 h. Bottom left: contours in the γbreakB2R plane. The black lines represent the expected value of γbreak in the presence of pure synchrotron cooling, assuming a value of βesc equal to (from left to right) 1/100, 1/50, and 1/6. Bottom right: contours in the ueuB plane. The black lines correspond to ue/uB values equal to (from bottom to top) 10, 50, and 100, respectively.

Open with DEXTER
In the text
thumbnail Fig. A.1

Comparison between the sampled values of the six observables considered in this study and their reconstructed values, expressed as a function of the SSC-model parameters, for the case of 1RXS J101015.9 - 311909. In a perfect fit, the points would follow a linear relation (thin solid line). The six subplots are in the order (from top to bottom, and left to right), the synchrotron peak frequency (νs), the synchrotron peak flux (expressed as νFν;s), the Fermi flux (measured at the decorrelation energy and expressed as νFν;GeV), the H.E.S.S. flux (measured at the decorrelation energy and expressed as νFν;TeV), the measured Fermi photon index (ΓGeV), and the measured H.E.S.S. photon index (ΓTeV).

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.