text:
# 1 Introduction

## 1 Introduction

Deep inelastic electron-proton scattering experiments at HERA have provided measurements of the inclusive structure function $`F_2(x,Q^2)`$ in very small values of the Bjorken variable $`x`$ ($`10^2`$ down to $`10^5`$). In these processes the proton target is analyzed by a hard probe with virtuality $`Q^2=q^2`$, where $`xQ^2/2p.q`$ and $`p,q`$ are the four-momenta of the incoming proton and the virtual photon probe. In that kinematical region, the gluon is the leading parton driving the small $`x`$ behavior of the deep inelastic observables. The small $`x`$ region is described in a formal way using the summation of gluon ladder diagrams, whose virtual contributions lead to the gluon reggeization. Such ladders diagrams are associated with the Pomeron, the leading Reggeon with the vacuum quantum numbers, which was introduced phenomenologically to describe the high energy behavior of the total and elastic cross-sections of the hadronic reactions and connected with the existence of large rapidity gaps in the produced final state .

At the leading order (LO) all powers of $`\alpha _s\mathrm{ln}(Q^2/\mu ^2)`$, with $`\mu ^2`$ the factorization scale, are summed by the DGLAP evolution equations , which take into account only the strongly ordered parton transverse momenta ($`k_T`$) ladders. At the present, the next-to-leading order (NLO) contribution is also considered, including non-ordered $`k_T`$ contributions in a covariant way. In the current HERA kinematical regime the DGLAP approach is quite successful, although the theoretical expectation of deviations due to the parton saturation phenomena could be present . Following this point of view, recently the structure function measured at HERA has been studied considering an analysis based on the first non-linear corrections to the DGLAP evolution , showing for the gluon distribution the nonlinear effects play an increasingly role at $`x10^3`$ and $`Q^210`$ GeV2. The procedure gives stable results for that distribution at both low $`x`$ and $`Q^2`$, instead of negative or constant values as in the recent NLO DGLAP analysis (see  for details).

On the other hand, at very small $`x`$ the leading logarithms $`\alpha _s\mathrm{ln}(1/x)`$ are shown to be important. In the leading logarithmic approximation (LLA) the QCD Pomeron corresponds to the sum of ladder diagrams with reggeized gluons along the chain, which are strongly ordered in momentum fraction $`x`$. Such sum is described by the Balitzkij-Fadin-Kuraev-Lipatov (BFKL) equation . The corrections at next-to-leading level (NLLA) are now known, leading to strong modifications in the LLA spectrum . At first glance, the results turned out conflicting and not useful for suitable phenomenological studies. Further, these features have been cured by canceling not desired singularities via a suitable resummation at all orders of the perturbative expansion, satisfying the renormalization group requirements . The main shortcoming in doing phenomenology with a BFKL evolution beyond the leading order is the complicated mismatch between QCD perturbative and non-perturbative inputs, as the corresponding factorization properties are more involved than the DGLAP evolution.

A very promising approach encoding all order $`\alpha _s\mathrm{ln}(1/x)`$ resummation is provided by the QCD dipole picture. It was proven that such approach reproduces the BFKL evolution . The main process is the onium-onium scattering, that is the reaction between two heavy quark-antiquark states, in such way that the process is perturbative since the onium radius provides the essential scale at which the strong coupling constant is evaluated. In the large $`N_c`$ limit, the heavy pair and the soft gluons are represented as a collection of color dipoles. The cross section is written as a convolution between the number of the dipoles in each onium state and the basic cross section for dipole-dipole scattering due to two-gluon exchange. The QCD dipole model can be applied to deep inelastic lepton-nucleon scattering, assuming that the virtual photon at high $`Q^2`$ can be described by an onium. On the other hand, the proton is approximately described by a collection of onia with an average onium radius to be determined from phenomenology. This model has produced a successful description of the old structure function data .

The QCD dipole picture allows a systematic framework for testing the resummed next-to-leading order BFKL evolution kernels, producing predictions for the proton structure function. A method for doing this was recently proposed , where the resummation schemes are tested through the Mellin transformed $`j`$-moments of $`F_2`$. Moreover, it has been shown that a geometric scaling for the photon-proton cross section and the symmetry between low and high $`Q^2`$ regions are associated to the symmetry of the two-gluon dipole-dipole cross section . There, the proton is assumed to be a collection of independent dipoles at the time the interaction whose sizes are distributed around $`1/Q_s(x)`$, with the latter being also the mean distance between the centers of the neighboring dipole. The quantity $`Q_s^2(x)=\mathrm{\Lambda }^2e^{\lambda \mathrm{log}(1/x)}`$ is the saturation scale . Further, it was shown that the local geometric scale can be testable experimentally . That is, the cross sections for exclusive processes only depend on the ratio of scales $`Q/Q_s(x,b)`$, where now the saturation scale $`Q_s^2(x,b)`$ has a impact parameter ($`b`$) dependence.

In this work we revisit the phenomenology using this successful approach. The analysis is organized as follows. In the next section we shortly review the main formulae for the QCD dipole picture and give the phenomenological expressions for the structure functions and for the gluon distribution. In Sec. 3 we analyze the fitting procedure to the recent structure function data, determining the quality of the fit procedure and the range of applicability of the model. At the last section we drawn the main conclusions and summarize the results.

## 2 The QCD color dipole picture and structure functions

Let us start by describing the main features and basic formulas of the QCD dipole picture at high energies. The forward amplitude, $`𝒜(s,t)`$, for the onium-onium scaterring, integrated over impact parameter is written as,

$`𝒜(s,t=0)`$ $`=`$ $`i{\displaystyle d^2r_1d^2r_2𝑑z_1𝑑z_2\mathrm{\Phi }^{(0)}(𝒓_1,z_1)\mathrm{\Phi }^{(0)}(𝒓_2,z_2)}`$ (1)
$`\times {\displaystyle }{\displaystyle \frac{d^2s_1}{𝒔_1}}{\displaystyle \frac{d^2s_2}{𝒔_2}}n(Y/2,𝒓_1,𝒔_1)n(Y/2,𝒓_2,𝒔_2)\sigma _{\mathrm{dd}}(𝒔_1,𝒔_2),`$
where $`\mathrm{\Phi }^{(0)}(𝒓_i,z_i)`$ is the squared wave function of the quark-antiquark part of the onium wave function, $`𝒓_i`$ being the transverse size of the quark-antiquark pair and $`z_i`$ the longitudinal momentum fraction of the antiquark. The rapidity is given by $`Y/2\mathrm{ln}(x_0/x)`$, where $`x_0`$ is a phenomenological constant and $`x`$ is the Bjorken variable labelling the softer end of the produced dipole. The dipole density is labeled as $`n(Y/2,𝒓_i,𝒔_i)`$, which is discussed below. The quantity $`\sigma _{\mathrm{dd}}`$ is the elementary dipole-dipole cross-section and reads as,

$`\sigma _{\mathrm{dd}}(𝒔_1,𝒔_2)=2\pi \alpha _s^2[\mathrm{min}(𝒔_1,𝒔_2)]^2\left[1\mathrm{ln}\left({\displaystyle \frac{\mathrm{max}(𝒔_1,𝒔_2)}{\mathrm{min}(𝒔_1,𝒔_2)}}\right)\right].`$ (2)
In the large $`N_c`$ limit and in the leading-logarithmic-approximation the radiative corrections are generated by emission of gluons with strongly ordered longitudinal momenta fractions $`z_i>>z_{i+1}`$. The onium wave function with $`n`$ soft gluons can be calculated using perturbative QCD. In the Coulomb gauge the soft radiation can be viewed as a cascade of color dipoles emitted from the initial dipole, since each gluon acts like a quark-antiquark pair. The dipole density $`n(Y,𝒓,𝒔)`$ was defined  in such way,

$`N(Y,𝒔)={\displaystyle 𝑑zd^2r\mathrm{\Phi }^{(0)}(𝒓,z)n(Y,𝒓,𝒔)}`$ (3)
is the number of dipoles of transverse size $`𝒔`$ with the smallest light-cone momentum in the pair greater than or equal to $`e^Yp_+`$, where $`p_+`$ is the light-cone momentum of the onium. The whole dipole cascade can be constructed from a repeated action of a kernel $`𝒦`$ on the initial density $`n_0(r,𝒔)`$ through the dipole evolution equation,

$`n(Y,𝒓,𝒔_1)=n_0(𝒓,𝒔_1)+{\displaystyle _0^Y}𝑑y{\displaystyle _0^{\mathrm{}}}𝑑𝒔_2𝒦(𝒔_1,𝒔_2)n(y,𝒓,𝒔_2).`$ (4)
The evolution kernel $`𝒦`$ is calculated in perturbative QCD. For fixed $`\alpha _s`$ and in the limit of large $`N_c`$ the kernel has the same spectrum as the BFKL kernel . Therefore, the two approaches lead to the same phenomenological results for inclusive observables. The solution of Eq. (4) is given by ,

$`n(Y,𝒓,𝒔)={\displaystyle \frac{1}{2}}{\displaystyle \frac{𝒓}{𝒔}}{\displaystyle \frac{\mathrm{exp}\left(\omega _{IP}Y\right)}{\sqrt{7\pi \overline{\alpha }_s\zeta (3)Y}}}\mathrm{exp}\left({\displaystyle \frac{\mathrm{ln}^2(𝒓/𝒔)}{28\overline{\alpha }_s\zeta (3)Y}}\right),`$ (5)
where $`\omega _{IP}=4\overline{\alpha }_s\mathrm{ln}2`$, with $`\overline{\alpha }_s=\alpha _sN_c/\pi `$, and $`\alpha _{IP}=1+\omega _{IP}`$ is the Pomeron intercept. Therefore, the cross section grows rapidly with the energy since the number of dipoles in the light cone wave function grows rapidly on energy. The result in Eq. (5) has been obtained considering a process having only one scale, that is the onium radius. In the electron-proton deep inelastic scattering two scales are present and this result should be affected by non-perturbative contributions. The hard scale is given by the photon virtuality and the soft one is associated to the proton typical size.

Now, let us review the obtention of the proton structure functions in the QCD dipole picture. In order to do so, an assumption is made that the proton can be approximately described by onium configurations. The basic assumption is given by,

$`\sigma _{L,T}^{\gamma p}(x,Q^2)=\sigma _{L,T}^{\gamma onium}(x,Q^2,𝒌^2;\mu _F^2)\stackrel{~}{n}(𝒌^2;\mu _F^2),`$ (6)
where $`\stackrel{~}{n}(k^2;\mu _F^2)`$ is the probability of finding an onium in the proton, with $`\mu _F^2`$ being a unknown perturbative scale characterizing the average onium size, as will be seen later on. We follow Refs.  to write down the virtual photon-onium cross section, where the $`k_T`$-factorization approach  was considered within the framework of the QCD dipole model. It reads as,

$`\sigma ^{\gamma onium}(x,Q^2)={\displaystyle d^2r𝑑z\mathrm{\Phi }^{(0)}(𝒓,z)\sigma ^{\gamma dipole}(x,Q^2,𝒓)},`$ (7)
where $`\mathrm{\Phi }^{(0)}(𝒓,z)`$ was already defined before. The virtual photon-dipole cross section is written as,

$`Q^2\sigma ^{\gamma dipole}(x,Q^2,𝒓)={\displaystyle ^{Q^2}}d^2k{\displaystyle _0^1}{\displaystyle \frac{dz}{z}}\sigma _{\gamma g}({\displaystyle \frac{x}{z}},{\displaystyle \frac{𝒌^2}{Q^2}})_{onium}(z,𝒌,𝒓).`$ (8)
where $`\sigma _{\gamma g}/Q^2`$ is the Born cross section of an off-shell gluon of transverse momentum $`𝒌`$ in the the subprocess $`\gamma ^{}gq\overline{q}`$ and $`_{onium}`$ is the unintegrated gluon distribution function of an onium state of size $`r`$, containing the hard Pomeron physics. The relation between this quantity and the dipole density is written as,

$`𝒌^2_{onium}(z,𝒌,𝒓)={\displaystyle \frac{d𝒔^2}{𝒔^2}_0^1\frac{dz^{}}{z^{}}n(z^{},𝒓,𝒔)\sigma ^{\gamma dipole}(\frac{z}{z^{}},(𝒔𝒌)^2)},`$ (9)
where $`n(z^{},𝒓,𝒔)`$ is the density of dipoles of transverse size $`𝒔`$ with the smallest light-cone momentum in the pair equal to $`z^{}p_+`$ in a dipole of transverse size $`𝒓`$, of total momentum $`p_+`$. This is given by the solution in Eq. (5).

The next step is to apply a Mellin transform in the $`x`$-space and further a second Mellin transform in the $`k^2`$-space, where the conjugated variables are $`\omega `$ and $`\gamma `$, respectively. After that, one obtains ,

$`Q^2\sigma ^{\gamma dipole}(x,Q^2,𝒓)={\displaystyle \frac{d\omega }{2i\pi }\frac{d\gamma }{2i\pi }\sigma _{\gamma g}(\omega ,1\gamma )_{onium}(\omega ,\gamma )(Q^2𝒓^2)^\gamma x^\omega }.`$ (10)
The unintegrated gluon distribution for an onium of size $`r`$ can be obtained from perturbative QCD and it is written as ,

$`_{onium}(\omega ,\gamma )=\overline{\alpha }_s{\displaystyle \frac{1}{\gamma }}{\displaystyle \frac{\nu (\gamma )}{\omega \overline{\alpha }_s\chi (\gamma )}},`$ (11)
where in the leading-order BFKL approach, $`\chi _{\mathrm{L0}}(\gamma )=2\psi (1)\psi (\gamma )\psi (1\gamma )`$ and $`\overline{\alpha }_s\chi (\gamma )`$ is the eigenvalue of the LO BFKL kernel. The Mellin transform of the gluon-dipole coupling, $`\nu (\gamma )`$, can be obtained from $`k_T`$-factorization approach as,

$`\nu (\gamma )={\displaystyle \frac{2^{\mathrm{\hspace{0.17em}2}\gamma 1}}{\gamma }}{\displaystyle \frac{\mathrm{\Gamma }(1\gamma )}{\mathrm{\Gamma }(1\gamma )}}.`$ (12)
The unintegrated gluon distribution, Eq. (11), has a pole at $`\omega _{IP}=\overline{\alpha }_s\chi (\gamma )`$ and then the integration over $`\omega `$ can be done analytically. Furthermore, the cross section $`\sigma _{\gamma g}`$ can be written in function of the LO photon impact factors $`h_{L,T}(\gamma )`$, corresponding to the perturbative coupling to the photon. Namely, $`\sigma _{\gamma g}4\pi ^2\alpha _{em}e_f^2h_{L,T}(\gamma )`$, with $`e_f^2`$ being the total charge of the quarks flavor contributing to the reaction and where,

$`\left(\begin{array}{c}h_T\\ h_L\end{array}\right)={\displaystyle \frac{\overline{\alpha }_s}{9\gamma }}{\displaystyle \frac{\left[\mathrm{\Gamma }(1\gamma )\mathrm{\Gamma }(1+\gamma )\right]^3}{\mathrm{\Gamma }(22\gamma )\mathrm{\Gamma }(2+2\gamma )}}{\displaystyle \frac{1}{1\frac{2}{3}\gamma }}\left(\begin{array}{c}(1+\gamma )(1\frac{\gamma }{2})\\ \gamma (1\gamma )\end{array}\right),`$ (13)
Putting Eqs. (13) and (11) in Eq. (10) and solving analytically on $`\omega `$ for the pole on $`\omega =\omega _{IP}`$, one obtains,

$`Q^2\sigma ^{\gamma dipole}(x,Q^2,𝒓)=8\pi ^2\alpha _{em}e_f^2\overline{\alpha }_s{\displaystyle \frac{d\gamma }{2\pi i}h_{L,T}(\gamma )\frac{\nu (\gamma )}{\gamma }\left(𝒓^2Q^2\right)^\gamma \mathrm{exp}\left[\overline{\alpha }_s\chi (\gamma )\mathrm{ln}\frac{x_0}{x}\right]}.`$ (14)
The virtual photon-onium cross section is obtained replacing Eq. (14) in Eq. (7). The result depends on the squared wave function $`\mathrm{\Phi }^{(0)}`$ of the onium state, which cannot be computed perturbatively. This is solved by eliminating that dependence through an averaging over the wave function of transverse size ,

$`(\mu _F^2)^\gamma =𝒓^2{\displaystyle 𝑑zd^2r(𝒓^2)^\gamma \mathrm{\Phi }^{(0)}(𝒓,z)},`$ (15)
where $`\mu _F^2`$ is a scale which is assumed to be perturbative. With this result, the photon-onium cross section, Eq. (7), takes the following form,

$`\sigma ^{\gamma onium}(x,Q^2)={\displaystyle \frac{4\pi ^2\alpha _{em}}{Q^2}}\mathrm{\hspace{0.17em}2}\overline{\alpha }_s{\displaystyle \frac{d\gamma }{2\pi i}h_{L,T}(\gamma )\frac{\nu (\gamma )}{\gamma }\left(\frac{Q^2}{\mu _F^2}\right)^\gamma \mathrm{exp}\left[\overline{\alpha }_s\chi (\gamma )\mathrm{ln}\frac{x_0}{x}\right]}.`$ (16)
In order to compute the photon-proton cross section one needs to know the probability of finding an onium in it, $`\stackrel{~}{n}(\gamma ;\mu _F^2)`$. Relying on renormalization group properties, namely the overall result should be independent of the scale perturbative $`\mu _F^2`$, an suitable ansatz is given by ,

$`\stackrel{~}{n}(\gamma ;\mu _F^2)=n_{\mathrm{eff}}(\gamma )\left({\displaystyle \frac{\mu _F^2}{Q_0^2}}\right)^\gamma ,`$ (17)
where $`Q_0^2`$ is interpreted as a characteristic non-perturbative scale of the hadron-onium scattering. The quantity $`n_{\mathrm{eff}}`$ means the average number of primary dipoles in the proton and $`r_02/Q_0`$ their average transverse diameter (see second reference on ). Placing the ansatz Eq. (17) in Eq. (6) and using the result on Eq. (16), the proton structure functions can be calculated. The integration over $`\gamma `$ can be performed using the steepest descent method, with the saddle point given by $`\chi ^{}(\gamma _s)=\mathrm{ln}(Q^2/Q_0^2)/\overline{\alpha }_s\mathrm{ln}(x_0/x)`$. Using the expansion of the BFKL kernel near $`\gamma =\frac{1}{2}`$, one obtains,

$`\gamma _s={\displaystyle \frac{1}{2}}\left[1\kappa (x)\mathrm{ln}\left({\displaystyle \frac{Q}{Q_0}}\right)\right],`$ (18)
where $`\kappa (x)=[\overline{\alpha }_s\mathrm{\hspace{0.17em}7}\zeta (3)\mathrm{ln}\frac{x_0}{x}]^1`$ is the diffusion coefficient at rapidity $`Y=\mathrm{ln}(x_0/x)`$. Putting all together, the structure functions are written in the simple form ,

$`F_{T,L}(x,Q^2)=H_{T,L}{\displaystyle \frac{\overline{\alpha }_s\pi ^3e_f^2n_{\mathrm{eff}}}{96}}\left({\displaystyle \frac{x_0}{x}}\right)^{\omega _{IP}}{\displaystyle \frac{Q}{Q_0}}\sqrt{2\kappa (x)/\pi }\mathrm{exp}\left[{\displaystyle \frac{\kappa (x)}{2}}\mathrm{ln}^2{\displaystyle \frac{Q}{Q_0}}\right],`$ (19)
where $`H_T=9/2`$ and $`H_L=1`$, defining $`n_{\mathrm{eff}}n_{\mathrm{eff}}(\gamma _s)`$ and $`r_0r_0(\gamma _s)`$. In the nest section we consider the equation above to determine the parameters of the model by fitting the recent data on the proton structure function $`F_2=F_T+F_L`$. The conditions to obtain the structure functions from the saddle-point method fix the conditions of its safe applicability. Hence, the following bound should be covered,

$`\kappa (x)\mathrm{ln}\left({\displaystyle \frac{Q}{Q_0}}\right){\displaystyle \frac{\mathrm{ln}\left(\frac{Q}{Q_0}\right)}{\mathrm{ln}\left(\frac{x_0}{x}\right)}}1,`$ (20)
which is realized for the region of moderate $`Q/Q_0`$ when compared to the range on $`x_0/x`$. It is worth to mention that at higher $`Q^2`$, the DLLA becomes valid, with the pole now at $`1/\gamma `$ in the kernel $`\chi (\gamma )`$.

## 3 Fitting results and discussions

In this section we will perform a fitting procedure using the recent HERA data on the proton structure function  using the QCD dipole picture given by Eq. (19). Having determined the parameters of the model from data, the longitudinal structure function and the proton gluon distribution are calculated without further adjustments. We made an simple attempt to describe phenomenologically the charm structure function, although the full expressions for it to be known in the literature (last reference in ).

For sake of simplicity, we have defined the overall normalization for the proton structure function, $`𝒩_p=(H_T+H_L)\overline{\alpha }_s\pi ^3e_f^2n_{\mathrm{eff}}/96`$. Also, for some fitting procedures we will include a non-singlet contribution and threshold (large $`x`$) effects. Namely, the hard Pomeron expression, Eq. (19), is multiplied by the factor $`(1x)^{n_s}`$, while the non-Pomeron contribution is multiplied by $`(1x)^{n_{ns}}`$. The expressions for the non-singlet piece and $`n_{s,ns}`$ read as,

$`F_2^{\mathrm{ns}}(x,Q^2)`$ $`=`$ $`A_Rx^{1\alpha _R}\left({\displaystyle \frac{Q^2}{Q^2+a_R}}\right)^{\alpha _R}(1x)^{n_{ns}(Q^2)}`$ (21)
$`n_{ns}(Q^2)`$ $`=`$ $`{\displaystyle \frac{3}{2}}\left(1+{\displaystyle \frac{Q^2}{Q^2+d}}\right),n_s(Q^2)={\displaystyle \frac{7}{2}}\left(1+{\displaystyle \frac{Q^2}{Q^2+b}}\right),`$ (22)
where $`\alpha _R=0.4`$ is the Reggeon intercept, considered fixed in this analysis. Other values can be considered, however the $`\chi ^2`$ quality of the fit is worst.

Here, we have considered the following fitting procedures:

* (I) Adjusting data in the whole range $`1.5Q^2150`$ GeV2 and all $`x`$, considering the QCD dipole, Eq. (19), plus the non-singlet contribution, Eq. (21). The motivation is to include the large $`x`$ data at high $`Q^2`$, where the fixed target data (E665 and NMC points ) were also added. The results for H1 and ZEUS data sets are presented in Tables (1) and (2). The number of points is equal to 230 (H1) and 259 (ZEUS), respectively. The total number of parameters is 7;
* (II) Adjusting data in the range $`1.5Q^260`$ GeV2 and all $`x`$, considering the QCD dipole plus the non-singlet contribution. We have fixed the parameter $`x_0=1`$ in this case, reducing the final number of parameters. The results for H1 and ZEUS data sets are presented in Tables (1) and (2);
* (III) Adjusting data in the whole range $`1.5Q^2150`$ GeV2 and restricting to small region $`x10^2`$, considering only the QCD dipole contribution. We have considered the fixed values $`x_0=1`$ and $`n_s=7`$.

Lets start discussing the procedure (I). The main difference between the two data sets are the values for the parameter $`x_0`$, being larger for ZEUS than for H1. On the other hand, the effective intercept for the hard QCD contributions is quite the same, that is $`1+\omega _{IP}=1.25(8)`$, with the quality of fit better for ZEUS data set. It is worth to mention that the low values obtained for the parameters $`b`$ and $`d`$ suggests that they can be fixed in the usual way, $`n_s=7`$ and $`n_{ns}=3`$, in the range considered here. This procedure would reduce strongly the number of free parameters. From the fit, it is also clear the non-singlet piece is $`Q^2`$ independent, since $`a_R`$ is quite small.

For the procedure (II), the picture is similar to the previous one, with the effective intercept $`1+\omega _{IP}=1.26(9)`$ coming out different for the two data set. The quality of fit is also slightly equal in the two cases (see Tables). The conclusion about $`n_{s,ns}`$ is the same as the procedure (I).

Finally, in the procedure (III) we have considered only the small $`x10^2`$ data (102 data points for H1 and 107 for ZEUS). In this case, $`x_0=1`$ and $`n_s=7`$ are taken as fixed, reducing the number of parameters equal to 3. Indeed, their values are in large extent similar with a better quality of fit for the H1 data set. This procedure is similar to the previous analysis on Refs. , with an even lower effective power ($`1+\omega _{IP}1.282`$ in ). The low value for the fixed coupling constant $`\alpha _s0.1`$ reveals the well known necessity of sizeable higher order corrections to the approach. Accurate analysis in this lines, considering resummed NLO BFKL kernels, has been proposed recently  producing a reasonable $`\alpha _s0.2`$ for the typical $`Q^2`$ range considered in phenomenology for structure functions.

In order to illustrate the analysis, in Figs. 1 and 2 we present the fit result (procedure II) for the inclusive structure function for H1  and ZEUS  data set, including the fixed target points  and extrapolated up to $`Q^2=150`$ GeV2. It is worth to mention that the new data are more precise than the previous analysis and the same quality of fit is obtained in a small range on $`Q^2`$, in contrast with the studies in . However, a rough data description can be obtained up to higher values. Here, our aim is discuss the model as the quality of fit. We have checked using the new H1 data that a $`\chi ^2/\mathrm{dof}=1.4`$ is obtained in the range $`1.5Q^2150`$ for all $`x`$. This in some extent shows the range on virtuality where a LO BFKL approach is suitable. That is, low $`x`$ and not so high $`Q^2`$. In virtualities of hundreds of GeV, the physics should be described by the DLA limit of the BFKL approach and which is a common limit also in the DGLAP evolution . An interesting further study is to realize a best fit test to the DLA limit, determining its range of applicability. After this discussion, having determined the parameters of the model from data on $`F_2`$, it is possible determine without further adjustments the structure function $`F_L`$ and the gluon distribution $`xG(x,Q^2)`$. From Eq. (19) and the definition of the overall normalization, we can see that $`F_L=(2/11)F_2`$, since $`H_T+H_L=11/2`$ and $`H_L=1`$. In Fig. (3) this result is shown using parameters from procedure II from the H1 data set against the recent $`F_L`$ data . The results are in good agreement with the recent data, obtained without further adjustments.

We made an attempt to describe phenomenologically the charm structure function, despite explicit expressions for it to be available. One suppose the following relation between the $`F_2`$ structure function and $`F_2^{c\overline{c}}`$,

$`F_2^{c\overline{c}}(x,Q^2)={\displaystyle \frac{2}{5}}\left(1+{\displaystyle \frac{4m_c^2}{Q^2}}\right)^1F_2(x,Q^2),`$ (23)
where we have used $`m_c=1.5`$ GeV and the factor $`2/5`$ corresponds to $`e_c^2/(e_f^2+e_c^2)`$ . The result is shown in Fig. (4), where the parameters from procedure II (ZEUS data set) were taken into account. The simple ansatz is reasonable only at low $`Q^2`$, overestimating the data on higher virtualities. It should be stressed that this function can be calculated using the complete expressions for the photon impact factors, $`h_{T,L}(\gamma ;m_f)`$, which consider also the dependence on the quark masses, $`m_f`$ (see last reference in ).

As a final analysis, in the QCD dipole approach the gluon distribution function can be calculated in a straightforward way from $`F_2`$. The result is also independent of the overall normalization, which contains part of the non-perturbative inputs of the model. The gluon function is given by,

$`xG(x,Q^2)`$ $`=`$ $`\left[h_T(\gamma =\gamma _s)+h_L(\gamma =\gamma _s)\right]^1F_2(x,Q^2),`$ (24)
$`=`$ $`{\displaystyle \frac{9\gamma _s}{\overline{\alpha }_s}}{\displaystyle \frac{1\frac{2}{3}\gamma _s}{1+\frac{3}{2}\gamma _s\frac{3}{2}\gamma _s^2}}{\displaystyle \frac{\mathrm{\Gamma }(22\gamma _s)\mathrm{\Gamma }(2+2\gamma _s)}{\left[\mathrm{\Gamma }(1\gamma _s)\mathrm{\Gamma }(1+\gamma _s)\right]^3}}F_2(x,Q^2).`$
In Fig. (5) we shown the gluon distribution function obtained using Eq. (24) and the parameters of the procedure (III) for the H1 and ZEUS data sets. The curves represent $`xG(x,Q^2)`$ as a function of $`x`$ for distinct virtualities, 1.5, 15 and 150 GeV2. It would be timely to compare those curves with the ones from using a resummed kernel. In order to compare these predictions with NLO QCD DGLAP fit results, in Fig. (6) one plots the QCD dipole predictions with the recent QCD analysis. In the upper plot, we show the results for the H1 QCD fit (central value, without bands) in Ref.  for $`Q^2=4`$ GeV2 and the earlier analysis in Ref. . The deviation between the recent H1 fit and the QCD dipole curves (parameters from procedure III) is sizeable, whereas it is consistent with the previous analysis. Deviations are found in the comparison with the ZEUS NLO QCD fit  at $`Q^2=7`$ GeV2, in the lower plot. However, it should be stressed that the gluon distribution is an indirect observable, and its determination is model-dependent. For sake of comparison with other determinations of the gluon distribution, we have plotted the result using the Regge-like hard Pomeron contribution , which it is shown in the lower plot for $`Q^2=7`$ GeV2. Such approach shows that the ratio of the gluon distribution to the hard-Pomeron part of the singlet quark distribution is about $`8\pm 1`$ in a large range on virtualities. Furthermore, the dependence on $`Q^2`$ of the hard-Pomeron contribution presents very good agreement with DGLAP evolution for $`Q^25`$ GeV2.

## 4 Conclusions

We have revisited the phenomenology to the structure functions in the QCD dipole picture, performing a fit to the recent data on the proton structure function. The model for the hard Pomeron provides a reduced number of free-parameters and the results produces a good quality fit. The adjusted constants are the overall normalization to $`F_2`$, the effective Pomeron intercept $`\alpha _{IP}=1+\omega _{IP}`$, the average transverse diameter of the primary dipoles in the proton $`r_0=2/Q_0`$ and $`x_0`$ scaling the $`\mathrm{ln}(1/x)`$ behavior. Moreover, the longitudinal structure function and the gluon distribution function can be determined without further adjustments. The results for $`F_L`$ are in good agreement with the recent measurements. The main shortcoming from the approach is a low value for the strong coupling constant, $`\alpha _s011`$, showing that a next-to-leading level analysis must be performed. The number of initial dipoles is either small, $`n_{eff}=3`$, and they have a large size $`r_01`$ fm. Resummed perturbative kernels should provide a consistent value for the strong coupling constant in the range of virtualities considered, as reported in Ref. . In this sense, the model presented here, although providing a remarkable data description, can be considered as effective and further phenomenology should consider higher order effects. Regarding the gluon distribution, the result is distinct from the recent NLO QCD fits, producing higher values at small $`x`$. We have also made a simple ansatz for describe the charm structure function, obtaining a very rough result. This is due we are disregarding the dependence of the photon impact factors on the quark mass, which modifies the final result in contrast with the assumption of avoiding the quark mass.

## Acknowledgments

M.V.T.M. thanks the support of the High Energy Physics Phenomenology Group (GFPAE, IF-UFRGS) at Institute of Physics, Porto Alegre.


format: text

text:
# A far infrared view of the Lockman Hole from ISO 95 𝜇m observations - I. A new data reduction method Based on observations obtained with the Infrared Space Observatory, an ESA science missions with instruments and contributions funded by ESA Member States and the USA (NASA).

## 1 Introduction

The star formation in local galaxies is univocally found to be associated with dense dust-obscured clouds which are optically thick in the UV and become transparent or even emissive at long infrared wavelengths. For this reason, the far infrared domain is expected to be quite instrumental for studying not only the physics of star formation in our and closeby galaxies, but also the early phases of galaxy evolution, when stellar formation was far enhanced compared to what happens in the local universe.

The IRAS satellite mission proved indeed the potential of IR observations to detect galaxies optically obscured by dust. Only a fraction of the 25.000 sources detected in the All Sky Survey were found to have bright optical counterparts (Soifer et al., 1987), and of these most are local late-type spirals.

The IRAS survey was devoted to investigate the properties of the IR emission by local galaxies, at redshift $`<`$ 0.2 (Ashby et al., 1996). Only few sources were detected by IRAS at higher redshifts, typically ULIRGs magnified by gravitational lenses, like F10214+4724 (z=2.28, Rowan-Robinson et al., 1991). IR galaxy counts based on the IRAS data (Rowan-Robinson et al., 1994; Soifer et al., 1984) showed some marginally significant excess of faint sources with respect to no evolution models (Hacking et al., 1987; Franceschini et al. 1988; Lonsdale et al., 1995; Gregorich et al., 1995; Bertin et al, 1997), but did not provided enough statistics and dynamic range in flux to discriminate between evolutionary scenarios.

The cosmological significance of far-IR studies was first emphasized by the COBE detection of an isotropic far-IR/submillimeter background, of extragalactic origin (CIRB), interpreted as the integrated emission by dust present in distant and primeval galaxies (Puget et al., 1996, Hauser et al., 1998), including an energy density larger than that in the UV/optical background (Lagache et al., 1999).

With the advent of the Infrared Space Observatory (ISO, Kessler et al., 1996) the improved resolution and sensitivities of its cameras made possible deeper IR surveys, allowing us for the first time to detect faint IR galaxies at cosmological distances, both in the mid and in the far infrared.

The deep 15 $`\mu `$m counts determined with ISOCAM (Cesarsky et al., 1996) have revealed a very significant departure from the euclidean slope (Elbaz et al., 1999; Gruppioni et al., 2002), which has been interpreted as evidence for a strongly evolving population of starburst galaxies (Franceschini et al., 2001; Chary & Elbaz, 2001; Xu et al., 2001; Oliver et al., 2001). ISO surveys at longer (far-IR) wavelengths with the photo-polarimeter ISOPHOT (Lemke at el., 1996) found some evidence of evolution in the 175 $`\mu `$m channel (the FIRBACK survey: Puget et al. 1999, Dole et al. 2001; ELAIS survey, Efstathiou et al., 2000; the Lockman Hole survey: Kawara et al. ,1998, Matsuhara et al., 2000; Juvela et al., 2000). Unfortunately, at such long wavelengths the ISO observatory was quite limited by source confusion to moderately faint fluxes ($`S_{175\mu }135`$ mJy).

In principle, the shorter wavelength 95 $`\mu `$m C100 channel of ISOPHOT could allow a substantial improvement, by a factor 2, in spatial resolution and a significantly lower confusion noise. Furthermore, the filter samples in an optimal way the dust emission peak in the Spectral Energy Distribution (SED) of star-forming galaxies around 60 to 100 $`\mu `$m. For luminous infrared galaxies, emitting more than 80% of the flux in the far-IR, this far-IR peak is the best measure of the bolometric luminosity of such galaxies, and the best estimator of their star formation rate.

The well-known problem with ISOPHOT C100 observations was the difficulty in the data reduction to account for all the instrumental effects, due to cosmic-ray impacts and transient effects in the detectors producing spurious detections which may contaminate the final source lists.

We present in this paper a new method for the reduction of ISOPHOT C100 data, which we developed along similar lines as the code designed for ISOCAM data reduction by Lari et al. (2001, hereafter L01). We illustrate the value of this method with application to a deep ISOPHOT C100 survey in the Lockman Hole, a region particulary suited for the detection of faint infrared sources due to its low cirrus emission. The good quality of these data and a careful reduction allow us to reach faint detection limits ($`20`$ mJy). We focus in particular on the implications for the evolutionary models as derived from galaxy counts. Our results favour a scenario dominated by a strongly evolving population, quite in agreement with the model discussed by Franceschini et al. (2001). In a forthcoming paper we will discuss the optical identifications of these Lockman 95 $`\mu `$m sources with radio and ISOCAM mid-IR counterparts (Rodighiero et al., in preparation) and we will explore the nature of our far-IR sources.

The present paper is organized as follows. In Section 2 we discuss our reduction technique. In Section 3 we comment on the simulations we used to compute the completeness of our sample, and the photometric corrections. Section 4 is devoted to the flux calibration. We then report in Section 5 on our application to the Lockman Hole 95 $`\mu `$m data and our results on source counts. Our conclusions are reported in Section 6. In the three Appendices we detail some technical aspects of our data reduction method.

We assume throughout this paper $`\mathrm{\Omega }_M`$=0.3, $`\mathrm{\Omega }_\mathrm{\Lambda }`$=0.7 and $`H_0`$=65 $`Kms^1Mpc^1`$.

## 2 A NEW TOOL FOR THE REDUCTION OF ISOPHOT-C DATA

The reduction of ISO data requires a careful treatment of various external and instrumental effects affecting the detectors. The results recently obtained by L01 in developing a data reduction technique for ISOCAM data prompted us to attempt a similar approach for the analysis of ISOPHOT data.

PHOT C100 is a $`3\times 3`$ array of Ge:Ga with 0.7$`\times `$0.7$`\times `$1mm elements. The effective size of the pixels on the sky is 43.5$`\times `$43.5 arcsec, the distance between the pixel centers is 46.0 arcsec. There are 6 filters available for C100, covering the wavelength range from $``$60 to $``$100 $`\mu `$m.

As for the case of ISOCAM (Long Wavelength) Si:Ga detectors, two main effects must be considered when dealing with ISOPHOT-C data, produced by cosmic ray impacts (glitches) and detector hysteresis (i.e. the slow response of the detector to flux variations). The method discussed by L01 was based on the assumption that the incoming flux of charged particles generates transient behaviours with two different time scales: a fast (breve) and a slow (lunga) one. The method basically consists in looking at the time history of each detector pixel and identifying the stabilization background level. Then it models the glitches, the background and the sources with all the transients over the whole pixel time history.

If the approach is similar to that for ISOCAM, some peculiarities of the far-infrared detectors need to be treated with specific care. We found, in any case, that the description of transients with the equations used in the case of ISOCAM pixels provides good fits also to the PHOT-C data (after adapting the temporal and charge parameters).

Let us mention that a model like that of Fouks & Schubert (1995) can fit, with suitable parameters, the brighest sources, both in the case of CAM (Coulais & Abergel, 2000) and of PHOT (Coulais et al., 2001), indicating very similar shapes of the short term transient. In the case of CAM the method is an alternative to this model. As far as PHOT is concrned, the glitches look very similar to those found in the time histories of CAM’s pixels. Moreover, long term transients have been observed also in the PHOT detectors. These considerations suggest that the Lari model is applicable also to PHOT data, using the appropriate parameters. We will show that, in spite of the numerous glitches recorded, it has been possible to define a reduction strategy with reliable results, including the stability of the temporal parameters.

As already mentioned, the method by Lari et al. (2001) describes the sequence of readouts, or time history, of each pixel of CAM/PHOT detectors in terms of a mathematical model for the charge release towards the contacts. Such a model is based on the assumption of the existence, in each pixel, of two charge reservoirs, a short-lived one $`Q_b`$ (breve) and a long-lived one $`Q_l`$ (lunga), evolving independently with a different time constant and fed by both the photon flux and the cosmic rays. Such a model is fully conservative, and thus the observed signal $`S`$ is related to the incident photon flux $`I`$ and to the accumulated charges $`Q_b`$ and $`Q_l`$ by the (see also Lari et al., 2002):

$$S=I\frac{\mathrm{d}Q_{tot}}{\mathrm{d}t}=I\frac{\mathrm{d}Q_b}{\mathrm{d}t}\frac{\mathrm{d}Q_l}{\mathrm{d}t},$$
(1)
where the evolution of these two quantities is governed by the same differential equation, albeit with a different efficiency $`e_i`$ and time constant $`a_i`$

$$\frac{\mathrm{d}Q_i}{\mathrm{d}t}=e_iIa_iQ_i^2\mathrm{where}i=b,l,$$
(2)
so that

$$S=(1e_be_l)I+a_bQ_b^2+a_lQ_l^2.$$
(3)
The values of the parameters $`e_i`$ and $`a_i`$ are estimated from the data and are constant for a given detector, apart from the scaling of $`a_i`$ for the exposure time and the average signal level along the pixel time history, which is governed by the

$$a_i=\frac{t}{t_0}\sqrt{\frac{S}{S_0}}a_{i,0},$$
(4)
where $`a_{i,0}`$ is the value of $`a_i`$ relative to a reference exposure time $`t_0`$ and average signal level $`S_0`$. The model for the charge release, however, is exactly the same for CAM and PHOT detectors.

The values for the model parameters of CAM and PHOT respectively are reported in Table 1, together with the times characteristic of the long and the short transients (dt/a\_i).

### 2.1 DATA REDUCTION

The Interactive Data Language (IDL) has been used to develop all the procedures needed for the reduction of PHOT-C data.

Before discussing in detail the reduction procedures, let us note the main difference between ISOCAM and PHOT. The ISOPHOT pixels can be considered as independent detectors, with entirely uncorrelated behaviours and responses. This is not the case for ISOCAM, where all pixels have a common electronics. An example is presented in Figure 1, where we report a statistically rich data-set from the ISOPHOT ELAIS surveys (Oliver et al., 2000) in the southern S1 field. For an ELAIS raster observation with ISOPHOT, we have taken all the ramps along the time history and computed the median over each ramp position (one ramp is composed of 64 readouts), excluding those readouts with a signal exceeding $`3\times \sigma `$ (where $`\sigma `$ is the standard deviation of the signal over the whole pixel history). In the figure we report for each PHOT C100 detector pixel the RMS of the median signal over the whole pixel history in every ramp position. This gives a quantitative idea of the intrinsic noise of each detector pixel. The values are similar, except for pixel (0,1), which clearly shows a noise $``$ 3 times higher than the others. A similar analysis is performed over all ELAIS rasters (both in northern and southern fields) and found very similar levels for each pixel. This consideration allowed us to consider the pixels’ responsivities stable as a function of time and of the orbital position. The different peculiarities of the nine ISOPHOT C100 pixels imply that any kind of correction must be computed as a function of the pixel.

### 2.2 From Raw Data through the fitting algorithm

The raw data (ERD level) are converted into a raster structure containing instrumental informations on the observation, and astrometric informations on every pointing. The ramps (in Volts) are corrected for the non-linear response of the detector using a new technique (Appendix A) and converted in ADU/gain/s. The standard PHOT Interactive Analysis (PIA, Gabriel & Acosta-Pulido, 1999) package process the data fitting the ramps (in units of Volts).

In our procedure, the data are corrected for short-time cosmic rays. Readouts affected by such events are masked and their positions stored before copying the “deglitched” data into a new structure (called ’liscio’, as for ISOCAM).

We evaluate the general background as the stabilization level along the whole time history of each pixel (for clarity, in the following we will refer to this quantity as the “stabilization background”).

We then apply a constant positive offset signal to the data in order to take into account the contribution of thermal dark current (which is not otherwise accounted for in the preliminary pipeline) when the latter is estimated to be important, i.e. when the deepest dipper’s depth exceeds 10% of the stabilization background.

The task computing the stabilization background also performs an initial guess of the fitting parameters, storing them in the ’liscio’ structure.

The signal as a function of time is finally processed, independently for every pixel. The fitting procedure models the transients along the time history, and the features on both short and long timescales produced by cosmic ray impacts. At this level the code estimates several quantities needed to build the final maps on which source extraction will be performed:

* the charges stored into the $`breve`$ and $`lunga`$ reservoirs at each readout.
* The local background, i.e. the signal to be expected on the basis of the previously accumulated charges if only the stabilization background were hitting the detector.
* The model signal produced by the incident flux coming from both the stabilization background and detected sources.
* The “reconstructed” signal, i.e. the model signal recovered not only from glitches, but also from the transients due to the pointing-to-pointing incident flux change.

Furthermore the code recognizes sources (above a given threshold level) and recovers all the time histories “reconstructing” the local background as it would appear in absence of glitches.

With the previously defined quantities we can define two different kinds of fluxes:

1) ”Unreconstructed” fluxes and maps are derived from the excess of the measured signal with respect to the ‘local background’, and represent the flux excess not recovered from transients but only from glitches.

2) ”Reconstructed” fluxes and maps are computed from the “reconstruted” signal, and take into account not only the glitches, but also the transients “on” the sources.

In other words, the main difference between these two flux estimates appears evident when a pixel ’sees’ a source. Reconstructed fluxes describe sources taking into account the real excess of signal with respect to the local background plus the transient modelling. The latter ’recovers’ for the flux loss due to the slow response of the pixel when it detects an intense prolonged signal, like a source, during a pointing exposure.

Unreconstructed fluxes do not take into account the effects of this source modelling, and thus represent the effective flux collected by the detector during the raster exposure. We will see with simulations that our code is not able to properly recover faint sources, and for this reason we will use unreconstructed fluxes in order to generate our final maps.

As an example, Figure 2 shows how the code fits and describes the background and the transients. In the top panel we plot an example of real and model data through the pixel history. The solid “noisy” line represents the observed data, with superimposed the best-fit model (continous line), which nicely follows the transients induced by cosmic rays. The dotted line is the data corrected for transients and deglitched (reconstructed signal). The dot-dashed horizontal line is the assumed stabilization background, while the local background corresponds to the three-dots dashed line. The bottom panel shows how the code works when it sees a bright source. In this case the real data have been smoothed. The ‘unreconstructed’ signal is computed as the difference of the observed data and the local background.

The fitting algorithm starts with the brightest glitches identified in the pixel time history, assumes discontinuities at these positions, and tries to find a fit to the time history that satisfies the mathematical model assumed to describe the solid-state physics of the detector. In the fit we use the same default parameter values for all the pixels (the physical parameters scaled only for the stabilization background), leaving as free parameters only the charge values at the beginning of the observations and at the ’peaks’ of glitches.

By successive iterations, the parameters and the background for each pixel are adjusted to better fit the data, until the rms of the difference between model and real data is smaller than a given amount (e.g. 15 ADU/gain/s ).

All the features described in Figure 2 are present in the ISOCAM LW observations as well, but sometimes (and not so rarely) we find some peculiar behaviour in our data that does not correspond to any usual transient. These drop-outs can be explained in terms of saturation of the detector: when a very strong glitch impacts on a pixel, this can reach the saturation regime (the top level of the instrument dynamical range) and then be reset. After the impact, the detector needs several readouts before loosing the memory of such a shocking event. Usually these drop-outs appear as isolated features. However we found that an intense repeated series of impacts can cause much more serious problems to the detector electronics, and many drop-outs can appear consequently in the data for a significant fraction of the pixel time history, as shown in Figure 3. It is useful to note here that such long-duration drop-outs are not resolved in ISOCAM, where they affect only one or two readouts, because of the different temporal resolution (1/32 seconds for PHOT against 2 seconds for ISOCAM, for each readout). All the readouts affected by this problem must be masked and excluded in the fitting procedure.

### 2.3 The Interactive Analysis

After the first run of the automatic fitting procedure, the next step is the interactive mending of fitting failures. This massive work of interactive analysis is carried out with an easy-to-use widget interface, which allows any kind of repairing which may be necessary. This stage strongly depends on the assumed level of the reduction, related to the “goodness” of the observational data set. For this reason observations characterized by different observing parameters (exposure time, raster step …) need a specific treatment during the interactive reduction procedure. In the case of PHOT we found more efficient to use smoothed data as input for the fitting procedure. For the Lockman Hole we chose a smoothing factor of 32, implying that the code reads and uses only 1 readout every 32 (1 second) to which it associates the median value of the previous 32 readouts. Thus short-term features and noise are reduced, allowing a better (and faster) modelling of the general trend of the signal.

The automatic detection of sources is checked through “eye-balling” on the time history of each pixel. If the code fails (because it finds a source where a source is not present, or when it fits improperly a real source), a local interactive fit is carried out. Furthermore, the limited number of pixels of the PHOT-C100 detector (a 3x3 array) allows to carefully check the presence of sources scanning the time history of each pixel where the median level of the signal in a pointing position peaks with respect to nearby values. In order to minimize the loss of sources, we also check all the raster positions with a signal greater than a given threshold level (15 ADU/gain/s for the Lockman Hole).

### 2.4 Map generation and source detection

Once a satisfying fit is obtained for all the pixels over the whole pixel history, our pipeline proceeds with the generation of sky maps. An image for each raster position is created, by averaging the signals of all the readouts relative to that pointing for each pixel. The signal is then converted to flux units (mJy pixel<sup>-1</sup>, see Appendix B and Section 4.1), glitches and bad data are masked and the images are then combined to create the final raster maps (one for each raster position). These images are projected onto a sky map (raster image) using the projection algorithm available for ISOCAM data in the CIA package (Cam Interactive Analysis, Ott et al., 2001). When projecting the signal on the sky, we make use of the nominal raster astrometry.

The redundancy of Lockman ISOPHOT observations allowed us to generate high resolution maps, rebinning the original data into a final map with pixel size of 15$`\times `$15 arcsec. The detector signal is distributed in a uniform way between the smaller pixels. This process allows a better determination of source positions.

The source detection is performed on the signal-to-noise maps, given by the ratio of the “unreconstructed” flux maps and the corresponding maps of the noise. For the source detection we do not need any calibrated map, a relative map is sufficient in order to find the positions of any positive brightness fluctuation (as discussed by Dole et al., 2001).

First, our task selects all pixels above a low flux threshold (0.6 mJy pixel<sup>-1</sup>) using the IDL Astronomy Users Library task called $`find`$ (based on DAOPHOT’s equivalent algorithm). This algorithm finds positive brightness perturbations in an image, returning centroids and shape parameters (roundness and sharpness). The algorithm has been carefully tuned in order to detected all sources and to miss only the spurious ones that could be found surrounding the brighest sources or close to the edges of the maps (the input parameters are in particular the FWHM of the instrument and the limits for the roundness and sharpness geometric acceptable galaxy values).

Then we extract from the selected list only those objects having a signal-to-noise ratio $`>`$ 3.

In the final stage of our reduction we use our simulation procedures to reproject the sources detected on the raster map onto the pixel time history (see Section 3). In this way we are able to check the different temporal positions supposed to contribute to the total flux of each source. This method allows to improve the fit of the data for all the sources that appear above the interactive check threshold, and to significantly reduce the flux defect of the detected sources.

## 3 SIMULATIONS

The only way to assess the capability of our data reduction method for source detection and flux estimate is through simulations. We use the ISOPHOT C100 PSF (stored in the PIA file PC1FOOTP.FITS), rebinned in order to have the same pixel size (15”x15”) of the final maps in which the source extraction has been performed. Through the projection task used to produce maps, we project the PSF scaled to any given input flux on the raster maps. This corresponds to generate a synthetic source in the real map. Furthermore, our code is able to modify the time histories in those raster positions (of every pixel) contributing to the total flux of a given source. This is done by converting the given input flux excess (in respect to the local background) from mJy to ADU/gain/s, and adding it to the real pixel histories (containing glitches and noise). When this excess is added to the underlying signal, it is modelled by the algorithm that takes into account the transient response of the detector pixels (see Appendix C for details) and finally reproduces a real source, as it usually appears along any original pixel time history.

In our approach we make use of this powerful instrument at two different but complementary levels:

1- Simulation of detected sources in the same positions as they are detected on raster maps. As already anticipated in Section 2.4, with this procedure we can improve the data reduction and the reliability of our final sample. Simulating all detected sources with signal-to-noise ratio $`>`$3, we can check all the pieces of pixel histories supposed to contribute to them. In this way we can reject all spurious detections, and perform a better fit when needed.

2- Simulation of a sample of randomly located sources, in order to study the completeness and reliability of our detections at different flux levels, and estimate the internal calibration of the source photometry. The strategy we have adopted for the Lockman Hole is here described. We added 40 randomly distributed point sources at five different total fluxes (37, 75, 150, 300, 600 mJy). To avoid confusion, we imposed a minimum distance of 135 arcsec ($``$ 3 pixels) between the simulated and real sources on the maps. For the same reason, at each flux level we have distributed the 40 sources between the four rasters (each covering an area of $``$ 0.13 deg<sup>2</sup>): 10 sources per raster, divided in 2 simulation runs each containing 5 sources.

We reduced the simulated data cubes exactly in the same way as we did for the original data, doing the same checks and repairs. We produced the simulated maps on which we extracted the simulated sources, following the procedures used for real rasters. For each detected simulated source we have measured positions and peak fluxes. As in L01, the peak fluxes measured on the maps will be referred to as $`f_s`$ and $`f_{sr}`$ (both for real and synthetic sources) respectively for ‘unreconstructed’ and ‘reconstructed’ maps. The corresponding theoretical peak fluxes associated to the excess flux maps, not reduced and containing neither glitches or noise, will be named $`f_0`$ and $`f_{0r}`$. The theoretical quantities are produced only by the Lockman Hole observational strategy and the ISOPHOT instrument, while the measured quantities are also affected by our reduction method. These simulated data cubes contain both real sources and simulated ones. They have also the same noise, the ‘glitches’, and background transients as the original data.

In figure 4 (top) we compare the output peak fluxes obtained for the simulated sources affected only by the mapping effects ($`f_0`$) with the output fluxes of the reduced simulated sources ($`f_s`$). There is a linear correlation, and we observe that the reduced fluxes are always slightly lower than the unreduced ones. We find a similar correlation in the corresponding ‘reconstructed’ peak fluxes ($`f_{0r}`$ and $`f_{sr}`$), although for faint sources our algorithm is not able to reconstruct correctly the fluxes (see Figure 4, $`bottom`$). In order to have a correct flux determination at each level (both bright and faint), we will always use the ‘unreconstructed’ fluxes.

Figure 5 shows the distribution of the ratio between $`f_s`$ and $`f_0`$, for the simulated sources detected above $`3\sigma `$, compared to the same ratio for the detections above 300 mJy only. The peak of the distribution at 0.88 indicates a general underestimation of the total fluxes (derived from $`f_s`$). Two combined effects can explain this failure of the code to correctly compute total fluxes: one is the fact that $`f_0`$ is computed on the measured positions of $`f_s`$; the other is the underestimate of the wings of faint sources by the reduction method. The predicted and unbiased distribution of $`f_s/f_0`$ for the brighter simulated sources (above 300 mJy) peaks at 0.88. This value was subsequently assumed to correct our measured fluxes.

### 3.1 Completeness

As mentioned in the previous section, with simulations in the Lockman Hole we have derived the distribution of the measured ($`f_s`$) to theoretical ($`f_0`$) peak flux ratio. This distribution is crucial in deriving the completeness of the catalogue and the internal flux calibration, as it allows to predict the number of detected sources at a given flux level (the $`gfunction`$ described in detail in Gruppioni et al. 2002).

Our data reduction method can introduce some additional incompleteness if a source is interpreted as a background transient, and lost from the final source catalogue. We can estimate the incompleteness of our method from simulations, by computing the ratio between the number of detections and the number of expected sources in different peak flux intervals.

The results of the simulations are reported in Table 2 and shown in Figure 6, where the resulting function describing the incompleteness of our survey is plotted as a function of the simulated input flux. The loss of bright sources happens only when they are located at the extreme edges of the raster maps. We need to consider also the areal coverage of our survey (i.e. the fraction of the survey area where a source of peak flux $`f_s`$ can be detected with $`S/N>3\sigma `$), which is reported in Figure 7.

The global correction to be applied to the observed source counts is then obtained by convolving the function describing the completeness of our method with the areal coverage function (see L01 and Gruppioni et al. 2002 for details).

## 4 FLUX DETERMINATION

The simulations performed in the Lockman field provided not only the completeness of our detections at different flux levels, but also the internal calibration of the source photometry and the distribution of the ratio between the measured and the theoretical peak fluxes.

Here we summarize some relevant definitions and the relations used to derive the final total fluxes of real detected sources, as in L01 and Gruppioni et al. (2002):

* $`f_s`$: is the peak flux measured on maps for both real and simulated sources. Its value depends both on data reduction method, on the Lockman observing strategy, and the ISOPHOT instrumental effects;
* $`f_0`$: is the ‘theoretical’ peak flux measured on simulated maps containing neither glitches or noise. Its value depends only on Lockman observing strategy and the ISOPHOT instrumental effects;
* $`qf_s/f_0`$;
* $`q_{med}`$: is the peak of the $`f_s/f_0`$ distribution (also called systematic flux bias) and is 0.88. This value is used to correct the measured flux densities;
* the flux density $`s`$ of a source is computed by applying a correction factor to the measured peak flux $`f_s`$ in order to have a measure of its ’total’ flux. This can be done using the output informations stored after the simulation process of real sources in the positions where they have been detected (see Section 3). This gives the ratio between the injected simulated total flux $`s_0`$ and the corresponding theoretical peak flux $`f_0`$ measured on the simulated maps. This ratio represents the correction needed to derive total fluxes from peak fluxes. In this way the flux density is:

$$s=\left(f_s\times \frac{s_0}{f_0}\right)/q_{med}$$
(5)
For simulations $`s_0`$ is the injected total flux, while for real data $`s_0`$ is derived through successive iterations starting from a rough estimated value:

$$s_0=f_s/f_s/s_0_{sim}$$
(6)
where $`f_s/s_0_{sim}=0.23`$ is the average value resulting from simulating a point source of unitary flux. We can consider $`s`$ as the measured flux density and $`s_0`$ as the ‘true’ flux density of a source;
* finally we need to increase the derived flux density of a factor 1.2. This correction is needed because the PSF we use to simulate sources is undersampled and it misses the flux presents in the external wings of an ISOPHOT source profile.

### 4.1 Calibration

We have computed a new estimate of the detector responsivities, that converts digital units \[ADU/gain/s\] to physical flux units \[mJy\]. The detailed description is reported in Appendix B. In order to check the general consistency of our calibration, we have reduced a set of external calibrators with our procedures, by following the same steps described in the previous sections.

It is very difficult to find good “standards” for far-infrared observations. An absolute calibration is hampered by the intrinsic uncertainties and by the low sensitivities of previous measurements (IRAS, COBE-DIRBE), especially at faint fluxes. Our set of “calibrators” includes four stars and two IRAS sources, as reported in Table 3. For IRAS sources we have a direct measure of the far-infrared flux (from the IRAS Point Source Catalogue). To get the 95 $`\mu `$m flux we made an interpolation between the fluxes at 60 and 100 $`\mu `$m. For stars we compare our fluxes with predictions from Spectral Energy Distribution (SED) models. In particular we select stars from the ISO Ground-Based Preparatory Programme (GBPP, Jourdain de Muizon, M. and Habing, H. J., 1992). Synthetic SEDs for these stars are available in the ISO calibration Web page ($`http://www.iso.vilspa.esa.es/users/expl\mathrm{\_}lib/ISO/wwwcal/`$). The optical and near-infrared photometry observed in the GBPP have been fitted with a Kurucz model to provide the flux densities at longer wavelengths, thus extending the SED out to 300 $`\mu `$m (Cohen et al., 1996). We have chosen to adopt these SEDs model fluxes for use in our calibration at 95 $`\mu `$m, after convolution with the PHOT C\_95 filter. A summary of the fluxes obtained with our reduction and the expected fluxes are reported in Table 3. Quite good agreement with the comparison fluxes is found over a wide range of fluxes. Only for source F10507+5723, the only IRAS galaxy in the Lockman Hole, our value is a factor 2 lower than the interpolated IRAS flux. Given the uncertainties of the IRAS 100 $`\mu `$m flux ($`1020`$%), we have chosen to assume our calibration as final estimate of fluxes, without applying any other scaling factor. Excluding this IRAS source, Table 3 shows that, over a wide range of fluxes, the error on the photometry is within 20%. A more secure constraint on the calibration will be available after the reduction of ELAIS 95 $`\mu `$m fields ($`15`$ deg<sup>2</sup>) that will provide a wider sample of IRAS galaxies (few tens).

### 4.2 Flux errors

To compute the errors associated to our flux estimates, we have taken into account the two major contributions to the uncertainties. The first one depends on the data reduction method, given by the distribution of the ratio between fluxes measured after the reduction ($`f_s`$), and those measured taking into account only mapping effects ($`f_0`$). If $`f_0`$ is the real flux and if we measure after the data reduction $`f_s`$, it means that the reduction has introduced an error in the measure (by modifing the real flux). We can estimate this error as the width of the distribution of $`f_s/f_0`$. The second error term is due to the noise of the map.

Combining these two quantities, we get the final errors on the photometry (as reported in Table 5). The median of the errors distribution is $`20`$%.

### 4.3 Positional accuracy

With the set of simulations used to derive the completeness and the photometric corrections, we can also estimate the uncertainties on the astrometric positions. For each simulated source we have the injected and the measured coordinates. The local background noise can affect the estimate of the position of each galaxy centroid. The resulting effect is that a source, simulated in a given position, will be detected by the extraction algorithm in a slightly but different location.

Figure 8 shows the distribution of the differences in RA (top) and DEC (bottom) between the injected and the estimated positions for the simulated sources. The mean positional accuracy is of the order of 20 arcsec. By fitting the histograms in Fig. 8 with a Gaussian function, the 3-$`\sigma `$ of the distributions are respectively 18 arcsec for the right ascension and 17 arcsec for the declination. Of the 150 simulated and detected sources, 90 (60%) lie within \[-10,+10\] arcsec.

## 5 OBSERVATIONS OF AN AREA IN THE LOCKMAN HOLE

The 95 $`\mu `$m observations of the Lockman Hole field (P.I. Y. Taniguchi) by the photo-polarimeter ISOPHOT represent one of the best dataset available in the ISO archive to test the performance of our reduction technique. The long elementary integration times ($`16`$ sec for each raster position) and the observing redundancy allow us an accurate evaluation and modelling of any transient effects, for both long and short time scales.

### 5.1 The ISOPHOT observation strategy

The Lockman Hole (Lockman et al., 1986) was selected for its high ecliptic latitude ($`|\beta |>50`$), to keep the Zodiacal dust emission at the minimum, and for the low cirrus emission. This region presents the lowest HI column density in the sky, hence is particularly suited for the detection of faint infrared extragalactic sources. This consideration triggered a number of multifrequency observing campaigns in the past several years. The spectral coverage includes the X-rays (e.g. Hasinger et al., 2001), the optical (e.g. Fadda et al. 2002, in preparation), the mid-infrared (Fadda et al., 2002), the far-infrared (Kawara et al. 1998, and this work), the submillimiter (Scott et al., 2002), and the radio bands (De Ruiter et al., 1997; Ciliegi et al. 2003, in preparation). The published spectroscopic information is on the contrary still sparse.

Two different regions in the Lockman Hole have been observed by ISOPHOT. Each one of the two fields, called LHEX and LHNW, covers an area of $``$ 44’x44’ and has been surveyed at two far-infrared wavelengths with the C100 and C200 detector (respectively at 95 and 175 $`\mu `$m, with the C\_90 and C\_160 filters), in the P22 survey raster mode (see Kawara et al., 1998). We focused our analysis on the LHEX field, which consists in a mosaic of four rasters, each one covering an area of $``$ 22’x22’. The ISOPHOT detector was moved across the sky describing a grid pattern, with about half detector steps (corresponding to 1.5 detector pixel, or 67 arcsec) in both directions. This strategy improves the reliability of source detections and the image quality, as each sky position is observed twice in successive pointings. Table 4 summarizes the observational parameters for the C\_90 filter. These data have been retrieved from the ISO data archive through the WEB interface $`http://isowww.estec.esa.nl/`$.

### 5.2 The 95 $`\mu `$m source catalogue in the Lockman Hole

The final catalogue obtained with our method contains 36 sources detected at 95 $`\mu `$m in the Lockman Hole LHEX, over an area of 0.5 deg<sup>2</sup>. All sources have a signal-to-noise ratio greater than 3 and a flux greater than 16 mJy.

In Figure 9 we show the final mosaiced map obtained by combining together the four rasters. The open circles, whose sizes are roughly proportional to the source fluxes, indicate our detected sources. IAU-conformal names, sky coordinates (right ascension and declination at Equinox J2000), the detection significance (signal-to-noise ratio), the 95 $`\mu `$m total fluxes (in mJy) and their uncertainties are reported in Table 5. As previously discussed, all sources have been extracted from the map and confirmed by visual inspection on the pixel history (by two independent people). This approach produces an highly reliable source catalogue.

In a forthcoming paper we will discuss the optical, radio, mid-IR identifications of the Lockman 95 $`\mu `$m sources with radio and ISOCAM counterparts (Rodighiero et al. 2003, in preparation; Fadda et al. 2003, in preparation; Aussel et al. 2003, in preparation) and we will analyse the nature of our far-IR sources and their redshift distribution where the spectroscopic information is available.

### 5.3 Source counts from the Lockman Hole LHEX 95 $`\mu `$m survey

We concentrate in the present paper on discussing the statistical properties of the sample, like the source counts, confusion, and the contribution of the detected sources to the cosmic far-IR background (CIRB).

In the small area covered by the present study ($``$ 0.5 deg<sup>2</sup>), we have computed the 95 $`\mu `$m source counts down to a flux level of 30 mJy (in total 32 sources have been taken into account for this analysis). The integral counts have been obtained by weighting each single source for the effective area $`A_{eff}`$ corresponding to that source flux (as derived in Sect. 3.1). The errors associated with the counts in each level have been computed as $`\sqrt{\mathrm{\Sigma }_i(1/A_{eff}^2(S_i))}`$ (Gruppioni et al., 2002), where the sum is for all the sources with flux density $`S_i`$ and $`A_{eff}(S_i)`$ is the effective area. So the contributions of each source to both the counts and the associated errors are weighted for the area within which the source is detectable. These errors represent in any case the Poissonian term of the uncertainties, and have to be considered as lower limits to the total errors.

Our estimated values of the integral counts at different flux levels are plotted in Figure 10a as starred symbols. Our results are compared here with those from other surveys: the preliminary analysis of the ISOPHOT ELAIS survey (Efsthatiou et al. 2000, open circles), and the counts derived from the IRAS 100 $`\mu `$m survey (open squares). Our data are in excellent agreement with these results in the flux range in common.

The slope of the counts is $`\alpha 1.6`$.

A comparison of these integral counts with those published by Linden-Vornle et al. (2000), Kawara et al. (1998), Juvela et al. (2000) is reported in Figure 10b. We see quite a substantial scatter in these data. We believe that our improved analysis and careful check of all systematic and noise terms have produced a most reliable outcome. Obviously, our results are limited by the small survey area of about 0.5 deg<sup>2</sup> and source statistics. However, we find encouraging our excellent agreement with the source counts by Efstathiou et al. (2000), which are based on a far larger survey area. We take this to indicate that our deeper survey should not be biased too much by pathological clustering effects.

In Figure 11 we report the differential 95 $`\mu `$m counts dN/dS normalized to the Euclidean law ($`NS^{2.5}`$), providing a statistically independent dataset to be compared with model predictions (the data are reported in Table 6). A comparison is performed with modellistic differential counts by Franceschini et al. (2001), Xu et al. (2001) and Lagache et al. (2002).

The spatial distribution of the sources in our map of Fig. 9 is clearly non-random. We see in particular significant clustering in the East sector of the map. The number of sources in the 4 quadrants are 7, 7, 13 and 9 going in a clockwise order from the upper-right corner. In this situation, the evaluation of the source confusion noise requires some care. Using the beam of the ISOPHOT C100 (see Sect. 3) detector that has a FWHM of 45 arcsec, we have estimated that at the flux limit of 20 mJy for each source in the map there are $``$90 independent cells. This value is above the formal confusion limit, which is classically reached for a source areal density of $``$1/(30 independent beams) assuming Euclidean number counts (see Franceschini 1982 and Franceschini 2000, section 8.3). Around 20 mJy the Lockman 95$`\mu `$m counts are still close to the Euclidean regime (see Fig. 11). In the map’s quadrant with the highest number of sources, the areal density is $``$1/(45 independent beams), still above confusion. We conclude that our 95 $`\mu `$m map achieves a sensitivity close to the confusion limit in its most crowded parts, but still should not be much affected by the confusion noise. This is partly due to the significant incompleteness that our survey suffers at the faintest flux limits. For an ideal complete survey, our best-fit model implies a (3$`\sigma `$) confusion limit of $``$20 mJy occurring at an areal density of 30 beams/source. Kiss et al. (2001) have computed the $`3\sigma `$ confusion noise for ISOPHOT C100 90$`\mu `$m around $`21\pm 2`$ mJy, while Matsuhara et al. (2000) report a value $`30`$ mJy . All these estimates confirm that our survey should not be confusion limited.

It may be instructive to compare these figures for the ISO 90 $`\mu `$m selection with the confusion noise at longer wavelengths for observations with ISOPHOT C200 170$`\mu `$m. From the analysis of the FIRBACK fields Kiss et al. (2001) report a $`3\sigma `$ value of $`45\pm 4`$ mJy, while Dole et al. (2001) estimate a $`3\sigma `$ confusion noise of 135 mJy. This discrepancy is partly reconciled by considering that, in the computation of the confusion limit, Kiss et al. have used the default ISOPHOT PSF, while Dole et al. have modeled the ISOPHOT 170$`\mu `$m beam, thus recovering the flux fraction stored in the external wings of an ideal source (which indeed is lost by using the standard PSF). This confirms that the C100 imager is much less affected by confusion noise with respect to the C200’s ($`20`$ versus $`135`$ mJy the respective limits). For this reason, in spite of the more severe problems related to the C100 data, the 90 $`\mu `$m ISOPHOT observations provide in principle a deeper view and smaller error-boxes compared with longer wavelength observations.

### 5.4 Source counts interpretation

We have compared in Figs. 10 and 11 our determined source counts with various modellistic estimates. The long dashed lines in Figs. 10b and 11, in particular, show a comparison with the predictions for a non-evolving source population: our observed counts reveal a significant excess above these curves at the faintest fluxes, whose significance is better appreciated in Fig. 11 in terms of the independent flux bins of the differential counts. These results then confirm and substantiate earlier claims for the existence of an evolving population of IR galaxies, as previously identified in deep ISOCAM mid-IR and ISOPHOT 175 $`\mu `$m counts.

A comparison is made in Fig. 10a and 11 with predicted counts by Lagache et al. (2002) (dot-dashed lines), Xu et al. (2001) (dashed line), and by Rowan-Robinson (2001) (dotted line in Fig. 10a). These models typically assume that the whole local galaxy population evolves back in cosmic time in source luminosity or number density (with the exception of the Lagache model which accounts for both luminosity and density evolution). All these curves fit well the bright IRAS counts, but tend to more or less exceed those at fainter fluxes.

A better fit is provided by the multi-wavelength evolution model of Franceschini et al. (2001, hereafted F01) (solid lines). This model was designed to reproduce in particular the observed statistics (counts, z-distributions, luminosity functions) of the ISOCAM mid-IR selected sources, but it also accounts for data at other IR and sub-millimetric wavelengths. This model assumes the existence of three basic populations of cosmic sources characterized by different physical and evolutionary properties (their separate contributions are detailed in Figs. 11 and 12). The main contributions come from non-evolving quiescent spirals (long dashed line in Figs. 11 and 12) and from a population of fast evolving sources (dotted line in Fig. 12), including starburst galaxies and type-II AGNs (a third component considered - but always statistically negligible - are type-I AGNs, dot-dashed line in Fig. 12). The fraction of the evolving starburst population in the local universe is assumed to be $``$10 percent of the total, consistent with the local observed fraction of interacting galaxies.

In this scenario, the active starbursts and the quiescent galaxies belong to the same population. Each galaxy is expected to spend most of its lifetime in the quiescent state, but occasionally interactions and mergers with other galaxies put it in a short-lived (few to several $`10^7`$ years) active starbursting phase. The inferred cosmological evolution for the latter may be interpreted as an increased chance to detect a galaxy during the active phase back in the past, possibly due to a simple geometrical effect increasing the probability of interactions during the past denser epochs.

F01 argued that this two-population model was needed to explain the particular shape observed for the 15 $`\mu `$m counts. These keep a roughly euclidean slope down to few mJy and quickly turn up at fainter fluxes (Elbaz et al. 1999; Gruppioni et al. 2002; Fadda et al. 2003 in preparation). The two-population scheme (one evolving, one not) proposed by F01 allowed a best reproduction of these data. Our new results about the galaxy counts at 95 $`\mu `$m seem to confirm this model against the predictions of single-population schemes.

If we now integrate our estimated 95 $`\mu `$m galaxy counts down to the limit $`S_{lim}=20`$ mJy, we get a contribution $`I`$ of our resolved sources to the CIRB of $`I=_{\mathrm{}}^{S_{lim}}S𝑑N/𝑑S𝑑S2.0nW/m^2/sr`$. The CIRB intensity is rather uncertain at 100 $`\mu `$m due to the uncertain Zodiacal and Galactic contributions, published values ranging from $``$11 (Lagache et al. 1999) to $``$22 $`nW/m^2/sr`$ (Hauser et al. 1998). Our estimated resolved contribution then corresponds, respectively, to $``$18% to $``$9% of the CIRB intensity close to its peak wavelength. By using the models of Rowan-Robinson et al. (2001) and Xu et al. (2001) we derive a contribution of $`I2.8nW/m^2/sr`$ and $`I5nW/m^2/sr`$, respectively, to the CIRB at the limit of $`S_{lim}=20`$ mJy. These values are consistent with the spread observed in the counts predictions for the different models at 20 mJy (a factor 3 between the Franceschini and the Xu models).

Deep confusion-limited maps with MIPS on SIRTF at 70 $`\mu `$m could not be compared with the CIRB intensity (not directly measurable at this wavelength), while at 160 $`\mu `$m they are expected to be limited by confusion at $`S_{160}80`$ mJy (Franceschini et al., 2001). A significant improvement will require the substantially better spatial resolution of the Herschel Space Observatory in 2007.

## 6 CONCLUSIONS

We have developed a new procedure to clean and reduce deep survey data obtained with the photo-polarimeter ISOPHOT C100 on board the Infrared Space Observatory. These deep imaging data would have in principle the advantage over longer wavelength ISOPHOT 175 $`\mu `$m observations of a better spatial resolution and a lower confusion noise, but there use was limited by highly unstable detector background and responsivities.

Our procedure consists of a parametric algorithm fitting the signal time history of each detector, and able to extract from it the background level, to identify the singularities induced by cosmic-rays impacts and by transient effects in the detectors, and to identify and extract real sky sources.

The source extraction process, completeness, the photometric and astrometric accuracies of the final catalogues have been tested by us with extensive sets of simulations, by inserting into the real image sources with known flux and position, and accounting for all the details of the procedure. The flux calibration was also verified by reducing with the same technique C100 observations of calibrating stars.

We have tested our procedure by re-analysing data from a deep imaging survey performed at 95 $`\mu `$m with ISOPHOT C100 over a 40$`{}_{}{}^{}\times 40^{}`$ area within the Lockman Hole. Within this area we detect thirty-six sources with S/N$`>3`$, making up a complete flux-limited sample for $`S_{95\mu m}100`$ mJy. Reliable sources are detected, with decreasing but well-controlled completeness, down to $`S_{95\mu m}20`$ mJy.

These results provide us with the currently deepest far-IR image of the extragalactic sky. We estimate from it source counts down to a flux of $`30`$ mJy, at which limit we evaluate that of the order of 10% to 20% of the cosmic IR background has been resolved into sources.

The 95 $`\mu `$m galaxy counts reveal a slope at $`S_{95\mu m}100`$mJy quite steeper than that expected for a non-evolving source population. These observed counts are consistent with those determined from ISO surveys at 15 and 175 $`\mu `$m (Elbaz et al. 1999; Gruppioni et al. 2002; Puget et al. 1999; Dole et al. 2001). The detailed shape of these counts constrains the evolving population to dominate only below $`100`$mJy, whereas at brighter fluxes the majority of the sources are expected to be massive spirals at moderate to low redshift. We will report on the identifications and physical analyses of the 95 $`\mu `$m sources in separate papers (Aussel et al. 2003, in preparation; Fadda et al. 2003, in preparation; Rodighiero et al. 2003, in preparation).

## Acknowledgments

We are grateful to M. Vaccari for his careful reading of a preliminary version of the paper. G.R. wants to thank F. Pozzi for her introduction to the problem of ISOPHOT calibration. We thank the anonymous referee for his helpful suggestions that improved the paper.

This work was partly supported by the ”POE” EC TMR Network Programme (HPRN-CT-2000-00138).

## Appendix A THE LINEARITY CORRECTION

One of the main aspects of the PHOT-C detector is its non-linear response to the incoming flux, which requires a suited correction that has to be carefully derived. We have tried to follow a different approach from that of the standard PIA pipeline. An accurate analysis of the good stability of the detectors with time, enable us to construct a linearity correction starting from the dataset available for the S1 ELAIS field (for a total area of $``$ 4.5 deg<sup>2</sup>) and compare the result to the PIA internal calibration.

In order to obtain such correction, we have divided the data in 121 Volts intervals (from -1.2 to 1.2, covering the full dynamical range of the detector). Inside each step of voltage we have computed the median of the differences in ADU/gain/s and normalized to the median over the whole dynamical range. This means that the zero point of our correction is referred to the voltage range where the vast majority of the observed data in the ELAIS S1 field is located ($`0.8<V<0.4`$). We have taken into account for the computation only data around the background (to neglect the influence of deep and high transients in the final correction), discarding every readout corresponding to distructive readings, to the first two points of each ramp, to the sources and to the saturated readings.

The result for an example pixel is shown in Figure 13, where PIA standard correction (dashed line) is compared with ours, as a function of voltage. The corrections reported are percentual. There is a quite good agreement at $`V<0`$ (the deviation is of the order of 10$`\%`$), considering the larger quantity of data used in the computation of the PIA correction.

To check the consistency and the quantitative effects of the two corrections, we have estimated the residual correction we get after the data have been processed with PIA. In other words: we have derived with the same prescriptions our correction from data already PIA corrected. Such ”residual” correction is shown if Figure 14 (continuous line) and compared with the usual PIA percentual correction (dashed line, as in Figure 13). It is clear that at $`V>0`$ the PIA calibration is underestimated, and that we need to apply a further second order correction for a better linearity at the higher voltages (that we took into account).

We have considered the possibility of a further dependence of the linearity on the ramp position. In a similar way to that of the voltage dependence, we have thus construct a second order ”geometrical” correction from the data, normalized to the central region of the ramp, that allows to correct any residual non-linearity (a mean correction for the dependence of the ramp index erases completely the effects on the background, but can leave residuals on the stronger sources).

## Appendix B FLUX CALIBRATION

In order to get a good flux calibration, we have tried to compute a new estimate of the detector pixels responsivities. We mean with responsivities those factors that convert digital units \[ADU/gain/s\] to physical flux units \[mJy\]. We made use of a wide set of internal lamps, the so called FCS (Fine Calibration Sources, ISOPHOT Handbook, Laureijs et al., 2002). In particular we have taken into account all FCS measurements associated to raster observations in the ELAIS fields, in the Lockman Hole and to other minor mini-rasters of point sources, in order to cover an extended range of voltage levels. For every raster here considered, observed in the P22 mode, there are two FCS measurements available: the first ($`FCS_1`$) taken before the science observation, and a second ($`FCS_2`$) at the end of the observation. As discussed in section 2.1, we are confident about the stability of the detectors as a function of time. This justifies our attempt to compute the responsivities from a set of observed data.

We have reduced the FCS data following the same procedures described in this paper for raster science observations. However, we have previously applied the dark-current correction and the reset-interval correction, both part of the standard PIA package (Gabriel & Acosta-Pulido, 1999; Schulz et al., 2002). Our good statistics of internal calibrators enable us to reject all bad observations (mostly due to the massive presence of cosmic rays), and take into account only those for which the fitting algoritm has not failed. The main problem when dealing with FCS data is their short total exposure time (usually 1024 readouts), that prevents to determine a correct value of the stabilization level and may translate in a understimate of the flux.

Once the lamps have been reduced, we proceed to the statistical analysis for the calculation of responsivities. Figure 15 shows that with our reduction $`FCS_1`$ and $`FCS_2`$ are well correlated over the whole range of voltages considered, and for all pixels. The values reported are in units of Volts/s and represent the median values of every pointing (each FCS is composed of a single pointing). Different symbols refer to different pixels. The solid line represent the 1 to 1 relation. There seems to be a slight tendency for $`FCS_2`$ values to be statistically higher than $`FCS_1`$ values, in particular at lower fluxes. This could be due to the long integration times that divides the two lamps in each observation. After the science observation the detector is still quite hot, and a residual bias level could rest in the slowly varing memory of each pixel, offseting upward the stabilization level of each FCS2. This effect might be erased, or at least reduced, with the avalaibility of longer FCS exposure times, that would improve the determination of the stabilization.

In order to compute the responsivities, we chose to use the average value between each pair of $`FCS_1`$ and $`FCS_2`$. In the following we will call these values $`V_{FCS}`$ \[Volts/s\]. $`V_{FCS}`$ are related to the responsivities via the equation

$$RESP(i)=V_{FCS}(i)K(i)/P_{FCS}(i)$$
(7)
where $`RESP`$ are the responsivities, $`K`$ is a costant which depends only on the pixel and accounts for the illumination matrix and other instrumental effects of the detector, $`P_{FCS}`$ are the Inband fluxes (which represent the expected fraction of energy filtered through the system detector+filter, and collected as output signal on the detector), expressed in units of Watts. The $`i`$ index indicates the dependence of each variable on the different nine pixels. The previous equation is a syntax simplified version of the formula reported by Schulz et al. (2002) and in the ISOPHOT Handbook.

In order to compute the inband fluxes we have used the standard PIA power curve (Schulz et al. 2002), derived from a careful and complete analysis of external calibrators (stars, planets and asteroids). This calibration curve allows to correct the given heating powers (the total energy emitted from the FCS), to inband fluxes.

We can rewrite equation 7 as:

$$V_{FCS}(i)=RESP(i)P_{FCS}(i)/K(i)$$
(8)
and compute for every pixel the best linear fit between the two independent variables $`V_{FCS}(i)`$ and $`P_{FCS}(i)/K(i)`$ (each observation of our selected calibration data-set gives a point in this linear relation). The slope of this relation ($`RESP(i)`$) represents the responsivities. It is clear that the responsivities we get with this procedure are only a function of the pixel, and are assumed to be constant through the temporal history of the satellite. The values we derived are consistent with that of the standard PIA values.

The responsivities function of the pixels represent our best flat-field.

Finally, we can use our responsivities to calibrate our maps (to convert the signal from ADU/gain/s to flux units of mJy/pixel).

In order to check our calibration and look for any further physical scale factor (the absolute calibration) we have studied and reduced a few external calibrators (see discussion in Section 4.1). The good agreement we found at different flux levels and the intrinsic uncertainties on previous 95 $`\mu `$m measurements, made us confident to use the responsivities without appling any offset correction. We have estimated the photometric errors of our calibration throughout simulations and found they are of the order of 20-30 percent.

## Appendix C DETAILS ON THE SIMULATION PROCESS

To correctly simulate a source on a raster map and to determine its profile along the pixel time history, we make use of the ISOPHOT C100 PSF and of the projection algorithm. PHOT readout frequency is higher with respect to ISOCAM (1/32 seconds versus 2 seconds), the pixel size is greater (45 x 45 arcsec) and the distance between adjacent raster pointings is quite large (of the order of half detector array, 69 arcsec, in the case of the Lockman Hole, greater for other surveys like ELAIS). The combination of these effects requires the simulation of an high resolution source profile. In the ISOCAM procedure, the predicted flux level of a source was projected on every raster pointing and assumed to be there constant. For PHOT we have created a second raster structure with astrometric informations on each readout along the time history, including those where the satellite moves from a raster position to the other (not on target readouts). The projection of the PSF on these more gridded series of pointings and the lecture of these levels projected on every readout, is able to reproduce the exact profile of how the detector should see a source (if its response would be linear to the incoming flux). The detail description of this profile is mainly important between two raster pointings, when the detector moves on the sky. Given the size of the PHOT pixels, the detector starts to see a source when it is still moving and not yet positioned in a raster pointing. If we do not take into account this flux recorded by the detector (when the source enters and leaves a pixel) we could underestimate the total flux of the source. This correction is crucial when simulations are used to compute the corrections to the fluxes of the detected sources. In the following step, the projection algorithm we use takes into account the transient behaviour of the detector, and ”models” the simulated source profile along the time history in order to introduce the instrumental effects and obtain sources very similar to the real ones.

Our ”simulator” is very efficient, as it can describe sources independently of their position inside the pixel.