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/00046361/201220963  
Published online  02 October 2013 
Constraining the parameter space of the onezone synchrotronselfCompton model for GeVTeV detected BL Lacertae objects
^{1}
LUTH, Observatoire de Paris, CNRS, Université Paris Diderot,
5 place Jules Janssen,
92190
Meudon,
France
^{2}
HarvardSmithsonian Center for Astrophysics, 60 Garden Street,
Cambridge
MA, 02138, USA
email:
matteo.cerruti@cfa.harvard.edu
Received:
20
December
2012
Accepted:
19
May
2013
Context. The onezone synchrotronselfCompton (SSC) model aims to describe the spectral energy distribution (SED) of BL Lac objects via synchrotron emission by a nonthermal population of electrons and positrons in a single homogeneous emission region, partially upscattered to γrays by the particles themselves.
Aims. The model is usually considered as degenerate, given that the number of free parameters is higher than the number of observables. It is thus common to model the SED by choosing a single set of values for the SSC model parameters that provide a good description of the data, without studying the entire parameter space. We present here a new numerical algorithm that permits us to find the complete set of solutions, using the information coming from the detection in the GeV and TeV energy bands.
Methods. The algorithm is composed of three separate steps: we first prepared a grid of simulated SEDs and extracted from each SED the values of the observables; we then parametrized each observable as a function of the SSC parameters; we finally solved the system for a given set of observables. We iteratively solved the system to take into account uncertainties in the values of the observables, producing a family of solutions.
Results. We present a first application of our algorithm to the typical highfrequencypeaked BL Lac object 1RXS J101015.9  311909, provide constraints on the SSC parameters, and discuss the result in terms of our understanding of the blazar emitting region.
Key words: radiation mechanisms: nonthermal / relativistic processes / methods: numerical / BL Lacertae objects: general / galaxies: individual: 1RXS J101015.9311909
© 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 nonthermal 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 radioloud 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 nonthermal emission from the jet, enhanced by relativistic effects. The blazar class is divided into the two subclasses of flatspectrum radio quasars (FSRQs) and BL Lacertae (BL Lac) objects according to the strength of the nonthermal continuum emission relative to the thermal emission from the accretion disc that is partially reprocessed in the broadline 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 infrared and Xrays, 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 highfrequencypeaked BL Lac (HBL) if the peak frequency is located in UV/Xrays, or as lowfrequencypeaked BL Lac (LBL) if the peak frequency is located in infrared/visible light (see Padovani & Giommi 1995). The FSRQs always show a lowfrequency peak. The position of the first peak frequency seems to be inversely correlated with the luminosity (the socalled 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 nonthermal 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 (synchrotronselfCompton model, SSC, Konigl 1981) or with an external photon field (externalinverseCompton, 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 onezone 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 bestfit 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 FermiLAT (Atwood et al. 2009) and groundbased imaging atmospheric Cherenkov telescopes (IACTs), respectively.
Before describing our algorithm, we recall here the basis of the onezone 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 nonthermal population of leptons (e^{±}). The particle energy distribution is described as a broken powerlaw 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 Xray 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 timescale 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 onezone 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):
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 O_{i}, to obtain the solution of our model p_{j}. If the observables O_{i} include an intrinsic uncertainty (σ_{i}), we iteratively solve the system for O_{i} ∈ [O_{i}−σ_{i},O_{i} + σ_{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 × 10^{6}), 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 ultraviolet and Xray 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 Xray 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 KleinNishina 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 extragalactic 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 nontrivial). 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 FermiLAT 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).
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 OpenMP^{2}. 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 Xray 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 powerlaw function, using as minimum and maximum energies the values adopted in the spectral fitting of the FermiLAT and IACT detection. This defines the two indexes Γ_{GeV;TeV}, and the two fluxes (of the fitted powerlaw, 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 5^{6} = 15 625 SEDs. As an indication, the computing time on our 16core 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 software^{3}, 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 highorder 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 software^{4}, 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 Γ_{GeV}^{5}. An additional inequality is added to the system, relating the size and the Doppler factor to the variability timescale τ_{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 O_{i} ∈ [O_{i} − σ_{i},O_{i} + σ_{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 × 12^{5} = 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 chisquare () with respect to the observational data. This check allows us to estimate the true bestfit 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 bestfit value).
Fig. 2 SED of 1RXS J101015.9  311909 (Abramowski et al. 2012b, ; the H.E.S.S. spectrum is represented by the green bowtie and the blue points, the FermiLAT spectrum by the orange bowtie and the red empty circles; SwiftXRT data are shown by the pink crosses, SwiftUVOT 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 bestfit solution with . It is characterised by an extreme value of δ = 96.83, B = 0.015 G, R = 1.3 × 10^{16} cm, α_{1} = 2.0, K′ = 8.94 × 10^{8} cm^{3}, and γ_{br} = 5.31 × 10^{4}. The three different families of solutions, which can be distinguished in the range between 10^{11} and 10^{14} 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 hostgalaxy 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 multiwavelength 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 Xrays, 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 Xray data (case B in Abramowski et al. 2012b), and the Fermi spectrum is evaluated using a 1 GeV lowenergy 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 Xray spectrum measured with SwiftXRT (Γ_{X} = 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 × 10^{4},1.3 × 10^{5}], K′ ∈ [2 × 10^{9},2 × 10^{6}] cm^{3}, and R ∈ [4 × 10^{15},10^{17}] 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 nonlinear, and the fitted expression of the O_{i} 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 onezone 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 secondorder 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 O_{i} − σ_{i} to O_{i} + σ_{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 timescale (τ_{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 Xray and γray data only (i.e. excluding optical/UV and archival data), and the bestfit 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}B^{2}R, and u_{e}u_{B} (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 (topleft 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 timescale; solutions with a lower Doppler factor exist, but would imply a variability timescale 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 timescale 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 ≈10^{17} cm (at the edge of the values used to sample the grid). In Fig. 4 (topright 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.
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: violet1.6, pink1.8, and yellow2.0. 

Open with DEXTER 
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 bestfit 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 timescale higher than 24 h. Bottom left: contours in the γ_{break}–B^{2}R 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 u_{e}–u_{B} plane. The black lines correspond to u_{e}/u_{B} 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 (u_{B} = B^{2}/8π in CGS units) and the particle energy density (u_{e} = m_{e}c^{2}^{∫}dγ γN(γ)). In Fig. 4 (bottomright plot) all the solutions we found are shown in the u_{e}u_{B} plane, indicating that the system is significantly out of equipartition, with a u_{e}/u_{B} ratio comprised between 10 and 200. The evaluation of u_{e} depends on our assumptions on the lowenergy 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 u_{e}, reducing the equipartitionratio. 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 2155304).
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 lowenergy data cannot be used to define a unique α_{1} value because of the uncertainty on the correction by absorption and by the hostgalaxy 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 Xray observations), which is clearly excluded by the data. This might indicate that the homogeneous onezone model is too simple, and that more complex effects should be taken into account, such as nonlinear inverse Compton cooling, energydependent acceleration and escape from the emitting region, nonhomogeneity 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 β_{esc}c^{6}, and lose their energy while emitting synchrotron radiation, the value of γ_{break} can be expressed as (see T98,Eq. (30)) (5)In Fig. 4 (bottomright plot) we show the values of γ_{break} that we found, as a function of B^{2}R, 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 onezone SSC model, that can be successfully applied to any HBL with simultaneous spectral measurements in the Xray, GeV, and TeV energyrange. 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 bestfit 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 bestfit 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 lowenergy 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 preselects 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 onezone 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.
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.
Acknowledgments
The work has been performed using the computation facilities of the Paris Observatory and the HarvardSmithsonian 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
 Abdo, A. A., Ackermann, M., Agudo, I., et al. 2010, ApJ, 716, 30 [NASA ADS] [CrossRef] [Google Scholar]
 Abdo, A. A., Ackermann, M., Ajello, M., et al. 2011a, ApJ, 727, 129 [NASA ADS] [CrossRef] [Google Scholar]
 Abdo, A. A., Ackermann, M., Ajello, M., et al. 2011b, ApJ, 736, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Abramowski, A., Acero, F., Aharonian, F., et al. 2012a, A&A, 539, A149 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Abramowski, A., Acero, F., et al. (H.E.S.S. Collaboration) 2012b, A&A, 542, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Aharonian, F. A., Khangulyan, D., & Costamante, L. 2008, MNRAS, 387, 1206 [NASA ADS] [CrossRef] [Google Scholar]
 Angel, J. R. P., & Stockman, H. S. 1980, ARA&A, 18, 321 [NASA ADS] [CrossRef] [Google Scholar]
 Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071 [NASA ADS] [CrossRef] [Google Scholar]
 Dermer, C. D., & Schlickeiser, R. 1993, ApJ, 416, 458 [NASA ADS] [CrossRef] [Google Scholar]
 Dermer, C. D., Schlickeiser, R., & Mastichiadis, A. 1992, A&A, 256, L27 [NASA ADS] [Google Scholar]
 Finke, J. D., Dermer, C. D., & Böttcher, M. 2008, ApJ, 686, 181 [NASA ADS] [CrossRef] [Google Scholar]
 Fossati, G., Maraschi, L., Celotti, A., Comastri, A., & Ghisellini, G. 1998, MNRAS, 299, 433 [NASA ADS] [CrossRef] [Google Scholar]
 Franceschini, A., Rodighiero, G., & Vaccari, M. 2008, A&A, 487, 837 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ghisellini, G., Tavecchio, F., Foschini, L., & Ghirlanda, G. 2011, MNRAS, 414, 2674 [NASA ADS] [CrossRef] [Google Scholar]
 Inoue, S., & Takahara, F. 1996, ApJ, 463, 555 [NASA ADS] [CrossRef] [Google Scholar]
 Katarzyński, K., Sol, H., & Kus, A. 2001, A&A, 367, 809 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Katarzyński, K., Ghisellini, G., Tavecchio, F., Gracia, J., & Maraschi, L. 2006, MNRAS, 368, L52 [NASA ADS] [CrossRef] [Google Scholar]
 Kaufmann, S., Wagner, S. J., Tibolla, O., & Hauser, M. 2011, A&A, 534, A130 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Konigl, A. 1981, ApJ, 243, 700 [NASA ADS] [CrossRef] [Google Scholar]
 Mankuzhiyil, N., Ansoldi, S., Persic, M., & Tavecchio, F. 2011, ApJ, 733, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Mannheim, K. 1993, A&A, 269, 67 [NASA ADS] [Google Scholar]
 Moore, R. L., & Stockman, H. S. 1981, ApJ, 243, 60 [NASA ADS] [CrossRef] [Google Scholar]
 Mücke, A., & Protheroe, R. J. 2001, APh, 15, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Nieppola, E., Valtaoja, E., Tornikoski, M., Hovatta, T., & Kotiranta, M. 2008, A&A, 488, 867 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Padovani, P., & Giommi, P. 1995, ApJ, 444, 567 [NASA ADS] [CrossRef] [Google Scholar]
 Padovani, P., Giommi, P., & Rau, A. 2012, MNRAS, L423 [Google Scholar]
 Piranomonte, S., Perri, M., Giommi, P., Landt, H., & Padovani, P. 2007, A&A, 470, 787 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Salamon, M. H., & Stecker, F. W. 1998, ApJ, 493, 547 [NASA ADS] [CrossRef] [Google Scholar]
 Sikora, M., Begelman, M. C., & Rees, M. J. 1994, ApJ, 421, 153 [NASA ADS] [CrossRef] [Google Scholar]
 Stein, W. A., Odell, S. L., & Strittmatter, P. A. 1976, ARA&A, 14, 173 [NASA ADS] [CrossRef] [Google Scholar]
 Tavecchio, F., Maraschi, L., & Ghisellini, G. 1998, ApJ, 509, 608 [NASA ADS] [CrossRef] [Google Scholar]
 Urry, C. M., & Padovani, P. 1995, PASP, 107, 803 [NASA ADS] [CrossRef] [Google Scholar]
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.
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 SSCmodel 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 
Coefficients of the fit performed to obtain an expression of ν_{s} (left) and νF_{ν;s} (right) for the case of 1RXS J101015.9  311909.
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.
Coefficients of the fit performed to obtain an expression of Γ_{GeV} for the case of 1RXS J101015.9  311909.
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 SSCmodel parameterspace. 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 × 10^{6} 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 filename, 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 bestfit 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
Summary of the constrained parameters for the SSC modelling of 1RXS J101015.9  311909.
Coefficients of the fit performed to obtain an expression of ν_{s} (left) and νF_{ν;s} (right) for the case of 1RXS J101015.9  311909.
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.
Coefficients of the fit performed to obtain an expression of Γ_{GeV} for the case of 1RXS J101015.9  311909.
Coefficients of the fit performed to obtain an expression of Γ_{TeV} for the case of 1RXS J101015.9  311909.
All Figures
Fig. 1 Flow diagram for the new numerical algorithm for constraining the SSC model parameters. 

Open with DEXTER  
In the text 
Fig. 2 SED of 1RXS J101015.9  311909 (Abramowski et al. 2012b, ; the H.E.S.S. spectrum is represented by the green bowtie and the blue points, the FermiLAT spectrum by the orange bowtie and the red empty circles; SwiftXRT data are shown by the pink crosses, SwiftUVOT 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 bestfit solution with . It is characterised by an extreme value of δ = 96.83, B = 0.015 G, R = 1.3 × 10^{16} cm, α_{1} = 2.0, K′ = 8.94 × 10^{8} cm^{3}, and γ_{br} = 5.31 × 10^{4}. The three different families of solutions, which can be distinguished in the range between 10^{11} and 10^{14} 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 hostgalaxy contribution. 

Open with DEXTER  
In the text 
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: violet1.6, pink1.8, and yellow2.0. 

Open with DEXTER  
In the text 
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 bestfit 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 timescale higher than 24 h. Bottom left: contours in the γ_{break}–B^{2}R 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 u_{e}–u_{B} plane. The black lines correspond to u_{e}/u_{B} values equal to (from bottom to top) 10, 50, and 100, respectively. 

Open with DEXTER  
In the text 
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 SSCmodel 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 (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.