J Delabrouille

M Betoule

J.-B Melin

M.-A Miville-Deschênes

J Gonzalez-Nuevo

M Le Jeune

G Castex

G De Zotti

S Basak

M Ashdown

J Aumont

C Baccigalupi

A J Banday

J.-P Bernard

F R Bouchet

D L Clements

A Dasilva

C Dickinson

F Dodu

K Dolag

F Elsner

L Fauvet

G Faÿ 18

G Giardino

S Leach

J Lesgourgues

M Liguori

J F Macías-Pérez

M Massardi

S Matarrese

P Mazzotta

L Montier

S Mottet

R Paladini

B Partridge

R Piffaretti

G Prezeau

S Prunet

S Ricciardi

M Roman

B Schaefer

The pre-launch Planck Sky Model: a model of sky emission at submillimetre to centimetre wavelengths

Keywords: cosmology, miscellaneous -radio continuum, general -submillimeter, general -cosmic background radiation -ISM, generalgalaxies, general

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The cosmic microwave background (CMB), which is relic radiation from the hot Big Bang, carries an image of the Universe at an age of about 380 000 years. CMB anisotropies reflect the inhomogeneities in the early Universe, and are thus an observable that is very important for constraining the parameters of the Big Bang model.

For this reason, a large number of experiments dedicated to CMB anisotropy detection and characterisation have observed the sky at wavelengths ranging from the sub-millimetre to a few centimetres (frequencies from a few to a few hundred GHz). Stimulated by CMB science, astronomical observations in this wavelength range have enriched many fields of scientific investigation, ranging from the understanding of emission from the Galactic interstellar medium (ISM), to the study of a large population of extragalactic objects, including radio galaxies, infrared galaxies, and clusters of galaxies.

The WMAP satellite, launched by NASA in June 2001, has surveyed the sky in five frequency bands, ranging from 22 to 94 GHz (Bennett et al. 2003a). The data from this mission have provided an all-sky map of CMB temperature and an unprecedented data set for the study of sky emission at millimetre wavelengths -both in total intensity and in polarisation (Jarosik et al. 2011). The Planck mission1 , launched by ESA in May 2009, has started to provide even better observations of the sky in nine frequency bands, ranging from 30 to 850 GHz (Tauber et al. 2010;Planck Collaboration 2011a). These observations, which will be released to the scientific community in early 2013, are expected to become the new reference in this frequency range.

It has long been recognised that planning future observations and developing methods for analysing CMB data sets both require realistic models and simulations of the sky emission and of its observations (Gawiser et al. 1998;Bouchet & Gispert 1999;Giardino et al. 2002;de Oliveira-Costa et al. 2008;Sehgal et al. 2007Sehgal et al. , 2010)). For preparing the analysis of CMB data sets such as those of the Planck mission, it has proven useful to put a model of sky emission together, the Planck Sky Model (PSM), which can be used to predict or simulate sky emission for various assumptions and various parameter sets. The pre-launch version of this model, based on publicly available data sets existing before Planck observations (see Sect. 2.3), is used in particular as a simulation tool in the context of Planck.

This paper summarises the outcome of this preparatory modelling and simulation effort. Software, and simulations are made available to the scientific community on a dedicated website2 . The current version of the PSM, described in the present paper, is v1.8. We plan to regularly update models and software tools on the basis of upcoming observations (in particular those of the Planck mission), and when additional sophistication becomes useful for upcoming experiments and scientific investigations. The PSM developers welcome suggestions for improving the models or implementing alternate options.

Although nothing prevents running the PSM to simulate sky emission at any frequency, our sky model is actually valid, at present, only for frequencies ranging from about 3 GHz to about 3 THz, i.e., wavelengths ranging between 100 µm and 10 cm. At lower frequencies, Faraday rotation is important and our polarised maps are not representative of the sky emission, as this effect is presently neglected. Intensity maps, however, should be reasonably accurate down to about 500 MHz. At frequencies higher than 3 THz, the complexity of the dust emission is not properly taken into account by the present model, for which the representation of dust emission is based on spectral fits valid below 3 THz.

The paper is organised as follows: in Sect. 2 we review the basic features of the model as well as the architecture of the package; in Sects. 3-5 we describe the simulation and implementation of CMB anisotropies, diffuse Galactic and compact object emissions; and finally, in Sect. 6, we summarise the status of our sky model at the time of this publication and draw some conclusions.

The model

Rationale and examples of application

The development of a sky model and sky-emission simulations has several complementary objectives. First, a prediction of signals is useful for the optimisation of planned instruments (in terms of resolution, number and characteristics of frequency channels, acceptable noise level, and in general for making the necessary compromises between instrumental performance and costs or feasibility). Secondly, an emission model is necessary to simulate representative "observed" data sets, i.e., simulate observations of the modelled sky, which can then be used as test cases for developing dedicated data analysis software and for testing data analysis pipelines. Finally, such a model is a tool for easy comparison of any observed data set (at the time it becomes available) with what is expected based on previous empirical or theoretical understanding. Because of these different objectives, two different philosophies co-exist in the development of the model.

For some applications, we address the problem of making an accurate prediction of sky emission, considering existing observations and present knowledge about the emission mechanisms. The objective is then to give the best possible representation of our sky as it has been observed. Such a prediction can be used as a basis for calibration of future data against a common model, for selecting sky areas to be observed or for cleaning observations from expected contamination by a particular component. A prediction of sky emission is strongly data-driven; when only upper limits are available on the basis of existing observations (e.g., CMB B-modes, cluster peculiar velocities, ...), we assume our knowledge to be non-existent, and a prediction is not possible.

For other applications, such as the development of data processing and analysis pipelines, or comparison of methods such as in Leach et al. (2008), we need a realistic, statistically representative, sky simulation, with the appropriate complexity and with some level of randomness, compatible with current uncertainties in the observational data and their interpretation in terms of a model.

Our model is designed to implement both sky emission predictions and simulations. For practical uses, the two are combined to achieve a realisation of the sky, at the desired complexity and resolution, in the context of a given experiment or general study. For each of the components, several options exist for the modelling, this feature allowing the user to investigate the impact of theoretical uncertainties in our interpretation and parametrisation of sky emission. The PSM can also be used to generate model data sets similar to those which could be produced by an instrument. Basic "sky observation" tools allow for the integration of the sky emission in frequency bands, smoothing with instrumental beams, and observing along scans. A comparison of "observed" sky predictions with actual data being collected by upcoming instruments allows for real-time assessment of the proper operation of experiments.

The sky simulation tool implemented with the PSM was originally developed for investigating component separation methods, building on early work by Bouchet & Gispert (1999) for the Phase A study of the Planck space mission. Sky observations targeting the measurement of CMB intensity and polarisation anisotropies in fact contain a superposition of emission from several other astrophysical sources: the diffuse ISM of our Galaxy; radiogalaxies and infrared galaxies; the Sunyaev Zeldovich (SZ) effect from the hot intra-cluster electron gas in massive clusters of galaxies. Separating these components into contributions from distinct processes is important for the interpretation of the observations. In the past 15 years, component separation methods developed by several authors have been making use of simulated sky maps to assess the performance of the proposed approaches, or to validate them before application to real data sets (see Leach et al. 2008;Delabrouille & Cardoso 2009;Dunkley et al. 2009, and references therein, for an overview).

It is likely that astrophysical confusion (i.e., contamination by other astrophysical emission), rather than by instrumental noise, will set the limit on the performance of upcoming CMB observations, such as those of the Planck mission. For this reason, component separation is a very active topic of research, for which accurate predictions or realistic simulations of the sky emission, such as those provided by the present model, are necessary.

Sky emission maps from various astrophysical components in the PSM can be generated for different models, all of which are regularly updated. For each model, the sky emission depends on a set of parameters. The generation of sky maps within the PSM allows for comparison of such variants with actual observations, and hence to distinguish between models, or evaluate the impact of their parameters. The C tool developed by Lewis & Bridle (2002) is an approach for constraining cosmological parameters from CMB observations. In the same spirit, the PSM could be used in Monte-Carlo simulations (Jarosik et al. 2011), or user-defined. CMB Gaussian (constrained or unconstrained), non-Gaussianity, lensing component. Synchrotron Haslam et al. (1982) normalised and extrapolated with a power law (Miville-Deschênes et al. 2008, model 4). Free-free Dickinson et al. (2003) with normalisation of Miville-Deschênes et al. (2008). Spinning dust Draine & Lazarian (1999) with normalisation of Miville-Deschênes et al. (2008). Thermal dust Two-component model 7 of Finkbeiner et al. (1999). CO CO J=1 → 0, J = 2 → 1, J = 3 → 2 using Dame et al. (2001) and line ratios. SZ Cluster surveys, cluster counts, semi-analytic simulation or N-body+hydrodynamical simulation.

Radio sources

Radio surveys extrapolated with power laws. Infrared sources IRAS survey modelled with modified blackbody emission laws. CIB Clustering properties of Negrello et al. (2004).

UCH  regions IRAS catalogue with radio power law and modified blackbody function.

to constrain parameters which govern the physics of other sky emission mechanisms.

Astrophysical components

The total emission of the sky in the 3 GHz-3 THz frequency range is modelled as the sum of emission from different processes, which we classify for convenience into three large categories:

1. cosmic microwave anisotropies (including the dipole); 2. Galactic diffuse emission; 3. emission from compact objects (external galaxies, clusters of galaxies, and Galactic compact sources).

Each identified process with a given emission law is considered as an independent astrophysical component. The current software implements the CMB dipole, CMB anisotropies, synchrotron, free-free, thermal dust, spinning dust, CO molecular lines, thermal SZ effect, kinetic SZ effect, radio sources, infrared sources, far infrared background and ultracompact H  regions.

Solar-system emission from the planets, their satellites, from a large number of small objects (asteroids), and from dust particles and grains in the ecliptic plane (which generate zodiacal light), are not modelled in the current version.

The emission of all components is represented using parametric forms. Some components are modelled with a fixed emission law which does not depend on the sky location. In the current implementation, this is the case for the CMB, the nonrelativistic SZ effects (thermal and kinetic), free-free and spinning dust emission. Synchrotron and thermal dust are modelled with emission laws which vary over the sky. Each point source is also given its own emission law.

A summary of the astrophysical components included in the model is given in Table 1.

Data sets used in the model

The models of sky emission implemented in the PSM rely on maps and object catalogues obtained from observations made with different instruments. The main data sets serving as a basis for this are full-sky (or near full-sky) observations of the Galactic diffuse emission at various frequencies, and catalogues of observed compact objects. The selection of the data used in the model is based on a practical compromise between availability and usefulness for the present implementation, and will evolve with future releases.

We currently use: maps from IRAS, DIRBE and WMAP, as detailed below; the 408 MHz survey of Haslam et al. (1982)(version from the LAMBDA server3 , reprocessed to remove point sources and striping); maps of CO emission from Dame et al. (2001); and the Hα maps of Finkbeiner (2003). Most of these data sets, in HEALP format (Górski et al. 2005), have been downloaded from the LAMBDA web site.

The modelling of compact sources makes use of several radio surveys (listed in Sect. 5.2,T a b l e2), of the IRAS catalogue of infrared sources [START_REF] Beichman | IRAS Catalogs and Atlases[END_REF], and of catalogues of galaxy clusters observed with ROSAT (compilation from Piffaretti et al. 2011), as well as with the Sloan Digital Sky Survey (SDSS, Koester et al. 2007). Details about the use of each of the above data sets are postponed to the description of individual components below.

Sky prediction

Rather than simply interpolating existing observations, sky-emission predictions with the PSM implement recipes for combining available data to produce a set of modelled astrophysical components. Model predictions hence decompose observations into components, each of which is described according to a parametric model. Such predictions are deterministic, always returning the same products once global parameters (frequency, pixelisation scheme, coordinate system, model used to interpret the data) are set.

The accuracy of predictions is limited by two factors: the quality of sky observations (resolution, noise levels and systematic effects); and the appropriateness of the model assumed to make the predictions.

Predicted maps of emission at a given frequency are designed to best match the real sky emission as would be observed at a target resolution specified by the user. The effective resolution of the map, however, depends on the available data used to make the prediction. Prediction makes no attempt to extrapolate observations to smaller scales, whereas simulation does.

Predicted synchrotron maps are at present determined by 408 MHz observations and by WMAP observations in the lowest frequency channels, and hence have a resolution of about 1 • . Predicted dust maps, based on Model 7 or 8 of Finkbeiner et al. (1999), have a resolution of about 9 ′ (except in the 4 percent of the sky not covered by IRAS, where the resolution is that of DIRBE, i.e., about 40 ′ ). The predicted CMB is obtained from WMAP 5-year data, and has a resolution of order 28 ′ ,a f t e r applying a Wiener filter to minimise the integrated error in the reconstructed CMB map (see Sect. 3.2.2). A predicted SZ effect is obtained from known ROSAT and SDSS clusters, for which a model of SZ emission is produced from the universal cluster pressure profile, as presented by Arnaud et al. (2010). The scaling relations for estimating the SZ flux are normalised on the basis of WMAP and Planck data (Melin et al. 2011;Planck Collaboration 2011d,e). This prediction contains more than 15 000 clusters (see details in Sect. 5.1.1).

Simulations

Simulations generated with our model are random (or partially random) realisations. Each time a simulation is generated, a new sky is produced. It is possible however to re-generate exactly the same sky by rerunning the same version of the PSM, with the same sky model options and the same simulation seed (for random number generation).

For given cosmological parameters, the power spectrum of the CMB can be predicted using, e.g., the  software4 (Lewis et al. 2000), so that a plausible CMB sky can be simulated by drawing at random the a ℓm s of a simulated CMB, according to normal distributions with variance C ℓ . An interface to the CLASS5 software (Blas et al. 2011) is also implemented in the PSM. Specific realisations of sky emission can be produced for other components as well. A model of cluster counts as a function of mass and redshift, dN(M, z)/dMdz, can be used to simulate at random a population of clusters, and correspondingly, a model of SZ emission (Delabrouille et al. 2002a). The same is true for point sources, on the basis of number counts as a function of flux density.

Completely stochastic models of Galactic emission are more problematic, as the underlying statistical distribution of the emitting material is not known (if this concept makes sense at all). Nonetheless, a prescription to simulate realistic Galactic foreground emission is needed in any model of millimetre sky emission. For such components, we adopt an approach which uses the predictive model based on existing observations, and complements it with stochastic corrections representative of current uncertainties in the data and the model.

The generation of such constrained realisations,i nw h i c h simulations match observations within observational errors, is implemented for most sky emission components in the present model (see, for example, paragraphs about CMB anisotropies in Sect. 3.2.3, galactic emission with small scales in Sect. 4.6,and constrained thermal SZ effect in Sect. 5.1.2).

Use of the PSM in the Planck collaboration

The PSM is the standard sky modelling tool used for simulations in the Planck collaboration. Maps and catalogues obtained from the PSM are either used directly in component separation separation challenges (Leach et al. 2008) or used as inputs to the Level-S software for simulating Planck timelines (Reinecke et al. 2006). The PSM software package includes an interface with descriptions of the Planck HFI and LFI instruments for the generation of band-integrated sky maps, and the calculation of colour-correction coefficients.

The model in practice

In practice, the present sky model consists of a collection of data sets (maps and catalogues) and a software package, which can be used to produce maps of estimated or simulated sky emission. The software presently available is written mostly in the IDL programming language6 , and amounts to about 50 000 lines of code. Input data sets are stored in the form of ASCII and FITS files, for a total of about 500 GB of data. The software is designed to be used on a modest computer (e.g., a standard laptop); its data files are downloaded automatically from the associated data repository during the execution itself, which avoids storing on local disks large data files which are rarely used (i.e., those that are used only for very specific values of input parameters). The HEALP sky pixelisation scheme (Górski et al. 2005), in Galactic coordinates, is used in the current implementation.

All the input parameters are passed using a single configuration file, which is read upon execution of the code. This configuration file is stored with the outputs of the simulation for traceability of the simulation. The software consistently uses one single master seed for the generation of all the necessary chains of random numbers, which guarantees the reproducibility of a set of simulations, and the independence of the various random numbers used for the generation of different components of sky emission and instrumental noise.

The simulations can be generated at various resolutions, for various values of HEALP pixel size, and for various values of maximum harmonic mode number ℓ max . In principle, there is no limitation to the resolution (which can be set to 0), or to the HEALP pixel size (which is limited only by the HEALP software itself). The maximum CMB harmonic mode is set by the  (or CLASS) software (reference CMB power spectra have been precomputed for ℓ max = 10 000). In practice, the model is constrained by observational data on scales of 6 ′ to 1 • for intensity, and few degrees for polarisation, depending on the component considered. The accuracy of extrapolations at smaller scales are expected to grow worse with increased resolutions, and the model is expected to be useful down to the resolution of the Planck HFI (about 4 ′ resolution). Test cases ha v ebeenrunuptoHEALP N side = 4096, and ℓ max = 6000. Technical information about the model, released software, data sets and documentation can be obtained from the PSM web site. PSM developers will make all possible efforts to maintain the code, document the consecutive versions, and give information about limitations and bugs on this website as well.

The code runs in three consecutive steps: generation of a sky model; generation of sky emission maps at specific frequencies or in specified frequency bands, and generation of simulated data sets as observed by simulated instruments.

Step 1: the sky model

The sky model is generated at a finite resolution that is set by the user and is common to all components. The sky model therefore implements maps of the sky smoothed with a Gaussian beam of full width at half maximum (FWHM) θ sky . This ensures that all maps are near band-limited, and hence can be properly sampled (with an appropriate HEALP pixel size, which must be smaller than about one third of θ sky ).

A particular sky model uses a common pixelisation scheme for all components (HEALP, in Galactic coordinates, ring ordering, and a given pixel size appropriate for sampling maps with no power above a specific harmonic mode ℓ max ). Polarised emission modelling is optional.

The sky model consists of a set of maps of physical or phenomenological parameters describing the components: CMB temperature and polarisation anisotropy; synchrotron amplitude, polarisation, and spectral index maps; dust amplitude, polarisation and temperature maps; catalogues of point sources (described by their coordinates and by their emission laws for intensity and polarisation); catalogues of clusters (described by their coordinates, mass, redshift, temperature), etc. The sky model generated with the software is saved at the end of this first step; i.e., all the data sets used to model the sky are part of the output of a given run. It is then possible to apply steps two and three to the model at any time later on, without the need to re-run step one.

Step 2: sky emission maps

Sky emission maps result from the calculation of the observed emission at given frequencies, or in given frequency bands, on the basis of the sky model generated in the first step. The sky emission maps are what an ideal noiseless instrument would observe with the specified frequency bands and with ideal Gaussian beams with FWHM equal to the resolution of the model sky θ sky , andaresampledusingthesameHEALP scheme as component maps. Note that this step is not fully independent of Step 3, as it is the list of instruments used to "observe" the sky which sets the list of frequency bands used in Step 2.

Step 3: observations

Finally, sky emission maps generated in Step 2 can be "observed" by simple model instruments. Each channel of the instrument is defined with a frequency band, a beam, polarisation properties, an observing scheme, noise properties, and a format for the generated data sets. Although the generation of instrumental noise excludes at the present slow drifts correlated between detectors, which generate signals that interfere with component separation (Delabrouille et al. 2002b;Patanchon et al. 2008), it is possible to include the generation of such effects in future developments of the code.

The - software (Reinecke et al. 2006) can also be used as a means for sophisticated instrument simulation.

Limitations

The pre-launch PSM is based on our understanding of the sky prior to Planck observations. It is intended to make available representative maps, to be used essentially in simulations. It should not be used for inferring the properties of the real sky. Limitations specific to some of the models of the individual components are specified at the end of each of the relevant sections.

Cosmic microwave background

CMB anisotropies are modelled in the form of temperature fluctuations of a perfect blackbody. The assumed default CMB temperature is T CMB = 2.725 K (Fixsen et al. 1996;Fixsen 2009). Brightness fluctuations are modelled on the basis of a first order Taylor expansion of a blackbody spectrum, where ∆T (θ, φ) is a map of temperature anisotropies, and where the emission law A(ν) is the derivative with respect to temperature of the blackbody spectrum:

I CMB (θ, φ) = A CMB (ν)∆T (θ, φ), (1)
A CMB (ν) = ∂B ν (T ) ∂T T = T CMB . ( 2
)
We ignore any effects, such as a non-zero chemical potential, which move the distribution of the CMB photons away from a perfect blackbody.

The same emission law is used for CMB polarisation, which is represented with polarisation maps Q(θ, φ)a n dU(θ, φ)o r , alternatively, by the maps of scalar and pseudo-scalar modes, E(θ, φ)andB(θ, φ) (Zaldarriaga & Seljak 1997;Kamionkowski et al. 1997).

The CMB dipole

The CMB dipole is mostly due to the motion of the solar system with respect to the background. It is modelled on the basis of the measurement obtained with WMAP seven-year data (Jarosik et al. 2011). The measured amplitude and direction are

A dip = (3.355 ± 0.008) mK; l dip = 263.99 • ± 0.14 • ; b dip = 48.26 • ± 0.03 • . (3)
A dipole prediction uses by default these measured values.

Alternatively, these values can be supplied by the user for settingupadifferent predicted dipole.

A dipole simulation generates amplitude and direction assuming a Gaussian distribution with default average values and standard deviations as given by the WMAP measurement. Alternatively, the simulation can use amplitude, direction and error as supplied in the input parameter list.

Upon generation of a model sky, the amplitude of the CMB dipole map is modified to take into account the smoothing by the sky beam (even if for small beams the effect is negligible in practice). The default predicted CMB dipole generated by the PSM (here for a θ sky = 1 • sky map resolution) is displayed in Fig. 1.

Note that in addition to the dipole, the peculiar motion of the observer generates higher order terms (Doppler quadrupole and higher-order corrections), which are ignored in the current version of the PSM.

Temperature and polarisation anisotropies

On smaller scales, the CMB map is modelled as the outcome of a random process (the random generation of perturbations in the early universe). In the simplest model, fluctuations in the CMB are assumed to be Gaussian and stationary, so that its statistical distribution is entirely determined by its power spectrum C ℓ . In addition to this generic model, the software implements a non-Gaussian CMB model, as well as predictions and constrained realisations matching WMAP observations.

Gaussian CMB anisotropies

Anisotropies on scales smaller than the dipole are imprinted on the CMB mostly at z ≃ 1100, when primordial nuclei captured free electrons to form neutral atoms. In the simplest model, these primary anisotropies are assumed Gaussian and stationary. Additional secondary perturbations arise due to electromagnetic and gravitational interactions, while photons propagate towards us. Reionisation of the universe at late times smoothes the small scale structure of the acoustic peaks, and generates "reionisation bumps" in the polarisation power spectra. Gravitational interaction generates large-scale anisotropies via the integrated Sachs-Wolfe (ISW) effect, and modifies small scale anisotropies via lensing by large-scale structure (galaxies, and clusters of galaxies).

The statistics of the temperature and polarisation fluctuations of a Gaussian CMB are fully described by the multivariate temperature and polarisation power spectrum C ℓ of one temperature map T and two polarisation maps E and B. For a standard cosmology (ignoring such effects as, e.g., cosmological birefringence), the only non-vanishing terms in the power spectrum are C TT ℓ , C TE ℓ , C EE ℓ , C BB ℓ , the other two (C TB ℓ and C EB ℓ ) vanishing for parity reasons. Fast and accurate software is available to compute these power spectra for a given cosmological model, and to randomly generate CMB temperature and polarisation maps for any given such spectrum: in the current version, interfaces to the CAMB (Lewis et al. 2000)andCLASS (Blas et al. 2011) software are implemented to compute CMB power spectra for a given cosmological model. Alternatively, a default CMB power spectrum can be used, which matches current observational data. Figure 2 displays the power spectra used in the code for the current cosmological best fit, with no primordial tensor modes (and hence C BB ℓ from lensing of E modes only).

CMB prediction

CMB anisotropies have been observed by a number of experiments. In particular, CMB temperature maps have been obtained from publicly released WMAP data by a number of authors (e.g. Tegmark et al. 2003;Eriksen et al. 2004;Delabrouille et al. 2009;Basak & Delabrouille 2012). Each of these maps can be used to predict the expected level of CMB anisotropy at any point of the sky. In the current version of the model, we use the Wiener-filtered Needlet ILC (NILC5) map obtained from WMAP 5-year data by Delabrouille et al. (2009) as the standard CMB map at the resolution of the WMAP W-band channel. The harmonic coefficients x T ℓm of that map are connected to the true CMB harmonic coefficients a T ℓm by

x T ℓm = w ℓ B W ℓ a T ℓm + n ℓm , (4)
where B W ℓ are the coefficients of the expansion of the symmetrised beam of the WMAP W-band channel (as released with the WMAP 5-year data), w ℓ the coefficients of the harmonic Wiener filter applied to the NILC map, and n ℓm is a residual additive term that accounts for all sources of error in the NILC5 map (in particular instrumental noise and residual foreground emission).

A predicted map at a different sky resolution (as specified in the list of input parameters of the code) is obtained from this one by multiplying the spherical harmonic coefficients by B ℓ /B W ℓ , where B ℓ are the coefficients of the expansion of the equivalent Gaussian beam corresponding to the target sky resolution

āT ℓm = B ℓ B W ℓ x T ℓm . (5)
Predicted CMB polarisation maps for Q and U are obtained as well, on the basis of the model TE correlation. The predicted E-mode harmonic coefficients are

āE ℓm = C TE ℓ C TT ℓ āT ℓm , (6)
where āT ℓm are the predicted CMB temperature harmonic coefficients computed from Eq. ( 5). This polarisation prediction, however, is quite uncertain, since the TE correlation is weak. The predicted B-type polarisation vanishes in the current model implementation. Figure 3 shows the predicted CMB temperature and polarisation at a target sky resolution of 1 • .

The predicted CMB maps are produced in such a way that the error term (difference between predicted map, and true CMB map at the target resolution) is minimised. They hence provide t h eb e s tp o s s i b l ei m a g eo ft h et r u eC M Bs k yg i v e nt h ep r e s e n t observations used in the model. With this criterion however, the power spectrum of the predicted CMB maps does not match the theoretical model. Indeed, given a noisy observation x ℓm = s ℓm + n ℓm , the Wiener filter as a function of ℓ is w ℓ = S ℓ /(S ℓ + N ℓ ), where S ℓ is the power spectrum of the signal, and N ℓ that of the noise. Then the power spectra of the predicted CMB temperature map is i.e., the power spectrum is lower than expected, by a factor w ℓ . Similarly, the power spectrum of the predicted E-mode polarisation map is

C TT ℓ = āT ℓm 2 = w ℓ B 2 ℓ C TT ℓ , (7)
C EE ℓ = w ℓ C TE ℓ 2 C TT ℓ C EE ℓ B 2 ℓ C EE ℓ , (8)
and the cross spectrum of T and E is

C TE ℓ = w ℓ B 2 ℓ C TE ℓ . (9)

Constrained Gaussian realisations

Constrained CMB realisations, compatible with these observations, are generated in the PSM on the basis of a theoretical model (power spectra of temperature and polarisation), and constraints (the observed CMB temperature map described in Sect. 3.2.2).

Recall that if (X, Y) t is a two-dimensional random vector with mean (0, 0) and covariance matrix

σ 2 X σ XY σ XY σ 2 Y , (10)
then the conditional probability of X given Y is Gaussian with mean m and variance σ 2 given by

m = σ XY σ 2 Y Y, ( 11
)
and

σ 2 = σ 2 X - σ 2 XY σ 2 Y • (12)
We simulate the constrained a ℓm by drawing random variables following the conditional distribution of a ℓm given some observation a ℓm . The observation is performed at finite resolution, charact e r i s e db ya ne ffective beam B ℓ . Each a ℓm of the CMB map is also measured with an error ǫ ℓm . The multipole coefficients of the observed map are

a ℓm = B ℓ a ℓm + ǫ ℓm , (13)
where a ℓm and ǫ ℓm are centred and independent Gaussian random variables, with variance C ℓ and N ℓ , respectively. B ℓ is the response of the beam in harmonic space (the beam is assumed to be symmetric).

As cov(a ℓm , a ℓm ) = B ℓ C ℓ ,v a r ( a ℓm ) = C ℓ and var(â ℓm ) = B 2 ℓ C ℓ + N ℓ (we ignore correlations between the errors for different (ℓ, m), so that each mode can be generated independently of the others), we obtain, from the previous result, harmonic coefficients of the constrained realisation that are Gaussian with mean

āℓm = B ℓ C ℓ B 2 ℓ C ℓ + N ℓ a ℓm , (14)
and variance

σ 2 ℓ = C ℓ N ℓ B 2 ℓ C ℓ + N ℓ • (15)
We note that as N ℓ →∞or B ℓ → 0, this law becomes N(0, C ℓ ) (which is an unconstrained realisation of the CMB field), whereas if N ℓ = 0 the conditional mean is a ℓm /B ℓ .T h e mean value of the distribution of each harmonic coefficient matches the Wiener-filtered prediction described in Sect. 3.2.2.

Non-Gaussian CMB anisotropies

Non-Gaussian corrections to the statistics of CMB primary anisotropies are expected in the early Universe (inflationary) scenarios. Software to generate non-Gaussian CMB maps has been developed by Smith & Zaldarriaga (2007), Liguori et al. (2007) and Elsner & Wandelt (2009).

The first simulations of CMB temperature maps with primordial non-Gaussianity (NG) of the "local type" have been based on the generation of the underlying primordial perturbation in Fourier space (Komatsu et al. 2003). A different method has then been proposed in Liguori et al. (2003Liguori et al. ( , 2007)), where the authors work with filter functions to introduce the proper spatial correlations of the primordial potential. The current version of the model uses maps generated according to the method of Elsner & Wandelt (2009), a recent improvement of the method of Liguori et al. (2003) with better numerical efficiency.

The non-Gaussian CMB model assumes a linear-plusquadratic model for Bardeen's gauge-invariant curvature potential Φ7 , namely

Φ(x) =Φ L (x) + f NL Φ L (x) 2 -Φ L (x) 2 , ( 16
)
where Φ L is a Gaussian random field and the dimensionless nonlinearity parameter f NL sets the level of NG in the primordial potential. The kind of NG described by this linear-plus-quadratic model is theoretically motivated by detailed second-order calculations of the NG arising during -or soon after -inflation in the early Universe (Bartolo et al. 2004); the assumption that the NG level is set by a constant parameter f NL is a fair approximation, as long as | f NL |≫1.

Once realisations of the NG potential Φ are obtained, the harmonic coefficients for T and E are calculated by radial integration of the linear and non-linear parts of the potential independently, yielding maps of harmonic coefficients a ℓm for T and E and for both the linear and the non-linear part. The corresponding input templates are used to generate simulated CMB skies for different values of f NL .

Because the computation of a non-Gaussian simulation at Planck resolution can take up to about 20 CPU hours, the code to make such simulations is not directly included in the package. Instead, we use a number of pre-computed simulations with ℓ max = 3500, and pick among those to generate a simulated non-Gaussian CMB sky.

Full-sky maps of one of the non-Gaussian simulations included in the model are shown in Fig. 4, where the top panel is the linear part and the bottom panel the non-linear correction (for f NL = 1).

Lensed CMB

One of the most important physical mechanisms that generate secondary anisotropies is weak gravitational lensing of the CMB, which arises from the distortions induced in the geodesics of CMB photons by gradients in the gravitational matter potential. As a result, the CMB temperature anisotropy that we observe at a particular point on the sky in a direction n is coming from some other point on the last scattering surface in a displaced direction, n′ , such that,

X(n) = X(n ′ ), ( 17
) n′ = n + d(n), ( 18
)
where X stands for any of the lensed CMB Stokes parameters I, Q or U and X is its unlensed counterpart. The directions n and n′ are connected by the deflection field d(n).

In the Born approximation, the deflection field is related to the line-of-sight projection of the gravitational potential, where d A (χ) is the comoving angular diameter distance corresponding to the comoving distance χ,a n dχ s is the comoving distance to the last scattering surface. Although lensing depends upon all of the large-scale structure between us and the last scattering surface, most of the effect comes from structures of comoving size of order few hundred Mpc at redshifts below 10. The typical deflection angle is around 2 ′ to 3 ′ but is correlated over several degrees.

Ψ(d A (χ)n,χ)asd(n) = ∇Φ(n), such that, Φ(n) = -2 χ s 0 dχ d A (χ s -χ) d A (χ) d A (χ s ) Ψ(d A (χ)n,χ), (19)
Weak gravitational lensing has several important effects on the CMB (see, e.g., Lewis & Challinor 2006,forareview).Inaddition to being an extremely robust prediction of the cosmological concordance model, it is a probe of the structure formation process in cosmology, overlapping with the onset of cosmic acceleration, mainly in the linear or quasi-linear regime. The power spectrum of the deflection field is therefore useful for measuring parameters like the total neutrino mass (Lesgourgues et al. 2006) or the dark energy equation of state (see, e.g. Acquaviva & Baccigalupi 2006).

Besides its intrinsic cosmological interest, the distortion of primary CMB anisotropies by lensing is a source of confusion for several scientific objectives of sensitive CMB experiments. Lensing modifies the CMB power spectrum, and generates, through non-linear mode coupling, non-Gaussianity competing with the primordial one (Lesgourgues et al. 2005). It also mixes primordial E and B modes, the main effect being the conversion of a fraction of the E-modes into B-modes, which is a nuisance for the precise measurement of CMB B-modes from primordial gravitational waves.

Weak lensing of CMB anisotropies by large-scale structure (LSS) has been measured on WMAP data by cross-correlation with the NRAO VLA Sky Survey (NVSS), used as a high-redshift radio galaxy catalogue (Smith et al. 2007), as well as in a combined analysis with SDDS and NVSS (Hirata et al. 2008). Although marginally detectable in WMAP data, CMB lensing should be measured with high signal to noise by Planck without need to rely on an external data set (Hu & Okamoto 2002). However, in order to carry out this measurement in practice, the impact of instrumental (anisotropic beams, missing data, correlated noise) and astrophysical (Galaxy contamination, point sources, etc.) systematic effects on the CMB lensing estimators has to be studied with great care. This calls for the implementation, in the present model, of fast and accurate methods to simulate the lensing of CMB maps.

In principle, such lensing should be compatible with the distribution of matter between the last scattering surface and the observer. This matter distribution is connected to the distribution of galaxies contributing to the anisotropies of the cosmic infrared background, to the radio source background, and to high-redshift galaxy clusters, which are simulated as well by the PSM. This connection between the different components of sky emission is ignored in the current PSM implementation, but will be the object of future improvements.

Use of Non-equispaced FFT

In order to compute the lensed CMB field on HEALP grid points, we need to resample the unlensed CMB at arbitrary positions on the sphere. The method implemented for doing so is based on the possibility to recast remapping on a sphere into remapping on a 2-d torus, and uses the non-equispaced fast Fourier transform (NFFT) to compute lensed CMB anisotropies at HEALP grid points. The method is based on the previous work of Basak et al. (2009), which we briefly review here.

The basic idea of the NFFT is to combine the standard fast Fourier transform and linear combinations of a window function that is well localised in both the spatial and frequency domains. Suppose we know a function f by means of N evaluations f k in the frequency domain. According to the NFFT, the Fourier transform of that function evaluated at M non-equispaced grid points x j in the spatial domain can be written as,

f (x j ) = 1 √ 2π m∈Z Φ(σx j -m) × N/2-1 k=-N/2 exp - 2πimk σ N f k φ(2π k/σ N) , (20)
where φ(ξ) is a window function in Fourier space, Φ(ξ) its counterpart in pixel space, and σ is the oversampling factor. The indices j of the grid points x j take the values j = 1, 2, 3,...,M.

The efficient evaluation of f (x) on irregularly spaced grid points requires a window function Φ(ξ) that is well localised both in space and frequency domain. The computational complexity is O σN log N + KM ,w h e r eK is the convolution length, σ is the oversampling factor, M is the number of real space samples, and N is the number of Fourier modes. Among a number of window functions (Gaussian, B-spline, sinc-power, Kaiser-Bessel), the Kaiser-Bessel function is found to provide the most accurate results (Fourmont 2003;Kunis & Potts 2008), and is used in our simulations by default.

Simulation procedure

We summarise here the main steps of the simulation procedure of lensed CMB maps:

1. Generate a realisation of harmonic coefficients of the CMB, both temperature and polarisation, and lensing potential from their respective theoretical angular power spectra obtained using CAMB; 2. Transform the harmonic coefficients of the unlensed CMB fields and of the displacement field into their 2-d torus Fourier counterparts using the symmetry and recursion relation of the Wigner rotation matrix [START_REF] Varshalovich | Quantum theory of angular momentum[END_REF]Challinor & Lewis 2005); 3. Sample the displacement field at HEALP pixel centres using the NFFT on the 2-d torus, apply this displacement field to obtain displaced positions on the sphere, and compute the additional rotation needed for the polarised fields using identities of the spherical triangle (Lewis 2005); 4. Resample the temperature and polarisation fields at the displaced positions using the NFFT to obtain the simulated lensed CMB fields, sampled at the HEALP pixel centres.

In Fig. 5, we show one realisation of unlensed CMB temperature anisotropies, together with the lensing correction, the amplitude of the deflection field and the lensing potential map. We show a small portion of the same set of maps in Fig. 6 to illustrate the lensing effect more clearly. The correlation between the deflection field and the difference between the lensed and unlensed CMB temperature anisotropies is clearly visible. We have found that, for unlensed CMB and deflection field maps with N side = 1024 and ℓ max = 2048, 10 -8 accuracy is easily reached when setting the (σ, K) parameters to (2, 4). The current implementation of the method to simulate lensed CMB with these parameters requires 15 min wall-clock computing time, using 7.9 GB of RAM, on an ordinary PC configuration. Figure 7 shows that the average angular power spectra recovered from 1000 realisations of lensed CMB maps (N side = 1024) are consistent with the corresponding theoretical ones. This agreement between our numerical power spectra estimates and the CAMB predictions is a validation of both our numerical method and of the CAMB estimates based on partially re-summed perturbative calculations.

The accuracy of the lensing simulation is set by default in the current implementation, and is not a parameter of the software. If necessary it can be improved by increasing both the oversampling factor and the convolution length, at the expense of extra memory consumption and CPU time. Since the full precomputation of the window function at each node in the spatial domain requires storage of large amount of real numbers ([2K +1] 2 per pixel in 2-d), we have used a tensor product form for the multivariate window function (the default method within the NFFT library), which requires only unidimensional precomputations. This method uses a medium amount of memory to store 2[2K+1] real numbers per pixel, but at the cost of extra multiplication operations to compute the multivariate window function from the factors. In addition to that, pre-computation of window function at harmonic domain requires storage of 2[2ℓ max +2] real numbers.

An increase in the convolution length not only increases the computational cost of the interpolation part of the NFFT, but also increases the cost of the pre-computation of window function and memory requirement as one has to compute and store the window function at a larger number of grid points in the spatial domain before applying the NFFT. On the other hand, increasing The simulated displacement fields are computed from the gradients of Gaussian potential fields, neglecting non-linearities produced by the growth of structures and rotation induced by departure from the Born approximation. For CMB studies on scales comparable to or larger than arc-minutes, this approximation is excellent (see e.g. Bartelmann & Schneider 2001), and enables us to compare our lensed power spectrum estimates to the CAMB predictions, for which analytical estimates exist. However, our method is valid for any given displacement field; hence we could relax the Born approximation if needed, replacing it with ray-tracing in dark matter N-body simulations (Carbone et al. 2008(Carbone et al. , 2009)). Ray-tracing is affected by similar problems as the simulation of the lens effect on CMB maps, i.e., difficulties in accurately resampling a vector field on the sphere. Current state-of-the-art ray-tracing algorithms, such as described by Teyssier et al. (2009), could be made more accurate by using the technique described here.

It should be noted that the displacement field used to generate lensed CMB maps is sampled at the centres of HEALP pixels. For large pixels, the displacement field sampled at the pixel centre will not necessarily correspond to the average of the displacement field integrated over the pixel area (in other words, lensing and smoothing do not commute). For numerical accuracy, it is recommended to generate simulations with resolution and pixel sizes smaller than the typical displacement due to lensing.

While in general for temperature and polarisation E-modes the main effect of lensing is visible on small scales, it is not the case for CMB B-modes, for which lensing generates an appreciable amount of large-scale B-modes, as compared to primordial B-modes due to tensor perturbations. For accurate lensing B-modes, it is important to make a simulation that properly samples the E-polarisation field, i.e. maximum harmonic mode ℓ max ≃ 2000 or more, and HEALP nside of 1024 or more.

Limitations of the CMB models

The various CMB models implemented in the PSM suffer from the following limitations:

1. The prediction model is based on a CMB map extracted from WMAP data. That map is contaminated by noise and residuals of foreground emission, particularly in the vicinity of the Galactic plane. It has been Wiener-filtered to minimise the total error. This however results in a filtered CMB map, the power spectrum of which hence is not that expected for the cosmological model considered. 2. The CMB lensing and ISW as currently modelled are not connected to the distribution of large-scale structures. This will be developed in a future version. 3. The non-Gaussian CMB is limited to ℓ max = 3500, although the Gaussian signal can be extended to higher harmonic modes. 4. Lensing of the CMB is currently implemented in the PSM pipeline only for Gaussian CMB models. It is possible, however, to use the ilens code distributed with the PSM to lens Non-Gaussian CMB maps in a post-processing step.

Diffuse Galactic emission

Diffuse Galactic emission originates from the ISM of the Milky Way. The ISM is made up of cold atomic and molecular clouds, of a warm inter-cloud medium that is partly ionised, and of a hot ionised medium. The ISM also comprises magnetic fields and cosmic rays. Galactic emission from the ISM is usually separated into distinct components according to the physical emission process: synchrotron radiation from relativistic free electrons spiralling in the Galactic magnetic field; free-free emission from the warm ionised medium, due to the interaction of free electrons with positively charged nuclei; thermal (vibrational) emission from interstellar dust grains heated by radiation from stars; spinning dust emission from the dipole moment of rotating dust grains; line emission from atoms and molecules. As the interstellar matter is concentrated in the Galactic plane, the intensity of these diffuse emissions decreases with Galactic latitude approximately according to a cosecant law (the optical depth of the emitting material scales proportionally to 1/sin |b|).

The Galactic diffuse emission carries important information on the cycle of interstellar matter in the Milky Way and its link with the star formation process. The precise mapping of the microwave sky by the COBE and WMAP experiments and the detailed analysis of the Galactic contamination to the CMB anisotropies have contribued significantly to our understanding of the Galactic ISM (Wright et al. 1991;[START_REF] Bennett | [END_REF]Boulanger et al. 1996;Finkbeiner et al. 1999;Banday et al. 2003;Bennett et al. 2003b;Finkbeiner 2004;Davies et al. 2006;Miville-Deschênes et al. 2008). This has even led to the discovery of a new Galactic emission component: the anomalous microwave emission (Kogut et al. 1996a). The Planck experiment will certainly expand our knowledge on the Galactic ISM; it already has with its first series of papers (Planck Collaboration 2011l,m,n,o).

In the context of the study of CMB anisotropies, Galactic emission is a nuisance. Its impact on the accuracy that can be reached on the cosmological parameters has been a concern for a long time and it is still paramount. For this reason, there have been significant efforts to model Galactic emission over the whole sky in the frequency range where the CMB fluctuations can be best observed. Many studies deal with specific Galactic components: synchrotron (Giardino et al. 2002;Platania et al. 2003), free-free (Dickinson et al. 2003;Finkbeiner 2003), thermal dust (Finkbeiner et al. 1999) and, more recently, polarised synchrotron and dust (Waelkens et al. 2009;Fauvet et al. 2011;O'Dea et al. 2012). Some address the question more globally (Bouchet et al. 1996;Bouchet & Gispert 1999;Tegmark et al. 2000), trying to estimate the foreground contamination as a function of frequency and scale. Because of the highly non-stationary statistical properties of the Galactic diffuse emission, this exercise has limitations that can be partly overcome by using a tool that predicts maps of foreground emission for specific regions of the sky.

In that spirit de Oliveira-Costa et al. ( 2008) developed a tool, based on a principal component analysis decomposition of radio survey maps, to produce maps of the total-power Galactic diffuse A96, page 11 of 35 Fig. 7. In all panels, the solid black line is the theoretical angular power spectrum of the lensed CMB, and red filled circles are the average angular power spectrum recovered from 1000 realisations of lensed CMB maps at HEALP N side = 1024. Lensed CMB maps are obtained using the NFFT for the oversampling factor σ = 2 and convolution length

K = 4.
emission in the 10 MHz to 100 GHz range with an angular resolution of 1 • to 5 • . The PSM is similar but also models the submillimetre and far-infrared sky dominated by thermal dust emission, a dominant foreground in much of the Planck frequency range, and the polarised Galactic emission. The PSM also has the capability to produce simulations down to arcminute resolution. Also, contrary to the model of de Oliveira-Costa et al. ( 2008), the PSM is strongly linked to physical models of the Galactic emission processes. It is intended not only to simulate foreground emission for CMB studies but also to serve as a tool to gather relevant knowledge and models of the Galactic ISM.

For each Galactic emission component, the emission across frequencies is modelled using template maps (either observations or model maps that are constructed on the basis of existing data), which are interpolated in frequency according to specific emission laws. Each emission law depends on a set of physical or empirical parameters, which can vary across the sky. Note that Galactic emission is actually implemented as a collection of different models using various input templates, and any of these can be easily called to produce different results. The PSM has been designed in such a way as to facilitate the integration of new models, in particular the ones coming out of the analysis of the Planck data themselves. At the time of writing, much has already been learned about the Galactic emission with the Planck data. This knowledge is not in the PSM yet but it will be after the Planck data become public.

In the following we describe the different models implemented for Galactic emission and in particular the combination of models which we believe, at the current time, is the most plausible and complete given the available data. The default Galactic emission model of the PSM is based on the combination of the work of Miville-Deschênes et al. ( 2008) who modelled the 23 GHz WMAP data by a sum of synchrotron, freefree and spinning dust emission, and the model of Finkbeiner et al. (1999) for thermal dust emission. In polarisation, the default PSM model includes 1) a synchrotron component based on the 7-year 23 GHz WMAP polarisation data (Gold et al. 2011) extrapolated in frequency with the same spectral index as for the total intensity; and 2) a dust component based on a new model presented in Sect. 4.7.2. This subset of models was chosen in order to reproduce the WMAP data, in intensity and in polarisation, while including an unpolarised spinning dust component.

Synchrotron

Galactic synchrotron emission arises from relativistic cosmic rays accelerated by the Galactic magnetic field (e.g., [START_REF] Rybicki | Radiative processes in astrophysics[END_REF]. It is the dominant contaminant of the polarised CMB signal at low frequencies (below about 80 GHz). The intensity of the synchrotron emission depends on the cosmic ray density n e , and on the strength of the magnetic field perpendicular to the line of sight. The spectral signature of the synchrotron emission is determined by the energy distribution of cosmic rays. In the simplest model, in the radio-mm frequency range, for electron density following a power law of index p,( n e (E) ∝ E -p ), the synchrotron frequency dependence is also a power law,

I sync (ν) ∝ ν β s +2 , ( 21
)
or equivalently, in units of antenna (Rayleigh-Jeans) temperature,

T sync (ν) ∝ ν β s , ( 22
)
where I(ν)andT (ν) are related via the equation:

I(ν) = 2kν 2 c 2 T (ν), (23)
and where the spectral index, β s = -(p + 3)/2, is equal to -3for a typical value p = 3. The synchrotron spectral index depends on cosmic ray properties. It varies with the direction on the sky, and possibly, with the frequency of observation (see, e.g., Strong et al. 2007,fora review of propagation and interaction processes of cosmic rays in the Galaxy).

Observations in the radio frequency range from 408 MHz to 10 GHz have shown spatial variations of the spectral index β from -2.8t o-3.2 ( Reich & Reich 1988;Davies et al. 1996;Platania et al. 2003;Bennett et al. 2003b), with a general steepening of the spectrum towards high Galactic latitudes. In a more recent analysis, Miville-Deschênes et al. ( 2008)fi n dal o w e r scatter of the spectral index by taking into account the presence of additional emission in the radio domain from spinning dust grains. The spectral index map obtained in this way is consistent with β s = -3.00±0.06. Thus, as shown in this latter work, the inferred synchrotron spectral index map depends on assumptions made about other components.

Due to energy loss of the cosmic rays, a break in the synchrotron spectrum is expected at frequency higher than ∼20 GHz. A detection of this effect has been reported by the WMAP team in the analysis of the WMAP first year data (Bennett et al. 2003b). As mentioned in the analysis of the 7-year data, this apparent steepening could also be explained on the basis of a contribution of spinning dust to the total radio emission (Gold et al. 2011).

Synchrotron emission is modelled on the basis of the template emission map observed at 408 MHz by Haslam et al. (1982). The original map, reprojected on a HEALP grid in Galactic coordinates at N side = 512, and reprocessed to suppress residual stripes and point sources, is obtained from the LAMBDA website. We further correct the map for an offset monopole of 8.33 K (Rayleigh-Jeans), which is subtracted to match the zero level of the template synchrotron of Giardino et al. (2002). This offset includes the CMB (2.725 K), an isotropic background produced by unresolved extra-galactic radio sources, the background monopole of synchrotron emission, and any zero level error in the data.

This template synchrotron map is extrapolated in frequency using a spectral index map. Three options exist in the PSM for the spectral index. This third model, without spectral curvature, is the default synchrotron model in the PSM. The template 408 MHz map used to model synchrotron emission is displayed in Fig. 8,a n dt h e spectral index map in Fig. 9. In the case of a simulation with resolution better than 1 • , additional small scale fluctuations are added to the synchrotron intensity and spectral index maps, as described below in 4.6.

Free-free

Free-free emission is produced by electron-ion interactions in the ionised phase of the ISM. It is in general fainter than either the synchrotron or the thermal emission from dust, except in active star-forming regions in the Galactic plane. The free-free spectral index depends only slightly on the local value of the electronic temperature T e and it is a slowly varying function of frequency [START_REF] Bennett | [END_REF](Bennett et al. , 2003b;;Dickinson et al. 2003). Between 10 and 100 GHz, the spectral index varies only from 2.12 to 2.20 for 4000 < T e < 10 000 K. In the PSM, the spectral dependance of the free-free emission is based on the model described by Dickinson et al. (2003), assuming a constant electronic temperature (T e = 7000 K is the default value). The default free-free emission law is plotted in Fig. 10.

Three options exists in the PSM to model the spatial variations of the free-free emission.

1. The map of Dickinson et al. (2003)thatusesaHα emission template (a combination of WHAM (Reynolds et al. 1998) and SHASSA (Gaustad et al. 2001) data) corrected for extinction using the E(B -V)mapofSchlegel et al. (1998); 2. The Maximum Entropy Method (MEM) free-free map of WMAP (Bennett et al. 2003b); 3. The free-free map obtained by Miville-Deschênes et al. (2008) that is a combination of the two previous maps. It uses the WMAP MEM decomposition of Bennett et al. (2003b)in regions where the extinction is In more diffuse regions, where Hα can be used more reliably as a proxy for free-free emission, it uses the estimate of Dickinson et al. (2003), unless the WMAP MEM free-free is lower than the Hα estimate.

E(B -V) ≥ 2( o rA V ≥ 6).
These three maps have a resolution of 1 • . The third one, displayed in Fig. 11, is the default free-free model in the PSM.

Molecular lines

Molecular line emission is seen towards dense molecular clouds in our Galaxy and external galaxies. The strongest lines in the useful frequency range are those of 12 CO at frequencies of 115.27 GHz, 230.54 GHz, and 345.80 GHz (for (J = 1 → 0), (J = 2 → 1) and (J = 3 → 2) respectively). Other emission lines of interest, although fainter than those of 12 CO, include: (2001). The original map is not full sky, but contains most of the regions of strong emission (about 45% of the sky is covered, i.e., essentially all of the sky at Galactic latitudes |b| < 30 • ).

As surveys of other 12 CO lines are missing, we choose to model the (J = 2 → 1) and (J = 3 → 2) lines by scaling the (J = 1 → 0) map using the following line ratios in units of K RJ km s -1 (Boulanger 2010):

(J = 2 → 1)/(J = 1 → 0) = 0.65; (J = 3 → 2)/(J = 2 → 1) = 0.50.
The resulting relative integrated line intensities are displayed in Fig. 10. Emission from other molecular lines ( 12 CO lines of higher order or involving isotopes, and molecules other than CO) is neglected in the current model.

Thermal dust emission

The thermal emission from heated dust grains is the dominant Galactic signal at frequencies above about 80 GHz. Contributions from a large range of grain sizes and compositions are required to explain the infrared spectrum of dust emission from 3 µmto1mm (Désert et al. 1990;Li & Draine 2001;Compiègne et al. 2011). At long wavelengths of interest for CMB observations (above 100 microns), the emission from big dust grains at equilibrium with the interstellar radiation field should be the dominant source of the observed radiation.

The thermal emission from interstellar dust is estimated on the basis of Model 7 or 8 (up to the user, model 7 being the default) of Finkbeiner et al. (1999). This model, fitted to the FIRAS data (7 • resolution), is based on the hypothesis that each line of sight can be represented by the sum of the emission of two dust populations, one cold and one warm. To limit the number of parameters, the model rests on three important assumptions:

1. that the two dust populations are well mixed spatially and therefore both heated by the same radiation field; 2. that the optical properties of the two populations are constant (i.e., fixed β and opacity ǫ); 3. that the ratio of their abundance (in mass) is constant.

These assumptions imply that the ratio of absorbed (and emitted) power of the two populations is constant, a limitation that will be alleviated in the future using the 5 ′ all-sky dust temperature derived from combining the IRAS and Planck data, such as computed in Planck Collaboration (2011l). At present, however, given these assumptions, Finkbeiner et al. (1999)sho wed that the spectral and spatial variations observed are only determined by the local strength of the interstellar radiation field and the total dust column density. One can estimate the dust temperature of the two populations using only the ratio of the emission at two wavelengths. Following Finkbeiner et al. (1999), we use their map of the ratio of 100 µm to 240 µm emission based on DIRBE data, which has been processed to have a constant signalto-noise ratio over the sky, resulting in a lower angular resolution (several degrees) at high Galactic latitude than in the plane (40 ′ ). This map is shown in Fig. 12. From the temperature of each dust population at a given position on the sky, the assumed ratio of emitted power, and the observed 100 µm brightness, one can estimate the dust column density N i of the two populations at each sky position. The extrapolation at any frequency is then given by the sum of two modified blackbodies,

I ν = 2 i=1 N i ǫ i ν β i B ν (T i ), ( 24
)
where B ν (T i ) is the Planck function at temperature T i , N i is the column density of species i,andǫ i ν β i accounts for the normalisation and frequency dependence of the emissivity. The 100 µm map used is the one of Schlegel et al. (1998) or a modified version of this map, displayed in Fig. 13, for which residual point sources, including a catalogue of ultra-compact H  regions (see Sect. 5.4), were removed.

The main uncertainty in this representation of thermal dust emission comes from the model itself, rather than from noise in the observations used. We recall that the model is a fit to the very large scale dust emission (7 • ). Variations of dust properties are observed everywhere in the ISM (Planck Collaboration 2011n,o). In addition, some models of interstellar dust predict significant variations of the modified blackbody spectrum for amorphous silicate grains (Boudet et al. 2005). The validity of the hypothesis of a constant ratio of absorbed power is therefore to be questioned. Recent work, modelling dust emission as a two-level system in a disordered medium (Meny et al. 2007;Paradis et al. 2011), predicts significant curvature in the dust spectral index at frequencies below 1 THz and is a counterexample to the assumption of a simple modified blackbody.

Under the hypothesis where the model used here is correct, the main uncertainty would come from the low angular resolution of the map of the ratio of 100 µm to 240 µm emission which does not trace the small scale variations of the dust temperature.

Spinning dust

In the 10-100 GHz range several observations have pointed to excess emission, the so-called anomalous microwave emission (Kogut et al. 1996b;Leitch et al. 1997;de Oliveira-Costa et al. 1999;Watson et al. 2005). At frequencies above ∼20 GHz this diffuse emission has a spectrum falling with frequency, similar to synchrotron, but has a spatial structure closer to the distribution of thermal dust. One model proposed by Draine & Lazarian (1998) attributes this emission to small spinning dust particles. A detailed analysis of the WMAP data by Davies et al. (2006) revealed the presence in selected areas of a component spatially correlated with thermal dust and with a spectral signature compatible with the model of Draine & Lazarian (1998). The spinning dust hypothesis is also compatible with the analysis of the WMAP 23 GHz polarisation data (Miville-Deschênes et al. 2008), with the ARCADE 2 measurements between 3 and 10 GHz (Kogut et al. 2011) and with the fact that the microwave anomalous emission is well correlated with the 12 µm emission divided by the radiation field strength (Ysard et al. 2010a), as predicted by models (Ysard & Verstraete 2010b). In addition, a recent joint analysis of data from the Planck mission, from WMAP, and from the COSMOSOMAS experiment shows unambiguous detection of a component in some specific regions, the spectrum of which peaks around 20 GHz and is compatible with spinning dust (Planck Collaboration 2011m).

Similarly to the free-free component, the spinning dust emission is modelled using a spinning dust template map extrapolated in frequency using a single emission law. Three template maps are available. The emission law can be parameterised following the model of Draine & Lazarian (1998); the default adopted in the PSM, plot-tedinFig.10, is a mixture of 96% warm neutral medium and 4% reflection nebula. This specific parametrisation of the Draine & Lazarian (1998) model has not been selected on any physical ground but because it fits the average spinning dust spectrum found in Miville-Deschênes et al. ( 2008).

Adding small-scale fluctuations

Small-scale fluctuations are automatically added to Galactic emission simulations when the resolution of the simulated map is greater than the resolution of the original template. The method used is described in Miville-Deschênes et al. (2007). A Gaussian random field G ss having a power spectrum defined as

C ℓ = ℓ γ exp -ℓ 2 σ 2 sim -exp -ℓ 2 σ 2 tem , (25)
is generated on the HEALP sphere at the proper N side .H e r e σ sim and σ tem are the Gaussian widths (in radians) of the simulation and of the template to which small-scale fluctuations are added. The zero mean Gaussian field G ss is then normalised and multiplied by the template map exponentiated to a power β,i n order to generate the proper amount of small-scale fluctuations as a function of local intensity. This is motivated by the fact that the power of the Galactic emission fluctuations at a given scale generally scales with the local average intensity at a given power between 2 and 3 (Miville-Deschênes et al. 2007). The resulting template map with small-scale fluctuation added is then

I ′ tem = I tem + α G ss I β tem , (26)
where α and β are estimated for each template in order to make sure that the power spectrum of the small-scale structure added is continuous with the large-scale part of the template. This approach for generating small scales is similar to what has been done by Giardino et al. (2002), although in our case the weighting by the low-resolution map is slightly different. The addition of small scales is illustrated for free-free in Fig. 14 and for spinning dust in Fig. 15.

ISM Polarisation

Synchrotron polarisation

If the electron density follows a power law of index p, the synchrotron polarisation fraction is

f s = 3(p + 1)/(3p + 7). ( 27
)
For p = 3, f s = 0.75, a polarisation fraction that varies slowly for small variations of p. Consequently, the intrinsic synchrotron Instead of relying directly on a model of the Galactic magnetic field, our model of polarised synchrotron emission is built on the basis that synchrotron is the main low-frequency polarised emitter. Hence we use the WMAP 23 GHz channel polarisation maps Q 23 and U 23 , extrapolated in frequency using the same synchrotron spectral index as that of intensity. Other options may be implemented in the future, using either different models of the Galactic magnetic field or interfacing the PSM with dedicated software for simulating Galactic synchrotron emission, such as the publicly available  code (Waelkens et al. 2009).

To extend the resolution of Q 23 and U 23 , small scales are added in a similar way as for the total intensity templates (see Eq. ( 26)), but because Q and U are not strictly positive, the scaling of the Gaussian random field cannot be done with the template itself. Instead the small scale field is scaled with the unpolarised intensity field. For

Q 23 that is Q ′ 23 = Q 23 + α G ss I β 23 .
Figure 17 displays synchrotron Q and U maps, as well as the polarisation fraction and angle as modelled here. The power spectrum of the polarised sky emission model at WMAP frequencies, using the P06 mask, is displayed in Fig. 16

Dust polarisation

Polarisation of starlight by dust grains indicates partial alignment of elongated grains with the Galactic magnetic field (see Lazarian 2007, for a review of possible alignment mechanisms). Partial alignment of grains also results in polarisation of their emission in the far-infrared to millimetre frequency range.

Very few observational data are currently available to model polarised dust emission, although things will drastically change soon with the observations of the Planck mission (Tauber et al. 2010). So far, dust polarisation measurements have mostly concentrated on specific regions of emission, with the exception of the Archeops balloon-borne experiment (Benoît et al. 2004), which has mapped the emission at 353 GHz over 17 percent of the sky, showing a polarisation fraction around 4-5% in the Galactic plane, and up to 10-20% in some clouds. Recently, polarisation of dust emission at 100 and 150 GHz has been observed by BICEP over a smaller region of the Galactic plane (Matsumura et al. 2010) and by QUaD (Culverhouse et al. 2010). A large-scale estimate of the polarised dust power spectrum has been obtained with WMAP 7-year data by Gold et al. (2011). The measured dust polarisation fraction is in rough agreement with what is expected from the polarisation of starlight (Fosalba et al. 2002). Draine & Fraisse (2009) have shown that for particular mixtures of dust grains, the intrinsic polarisation of the dust emission could vary significantly with frequency in the 100-800 GHz range. As for synchrotron, geometrical depolarisation caused by integration along the line of sight also lowers the observed polarisation fraction; the observed polarisation fraction of dust depends on its three-dimensional distribution, and on the geometry of the Galactic magnetic field.

By making use of presently available data, we model polarised thermal dust emission by extrapolating dust intensity, I ν (p) (see Sect. 4.4), to polarisation intensity, assuming intrinsic polarisation fractions, f d , for the two dust populations and a model for the dust geometric depolarisation, g d , and dust polarisation angle, γ d . The Stokes Q and U parameters for dust are

Q ν (p) = f d g d (p) I ν (p)c o s ( 2 γ d (p)), (28) U ν (p) = f d g d (p) I ν (p)s i n ( 2 γ d (p)). ( 29
)
The value of f d is set to 0.15 for the two dust populations, consistent with maximum values observed by Archeops (Benoît et al. 2004) and in good agreement with the WMAP 94 GHz measurement. It should be noted that f d does not correspond to the average observed polarisation fraction P ν /I ν ,w h e r eP ν = U 2 ν + Q 2 ν .T h i si s f d g d (p), which has a mean of about 5% by default in the current model. Setting the same value of f d for the two dust populations implies that the polarisation fraction in the PSM does not change with frequency. The model has been designed to allow for different values of f d for the two dust populations, but in the absence of observational evidence, both dust populations have been given the same value for now.

Because of the lack of observations, modelling polarised dust emission comes down to estimating maps of the geometric depolarisation, g d , and of the polarisation angle, γ d .Oneapproach is to use the large-scale structure of the Galactic magnetic field with a prescription for the turbulent depolarisation along the line-of-sight as O'Dea et al. ( 2012) did. This produces very smooth maps of g d and γ d as it does not reproduce the structure projected on the sky due to the Galactic magnetic field associated with the ISM turbulent dynamics. Such structure can be included by using a code like H (Waelkens et al. 2009)t h a t simulates the turbulent part of the magnetic field in three dimensions. This was done by Fauvet et al. (2011) and is implemented as an option in the PSM but, even though the geometric depolarisation and polarisation angle maps produced that way have sensible morphology and statistical properties, their detailed structure does not match the observed sky. Another option would be to use polarisation angle derived from optical starlight measurements (Heiles 2000;Berdyugin et al. 2001Berdyugin et al. , 2004;;Berdyugin & Teerikorpi 2002). However, there are several drawbacks to using optical measurements. These data are under sampled on the sky, with very few measurements at high Galactic latitudes, which makes it difficult to use them to generate full-sky polarisation angles. Also, polarisation in absorption in the optical does trace only the dust located in front of the star on which the absorption is measured, unlike the dust polarisation in the sub-millimetre that samples the whole line of sight.

Hence, in order to simulate dust polarisation maps that stay as close as possible to the data, and to be coherent with the polarised synchrotron model, we make use of the 23 GHz WMAP polarisation data to estimate the geometric depolarisation and polarisation angle maps for dust. Here we explicitely make the assumption that, on average, synchrotron and thermal dust polarised emission are influenced by the same magnetic field and that their maps of g and γ should be correlated. Therefore, in our model g d and γ d are based on equivalent maps estimated for the synchrotron emission using the WMAP 23 GHz and the Haslam 408 MHz data, the synchrotron polarisation angle

γ s = 1 2 tan -1 (-U 23 , Q 23 ) , (30)
and geometric depolarisation factor

g s = U 2 23 + Q 2 23 1/2 f s I 408 (23/0.408) -3 , (31)
where all templates are smoothed to 3 • . In order to imprint a specific signature of dust on g s and γ s , we have modelled what would be the large-scale structure of those two quantities for dust and synchrotron. In order to do that, we have built maps of I, Q and U for dust and synchrotron given a model of the Galactic magnetic field and a simplistic model of the distribution of dust and electrons in the Milky Way. The maps of g and γ are then obtained by inverting Eqs. ( 28)and( 29). Assuming only one dust population of volume density n d , the total intensity of dust is

I d (ν) = z n d ǫ d ν β d B ν (T d )dz, ( 32
)
where the integral is along the line of sight z.F orQ and U the integrals include terms related to the angle between the magnetic field direction and the line of sight:

Q d (ν) = f d z n d ǫ d ν β d B ν (T d )cos2φ sin α dz, (33)
and

U d (ν) = f d z n d ǫ d ν β d B ν (T d )sin2φ sin α dz, (34)
where

cos 2φ = B 2 x -B 2 y B 2 ⊥ , ( 35
) sin 2φ = -2B x B y B 2 ⊥ , ( 36
) sin α = 1 -B 2 z /B 2 , (37)
and

B ⊥ = B 2 x + B 2 y , (38)
with the x-y plane corresponding to the plane of the sky. For synchrotron, the unpolarised intensity I is given by

I s (ν) = z ǫ s ν β s n e B (1+p)/2 ⊥ dz, (39)
A96, page 18 of 35 where β s is the synchrotron spectral index defined in Eq. ( 22). As in the dust case, the Stokes parameters of polarised synchrotron emission are

Q s (ν) = f s z ǫ s ν β s n e B (1+p)/2 ⊥ cos 2φ sin α dz, (40)
and

U s (ν) = f s z ǫ s ν β s n e B (1+p)/2 ⊥ sin 2φ sin α dz. ( 41
)
The polarisation fraction f s is related to the cosmic ray energy distribution according to Eq. ( 27). 2008), which has a regular component following a bi-symmetrical spiral with a pitch angle of -8.5 • and a turbulent component with an amplitude of 1.7 µG( i . e . , 0.57 times the local value of the magnetic field B 0 = 3 µG). The turbulent part of the field used here cannot reproduce the multiscale properties of the interstellar emission on the plane of the sky; it is included to take into account the fact that longer lines of sight in the Galactic plane have a stronger depolarisation effect, due to a larger number of turbulent fluctuations of the magnetic field than at high Galactic latitude. The turbulent part of each line of sight is simulated independently of its neighbours, so this method produces no turbulent structures on the plane of the sky.

The maps of g and γ obtained using this model reflect polarisation effects on large scales (>20 • ). These maps, shown in Fig. 18, are similar to the ones used in the model of O'Dea et al. (2012). As mentioned previously, local variations of the dust and electron densities and of the magnetic field properties caused by the complex physics of the ISM will produce specific structures in the real maps that cannot be reproduced in this simple model. The objective here is to capture large-scale differences between the synchrotron and dust polarised emission due to the fact that 1) dust grains and energetic electrons have different scale-heights and lengths in the Galaxy; and that 2) their respective polarised emissions do not depend in the same way on the magnetic field (unlike the dust, the synchrotron emissivity depends only on B ⊥ -seeEq.(39)). Figure 18 shows indeed that the geometrical depolarisation and polarisation angle maps of these two processes are slightly different.

To produce the dust geometrical depolarisation (g d ) and polarisation angle (γ d ) maps we combine this large-scale model (that provides the structure at scales larger than 20 • ) with the synchrotron g s and γ s estimated with the 23 GHz and 408 MHz templates (see Eqs. ( 30)and( 31)) to add structures between 3 • and 20 • . The dust geometrical depolarisation and polarisation angle maps used in the model are then defined as 43)and( 44)).

g d = g ′ d + (g s -g ′ s ), (43)
and

γ d = γ ′ d + (γ s -γ ′ s ), (44)
where g ′ i and γ ′ i are the maps built with the large-scale model and shown in Fig. 18. The addition of fluctuations at scales smaller than 3 • is performed in a similar way as for the synchrotron polarisation: they are added in Q and U directly with a scaling following the intensity I.

The dust polarisation maps estimated with our model at 200 GHz (in the frequency range of interest for high resolution CMB observations) are shown in Fig. 19. As expected, the polarisation fraction and polarisation angle maps are very similar to the ones for the synchrotron shown in Fig. 17 hand, the maps of Q and U are significantly different from those for the synchrotron, because of the distinct spatial structure of the thermal dust and synchrotron unpolarised emission. We expect the modelling of polarised dust to be significantly improved in the future, with the use of the Planck polarised data between 100 and 353 GHz.

Polarisation of other ISM emission

Free-free emission is intrinsically unpolarised. At second order, one expects some polarisation towards ionised regions by a mechanism similar to that giving rise to CMB if incident light, scattered by free electrons in the ionised gas has a quadrupole moment perpendicular to the line of sight, then scattered light is preferably polarised perpendicular to the direction of incidence of the strongest incoming radiation, by reason of the polarisation dependence of the Thomson cross section. This is neglected in the current version of the model. The Thomson scattered free-free light is expected to be small (a few percent or likely less, Page et al. 2007;Gold et al. 2011;Macellari et al. 2011), with a maximum of about 10% on the edge of bright H  regions (e.g., [START_REF] Rybicki | Radiative processes in astrophysics[END_REF].

Spinning dust emission is observed to be at most weakly polarised, of order a few percent (Battistelli et al. 2006;López-Caraballo et al. 2011;Dickinson et al. 2011). On the theoretical side, Lazarian & Draine (2000) predict very low (but possibly measurable) levels of polarisation for spinning dust emission. Even if anomalous dust emission is due to magnetic dust (Draine & Lazarian 1999), the polarisation is expected to be low (∼0.5%) if the iron inclusions are random (Lazarian & Draine 2000). In the current software implementation, given the absence of observational evidence or theoretical predictions of significant spinning dust polarisation, its emission is assumed unpolarised.

Limitations of the models of galactic emission

There is no theoretical statistical model for ISM emission. Models for the emission of all the galactic components are phenomenological descriptions based on observed emission. As such, the models are at least as uncertain as sky observations.

In particular, free-free emission, which does not dominate at any frequency, is poorly measured so far. Maps of free-free emission used in the PSM are derived from either or both of an Hα emission composite map that combines two different surveys and a map based on the WMAP MEM decomposition. The former, using Hα as a proxy for free-free emission, is subject to various uncertainties (dust absorption, scattering of Hα by dust grains, variations in electron temperature). The latter strongly depends on the model assumed in the WMAP analysis, and does not seem to match well the radio recombination line measurements in Alves et al. (2010) (albeit over a limited sky region). Hence none of the two basic free-free templates used in the PSM is very reliable, and the modelled PSM free-free emission should be considered at best as indicative.

Thermal dust emission is modelled following the study of Finkbeiner et al. (1999) based on an analysis of the FIRAS data. These authors showed that the data can be decomposed by the sum of the emission from two dust populations. First, it is important to note that, even though it is thought that interstellar dust is made of carbonaceous and silicate grains, the dust spectral indices and equilibrium temperatures used in this model are not compatible with the current knowledge of thermal emission from big grains (Li & Draine 2001;Compiègne et al. 2011). Combined with the limited signal-to-noise ratio of the FIRAS low frequency data used to develop this model, the spectral shape of the dust emission at frequencies lower than 300 GHz is not expected to be very accurate. Second, the template of the dust temperature used in this model is based on a version of the 100/240 µm ratio map from the DIRBE data, smoothed at high Galactic latitude to increase the signal-to-noise ratio. Thus the effective resolution of the dust temperature template is varying across the sky, from 30 arcmin in the plane to only several degrees at high Galactic latitude. That has to be kept in mind when investigating the variations of the dust emission properties at high Galactic latitude. Similarly, the template and emission laws of spinning dust emission are representative, but not exact. The models for thermal and spinning dust emissions will be adjusted in future releases of the PSM, after the analysis of the Planck observations has been completed.

As discussed by Sun et al. (2008), the synchrotron emission is absorbed at low frequency by foreground thermal electrons that produce free-free emission. This absorption is less important at 23 GHz than at 408 MHz, which produces an artificial flattening of the synchrotron spectral index. Therefore, Sun et al. (2008) showed that it is likely that in the inner portion of the Galactic plane, where the free-free absorption can be significant, the 408 MHz map underestimates the synchrotron emission if extrapolated at higher frequencies with a constant spectral index. At this point, given the significant uncertainty on the freefree emission in the Galactic plane and on the filling factor of the diffuse ionised gas (the opacity is proportional to the square of the electron density), it is difficult to quantify this effect on the 408 MHz map. In addition, the fiducial synchrotron model of the PSM is not based on an extrapolation of the 408 MHz map; this synchrotron template has been derived from the WMAP 23 GHz polarisation data, corrected for a model of the synchrotron depolarisation Miville-Deschênes et al. (2008). In this case the uncertainty on the synchrotron map is coming from the model of the depolarisation itself.

Emission from compact sources

Numerous galactic and extragalactic compact sources emit radiation in the frequency range modelled in the PSM. Extragalactic emission arises from a large number of radio and far-infrared galaxies, as well as clusters of galaxies. Compact regions in our own galaxy include molecular clouds, supernovae remnants, compact H  regions, as well as galactic radio sources.

The thermal and kinetic SZ effects, due to the inverse Compton scattering of CMB photons off free electrons in ionised media, are of special interest for cosmology. These effects occur, in particular, towards clusters of galaxies, which are known to contain hot (few keV) electron gas. SZ emission is modelled here both on the basis of known existing clusters, observed in the optical and/or in X-rays, and on the basis of cluster number counts predicted in the framework of a particular cosmological model.

Distant radio and far-infrared sources are not resolved by CMB instruments, and are represented in our model as a population of point sources, each of which is described by a number of parameters (type of the source, position on the sky, flux and polarisation as a function of frequency). Although unresolved, the brightest objects are detected individually in observations, and represented in the form of a catalogue (rather than a map).

In addition to those brightest sources, a diffuse source background arises from the integrated emission of a large number of objects that are too faint to be detected individually, but which contribute to sky background inhomogeneities. Background sources are represented using maps of resulting brightness fluctuations.

Finally, our Galaxy itself contains a number of compact regions of emission that can be modelled as point sources. Some of those are modelled as part of the catalogues of radio and far-infrared sources (without distinction between galactic and extragalactic objects). A specific population, identified as ultra-compact galactic H  regions, is singled out and treated separately.

Sunyaev-Zeldovich effects

Hot electron gas, as found in the largest scale structures in the universe, interacts with incoming photons by inverse Compton scattering. In particular, background photons are scattered by hot gas in clusters of galaxies, giving rise to secondary anisotropies of the CMB by shifting the energy of the interacting photons. The general properties of the effect have been described first by Sunyaev & Zeldovich (1972)( s e ea l s oBirkinshaw 1999; Carlstrom et al. 2002, for recent comprehensive reviews).

Classically, a distinction is made between two SZ effects. The thermal SZ effect (tSZ) corresponds to the interaction of the CMB with a hot, thermalised electron gas, while kinetic SZ (kSZ) corresponds to CMB interaction with electrons having a net ensemble peculiar velocity along the line of sight.

In the thermal case, the electron population is characterised by a Fermi distribution at a temperature T e much larger than the CMB temperature (in typical clusters of galaxies, T e is of order few keV). For non-relativistic electrons, the change in sky brightness is:

δI ν = y f (ν)B ν (T CMB ), ( 45
)
where B ν (T CMB ) is the CMB blackbody spectrum, f (ν) is a universal function of frequency that does not depend on the physical parameters of the electron population (and hence on the galaxy cluster they are associated with), given by

f (ν) = xe x e x -1
x (e x +1)

e x -1 -4 , (46)
where x = hν/kT CMB . The Compton parameter y is proportional to the integral along the line of sight of the electronic density n e multiplied by the temperature of the electron gas:

y = kT e m e c 2 n e σ T dl. ( 47
)
For a standard cluster detectable at millimetre wavelengths from its thermal SZ effect with current instruments, T e ≃ 5k e V ,t h e optical depth of the gas, τ = n e σ T dl,isoforder10 -2 , y ≃ 10 -4 and the Compton parameter integrated over the cluster angular extent, Y SZ ,isoforder10 -4 arcmin 2 for a typical cluster angular size 1 ′ . In the kinetic case, the interaction of CMB photons with a population of moving electrons results in a shift of the photon distribution seen by an observer on Earth:

δI ν = -β r τ ∂B ν (T ) ∂T T =T CMB , (48)
where β r = v r /c is the dimensionless cluster velocity along the line of sight. Hence, the kinetic SZ effect has the same frequency dependence as CMB primary anisotropies themselves. For typical radial velocities of 300 km s -1 (β r ≃ 10 -3 ), and for τ ≃ 10 -2 the net kinetic SZ effect is ∆T/T ≃ 10 -5 . At present, only limited SZ observations are available. Those observations aim both at studying previously known clusters (Basu et al. 2010;Zwart et al. 2011;Hincks et al. 2010)a n d at discovering new objects blindly (Vanderlinde et al. 2010;Williamson et al. 2011;Planck Collaboration 2011c;Marriage et al. 2011a;Reichardt et al. 2013). They have not only confirmed the existence of the effect but they have also made possible a variety of studies which allow predictions of the general properties of the SZ sky, such as the shape and normalisation of its angular power spectrum (Shaw et al. 2009), and the connection between SZ observables and cluster masses or the X-ray luminosity from direct free-free emission of the hot ionised gas (Melin et al. 2011;Planck Collaboration 2011d,e).

Several options are available to simulate SZ emission with the present sky model. A prediction model generates thermal SZ effect as expected from catalogues of known clusters, observed in X-rays with ROSAT or in the visible with SDSS (Sect. 5. 1.1). A simulation of both thermal and kinetic SZ emission is possible on the basis of cluster number counts from cluster mass functions (Sect. 5.1.2), or using N-body and hydrodynamical simulations of large-scale structures (Sect. 5.1.3), or making use of hydrodynamical simulations at low redshift (z < 0.25), complemented by cluster counts at higher redshift (Sect. 5.1.4).

SZ emission prediction

For any given cluster of known mass and redshift, scaling relations exist that allow one to predict the level of thermal SZ emission originating from that particular cluster. The prediction of the thermal SZ emission of the sky in our model is based on the observations of 1743 galaxy clusters with the ROSAT X-ray satellite (NORAS Böhringer et al. 2000;REFLEX Böhringer et al. 2004, serendipitous surveys), and of 13 823 optically selected clusters extracted from the SDSS galaxy survey. 215 clusters are common to these two independent catalogues.

The ROSAT catalogue used in the model is the Meta-Catalogue of X-ray detected Clusters of galaxies (MCXC) by Piffaretti et al. (2011). It covers the full sky, although with uneven depth. These clusters are located at redshifts ranging from 0.003 to 1.26, the median redshift of the sample being 0.14. For each cluster observed in X-rays, the predicted SZ flux Y 500 is given by the combination of the L X -M relation of P r a t te ta l . ( 2009)andtheM-Y relation of Arnaud et al. (2010). The cluster model adapted here and its normalisation on the basis of cluster observations is detailed in Planck Collaboration (2011d). We simply cite here the two key scaling laws used in the model. The first one connects the X-ray luminosity and the mass:

E(z) -7/3 L 500 10 44 erg s -1 = C LM M 500 3 × 10 14 M ⊙ α LM , (49)
where E(z) is the Hubble parameter normalised to its present value, L 500 is the X-ray luminosity emitted, in the [0.1-2.4] keV band, from within the radius R 500 at which the mean mass density is 500 times the critical density of the universe at the cluster redshift, M 500 is the cluster mass within R 500 , and the numerical values of the constants are log (C LM ) = 0.274 and α LM = 1.64.

The second one connects the SZ flux and the mass:

Y 500 = 1.383 × 10 -3 I(1) M 500 3 × 10 14 M ⊙ 1 α MY X × E(z) 2/3 D A (z) 500 Mpc -2 arcmin 2 , ( 50
)
where Y 500 is the SZ flux in R 500 , I(1) = 0.6145 is a numerical factor arising from the volume integral of the pressure profile, D A (z) is the angular distance and α MY X = 0.561. MaxBCG clusters (Koester et al. 2007) are detected mostly in the northern Galactic hemisphere (see Fig. 20). These clusters are located at redshifts ranging from 0.1 to 0.3, the median redshift of the sample being 0.23. We compute the SZ flux Y 500 using the Y 500 -N 200 relation derived by the Planck Collaboration (2011f), where N 200 is the cluster richness, Y 20 = 7.4 × 10 -5 arcmin 2 and α 20 = 2.03.

Y 500 = Y 20 N 200 20 α 20 E 2/3 (z) D A (z) 500 Mpc -2 arcmin 2 , (51)
The SZ prediction implemented in the current version of our model does not include kinetic or polarised SZ effects.

SZ emission from cluster mass functions

An approach for simulating the SZ effect from galaxy clusters consists in assuming a cluster mass function N(M, z). Following the simulation method described in Delabrouille et al. (2002a), we start from a set of cosmological parameters (Hubble parameter h; density parameters for matter Ω m , vacuum energy Ω Λ , and baryons Ω b ; amplitude σ 8 of fluctuations on 8 h -1 Mpc scales). This set of cosmological parameters matches what is used for simulating CMB anisotropies (and hence is defined uniquely for both). The cluster distribution in the mass-redshift plane (M, z) is then drawn from a Poisson law whose mean is given by a mass function: Press-Schechter (Press & Schechter 1974), Sheth-Tormen (Sheth & Tormen 1999), Jenkins (Jenkins et al. 2001), Evrard (Evrard et al. 2002), or Tinker (Tinker et al. 2008). Cluster Galactic coordinates (l, b) are then drawn at random on the sphere, with a uniform probability distribution function (which here neglects correlations in the large-scale structures traced by galaxy clusters). Baryons are pasted on clusters assuming an electron density profile (either a Navarro et al. 1996 profile, modified according to Nagai et al. 2007 andArnaud et al. 2010 or a β-model) and using a mass-temperature relation (Pierpaoli et al. 2003), the normalisation of which is a free parameter, or a Y SZ -M relation from Chandra or XMM observations (Nagai et al. 2007;Arnaud et al. 2010).

The three-dimensional velocities of the clusters are drawn from a Gaussian distribution with zero mean and standard deviation given by the power spectrum of density fluctuations (Peebles 1993). This, again, neglects correlations between cluster motions, i.e., largest-scale bulk flows. The kinetic SZ map is obtained using this velocity field and the electron density of the clusters.

This semi-analytic method of simulating SZ clusters requires relatively small amount of CPU time. It is hence fast enough to allow for a new generation of the cluster catalogue with each run of the code, and thus makes it possible to produce many different sky maps with different input parameters as, for example, a varying value of σ 8 . The gas properties and scaling relations can be easily tuned to the observed ones (and match those measured recently from WMAP and Planck data). It is possible to include the contribution from clusters in an arbitrarily large mass A96, page 22 of 35 range, including groups of galaxies. Simulations can be repeated many times with different seeds, which permits us to characterise the selection function of a data processing chain (Melin et al. 2005) and to evaluate errors in parameter estimation by Monte-Carlo methods. Current drawbacks are the lack of proper halo-halo correlations and of large-scale velocity flows, and the absence of halo substructure for the resolved clusters. Note also that the approach based on a cluster mass function does not include the faint SZ emission of the filamentary structures of the cosmic web.

Comparison of cluster counts

Figure 21 shows the cumulative counts for clusters of mass M 500 > 10 14 M ⊙ , for a subset of the mass functions available to simulate the SZ effect and for different values of σ 8 . Counts are quite different between the various options. The dependence of counts on σ 8 is important and its precise value still uncertain, with a typical allowed range extending from about 0.81 as obtained from primary CMB constraints (Komatsu et al. 2011)to about 0.87 as obtained from weak lensing (Massey et al. 2007). Moreover, the dependence on the mass function is weaker but not negligible: for example, for the same set of WMAP 7-year cosmological parameters, the Jenkins and Sheth-Tormen mass functions predict about 50% more cluster above M 500 > 10 14 M ⊙ than the Tinker mass function. The sky model software can be used to produce maps with different mass functions and values of σ 8 .

Constrained SZ simulations

The simulation of cluster catalogues on the basis of a mass function can be refined to constrain the catalogue to match available observational data. We then replace a subset of clusters of the randomly drawn catalogue by the clusters of the ROSAT all-sky survey and SDSS. The subset of replaced clusters in our catalogue is then chosen to match best the redshift-mass distribution of the MCXC and MaxBCG catalogues. The resulting maps can be computed quickly and are fully consistent with the X-ray and optical characteristics of the known clusters.

SZ emission from large-scale structure simulations

A complete and accurate simulation of the SZ effect requires a simulation of the distribution of hot electron gas in the universe. As full scale simulations, including all the physics of the baryons in any given cosmological scenario, are not available at present, we adopt a simplified approach which relies both on N-body+hydrodynamical simulations of the distribution of baryons for redshifts z < 0.025 (the local universe), and on pure N-body simulations of dark matter structures in a Hubble volume.

The local simulation of baryonic matter is obtained from the N-body+hydrodynamical simulation of Dolag et al. (2005). The distribution is constrained by the observed IRAS 1.2-Jy all-sky redshift survey, so that the positions and masses of the closest galaxy clusters are constrained to match quite closely those of real clusters. The distribution and peculiar motion of the baryonic gas as obtained in the simulation permits us to simulate directly the corresponding thermal and kinetic SZ effect.

The distribution of matter at higher redshift is obtained starting from a pure N-body simulation of a Hubble-volume (a cube of 1 Gpc size). The SZ map for this latter simulation is obtained by putting, in the position of each dark matter concentration in the Hubble-volume simulation, a suitable cluster extracted from a template obtained using a different N-body+hydrodynamical simulation of a smaller volume (Schäfer et al. 2006).

Merging the two simulations, thermal and kinetic SZ maps are obtained by co-adding the contributions of the high redshift and low redshift structures. It is worth noting that while the higher redshift map contains only clusters, the hydrodynamical map used to describe the SZ emission from the local universe features the SZ from all the filamentary structure of the cosmic web. One pre-generated SZ simulation obtained in this way is currently available in the model.

The N-body+hydrodynamical SZ simulation has several advantages. First of all, it is the only model presently implemented in the code that properly takes into account halo-halo correlations, models halo substructure and adiabatic gas physics, and provides cluster peculiar velocities in agreement with what expected from the actual density contrast and its evolution with time. The simulation also models properly the filamentary structure of the cosmic web for the low redshift universe. It also matches the observed low redshift matter distribution. The main drawbacks is the lack of flexibility: one single simulation is available, generated for fixed cosmological parameters (h = 0.7, Ω m = 0.3, Ω Λ = 0.7, Ω b = 0.04, σ 8 = 0.9, n S = 1). In addition, the SZ map is generated on the basis of the extraction of a number of halos from the output of the LSS simulation. Hence, the mass range of the clusters included in the map is limited by the cluster templates used: many low-mass systems, and even a few high-mass ones, may be missing in the final map. Finally, the cluster gas properties and scaling relations strongly depend on the physics included in the N-body+hydrodynamical code used to generate the cluster templates. If some of these properties are not in complete agreement with X-ray cluster observations, this approach does not easily allow us to modify them to correct for this effect.

In order to cover redshifts out to the anticipated limit for Planck, several light-cone outputs were combined in the Hubble volume simulation: first, a sphere covering the full solid angle of 4π was used for redshift radii up to z = 0.58. For redshifts exceeding z = 0.58, octant data sets, spanning a solid angle of π/2, were replicated by rotation in order to cover the full sphere. Note that this generates fake replications of the farthest simulated clusters, the effect of which can be seen by careful inspection of the maps displayed in the left column of Fig. 22.

SZ from LSS simulations and high z cluster counts

In the previous approach, full hydrodynamical simulations at low redshift are complemented with a single high-redshift simulation of large-scale structures. An alternate solution consists in merging a low redshift full hydrodynamical simulation, containing the constrained local SZ map, with a high redshift model based on cluster number counts.

In this model, the low redshift SZ signal is simulated in concentric layers around the observer using snapshots from full hydrodynamical simulations up to z ≃ 0.25. Seven layers, each with a comoving thickness of 100 h -1 Mpc, are used to generate the required volume. The thermal and kinetic SZ signals are then integrated using the method in da Silva et al. (2000,2001)and the HEALP sky pixelisation scheme. For the innermost layer we used the local constrained simulation in Dolag et al. (2005), i.e., the baryon distribution and the resulting local constrained S Zm a pu s e di nS e c t .5.1.3. For the remaining layers up to z = 0.25 we use the ΛCDM simulation in De Boni et al. (2010). Both hydrodynamical simulations feature full baryonic gas dynamics and state-of-the-art modelling for gas cooling, heating by UV, star formation and feedback processes.

We then draw at random additional clusters modelled as in Sect. 5.1.2 for z > 0.25. As in the previous approach, the simulation is generated for fixed cosmological parameters (h = 0.70, Ω m = 0.27, Ω Λ = 0.73, Ω b = 0.044, σ 8 = 0.78, n S = 0.95). The main advantage of this approach with respect to the previous one is that the simulation features a model for the local SZ diffuse component and filamentary structures up to z = 0.25. Note that the cosmological parameters used are different from those of Sect. 5.1.3.

Polarised SZ effect

The polarised SZ effect arises from local quadrupoles in the incoming radiation pattern as seen from the cluster (Audit & Simmons 1999;Sazonov & Sunyaev 1999;Liu et al. 2005). Such quadrupoles can be either cosmological (the CMB quadrupole as seen from the location of the cluster) or kinetic (due to a A96, page 24 of 35 (2010). Lower curves are obtained when the value of σ 8 is set to 0.75 instead of 0.809, and assuming an hydrostatic bias of 15% in the Y-M scaling law, i.e., the Y parameter is 15% lower than in Arnaud et al. (2010) for a given cluster mass. Red (resp. blue) curves are tSZ (resp. kSZ) power spectra for the model described in 5. 1.4, where the high redshift part is generated with the same two modelling alternatives (upper curves for σ 8 = 0.809 and no hydrostatic bias, lower curves for σ 8 = 0.75 and 15% hydrostatic bias). As the low-redshift part of the SZ maps uses maps at HEALP N side = 1024, power spectra for this case are computed only up to ℓ = 2000. The red diamond at ℓ = 3000 is the tSZ measurement of D tS Z 3000 = 3.65 ± 0.69 µK 2 obtained recently by SPT (Reichardt et al. 2012), improving on previous measurements from Lueker et al. (2010), Dunkley et al. (2011), andShirokoff et al. (2011). The blue arrow is the corresponding upper limit for kSZ. relativistic effect of the transverse motion of the cluster). This effect is neglected in the current implementation.

Limitations of the SZ models

The SZ maps derived from N-body+hydrodynamical simulations are stored in HEALP format, with N side = 1024 (pixel size of 3.44 ′ )a n dN side = 2048 (pixel size of 1.72 ′ ) respectively. Running the PSM at smaller pixel size would not upgrade the resolution of the maps, i.e., in spite of the larger number of pixels, no cluster substructures at scales below 1.72 ′ would be present.

The semi-analytic model does not suffer from this limitation, as clusters are modelled using analytic profiles which are scaled and normalised using scaling laws between flux, size and redshift, mass. This model, however, does not include cluster substructure, nor scatter in the scaling relations and profiles. Hence, two distinct clusters that have the same mass and redshift would be strictly identical on the modelled SZ map. In addition, the diffuse SZ effect is not modelled in the current version, as the SZ effect is generated by spherically symmetric clusters distributed on the sky with uniform probability and no correlation.

Figure 23 illustrates the large fluctuations these different options yield in terms of the thermal and kinetic SZ power spectra. The semi-analytic model can be adjusted to fit current measurements at high ℓ, adjusting the value of σ 8 and the Y-M scaling law, but has significantly less power on large scales than simulated SZ maps based on hydrodynamical simulations at low redshift.

The current implementation of the thermal SZ effect is based on maps of the Compton y parameter. The relativistic SZ effect is included only at some level for the models that are based on the generation of a cluster catalogue.

Radio sources

Radio sources in our model are mainly modelled on the basis of surveys of radio sources at frequencies ranging from 0.85 GHz to 4.85 GHz. A summary of these surveys is given in Table 2. Figure 24 illustrates their sky coverage. For the present purpose it is useful to distinguish four cases: Extrapolation of the flux of each source to all useful frequencies requires the knowledge of the spectral behaviour of the sources which, for any individual object frequently is quite complex (see e.g. [START_REF] Sadler | [END_REF]. For the present work, we resort to a power law approximation of the spectra of the sources, of the form S ν ∝ ν -α . Spectral indices, however, are different in a few pre-selected intervals of electromagnetic frequency.

For sources of the group (1) above, which have measurements at two frequencies, the individual spectral indices α below 5 GHz can be directly estimated as α =log (S 4.85 GHz /S low )/log (4.85/ν low ), (52) where ν low is either 1.4 GHz (if we are in the region covered by the NVSS) or 0.843 GHz (if we are in the lower declination region covered by the SUMSS). However, the calculation is not as straightforward as it may appear, because the surveys at different frequencies have different resolutions, implying that a single source in a low resolution catalogue can correspond to multiple sources at higher resolution. The spectral index estimates were carried out by degrading the higher resolution survey to the resolution of the other. In practice, whenever the higher resolution (NVSS or SUMSS) catalogue contains more than one source within a resolution element of the 4.85 GHz one, we have summed the NVSS or SUMSS fluxes, weighted with a Gaussian beam centred on the nominal position of the 4.85 GHz source, and with FWHM equal to the resolution of the 4.85 GHz survey.

On the other hand, the low frequency surveys, and especially the NVSS, are much deeper than the 4.85 GHz ones, which also have inhomogeneous depths. Simply summing all the lower frequency sources within a resolution element has the effect of including a variable fraction of the background, unresolved at 4.85 GHz, thus biasing the spectral index estimates. To correct for this, we have selected 159 195 control fields, free of have been subtracted from the summed NVSS or SUMSS fluxes associated with 4.85 GHz sources. In this way we obtained spectral indices from 1 to 5 GHz for a combination of complete 5 GHz selected samples with somewhat different depths, summing up to 109 152 sources over about 95% of the sky.

Although this is the best we can do with the available data, it is clear that the derived individual spectral indices are uncertain, due to a combination of several factors: measurement errors, uncertainties associated to the corrections applied, that are of statistical nature, and variability (the surveys have been carried out at different epochs). As a result, the absolute values of several individual spectral index estimates turned out to be unrealistically large, and the global distribution was found to be much broader than indicated by accurate studies of smaller samples. Furthermore, we obtain an estimate of the 5 GHz counts exceeding by an average factor of roughly 1.8 those directly observed, again indicating that the spectral index distribution is spuriously broadened. Therefore, we use these spectral index estimates only to assign the sources either to the steep-or to the flat-spectrum class (the boundary value being α = 0.5) and to determine the mean spectral index of each population. We find α steep = 1.18, α flat = 0.16. The spectral index distributions are approximated by Gaussians with variances σ steep,flat = 0.3, consistent with the results by Ricci et al. (2006). We then extrapolate the 5 GHz fluxes to 20 GHz by assigning to each source a spectral index randomly drawn from the Gaussian distribution for its class.

Sources with flux measurements at a single frequency have been randomly assigned to either the steep or to the flat-spectrum class in the proportions observationally determined by Fomalont et al. (1991) for various flux intervals, and assigned a spectral index randomly drawn from the corresponding distribution. The holes in the sky coverage have been filled by randomly copying sources from other regions in proportion to the surface density appropriate for each flux interval. The same procedure was adopted to add fainter sources in the regions where the existing surveys are shallower, until a coverage down to at least ≃20 mJy at 5 GHz over the full sky was achieved. We have checked that still fainter sources do not appreciably contribute to fluctuations in Planck channels for detection limits in the estimated range (approximately 200 to 500 mJy; López-Caniego et al. 2006), as expected since fluctuations are dominated by sources just below the detection limit. This check is carried out computing the power spectra of fluctuations due to sources below such limits in regions covered by the NVSS (the deepest survey) and in regions covered to shallower limits: the results are indistinguishable. This, however, holds for simulations of Planck observations. For other experimental setups, this check should be repeated considering a correspondingly different combination of resolution and sensitivity.

The regions less covered by the surveys used in the model, where the fraction of simulated sources is larger, are mostly around the Galactic plane, where they have a small effect compared to free-free and synchrotron emissions. At Galactic latitudes |b| > 10 • the fraction of real sources (at least as far as positions are concerned) is about 97%; over the full sky it is about 95%. Therefore we expect that the simulated maps faithfully reflect also the clustering properties of radio sources.

In Fig. 25 the source counts at 5 and 20 GHz obtained from our map are compared with observed counts, with the model by Toffolatti et al. (1998), and with an updated version of the model by de Zotti et al. (2005), allowing for a high-redshift decline of the space density of both flat-spectrum quasars (FSQs) and steep-spectrum sources (not only for FSQs as in the original model). The model adopts luminosity functions (in units of Mpc -3 (d log L) -1 ) of the form

Φ(L, z) = n 0 (L/L * ) a + (L/L * ) b , (53)
and lets those of steep-spectrum sources and of FSQs evolve in luminosity as

L * ,FSRQ (z) = L * (0)10 k ev z(2z top -z) , (54)
while for BL Lac objects a simpler evolutionary law is used

L * (z) = L * (0) exp[k ev τ(z)], (55)
where τ(z) is the look-back time in units of the Hubble time, H -1 0 , and where k ev and z top parametrise the luminosity evolution. The Fomalont et al. (1991), andHaarsma et al. (2000). Data in the 20 GHz panel are from the 9C survey (Waldram et al. 2003) at 15 GHz and from the ATCA survey at 18 GHz (Ricci et al. 2004a); no correction for the difference in frequency was applied.

new values of the parameters are given in Table 3. Figure 26 similarly shows a comparison of the radio source number count in the present model with predicted number counts from the Degree Angular Scale Interferometer (DASI) and the Very Small Array (VSA), and observed number counts from Planck, Planck ACTA Coeval Observations (PACO), the South Pole Telescope (SPT) and the Atacama Cosmology Telescope (ACT).

Note that we do not distinguish, for the moment, Galactic radio sources from extragalactic ones. While comparisons of number counts are made with cosmological evolution models for extragalactic radio sources, Galactic sources do affect only number counts at low galactic latitudes, typically |b| < 5 • .

Due to the complex spectral shape of radio sources, the power-law approximation holds only for a limited frequency range. To extrapolate the fluxes beyond 20 GHz we use the multifrequency WMAP data (Bennett et al. 2003) to derive the distributions of differences, δα, between spectral indices above and below that frequency. Such distributions can be approximated by Gaussians with mean 0.35 and dispersion 0.3. To each source we associate a spectral index change drawn at random from the distribution. Finally, for a small number of sources with exceedingly low negative values of α, yet another spectral index is assigned at frequencies higher than 100 GHz. In summary, each radio source in the model is simulated on the basis of four different power laws: One at frequencies below 5 GHz, another one between 5 and 20 GHz, another one between 20 and 100, and a last one at ν>100 GHz (which is, for most sources, the same as between 20 and 100 GHz). Our method of assigning to each source a different spectral index in these three regions of the electromagnetic spectrum naturally results in a source population that includes some sources peaked at 5 or at 20 GHz, and some sources with inverted spectra up to ν = 100 GHz.

We also produce maps and catalogues of polarised radio source emission attributing to each source a polarisation degree randomly drawn from the observed distributions for flat-and steep-spectrum sources at 20 GHz (Ricci et al. 2004b), and a polarisation angle randomly drawn from a uniform distribution.

The special case of WMAP sources

Sources available in the WMAP Point Source catalogue (Gold et al. 2011) have been included in the model as follows:

Given a WMAP source, we search for a match in the available radio sources in the catalogue built from low frequency surveys. The identification is made with the brightest radio source at 4.85 GHz within 11 ′ from the WMAP source location. This method provides radio counterparts for all the WMAP sources. The value of 11 ′ coincides with the one used by WMAP authors for source identification at 5 GHz.

The matching radio sources are removed from the radio source catalogue described in Sect. 5.2, and are included in a separate catalogue which merges measured low frequency and WMAP fluxes.

These WMAP sources are simulated with seven power laws: The first at frequencies below 5 GHz, and then between 5 and 23 GHz (WMAP K band), between 23 and 33 GHz (Ka band), between 33 and 41 GHz (Q band), 41 and 61 GHz (V band), 61 and 94 GHz (W band), and the final emission law for frequencies over 94 GHz. In this last frequency range, the spectral index is taken to be the same as between 61 and 94 GHz, unless it is positive (flux increasing with frequency) in which case it is set to 0.

Note that most of WMAP radio sources are highly variable (González-Nuevo et al. 2008). Including them in the simulations hence makes sense for simulating WMAP sky observations, but not for all types of simulations.

If the simulation is polarised, then the WMAP sources polarisation is set at the (simulated) value of their radio counterpart.

Limitations of the radio source model

Extrapolations of radio source counts to frequencies ≥30 GHz show evidence of substantial incompleteness below about 0.1 Jy (see Fig. 28). This is not a serious problem for WMAP or Planck-related simulations because the fluctuations due to unresolved radio sources are dominated by sources brighter than this limit. On the other hand, the current version of the PSM underestimates the amplitude of radio shot noise for experiments with higher spatial resolution and correspondingly lower confusion noise levels. Note, however, that at frequencies greater than 200 GHz, the radio shot noise becomes quickly negligible compared to that due to far-IR galaxies. In the next release of the PSM this problem will be cured by taking into account the data from the AT20G survey (Hancock et al. 2011;[START_REF] Murphy | [END_REF], that allows a much better determination of the fluxdensity dependent distribution of spectral indices up to 20 GHz, as well as the high frequency counts and spectral information from the SPT survey at 150 and 220 GHz [START_REF] Vieira | [END_REF] and from the ACT survey at 148 GHz (Marriage et al. 2011b). Further information on high-frequency source spectra is provided by the quasi-simultaneous multifrequency observations reported by Bonavera et al. (2011);Massardietal.(2011); Planck Collaboration (2011h,i); Procopio et al. (2011).

Far-infrared sources

IRAS sources

The software uses a compilation including all FIR sources taken from the IRAS Point Source Catalog (PSC, [START_REF] Beichman | IRAS Catalogs and Atlases[END_REF]) and the Faint Source Catalog (FSC, [START_REF] Moshir | IRAS Faint Source Survey[END_REF]).

A96, page 27 of 35 Sources identified as ultra-compact H  regions are removed from the catalogue, and given a special treatment (see Sect. 5.4).

For the remaining sources, fluxes are extrapolated to Planck frequencies, adopting modified blackbody spectra ν b B(ν, T ), B(ν, T ) being the blackbody function. For sources detected at only one IRAS frequency, the values of b and T are taken to be those of the average spectral energy distribution of the sample by Dunne et al. (2000), i.e., b = 1.3andT = 35 K. If the source is detected at 60 and 100 µm, then b = 1.3 is still assumed but the temperature is obtained from a fit to the data.

Since the PSC does not contain objects where the confusion from Galactic sources is high, and thus does not penetrate the Galactic centre well, and since the FSC is restricted to regions away from the Galactic plane, the source density is a function of Galactic latitude. For each PSM simulation, we add randomly distributed sources until the mean surface density as a function of flux matched everywhere the mean of well covered regions down to S 857 GHz ∼ 80 mJy. The coverage gaps of the catalogue (the IRAS survey missed about 4% of the sky) are filled by the same procedure.

To each source we assign a polarisation fraction drawn from a χ 2 distribution with one degree of freedom, the global level of which can be adjusted as an input parameter (and is 1% by default). The polarisation angle is randomly drawn from a uniform distribution.

In the near future, the catalogue of observed far-infrared sources implemented in the model will be complemented with

Cosmic Infrared Background anisotropies

An important, and possibly dominant, contribution to submillimetre small-scale anisotropies comes from galaxies selected by SCUBA and MAMBO surveys (see, e.g., Scott et al. 2006;Coppin et al. 2006), that are probably strongly clustered (Negrello et al. 2004(Negrello et al. , 2007)). These galaxies are interpreted as massive proto-spheroidal galaxies in the process of forming most of their stars in a gigantic starburst. We adopt the counts predicted by the model of Lapi et al. (2006), which successfully accounts for a broad variety of data including the SCUBA and MAMBO counts and the preliminary redshift distribution (Chapman et al. 2005).

The clustering properties of these sources are modelled as in Negrello et al. ( 2004 As first pointed out by Blain (1996), the combination of the extreme steepness of the counts determined by SCUBA surveys and of the relatively large lensing optical depth corresponding to the substantial redshifts of these sources maximises the fraction of strongly lensed sources at (sub)-mm wavelengths. We include such sources in our simulation by randomly distributing them with flux-dependent areal densities given by Perrotta et al. (2003). The frequency extrapolations are made via the spectral indices obtained in the same way as for the proto-spheroidal galaxies.

A96, page 28 of 35 Figure 29 shows the power spectrum of cosmic infrared background anisotropies as computed on a simulated sky, and as measured by Planck (Planck Collaboration 2011k), and ACT (Dunkley et al. 2011). The agreement between the model and the measured power spectra is reasonable, but not perfect, and discrepancies will be addressed with improved modelling in future releases, in the light of recent advances on number counts and on correlations from observations with the BLAST balloon-borne experiment (Patanchon et al. 2009;Viero et al. 2009), SPT (Hall et al. 2010), Herschel (Maddox et al. 2010;Cooray et al. 2010), and from a combined analysis of ACT and BLAST (Hajian et al. 2012). Also, the correlation between the CIB maps across the frequency range is almost unity in the current model, a feature that will have to be improved in the future using a more detailed model constrained by upcoming additional Planck observations. Cosmic infrared background anisotropies as implemented in the current model are assumed unpolarised.

Limitations of the far-infrared source model

The main uncertainties in the model for relatively nearby farinfrared galaxies arise from the extrapolations of IRAS flux densities. The Planck data (Planck Collaboration 2011j) show evidence for colder dust than has previously been found in external galaxies so that our model likely underestimates the (sub)-mm flux densities of at least a fraction of IRAS galaxies and, consequently, the source counts. Moreover, the CO emission lines are not included in the simulations. The CO(J = 1 → 0), CO(J = 2 → 1), CO(J = 3 → 2), that, for low-z galaxies, fall within the Planck 100, 217 and 353 GHz passbands, respectively, may contribute significantly to the observed fluxes, while higher order CO transitions are probably negligible. A careful investigation of the possible contamination by these lines is in progress in preparation for future releases of the PSM.

A wealth of data on faint far-IR sources, those that dominate the small-scale fluctuation at frequencies higher than 200 GHz, have become available after the extragalactic source model was worked out. These include source counts (Clements et al. 2010;Oliver et al. 2010;Béthermin et al. 2010;Glenn et al. 2010) and estimates of the power spectrum of the cosmic infrared background anisotropies (Planck Collaboration 2011c; Amblard et al. 2011). The models implemented in the PSM fare reasonably well with the new data, but the steep (sub-)mm spectrum of dusty galaxies, coupled with the fast cosmological evolution, make both the faint counts and the CIB power spectrum exceedingly sensitive to the detailed spectral energy distribution (SED) of model galaxies. For example, the amplitude of the (sub-)mm power spectrum scales as the frequency to a power of 7-8 so that relatively small differences in the SED are enough to yield large discrepancies between the model and the data. Much improved models, successfully reproducing the recent data, are now available both for the epoch-dependent luminosity function (hence for the source counts; e.g. Béthermin et al. 2011;Lapi et al. 2011) and for the clustering properties (e.g. Xia et al. 2012;Pénin et al. 2012). These models will be exploited in next releases of the PSM.

Galactic point sources

Most Galactic compact sources are currently treated in the model as part of the diffuse emission. IR sources present in the IRAS catalogue are treated as extragalactic sources (dusty galaxies). A special treatment has been given to simulate emission from ultra-compact H  (UCH ) regions, however.

A list of UCH  regions and UCH  region-candidates has been extracted from the IRAS Point Source Catalogue with fluxes from radio counterparts. A total of 864 sources were selected from IRAS PSC according to the following criteria.

1. IRAS colours as specified in Kurtz et al. (1994). 2. Fluxes in the IRAS channels 100 µm, 60 µm and 25 µmare not upper limits (i.e., they are at least of "medium" quality, according to definitions given in [START_REF] Beichman | IRAS Catalogs and Atlases[END_REF]). 3. The flux at 100 µm was required to be greater than or equal to 100 Jy. 4. The Galactic latitude is below 10 • . This list includes all confirmed UCH  regions from the original sample of Kurtz et al. (1994) that have an IRAS counterpart (48 objects) and 38 UCH  regions out of 53 of the sources in Wood & Churchwell (1989).

Note that as some of these sources are present in the dust map used for modelling thermal dust emission, that map has been processed to subtract the UCH  regions, to avoid doublecounting.

Radio counterparts have been searched for in the GB6 and NVSS catalogues. The search radius is 1 ′ ; when more than one counterpart was found within the search radius (within that catalogue), the brightest was selected. None of these radio catalogues covers the Southern hemisphere: search for radio counterparts in catalogues of the southern hemisphere is in progress.

For all 864 sources, the fluxes in the useful frequency range are estimated by performing a modified blackbody fit to the IRAS fluxes at 100 and 60 µm. The emission is parametrised with the form

S = ν β B(ν) = ν β 2hν 3 c 2 1 exp(hν/kT ) -1 , ( 56
)
with β = 1.5(seee.g. Gear et al. 1988;Hoare et al. 1991;Maxia et al. 2001). When a radio counterpart to the IRAS source is found, fluxes at low frequencies are corrected by adding to the modified blackbody fit an estimated low-frequency flux extrapolating the flux at the frequency of the radio catalogue (where the counterpart was found) with a free-free spectral index. When a radio counterpart to the IRAS source is found in more than one of the three radio catalogues priority was given to GB6, then to Giveon et al. (2005) and then NVSS.

In case of an erroneous radio counterpart this procedure may lead to slightly overestimated fluxes at the lowest HFI frequencies, but this has little impact for frequencies greater then 200 GHz, where we expect these sources to be very relevant for Planck.

UCH  regions are free-free emitters and lack magnetic fields to cause grains to align. Hence, they are not intrinsically polarised, and are thus modelled as unpolarised in the current version of the model.

Conclusions

We have developed a complete, flexible model of multicomponent sky emission, which can be used to predict or simulate astrophysical and cosmological signals at frequencies ranging from about 3 GHz to 3 THz. The model, which actually implements several options for each of the components, has been developed for testing component separation methods in preparation for the analysis of Planck data, but has been used also for simulating data in the context of analyses of WMAP observations, or for planning future missions such as COrE (The COrE Collaboration et al. 2011). Given the high sensitivity of these instruments, the quality of the reconstructed CMB temperature map depends strongly on our ability to remove the contamination by astrophysical foregrounds. It is therefore necessary to build simulations as close as possible to the complexity of the real sky over a large frequency range. We present here what has been achieved through a series of improvements accomplished over several years and the contributions of experts in different fields: CMB, Galactic emission and compact sources, extragalactic radio and far-IR sources, Sunyaev-Zeldovich signals from clusters of galaxies and the intergalactic medium. Note however that the current version of the model cannot be expected to match perfectly data sets that have not been used in the model, in particular those of the Planck mission and of other upcoming observations. Updates will be proposed as such data sets become available.

The maps of each astrophysical component have been repeatedly updated by requiring a close match with the constantly increasing amount of observational data. Figure 30 shows the spectrum of root mean square (rms) fluctuations in simulated maps at WMAP frequencies for |b| > 20 • and 10 • resolution. The various diffuse components have different spectral characteristics while the total signal is a good match to the WMAP 7-year data; the rms values agree to within better than 5 percent. The sky emission maps observed with WMAP, smoothed to 1 • resolution, are compared in Fig. 31 to PSM predicted emission maps at the same frequencies and resolution. The agreement is excellent over most of the sky, except in the galactic ridge and around a few compact regions of strong emission, for which the model does not exactly reproduce the observations. Discrepancies are due to a mixture of uncertainties in the exact resolution of the model (and, in particular, a lack of resolution of the map of synchrotron spectral index), errors in extrapolation of the dust emission, residuals of galactic emission in the predicted CMB map, and insufficient number of parameters used in the current PSM in regions of strong emission (in particular when emission comes from several distinct regions along the line of sight). The model will be refined in future version, in particular using Planck observations.

The CMB temperature and polarisation maps, currently based on WMAP observations (maps and best-fit cosmological model), are complemented with simulations to allow for the weak gravitational lensing by gradients in the large-scale gravitational field. Weak gravitational lensing has been dealt with in the Born approximation. Simulations of CMB temperature maps with primordial non-Gaussianity of local type are also accessible via our sky model.

The Galactic diffuse emission model includes five components: synchrotron, free-free, spinning dust and thermal dust radiation, and 12 CO (J = 1-0), (J = 2-1), and (J = 3-2) molecular lines at 115.27, 230.54, and 345.80 GHz, respectively. For each component, alternative models are available but we propose a combination of models that reproduces best the available data. In this default model the synchrotron emission is based on the 408 MHz all-sky map by Haslam et al. (1982) Fig. 31. Comparison of sky emission as observed by WMAP (7-year data) and as predicted by the PSM in the same frequency bands, at a resolution of 1 • . For each frequency channel, the colour scale is the same for WMAP (left column) and PSM prediction (middle column). An histogram equalised colour scale is used for the K, Ka, and Q channels, and a linear scale for the V and W channels. Maps are saturated to highlight common features away from the galactic plane. Maps of difference between PSM prediction and WMAP observation are displayed in the right column (note that the colour scale is different from that used to display the K, KaandQ maps), highlighting discrepancies in the galactic plane specifically and at the location of a few regions of compact emission. The agreement is excellent over most of the sky away from the galactic ridge and a few compact regions. extrapolated in frequency exploiting the spectral index map by Miville-Deschênes et al. (2008;model 4). The free-free template is obtained from the WMAP MEM map, complemented with the Hα all-sky template by Finkbeiner (2003) corrected for extinction as in Dickinson et al. (2003). For thermal dust we adopted model 7 of Finkbeiner et al. (1999), which features spectral variations across the sky. For spinning dust we adopted the template produced by Miville-Deschênes et al. (2008) from an analysis of WMAP data. The CO emission map is based on the 12 CO(1-0) survey by Dame et al. (2001). Standard intensity ratios with the (J = 1-0) line have been used for the (J = 2-1) and (J = 3-2) transitions. Since all these templates have a resolution insufficient for our purposes, the code allows the possibility to add small-scale fluctuations following Miville-Deschênes et al. (2007).

In our model as currently implemented, the only diffuse emissions that are polarised (apart from the CMB) are synchrotron and thermal dust. For synchrotron we rely on the WMAP 23 GHz polarisation maps, extrapolated in frequency using the same spectral index map used for intensity. Modelling dust polarisation is made difficult because of the paucity of data. In our model the dust Q and U maps are built assuming a constant intrinsic polarisation fraction, geometrical depolarisation and polarisation angle maps constructed using a Galactic A96, page 31 of 35 magnetic field model for scales larger than 20 • and constraints from 23 GHz WMAP polarisation data to reproduce the projected structure at smaller scales due to the turbulent magnetic field of the ISM.

While generally compact Galactic sources are part of the diffuse emission map and Galactic point sources are treated in the same way as extragalactic point sources (catalogues do not distinguish between Galactic and extragalactic point sources), ultra-compact H  regions have been included in the model as a separate population. They have been selected from the IRAS catalogue and the extrapolation in frequency of their flux densities is made adopting a grey-body spectrum. Radio counterparts have been searched for in the NVSS and GB6 catalogues.

The Sunyaev-Zeldovich effect from galaxy clusters is simulated in two ways. One can start from the epoch-dependent cluster mass function, for the preferred choice of cosmological parameters, plus some recipes to model the density and temperature distribution of hot electrons. A fraction of simulated clusters can be replaced by real clusters drawn from the ROSAT all-sky or the SDSS catalogue. Alternatively, we can resort to N-body+hydrodynamical simulations that also contain the SZ emission from the web of intergalactic hot gas.

Most radio sources in the model are real, taken from the relatively deep all-sky low-frequency catalogues. Each source has its own spectral properties, either determined directly from multifrequency data or randomly extracted from the observed distributions of spectral indices. Extrapolations to higher frequencies are made using different spectral indices for different frequency intervals, in order to ensure consistency with multi-frequency source counts. This approach turned out to be remarkably successful in predicting the counts that have been later measured by Planck. The polarisation degree of each source is randomly drawn from the distributions for flat and steep-spectrum sources determined by Ricci et al. (2004b) at 20 GHz, while the polarisation angle is drawn randomly from a uniform distribution.

Bright far-IR sources include those taken from the IRAS Point Source Catalog, with flux densities extrapolated in frequency using a grey-body spectrum, plus high-z galaxies strongly gravitationally lensed, based on the model by Perrotta et al. (2003). For fainter galaxies, which make up most of the cosmic infrared background, we adopted the model by Granato et al. (2004). Their clustering properties are modelled following Negrello et al. (2004). The implementation of their spatial distribution is made using the algorithm by Gonzalez-Nuevo et al. (2005). The resulting power spectrum of CIB fluctuations turns out to be quite close to that measured by the Planck collaboration (Planck Collaboration 2011k). The mean polarisation degree of individual infrared galaxies can be adjusted as an input parameter, the default value being 1%. The polarisation angle is randomly chosen from a uniform distribution.

On the whole, these simulations should provide a useful tool for several purposes: to test data analysis pipelines for CMB experiments, identify the most convenient areas of the sky for specific purposes (i.e., areas least contaminated by Galactic emissions), identify the optimal set of frequency channels for CMB temperature and polarisation experiments or for other studies (i.e., for studies of the cosmic infrared background), predict the levels of confusion noise, including the effect of clustering for a given frequency and angular resolution, and much more. Access to documented versions of the package and to reference simulations is available from a dedicated website9 .

Fig. 1 .

 1 Fig. 1. CMB dipole prediction as generated by the code (here for a 1 • "sky resolution", in Galactic coordinates).

Fig. 2 .

 2 Fig. 2. CMB intensity and polarisation power spectra used by default (for Gaussian CMB simulations) -black: C TT ℓ ;red :C TE ℓ (dashed parts are negative); green: C EE ℓ ;blue:C BB ℓ (from scalar modes only).

Fig. 3 .

 3 Fig.3. CMB temperature and polarisation anisotropy predictions as generated for a 1 • beam, N side = 256, in Galactic coordinates. This prediction is based on the CMB extracted from WMAP 5-year data using an ILC in needlet space. Note that the maps are not fully exempt from contamination by foregrounds and noise (seeDelabrouille et al. 2009, for a complete discussion of the component separation method used to obtain the CMB temperature map).

Fig. 4 .

 4 Fig. 4. Linear part (top panel) and non-linear correction (bottom panel) for a non-Gaussian CMB realisation, at 1 • resolution. Note the different colour scales, and that f NL ∼ 1 (or a few) is a typical expected level. The cosmological parameters assumed for this simulation are from the fit of WMAP 7-year + BAO + SN observations, as made available on the LAMBDA web site.

Fig. 5 .

 5 Fig. 5. Simulated lensed CMB temperature anisotropy map (top left), the amplitude of the difference of a simulated lensed and unlensed CMB temperature anisotropy map (top right), the amplitude of the simulated deflection field map (bottom left) and the lensing potential map (bottom right) with HEALP pixelisation parameter N side = 1024. These maps are obtained using the NFFT for the oversampling factor σ = 2a n d convolution length K = 4.

Fig. 6 .

 6 Fig. 6. Portion of a simulated lensed CMB temperature anisotropy map (top left), the amplitude of the difference of a simulated lensed and unlensed CMB temperature anisotropy map (top right), the amplitude of the simulated deflection field map (bottom left) and the lensing potential map (bottom right) with HEALP pixelisation parameter N side = 1024. These maps are obtained using the NFFT for the oversampling factor σ = 2 and convolution length K = 4. The size of the maps displayed here is about 16 • on a side.

  1. A constant value on the sky (set to -3.0 by default); 2. the spectral index map of Giardino et al. (2002) obtained using a combination of 408, 1420 and 2326 MHz data; 3. the spectral index map provided by model 4 of Miville-Deschênes et al. (2008) obtained using a combination of the 408 MHz and WMAP 23 GHz data. This model can be used with or without steepening at frequencies higher than 23 GHz.

Fig. 8 .

 8 Fig. 8. Template map at 408 MHz used to model the intensity of synchrotron emission in the PSM. The colour scale of the bottom panel is histogram-equalised to increase the dynamic range.

Fig. 9 .

 9 Fig. 9. Synchrotron spectral index map used by default in the PSM.

Fig. 10 .

 10 Fig. 10. Emission law of the free-free, spinning dust, and main CO molecular line emission (in units of spectral brightness I ν , e.g., ∝MJy sr -1 ). The spinning dust emission law is here normalised to unity at ν = 23 GHz and the free-free emission law to I ν = 2 at the same reference frequency. The integrated amplitude of the first 12 CO emission lines (transition (J = 1-0)) is normalised to unity, and the other lines (transitions (J = 2-1) and (J = 3-2)) are displayed in proportion to their integrated intensity relative to the first one.

Fig. 11 .

 11 Fig. 11. Template map at 23 GHz used to model the intensity of freefree emission in the PSM. As seen in the top panel, free-free emission is strongly concentrated in compact regions of the Galactic plane. The colour scale of the bottom panel is histogram-equalised to increase the dynamic range, and to show more extended, diffuse structures.

1 .

 1 those of CO molecules involving other isotopes of carbon and oxygen (the most important being the 13 CO emission at multiples of 110.20 GHz); 2. a set of HCN lines around multiples of 88.63 GHz; 3. the HCO + lines at multiples of 89.19 GHz 8 . Many other lines, such as those of CN, C 2 H and HNC, which have also been recently detected in high density molecular clouds in M82 (Naylor et al. 2010), are known to exist in the frequency range covered by the PSM. However, their contribution to the observed broad-band emission in the frequency range covered by the PSM is small. The 12 CO lines are of special interest, as they are the strongest ones located inside three of the Planck HFI frequency bands (the 100, 217, and 353 GHz channels; see Planck HFI Core Team 2011). They are thus important in the global model of sky emission. The 12 CO emission is implemented in our model on the basis of the 12 CO (J = 1 → 0) survey map of Dame et al.

Fig. 12 .

 12 Fig. 12. Map of the ratio of 100 µm to 240 µm emission, as obtained by Finkbeiner et al. (1999) and used in the present model.

Fig. 13 .

 13 Fig. 13. Map of thermal dust emission 100 µm. Colours are saturated at 100 MJy sr -1 , and the colour scale is histogram-equalised to increase the dynamic range and make both regions of low and high intensity visible.

1 .Fig. 14 .

 114 Fig.14. Top: small patch of free-free emission before and after adding random small scale fluctuations. Middle and bottom: power spectrum of free-free emission before and after adding small scales.

Fig. 15 .

 15 Fig. 15. Top: small patch of spinning dust emission before and after adding random small scale fluctuations. Middle and bottom:p o w e r spectrum of spinning dust emission before and after adding small scales.

Fig. 16 .

 16 Fig. 16. E and B power spectra of diffuse Galactic emission simulated with the model at WMAP central frequencies (solid and dashed thick lines respectively). The P06 Galactic mask is used. The WMAP derived foreground levels from Gold et al. (2011) are also shown (thin lines).

Fig. 17 .

 17 Fig. 17. Synchrotron Q and U maps at 23 GHz for a sky resolution of 3 • . At this frequency the synchrotron Q and U maps of the sky model are exactly the WMAP 7-year data. The two bottom panels show the polarisation fraction and the polarisation angle as implemented in our model. The average polarisation fraction is about 18%.

  Equations from (32)t o( 41) show the difficulties of modelling polarised emission. This requires a three-dimensional model of all dust properties (n d , ǫ d , β d and T d ), of the energetic electrons, and of the magnetic field geometry. To simplify this task we assume that ǫ d , β d and T d are constant along a given line of sight. Then, following Page et al. (2007), Miville-Deschênes et al. (2008)andFauvet et al. (2011), we model the density distribution of dust and electrons with a galacto-centric radius and height function n = n 0 exp(-r/h r ) sech 2 (z/h z ), (42) with scaling parameters for the galacto-centric radius and height dependence equal to h r = 3 kpc and h z = 0.1 kpc for dust, and h r = 5 kpc and h z = 1 kpc for electrons. Finally we use the model of the magnetic field described in Miville-Deschênes et al. (

Fig. 18 .

 18 Fig.18. Maps of the large-scale geometrical depolarisation g ′ and polarisation angle γ ′ for the synchrotron and dust, based on a model of the Galactic magnetic field and density distributions of the energetic electrons and dust grains. These maps are used to correct the dust depolarisation and polarisation angle maps deduced from the 23 GHz polarisation data for the difference in the large-scale structure between synchrotron and dust (see Eqs. (43)and(44)).

Fig. 19 .

 19 Fig. 19. Model of the dust polarisation at 200 GHz and at 3 • resolution. The polarisation fraction and polarisation angle maps are slightly modified versions of the same maps obtained for synchrotron at 23 GHz using the WMAP and 408 MHz data.

Fig. 20 .

 20 Fig. 20. Density of MaxBCG clusters. The clusters are mainly concentrated in the northern Galactic hemisphere.

Fig. 21 .

 21 Fig. 21. Top panel: cluster number counts for different mass functions. Bottom panel: cluster number counts for different values of σ 8 ,u s i n g the Tinker mass function.

Fig. 22 .

 22 Fig. 22. Maps of thermal and kinetic SZ effect from the Hubble volume simulation (left column), from the local universe (middle column), and total thermal and kinetic SZ effects from both simulations together, smoothed here to a sky resolution of 5 arc-min (right column). All the displayed maps are small 5 • × 5 • patches centred on the North Galactic Pole. The original full sky maps are at HEALP N side = 2048 for the Hubble volume, and at N side = 1024 for the local hydrodynamical simulations. The contribution from the local universe comprises a large cluster at sky coordinates close to those of the Coma cluster, which dominates the total SZ emission in the maps displayed here.

Fig. 23 .

 23 Fig. 23. Power spectra D ℓ = ℓ(ℓ+1)C ℓ /2π of thermal and kinetic SZ effect, as modelled by the PSM for various options. Black (resp. green) curves give the spectrum of the tSZ (resp. kSZ) maps as modelled on the basis of a population of clusters as described in 5.1.2. Upper curves are obtained using the mass function of Tinker et al. (2008) with WMAP 7year cosmological parameters and the Y-M scaling law of Arnaud et al. (2010). Lower curves are obtained when the value of σ 8 is set to 0.75 instead of 0.809, and assuming an hydrostatic bias of 15% in the Y-M scaling law, i.e., the Y parameter is 15% lower than in Arnaud et al. (2010) for a given cluster mass. Red (resp. blue) curves are tSZ (resp. kSZ) power spectra for the model described in 5.1.4, where the high redshift part is generated with the same two modelling alternatives (upper curves for σ 8 = 0.809 and no hydrostatic bias, lower curves for σ 8 = 0.75 and 15% hydrostatic bias). As the low-redshift part of the SZ maps uses maps at HEALP N side = 1024, power spectra for this case are computed only up to ℓ = 2000. The red diamond at ℓ = 3000 is the tSZ measurement of D tS Z 3000 = 3.65 ± 0.69 µK 2 obtained recently by SPT (Reichardt et al. 2012), improving on previous measurements from Lueker et al. (2010), Dunkley et al. (2011), and Shirokoff et al. (2011).The blue arrow is the corresponding upper limit for kSZ.

1 .

 1 green points are sources with fluxes measured at both approximately 1 GHz (NVSS or SUMSS) and 4.85 GHz (GB6 or PMN): 109 152 sources; 2. blue (yellow) points are sources present only in the NVSS (SUMSS) and thus with fluxes only at 1.4 (0.843) GHz: 1 813 551 (158 284) sources; 3. red points are sources with only PMN (4.85 GHz) fluxes: 8656 sources; 4. the small white regions are those not covered by any survey.

  4.85 GHz sources, and computed the average flux of NVSS or SUMSS sources within a 4.85 GHz beam pointing on the field centre, again taking into account the 4.85 GHz response function. The average fluxes of control fields, ranging from 1.67 to 3.16 mJy for the different combinations of low and high frequency catalogues (NVSS/GB6, NVSS/PMN, SUMSS/PMN), A96, page 25 of 35A&A 553, A96 (2013)

Fig. 24 .

 24 Fig. 24. Sky coverage of the surveys listed in Table 2, in Galactic coordinates. Green points: sources present in both ≃1 GHz (NVSS or SUMSS) and 4.85 GHz (GB6 or PMN) catalogues; blue points: sources in the NVSS catalogue only; yellow points: sources in the SUMSS catalogue only; red points: sources in the PMN catalogue only; white regions: not covered by any survey.

Fig. 25 .

 25 Fig. 25. Modelled source number counts at 5 and 20 GHz for one sky realisation, normalised to ∆N 0 = S (Jy) -2.5 , compared with models and observational data. Data at 5 GHz are from Kellermann et al. (1986),Fomalont et al. (1991), andHaarsma et al. (2000). Data in the 20 GHz panel are from the 9C survey(Waldram et al. 2003) at 15 GHz and from the ATCA survey at 18 GHz(Ricci et al. 2004a); no correction for the difference in frequency was applied.

Fig. 26 .

 26 Fig. 26. Comparison of modelled radio source counts (red points) and the Planck radio counts (blue points, Planck Collaboration 2011g). Also shown are: the counts estimated at 31 GHz from DASI (grey dashed box, Kovac et al. 2002) and PACO (grey diamonds, Bonavera et al. 2011), at 33 GHz from the VSA data (grey box, Cleary et al. 2005), and the SPT (grey squares, Vieira et al. 2010) and ACT (grey stars, Marriage et al. 2011b) counts of radio sources.

Fig. 27 .

 27 Fig. 27. Emission law for a selected sample of typical radio sources produced in a realisation of the sky emission with our model.

  ), using their more physical model 2. The simulation of their spatial distribution is produced using the method by González-Nuevo et al. (2005) .O n c eam a p of the source distribution is obtained at the reference frequency ν ref , extrapolations to any other frequency ν i are obtained via the flux-dependent effective spectral indices α = log (S ref /S i )/log (ν ref /ν i ), where S i is defined by n(>S i ; ν i ) = n(>S ref ; ν ref ). The spectral indices are computed in logarithmic steps of ∆ log (S ref ) = 0.1. We check that the counts computed from the extrapolated maps accurately match those yielded, at each frequency, by the model.

Fig. 28 .

 28 Fig. 28. Similar to 26, but with additional green points showing the counts of WMAP radio sources as represented in our model. The solid black curve shows the total number counts of extragalactic radio sources as predicted by the updated model of de Zotti et al. (2005).

Fig. 29 .

 29 Fig. 29. Cosmic infrared background power spectrum. Solid lines are obtained from a simulation. Data points are from Planck observations (Planck Collaboration 2011k). Dashed lines at high multipoles for 143 and 217 GHz are from the Dunkley et al. (2011) best-fit model for IR sources (at 148 and 218 GHz, extracted from Fig. 2 of their paper). The CMB power spectrum (in black) is a theoretical model fitting WMAP observations.

Fig. 30 .

 30 Fig. 30. Temperature rms fluctuations at WMAP frequencies for |b| > 20 • at a resolution of 10 • . The symbols represent the fluctuations in the various diffuse components of the sky model, the total simulated fluctuations, and WMAP 7-year maps. The total signal is a good match to the WMAP data.

Table 1 .

 1 Summary of the main astrophysical components included in the current version of the model. See text for detailed description of each component.

	Component	Description
	Dipole	WMAP dipole

Table 2 .

 2 Summary of the large-area surveys of point sources used in this work.

	Frequency Catalogue S lim (mJy)	Dec range	Angular resolution
		GB6	18	0-+75	3.5 ′	Gregory et al. (1996)
	4.85 GHz 1.4 GHz 0.843 GHz	PMNE PMNT PMNZ PMNS NVSS SUMSS	40 42 72 20 2.5 18 8	-9.5-+10 -29--9.5 -37--29 -87.5--37 -40-+90 -50--30 -90--50	4.2 ′ 4.2 ′ 4.2 ′ 4.2 ′ 45 ′′ 45 ′′ cosec(|δ|) 45 ′′ cosec(|δ|)	Griffith et al. (1995) Griffith et al. (1994) Wright et al. (1996) Wright et al. (1994) Condon et al. (1998) Mauch et al. (2003)

Table 3 .

 3 Best-fit values of the parameters of the evolutionary models for canonical radio sources.

	Source type		Luminosity function	Evolution
		log n 0 (Mpc -3 )	ab	log L k ev	z top
	FSQ	-8.989	0.658 2.938	34.043	0.224 2.254
	BL Lac	-7.956	0.975 1.264	32.831	1.341
	Steep	-7.389	0.729 2.770	33.177	0.262 2.390

* (erg s -1 Hz -1 at 5 GHz, z = 0)

Planck (http://www.esa.int/Planck) is a project of the European Space Agency -ESA -with instruments provided by two scientific Consortia funded by ESA member states (in particular the lead countries: France and Italy) with contributions from NASA (USA), and telescope reflectors provided in a collaboration between ESA and a scientific Consortium led and funded by Denmark.

http://www.apc.univ-paris7.fr/~delabrou/PSM/psm. html

http://lambda.gsfc.nasa.gov/

http://camb.info

http://class-code.net

http://en.wikipedia.org/wiki/IDL_(programming_ language)

Bardeen's gauge-invariant potential Φ is defined in such a way that the temperature anisotropy reduces to ∆T/T = -Φ/3 in the pure Sachs-Wolfe limit.

http://www.splatalogue.net/

http://www.apc.univ-paris7.fr/~delabrou/PSM/psm. html

Acknowledgements. We wish to thank our colleagues from the Planck Collaboration, and in particular members of the component separation working group (WG2), for useful suggestions and discussions, as well as for beta-testing the successive versions of the software. We thank Karim Benabed, Frédéric Guilloux, Frode Hansen and Benjamin Wandelt for useful discussions concerning some of the CMB models. François Boulanger, Rodney Davies and François-Xavier Désert have provided expertise for developing the model of Galactic emission used in this work. We thank Francesca Perrotta, Grazia Umana and Stephen Serjeant for useful discussions concerning the point source model, Benjamin Walter for his help in the editing of the present paper, Hans-Kristian Eriksen, Torsten Ensslin, Andrew Jaffe, Lloyd Knox and Douglas Scott for useful comments on the PSM paper draft, and Pei Yu for her help setting-up a data repository for PSM activities. This work has been developed largely within the Planck Collaboration. A description of the Planck Collaboration and a list of its members, indicating which technical or scientific activities they have been involved in, can be found at (http://www.rssd.esa.int/ index.php?project=PLANCK\&page=Planck_Collaboration). The development of Planck has been supported by: ESA; CNES and CNRS/INSU-IN2P3-INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MICINN and JA (Spain); Tekes, AoF and CSC (Finland); DLR and MPG (Germany); CSA (Canada); DTU Space (Denmark); SER/SSO (Switzerland); RCN (Norway); SFI (Ireland); FCT/MCTES (Portugal); and PRACE (EU). SB has been supported by a postdoctoral grant from the "Physique des deux Infinis" (P2I) Consortium. CD acknowledges an STFC Advanced Fellowship and an ERC IRG grant under the FP7.