3 Searching for Continuous Gravitational Wave signal with PTA RÉSUMÉ La Relativité Générale (RG) prédisait l'existence des Ondes Gravitationnelles (OG) il y a 100 ans de cela. Dans cette théorie, l'espace-temps n'est plus une entité rigide et il peut être courbé par la présence de matière et d'énergie. Les OGs sont précisément le résultat de la propagation de ces courbures dans l'espace-temps. Elles sont de très faible amplitude et leur détection représente un véritable défi technologique. Seules des sources très massives peuvent émettre des OGs suffisamment puissantes pour être mesurées. En 2016, la collaboration LIGO/Virgo détectait pour la première fois une OG provenant de la fusion de deux trous noirs. Cette découverte a révélé un tout nouveau spectre d'observations astrophysiques qui est aujourd'hui un domaine de recherche très actif. Pour détecter une OG, il faut pouvoir mesurer les courbures de l'espace-temps induites par le passage de celle-ci. Une méthode, basée sur un réseau de pulsars chronométrés ou Pulsar Timing Array (PTA), propose d'utiliser les pulsars (des étoiles à neutron en rotation rapide qui sont perçues sur Terre comme des séries de pulses de lumière très régulièrement espacés dans le temps) comme des "horloges cosmiques" pour mesurer les perturbations dues à l'OG. Ce manuscript se concentre sur la détection d'OGs basse fréquence (nanohertz) en utilisant les données de chronométrage haute précision d'un réseau de pulsars milliseconde obtenu à l'aide de radiotéléscopes.

Le premier chapitre donne un aperçu du bagage théorique nécessaire pour comprendre l'origine des OGs. Pour comprendre la RG, il est essentiel de comprendre les bases de la géométrie différentielle et Riemanienne. Elles donnent un cadre d'étude des espaces courbes qui est le coeur de la description mathématique de la RG. L'équation centrale de la RG est l'équation de champs d'Einstein. Elle décrit la relation qui existe entre la géométrie de l'espace-temps d'une part et la présence de matière d'autre part. Pour des très petites perturbations de la métrique de l'espace-temps, il est possible de montrer que cette équation s'apparente à une équation de propagation d'une onde. Cette dérivation est détaillée dans ce chapitre car elle correspond à la description mathémathique d'une OG qui nous sera utile pour la suite. Typiquement, les sources d'OG seront des systèmes binaires d'objets très massifs (trous noirs, étoiles à neutrons, ...) comme il a été le cas pour la première OG détectée en 2016. Nous présentons ici la forme des OGs émises par des systèmes binaires en orbite circulaire l'un autour de l'autre. Il apparaît que l'amplitude de ces ondes est proportionnelle à la masse des objets du système puis inversément proportionnelle à la distance de la source. Il en vient alors que plus le système binaire est massif ou plus il est proche, plus les effets dues à l'OG émise se feront sentir. Les objets les plus massifs de l'Univers que l'on connaît à ce jour sont les trous noirs supermassifs. Leur origine reste encore un mystère car il est difficile d'établir un modèle de formation de ces astres en un temps aussi "court" que celui de l'âge de l'Univers. Notamment, l'observation de noyaux de galaxies actifs très anciens compor-tant des trous noirs extrêmement massifs en leur centre indique que ces objets se sont formés rapidement, probablement par fusion d'autres trous noirs supermassifs. Cependant, nous ne disposons d'aucune observation de fusion de trous noirs supermassifs à ce jour. Nous argumentons qu'un réseau de pulsars chronométrés (Pulsar Timing Array ou PTA) pourrait détecter des OGs générées par des systèmes de trous noirs super massifs binaires. Comme indiqué précédemment, les pulsars peuvent êtres utilisés comme des "horloges cosmiques" donnant accès à une mesure des propriétés de l'espace-temps dans notre Galaxie. Nous dérivons une expression mathématique pour les résidus temporels induits par les OGs dans la donnée chronométrée des pulsars. Lorsqu'un tel signal est stable et monochromatique, on le qualifie d'Onde Gravitationnelle Continue (OGC).

Ce manuscript est centré sur l'étude des OGCs émises par des systèmes binaires de trous noirs supermassifs.

Dans le second chapitre, nous présentons les méthodes d'analyse de donnée pour la détection d'OGs dans PTA. Pour cela, nous expliquons brièvement l'inférence Bayesienne qui est au coeur de ces méthodes. Elle nous permet d'évaluer la significativité statistique de nos modèles ainsi que les valeurs des paramètres de signaux que l'on cherche à détecter en fonction de la donnée mesurée (dans notre cas, la donnée chronométrée des pulsars). Une quantité centrale de l'inférence Bayesienne est le Facteur de Bayes. Il donne une mesure comparative des "poids" statistiques de différents modèles et est crucial pour sélectionner ceux qui reproduisent le mieux la donnée observée. Nos modèles et méthodes d'analyse doivent être capables de : (i) quantifier tous les bruits présents dans la donnée pour pouvoir isoler le signal d'OG désiré (ii) estimer un facteur de Bayes pour les OGs que l'on cherche à détecter pour affirmer leur présence dans la donnée (iii) établir une limite haute d'amplitude de l'OG pour sa détection dans le cas où on ne détecte rien.

Les principaux bruits, outre les erreurs systématiques et bruits des appareils de mesures (ici les téléscopes radio), présents dans la donnée de chronométrage sont intrinsèques aux pulsars et corrélés dans le temps. Nous avons des bruits dits "rouges", dominants en basse fréquence, dus à une irrégularité dans la rotation même des pulsars et induisant de légères variations dans les temps d'arrivée des pulses chronométrés sur Terre. Nous avons aussi un bruit de mesure de dispersion, du à la propagation des photons composants les pulses dans le milieu interstellaire. Durant son voyage depuis le pulsar jusqu'à la Terre, un pulse de lumière interagit avec les électrons interstellaire, ce qui a pour effet de retarder son temps d'arrivée. Lorsque la densité d'électrons varie sur l'axe pulsar-Terre, cela induit une perturbation dans le temps d'arrivée des pulses chronométrés. La présence de ces bruits détériore notre sensibilité aux OGs et ils doivent être modélisés avec soin. Enfin, nous présentons comment ces modèles sont testés grâce à des méthodes de Monte-Carlo par chaînes de Markov (MCMC) dont le fonctionnement est expliqué.

Le troisième chapitre présente les résultats obtenus avec le second ensemble de données de la collaboration PTA internationale (IPTA DR2) qui combine les données de : l'observatoire d'OG nanohertz d'Amérique du Nord (NANOGrav), le PTA Européen (EPTA) et le Parkes PTA (PPTA) en Australie. Cet ensemble de donnée étant assez conséquent, les ressources en calcul nécessaires pour l'analyser seront elles aussi condisérables. Pour palier ce problème, il est possible de classer les pulsars en fonction de leur sensibilité moyenne à une OGC de fréquence donnée. Pour ce faire, de fausses OGCs sont simulées dans la donnée de chronométrage pour estimer le rapport signal sur bruit de chaque pulsar. Ce rapport signal sur bruit quantifie la réponse individuelle des pulsars à l'OGC injectée et leur classement se fait suivant les valeurs obtenues. Il apparaît alors que certains pulsars dominent complètement le réseau de par leur sensibilité supérieure. En sélectionnant un sous-ensemble de pulsars qui complètent au moins 95% de la sensibilité totale du réseau, nous réduisons la donnée nécessaire aux analyses, et donc le temps de calcul, tout en conservant une bonne sensibilité aux OGCs. Après avoir classé les pulsars, nous appliquons les méthodes présentées dans le chapitre précédent pour tester la présence d'OGCs dans IPTA DR2. Nous n'avons trouvé aucune preuve de leur présence.

En conséquence, nous avons établi une limite haute de détection pour l'amplitude des OGCs et avons trouvé que IPTA DR2 donne la meilleure sensibilité à ce jour pour des signaux dont la fréquence est plus haute que 11 nHz. Dernièrement, les collaborations PTA ont mis en évidence la présence d'un bruit "rouge" commun à tous les pulsars (ayant la même densité spectrale de puissance entre chaque pulsar). L'origine de ce bruit est encore à déterminer. Son inclusion dans notre modèle détériore la sensibilité aux OGCs du réseau en basse fréquence. De plus, nous montrons qu'une mauvaise modélisation des bruits intrisèques aux pulsars peut mener à une fausse détection. En général, les bruits sont décrits de la même façon pour tous les pulsars. Cependant, des études récentes ont montré qu'adapter nos modèles de bruits individuellement pour chaque pulsar donnait une meilleure représentation de la donnée. Les meilleurs modèles de bruits sont sélectionnés par calcul du facteur de Bayes (introduit dans le second chapitre). Avant l'utilisation des modèles "sur-mesure", nous détections un signal similaire à une OGC de 51nHz, supporté par le facteur de Bayes. L'introduction des nouveaux modèles a complètement absorbé cette composante. L'origine de ce signal est encore à détérminer mais il est probable qu'il résulte d'interactions complexes entre différentes éléments de notre modèle. C'est pourquoi une bonne modélisation des modèles de bruits dans PTA est cruciale à la recherche d'OGs car des erreurs dans celui-ci peuvent mener à des fausses dé- En résumé, ce manuscript explore différents aspects de l'analyse de donnée pour la recherche de signaux gravitationnels. Il est principalement axé sur la donnée PTA (en particulier la donnée de la collaboration internationale IPTA) et la recherche d'OGC. Ces travaux ont mené à l'écriture de 3 articles scientifiques (dont 2 en tant que premier auteur) qui sont, alors que j'écris ces lignes, en cours de soumission à des revues scientifiques. Pour le moment, aucune détection d'OGC n'est à signaler dans la donnée IPTA. L'observation des pulsars et leur chronométrage est constamment en cours, donnant accès à toujours plus de donnée. De plus, la construction de nouveaux téléscopes radio haute précision (comme MeerKAT ou SKA) devrait augmenter notre sensibilité aux OGs, faisant de la première détection d'un signal gravitationnel en provenance d'un système binaire de trous noirs supermassifs un espoir bien concret pour les années à venir. The first chapter gives an overview of the theoretical background necessary to understand the origins of GWs. After presenting the basics of General Relativity, we argue that a Pulsar Timing Array (PTA) could detect GWs generated by Super Massive Black Hole Binaries (SMBHBs) and derive a mathematical expression for the expected timing residuals induced in the pulsar timing data. When such a signal is stable and monochromatic, we call it a Continuous Gravitational Wave (CGW). In the second chapter, we present data analysis methods for the detection of GW signals in PTA. We briefly explain Bayesian inference and detail how the noises intrinsic to pulsars can be modeled. We show that to assess the significance of our noise and signal models, we compute the Bayes Factor (BF) using hypermodels and Markov Chain Monte-Carlo (MCMC). In chapter three, we discuss the results obtained with the international PTA second data release (IPTA DR2) combining data from : the North American Nanohertz Observatory for Gravitational Waves (NANOGrav), the European Pulsar Timing Array (EPTA) and the Parkes Pulsar Timing Array (PPTA). We have searched for CGW signals and found no evidence of their presence in the data. We have set an upper limit on the amplitude of CGW and found that IPTA DR2 gives the best sensitivity to date for signals with frequencies above 11 nHz. It has been recently discovered that a low frequency noise common to all pulsars was present in the data. The inclusion of the latter in the model deteriorates the sensitivity of the array to CGW. Moreover, we show that a bad modeling of the intrinsic noise properties of pulsars can lead to false detection. The last chapter presents a new method based on Kernel Density Estimation (KDE) that was developed during my PhD to optimize the data analysis pipelines for GW search. The method is used with MCMC to improve the effective sampling time. It was tested with IPTA DR2 for CGW search and the simulated Laser Interferometer Space Antenna (LISA) dataset "Sangria" for Galactic Binaries (GBs) search. We show that the method is efficient for IPTA DR2 but is less performing with LISA.

Keywords

General relativity

Albert Einstein's General relativity (GR) came out in 1915 as a new theory of gravity and a generalization of special relativity. It no longer describes gravity as an attractive force but rather as a manifestation of the curvature of space-time that is modified in the presence of masses, where free-falling objects are actually inertial and following geodesic trajectories. The latter is a consequence of the equivalence principle saying that no distinction can be made between the motion of a falling object on a planet or in an accelerated frame of reference. The gravitational mass is the same as the inertial mass. Therefore, a free-falling object experiencing weightlessness is in an inertial frame of reference. Special relativity already tells us that a particle in motion will experience time dilation and length contraction in the tangent direction to its trajectory, which was raising contradictions in the case of circular motions and orbits. One famous Gedankenexperiment that is behind the intuition of curvature in GR is called the Ehrenfest paradox : a rotating disk of radius R would, due to length contraction, see its circumference smaller than 2πR.

C = 2πr < 2πR ⃗ ω R C = 2πr < 2πR r R ⃗ ω
But its radius R, in the direction perpendicular to its trajectory, would remain unchanged. This is only possible in curved space (see figure 1.1). Planets orbiting stars should then be experiencing a curvature of space while following an inertial trajectory. It is possible if heavy masses curve space-time and orbits are geodesics in a curved space.

After years of development, the so called Einstein field equations were derived to express the equivalence between curvature of space-time and distribution of energy. To quote John Archibald Wheeler : "Spacetime tells matter how to move; matter tells spacetime how to curve". The introduction of curvature in the equations of physics led to countless predictions that were all verified experimentally like the anomalous precession of Mercury, gravitational lensing, the expansion of the universe, the existence of black holes, and more recently, gravitational waves, the main topic of this thesis.

The notion of curved space

From simple 1D curve to more elaborate n-D hypersurfaces, the curvature is a measurement of how much we deviate from a straight line or a flat surface. In curved space, many of the general results of flat euclidian space geometry are no longer true. For instance, the sum of the 3 angles of a triangle may not equal 180 o or two locally parallel straight lines (geodesics) may eventually cross.

Carl Friedrich Gauss worked extensively on the matter and derived a scalar quantity, the Gaussian curvature, giving a mathematical measurement of how curved is a surface. He realized that for 2D surfaces, this quantity was invariant under isometric transformations and independent of the 3D space in which it was embedded. This observation is stated in the Theorema Egregium.

The curvature is an intrinsic property of the space. A curved space cannot be folded back into a flat space or the other way around. For example, a world map, as a 2D flat representation of the spherical surface of the Earth, does not accurately reproduce distances and areas. There exists no direct mapping of a curved surface on a flat one that comes without scale distortion issues. The sphericity of the Earth implies positive curvature that is not shared by a flat sheet of paper of null curvature. To that extent, length and angle measurements are fundamentally different in the two.

A pupil of Gauss, Bernhard Riemann, continuing the work previously undertaken by his teacher, generalized the mathematical description of the curvature for N-dimensional spaces. In GR, we make use of the Riemannian geometry to characterize the curvature of space-time.

Mathematical description of a curved space

Historically, different approaches were used to mathematically describe the curvature.

The first and maybe simplest is to define a tangent osculating circle at any point of a given parametrized curve s(t ) (see figure). The curvature κ is defined as the inverse of the radius of the circle. A straight line has an osculating circle of infinite radius, hence no curvature.

This concept can easily be generalized in differential geometry for oriented surfaces of higher dimensions. However it requires the d -dimensional surface to be embedded in a

d + 1-dimensional space.
Riemannian geometry allows the description of the curvature without immersing the surface in a d + 1 space, following the Theorema Egregium defining it as an intrisic property of the surface. It is done using the concepts of geodesic deviation (see figure 1.5).

Here we will present the standard "toolbox" that is used in GR for calculation in curved space. One can refer to references for more details (see [START_REF] Thorne | Gravitation Freeman[END_REF][START_REF] Schutz | A first course in general relativity[END_REF]; Hausner 2016).

Basis vectors

Basis vectors are tangent vectors at each point of the parameterized surface f (x µ ), we have :

e µ = ∂ µ f . (1.1)
They define the basis of a local tangent (and flat) space at point x µ . We use the notation :

∂ µ ≡ ∂ ∂x µ , (1.2)
to express the derivative with respect to coordinate x µ .

Metric tensor

The metric tensor is the mathematical object that redefines scalar product hence distances/angle measurement. In Euclidian flat space, it is just the identity matrix, but more generally we have :

g µν = e µ • e ν , (1.3)
that is not necessarily diagonal or the identity. It quantifies length and angle distortions at every point in space.

Length measurement Length measurement depends on the metric tensor. We define the element on length of a parametrized curve with components X µ as :

d s 2 = g µν d X µ d X ν , (1.4)
which defines the scalar product of infinitesimal vector d X with itself.

Christoffel symbols

The Christoffel symbols Γ β µα are related to the derivative in the µ direction of basis vector e α as :

∂ µ e α = Γ β µα e β , (1.5)
They are the coefficients that quantify how basis e α varies with respect to coordinate µ. They can be expressed in terms of the derivatives of the metrix tensor g µν .

Γ β µα = 1 2 g βλ ∂ α g λµ + ∂ µ g λα -∂ λ g µα (1.6)

Covariant derivative

The covariant derivative is the derivative of a vector v taking into account the variations of the basis vectors e α :

∇ µ v = ∇ µ v α e α = ∂ µ v β + v α Γ β µα e β .
(1.7)

Geodesic equation

The geodesic equation is the equation of a straight line1 in curved space. Physically, it corresponds to the equation of motion of a non-accelerated object in a curved space.

d 2 X µ d λ 2 + Γ µ αβ d X α d λ d X β d λ = 0, (1.8)
with λ the parameter of the curve (corresponding to the proper time τ for massive objects).

Riemann tensor

The Riemann tensor is a measurement of geodesic deviation. Indeed, in flat space, two neighbouring geodesics are parallel straight lines but this is no longer true in curved space. We can quantify this effect by measuring the difference in variation of the basis vector e α on a closed loop in two directions µ and ν (figure 1.5). Using equation 1.5 we get :

R β µνα e β = ∂ µ ∂ ν -∂ ν ∂ µ e α = ∂ µ Γ β να -∂ ν Γ β µα + Γ λ να Γ β µλ -Γ λ µα Γ β νλ e β . (1.9)
This definition can be generalized to any vector using parallel transport and covariant derivatives [START_REF] Hausner | General Relativity and Differential Geometry[END_REF]). The Riemann tensor appears when calculating the evolution of the separation vector between two objects free-falling side by side. Physically, it corresponds to the tidal forces.

Ricci tensor

The Ricci tensor is obtained from the Riemann tensor R β µνα by contraction of first index β and third index ν2 . It is analog to taking the trace of the matrix containing all possible deviation vectors R associated with e α being first transported in the µ direction (see figure 1.5), quantifying the strain felt in the directions perpendicular to µ.

R µα = R β µβα (1.10)
It contains information about volume deformation along geodesics. Component R µα is a measurement of volume change when transporting basis vector e α along direction µ. A null Ricci tensor does not necessarily mean that there is no curvature. It just corresponds to volume conservation along geodesics (i.e. experienced contraction and dilation in directions perpendicular to µ compensate). This tensor is symmetric and is more often noted R µν .

Einstein field equations

The Einstein field equations are the central equations of GR. They tell how space-time geometry is affected by the distribution of energy within it. Finding a solution to this equation gives the metric tensor g µν that we will use to study the dynamics of objects using the geodesic equation. The distribution of energy and matter is described by the stress energy tensor T µν . It follows the continuity relation :

∇ µ T µν = 0, (1.11)
expressing the local conservation of energy and momentum.

At first sight, it could be tempting to write R µν = κT µν to map the curvature of spacetime to the presence of energy. But this cannot work because R µν does not have a null covariant derivative. Nonetheless, it possible to build a curvature tensor from the Ricci tensor R µν that also has null covariant derivative. We call it the Einstein tensor G µν :

G µν = R µν - 1 2 R g µν , (1.12)
∇ µ G µν = 0, (1.13)
We can now empirically postulate an equation that relates the stress-energy tensor T µν and the Einstein tensor G µν as :

G µν = κT µν . (1.14)
The value of the constant κ is set in such a way that equation 1.14 reduces to Poisson's equation for the gravitational potential in the Newtonian limit. Including the cosmological constant Λ, we have the final expression for the Einstein field equations :

G µν + Λg µν = 8πG c 4 T µν . (1.15)
For a given T µν , this equation gives a system of non linear second order differential equations for the metric tensor g µν . It is possible to linearize this equation by introducing small perturbations around the flat-space metric.

Gravitational waves

The following section is a mathematical introduction to gravitational waves inspired from the present references : Maggiore 2007; Babak 2020; Carroll 1997.

Linearized Einstein field equations

We define a metric tensor that is :

g µν = η µν + h µν , (1.16)
with inverse :

g µν = η µν -h µν , (1.17)
where h µν << η µν , is a small perturbation we add to the flat-space Minkowski metric η µν .

To express the Einstein field equations using this approximation, we first need to write the Ricci tensor R µν . Keeping only first order terms in h µν and getting rid of higher order terms. We have for the Christoffel symbols :

Γ β µα = 1 2 η βλ ∂ α h λµ + ∂ µ h λα -∂ λ h µα + O(h 2 ), (1.18)
which gives for the Ricci tensor :

R µν = ∂ λ Γ λ µν -∂ µ Γ λ λν + O(Γ 2 ) = 1 2 ∂ β ∂ ν h βµ + ∂ β ∂ µ h βν -∂ β ∂ β h µν -∂ µ ∂ ν h λ λ , (1.19)
and the Ricci scalar :

R = ∂ β ∂ λ h βλ -∂ β ∂ β h λ λ (1.20)
We can now write the Einstein tensor :

G µν = 1 2 ∂ β ∂ ν h βµ + ∂ β ∂ µ h βν -h µν -∂ µ ∂ ν h -η µν ∂ β ∂ λ h βλ -h . (1.21)
Where we have used to simplify notation convention :

h λ λ = h and ∂ β ∂ β = (1.22)
The physics of GR do not depend on the choice of coordinate system. That is to say, for any local coordinate transformation :

x µ → x µ + ξ µ (x), (1.23)
the Einstein field equations will remain unchanged. For such a coordinate transformation, the metric perturbation h µν rewrites as :

h µν → h µν + ∂ ν ξ µ + ∂ µ ξ ν . (1.24)
One can verify that the additional terms actually cancel out in equation 1. [START_REF] Corbin | Pulsar Timing Array Observations of Massive Black Hole Binaries[END_REF]. That means we have a gauge (coordinate) freedom which we can use to simplify the expression.

We introduce the trace reversed metric 3 that is

hµν = h µν - 1 2 η µν h, (1.25)
for which the gauge condition becomes :

hµν → hµν + ∂ ν ξ µ + ∂ µ ξ ν -η µν ∂ β ξ β , (1.26)
We can rewrite the linearised Einstein tensor in term of this trace reversed metric, yielding

G µν = 1 2 -hµν + ∂ β ∂ ν hµβ + ∂ β ∂ µ hνβ - 1 2 η µν ∂ β ∂ λ hβλ . (1.27)
Using the coordinate freedom expressed in (1.26) we can now derive the divergence of hµν as

∂ µ hµν → ∂ µ hµν + ξ ν , (1.28)
then choose the ξ µ in such a way so that they obey ξ ν = -∂ µ hµν to give the Lorenz gauge condition :

∂ µ hµν = 0. (1.29)
Using that gauge, the expression of the Einstein tensor enormously simplifies to become :

G µν = - 1 2 hµν (1.30)
giving the final expression of the linearized Einstein field equations for the trace re-

versed metric : hµν = - 16πG c 4 T µν (1.31)
We can already see that this equation corresponds to the equation of propagation of a wave with a source term. Finding solutions to (1.31) give periodic first order perturbations of the metric tensor that we identify as gravitational waves.

Polarization tensor

When considering gravitational waves propagating in the vacuum, we can set the stressenergy tensor to be T µν = 0. In that case, equation (1.31) takes the simple form :

hµν = 0. (1.32)
Solutions to the latter can be expressed as a linear combination of plane waves with polarization tensor µν : hµν = µν e i k σ x σ , (

with k µ the wave vector.

This expression of the linearized Einstein field equations in the vacuum offer additional coordinate freedom. Indeed, plugging equation (1.26) in (1.32) we have :

hµν + ∂ ν ξ µ + ∂ µ ξ ν -η µν ∂ β ξ β = 0, (1.34)
that is satisfied if the ξ µ follow the wave equation ξ µ = 0. Solving ξ µ = C µ e i k σ x σ and using the plane wave solution of hµν , the gauge condition (1.26) becomes :

hµν = µν e i k σ x σ → µν + i k ν C µ + k µ C ν -η µν k β C β e i k σ x σ , (1.35)
and the trace of the transformed polarization tensor gives :

µ µ → µ µ -2i k µ C µ (1.36)
There are infinite number of ways to pick the C µ without affecting the Lorenz gauge 4 .

We can fix the C µ (hence the ξ µ ) so they cancel the trace of hµν . This choice of coordinate system is called the Transverse Traceless (TT) gauge. If the trace is null, one immediately sees that we have :

hµν = h µν (1.37)
So plane wave solution of hµν are the same for the initial metric perturbation h µν .

Given equation (1.32), the harmonic gauge and the TT gauge, we have :

k µ k µ = 0 k µ µν = 0 µ µ = 0 (1.38)
The first condition shows that k µ is light-like so the wave is propagating at the speed of light c, the second condition expresses the transversality of the wave and the last condition gives the null trace of the polarization tensor. Therefore, for a wave propagating along k 3 with wave vector k µ = (k 0 , 0, 0, k 3 ), we have :

(k 0 ) 2 -(k 3 ) 2 = 0 0ν = 3ν = 0 11 + 22 = 0.
(1.39)

Because h µν is symmetric, we also have for the two last unconstrained components of µν : 12 = 21 .

(1.40)

Finally, we define 11 = h + and 12 = h × to get the polarization tensor of a GW :

µν =       0 0 0 0 0 h + h × 0 0 h × -h + 0 0 0 0 0       (1.41)
This tensor can be decomposed in two independent components as :

µν = + µν + × µν =       0 0 0 0 0 h + 0 0 0 0 -h + 0 0 0 0 0       +       0 0 0 0 0 0 h × 0 0 h × 0 0 0 0 0 0       , (1.42)
that we call the "plus" and "cross" polarizations.

Quadrupole formula

In the previous section we have presented a solution to the linearized Einstein field equations in the vacuum but we did not explain how GWs are generated. Just like the electromagnetic waves originate from electric charges in motion, we expect GWs to be produced by moving masses. However, there is a fundamental difference with GWs. Take a spherically symmetric gravitational potential Φ(r ) ∝ M /r with M = i m i the total mass of the source object and write its multipole expansion :

Φ(r -R) ∝ 1 R i m i - 1 R 3 i m i r i • R + ..., (1.43)
In virtue of the momentum p conservation principle, the massive dipole i m i r i is a constant. Indeed, we have in the center of mass frame :

i p i = i d d t m i r i = 0, (1.44)
Therefore, it is impossible to generate time-dependent GWs with a massive dipole.

The source must have asymmetric (at least quadrupolar) features to be able to radiate GWs. In Einstein field equations, this condition is satisfied by ∇ µ T µν = 0.

To solve equation (1.32) we make use of Green's function formalism. The solution for the wave equation has already been studied and we can find in the literature the expression of the Green function for the wave operator : .45) Using this result it comes naturally that the h µν will be expressed in terms of the Green function as :

G(t , x, t , x') = - δ(t -[t -[x -x']]) 4π|x -x'| . ( 1
hµν = 4G c 4 dx' T µν (t -[t -[x -x']]) |x -x'| . (1.46)
The distance to the source is usually very large. To that extent, |xx'| becomes the cosmological luminosity distance d L , taking into account the expansion of the universe for high redshift objects. Moreover, in the TT gauge, we only have 2 independent spatial components of h i j . Using ∂ µ T µν = 0 and T 00 >> T i j in the non-relativistic limit, it is possible to show after some calculation that (see Babak 2020) :

∂ 2 0 dx(T 00 x i x j ) = 2 dxT i j (1.47)
which is only valid for isolated sources and leads to :

hi j = 2G d L c 4 ∂ 2 0 dx'T 00 x i x j , (1.48)
where T 00 = ρc 2 is the rest energy density. Defining the quadrupole moment tensor as :

I i j = dx'ρx i x j , (1.49)
and using its expression in the TT gauge :

I i j = I i j - 1 3 δ i j trace(I ) (1.50)
we have the final expression of our quadrupole formula for the generation of GWs :

h i j = 2G d L c 4 Ïi j . (1.51)
The prefactor being extremely small, GWs have very low amplitude. The only way to produce loud and detectable GWs is to have a very massive source. Furthermore, as expected according to the initial comment in this section, the lowest possible order of moment is the quadrupole because of the conservation of energy and momentum. Gravitational dipoles are forbidden by the equations of GR. It means that for GWs to be produced, a system must present asymmetries in the distribution of its mass. We can find such asymmetries in orbiting binary systems.

Circular binary system

Consider a binary system of reduced mass µ and total mass M = m a + m b in circular orbit in the x-y plane at orbital frequency ω with separation r . We have for the positions of the two objects x a and x b in the center of mass frame :

x a (t ) = r µ m a cos(ωt ), sin(ωt ), 0

x b (t ) = - m a m b
x a (t ).

(1.52)

For two point objects, the mass distribution ρ( x) is :

ρ( x) = m a δ( x -x a ) + m b δ( x -x b ). (1.53)
We can now calculate the components of quadrupole moment I i j using equations (1.49) and (1.50) which yield :

I i j = 1 2 µr 2      1/3 + cos(2ωt ) sin(2ωt ) 0 sin(2ωt ) 1/3 -cos(2ωt ) 0 0 0 0      , (1.54)
and then take the second time derivative of I i j :

Ïi j = -2µr 2 ω 2     cos(2ωt ) sin(2ωt ) 0 sin(2ωt ) -cos(2ωt ) 0 0 0 0     .
(1.55)

We use Kepler's third law to constrain the value of r :

r 3 ω 2 = G M ⇔r = G M ω 2 1 3 , (1.56)
now using (1.51) we have for the two independent components h + and h × :

h + = - 4M 5/3 c G 5/3 ω 2/3 d L c 4 cos(2ωt ) h × = - 4M 5/3 c G 5/3 ω 2/3 d L c 4 sin(2ωt ), (1.57)
where we have introduced the chirp mass M c = µ 3/5 M 2/5 . We can see that the GW signal oscillates at twice the orbital frequency. To that extent, we have :

ω = π f g w (1.58)
with f g w the GW signal frequency. Like the luminosity distance d L , the frequency of the signal can be subject to redshift when the binary system is far away. Therefore, for redhsift z, it transforms as :

f g w → f g w 1 + z (1.59)
A gravitational binary system will radiate energy in the form of GWs. As a consequence, the orbital frequency ω will increase with time and inversely, the separation r will decrease until the two objects eventually merge. We can estimate this frequency and separation evolution if we know the amount of energy radiated by the system per unit time d E /d t given by [START_REF] Carroll | Lecture Notes on General Relativity[END_REF][START_REF] Schutz | A first course in general relativity[END_REF] :

d E d t = G 5c 5 〈 ... I i j ... I i j 〉, (1.60)
where ... I i j is the third time derivative of the quadrupole and 〈.〉 denote the averaging over some periods of the wave. After some straightforward calculation, we have :

d E d t = - 32 5c 5 G 7/3 M 4/3 µ 2 ω 10/3 = - 32G 4 µ 2 M 3 5c 5 r 5 , (1.61)
writing ω in terms of r using Kepler's third law. Making use of the Virial theorem for a gravitational system, we also have :

E = - 1 2 G M µ r ⇔ d E d t = 1 2 G M µ r 2 d r d t . (1.62)
Equating the two last expressions, we get an equation for r that we can solve. Once solved, it gives :

r (t ) = r 4 0 - 256 5 G 3 M 2 µ c 5 t 1/4 = r 0 1 - t τ c 1/4 , (1.63)
where we have defined r 0 the initial separation and τ c the coalescence time. At t = τ c , separation r (τ c ) = 0, thus the two orbiting objects come into contact and merge 5 . Again, using Kepler's third law and the chirp mass M c , we can write the evolution of the orbital frequency ω(t ) :

ω(t ) = ω 0 1 - 256 5
G 5/3 M 5/3 c c 5 ω 8/3 0 t -3/8 = ω 0 1 - t τ c -3/8 , (1.64)
where ω 0 is the initial orbital frequency, showing that indeed, ω(t ) goes up with time, hence, so does the GW frequency f g w (t ). It means that the phase Φ(t ) of the GW will evolve with time like :

Φ(t ) = Φ 0 + t t 0 d t ω(t ) = Φ 0 + c 5 32G 5/3 M 5/3 ω -5/3 0 -ω(t ) -5/3 (1.65)
In equation (1.57) we give the h + and h × components of h µν for a binary system orbiting in the x-y plane with angular momentum L, emitting GW in the z direction. For a more general description, we need to rotate around the x-axis the transverse quadrupole moment by an inclination angle ι giving :

h + = 4M 5/3 c G 5/3 ω 2/3 d L c 4 (1 + cos 2 ι) cos(2Φ(t )) h × = 4M 5/3 c G 5/3 ω 2/3 d L c 4 cos ι sin(2Φ(t )).
(1.66)

Finally, we can perform a rotation in the polarization plane by an angle 2ψ, introducing the polarization basis following ê1 • ê2 = 0 and ê1 • k = ê2 • k = 0. We have for e A i j the + and × GW polarization tensor :

+ i j = ê1 ⊗ ê1 -ê2 ⊗ ê2 × i j = ê1 ⊗ ê2 + ê1 ⊗ ê2 , (1.67)
h + h × → cos 2ψ -sin 2ψ sin 2ψ cos 2ψ h + h × , (1.68)
and we have h µν in its final form :

h µν = A=+,× A µν h A = + µν h + + × µν h × h + = 2M 5/3 c G 5/3 ω 2/3 d L c 4 (1 + cos 2 ι) cos 2ψ cos(2Φ(t )) -2 cos ι sin 2ψ sin(2Φ(t )) h × = 2M 5/3 c G 5/3 ω 2/3 d L c 4 (1 + cos 2 ι) sin 2ψ cos(2Φ(t )) + 2 cos ι cos 2ψ sin(2Φ(t )) Φ(t ) = Φ 0 + c 5 32G 5/3 M 5/3 ω -5/3 0 -ω(t ) -5/3 ω(t ) = π f g w (t ) = ω 0 1 - 256 5
G 5/3 M 5/3 c c 5 ω 8/3 0 t -3/8
(1.69)

In summary, we have 5 parameters that characterize our GW :

Supermassive black holes

In 1916, Karl Schwarzschild came up with a metric tensor g µν solution to the Einstein field equations in the vaccum (see Schwarzschild 1999, a translation to English of the original 1916 paper), describing the geometry of spacetime in spherical coordinates when a nonrotating point mass M is present in the center of the frame. This solution is consistent with classical gravity in the Newtonian limit but it predicts peculiar behaviour when approaching the center. Using that metric the length element d s 2 writes as :

d s 2 = 1 - 2G M r c 2 c 2 d t 2 -1 - 2G M r c 2 -1 d r 2 -r 2 d Ω 2 , (1.70)
where r is the distance to the center for an external (far away) observer and d Ω = sin 2 θd θd φ is the element of solid angle. One immediately sees that for r = 2G M /c 2 , the coordinate system diverges. This radius of singularity6 is called the Schwarzschild radius and it defines the boundaries of what we call the event horizon. Beyond this limit, nothing can escape and this region of spacetime appears completely black to external observers :

we have a Black Hole (BH).

These astrophysical objects are completely characterized by their mass that can range from a few solar masses M to billions of solar masses. If it is greater than 10 6 M , we call it a supermassive black hole (SMBH). They are the largest category of black holes and we expect to find them at the center of galaxies. How they form is still an open question in astrophysics. The observation of Active Galactic Nuclei (AGN) at very high redshift suggests that they were already present in the early Universe. The proposed scenarios for their formation rely on the production of black hole "seeds" that already have masses around 10 The most accepted scenario to explain the birth of such SMBH pairs is by means of galaxy merger. When their nuclei get closer, the stars and gas making up the galaxies steal kinetic energy from the SMBHs through gravitational interaction. The stars get ejected by slingshot effect and the gas generate dynamical friction. This loss of energy will drive the SMBHs closer and closer until they form a stable binary, separated by a distance that is only in the order of some parsecs. At that stage, there might not be enough matter present between the two SMBHs to significantly extract energy from the system. The only way it can give away energy is through GW emission that are extremely weak. It is estimated that such an inspiralling phase should be as long as the age of the Universe before entering the merging phase where the GW emissions become efficient. The observations seem to contradict this prediction because we see very massive AGNs at high redshift whose formation we can explain only with SMBHB merger. This is the "final parsec" problem.

The mechanisms allowing the SMBHs to merge within a reasonable amount of time are still unclear but some solutions have been proposed. In Milosavljevic 2003 they describe two of the suggested explanations : (i) they show that taking into account the shape of the galaxies and the inhomogeneous distribution of matter within them could make the dynamical friction more efficient, allowing the SMBHs to enter merging phase, (ii) they explain that a N -SMBH (galaxy) merger could solve to the issue, where the inspiralling phase is described as a N -body problem which could allow the SMBHs to get closer under certain conditions.

In practice, it is possible to detect such GWs thanks to the huge mass of the SMBHs (i.e.

huge chirp mass M c of the binary system) that make up for loud GWs. A large population of SMBHB would produce a superposition of many GWs filling-up the entire Universe 7 .

The resulting signal can de divided into two main categories : the stochastic Gravitational Wave Background (GWB) and Continuous Waves (CW).

Gravitational wave background

The vast population of observable galaxies indicates that we probably have a large population of SMBH in the Universe. The superposition of many individual SMBHB GW signals would result in a low frequency noise spread on a broad band. For a population of circular SMBHB, given the formula for GW emissions, we can show that such a noise signal has a frequency spectrum following a powerlaw. We define the characteristic strain of the GWB as (Taylor 2021):

h 2 c ( f ) ≡ f S h ( f ) ∝ f -2 d ρ d ln f , ( 1.71)
with S h ( f ) the power spectral density for the considered population of sources (here circular SMBHBs) and d ρ/d ln f the energy density of sources in a logarithmic interval of 7 We can only see the local Universe. Considering that we have N 0 emitting binaries and ignoring the effects due to redshift 8 , we can approximate the density d ρ/d ln f for circular binaries using equations (1.61) and (1.76) :

d ρ d ln f ∼ N 0 d E d ln f ∝ f d E d t d t d f ∝ f × f 10/3 × f -11/3 ∝ f 2/3 , (1.72)
yielding a frequency dependence following a powerlaw with spectral index 2/3.

Given equation 1.71, we have a h c that follows a powerlaw as :

h c ( f ) ∝ f -2/3 . (1.73)
On figure 1.9, we see a simulated population of individual SMBHBs with their corresponding GW strain amplitude h and frequency f . The higher frequency range appears less crowded because the frequency evolution rate of the circular SMBHB increases as

f 11/3 (equation 1.76
) so the binaries spend less time there. As a consequence, the characteristic strain h c decreases with frequency, like the GW energy density, and the signal looses its stochasticity.

When we search for such stochastic signal, we model the characteristic strain as a powerlaw with amplitude A GW B and spectral index γ GW B within our prior assumptions on their values :

h c = A GW B f 1yr -1 γ GW B
.

(1.74)

The detection of the GWB and the measurement of (A GW B , γ GW B ) would give crucial information about the population and distribution of SMBHBs across the Universe. It is the main target of PTA collaborations.

Continuous waves

The amplitude of a GW decreases with distance as 1/d L . Some of the SMBHB responsible for the stochastic GWB may be closer to us and appear brighter in their GW emission, resulting in GW signal resolved at single frequencies on top of the GWB (see figure 1.9).

These GWs are called Continuous Gravitational Waves (CGW) and are described as deterministic signals. According to equation (1.64) showing the evolution of GW frequency with time, CGWs should see their frequency increase. However, for supermassive systems, this evolution can be neglected over some orbital periods (for frequencies of ≤ 10 (red) by fitting a timing model. On the right panels, the same thing but in the presence of GW signal. The timing model does not account for the delay induced by the GW and timing residuals δ t are measured.

Hz only, see section 1.4.1 for more details). This manuscript focuses on this type of signal and we will present in the next section how it is possible to detect them using pulsars.

Pulsar Timing Array

A Pulsar Timing Array (PTA) is an array of Milliecond Pulsars (MSPs) observed from the Earth using radiotelescopes to monitor the time of arrivals (TOAs) t of their pulses. Among these pulsars, some have been observed for long periods of time T >20 years giving us robust datasets to work with. This timing data allows to use pulsars as "cosmic clocks". We can predict the TOAs by characterizing the properties of the pulsar and modeling perturbations of physical origin. The difference between the predicted TOAs t pr ed and the measured TOAs t obs is called the timing residuals δ t . If our PTA is sensitive enough, the effect of GWs should be present and detectable in the timing residuals δ t .

MSPs usually have a characteristic and well defined pulse profile. However, the observation of a single pulse can be challenging because of its low amplitude and the ambient noise levels. For that reason, we measure series of consecutive pulses that we stack to-gether to obtain an averaged pulse profile with a significantly better signal to noise ratio.

We then cross-correlate a previously determined template of the pulse profile 9 to extract a value of the TOA. One new TOA is obtained with a cadence of about 1 week. We can estimate the uncertainty of a TOA measurement as (Lorimer and Kramer 2012) :

σ T O A W S/N , (1.75)
with W the width of the pulse and S/N the signal to noise ratio of the measured (averaged) pulse profile.

Typically, we have σ T O A 10 -6 s so we are accurate to the micro-second, which is remarkable for astrophysical observations. This timing precision led physicists to believe that the effect of GWs could be noticeable in the TOAs. This idea was first proposed in [START_REF] Sazhin | Opportunities for detecting ultralong gravitational waves[END_REF][START_REF] Detweiler | Pulsar timing measurements and the search for gravitational waves[END_REF] where it is argued that only GW sources with masses around 10 8-10 M could be detectable within the sensitivity σ T O A of our MSP timing measurements. The natural candidates for such GW sources are SMBHBs.

Sensitivity of the array

The long period of observation T of the PTA makes it sensitive to GWs down to the nanohertz band (1/T ∼ 10 -9 Hz). The cadence of observation is about 1 week 10 so signals with frequencies above 1/7days ∼ 10 -6 Hz are under sampled. Moreover, many observation gaps are present in the dataset, sometimes of several months due to maintenance of the radiotelescopes. As a consequence, the TOAs are not evenly sampled and it is not possible to directly visualize the Fourier spectrum of the pulsars. Instead, we use Gaussian processes to model the time-correlated features intrinsic to pulsars and interpolate between gaps (see section 2.2.2).

In that window of frequency, the GWB should give rise to measurable effects thanks to the high precision timing measurements of MSPs (see figure 1.9). For CGWs, we can argue whether the frequency evolution of the SMBHBs is going to introduce a bias in the measurement of the signal (Sesana and Vecchio 2010). We derive from equation (1.64) :

d ω d t = 3 8 1 τ c ω 11/3 = 96 5 G 5/3 M 5/3 c c 5 ω 11/3 , (1.76)
the expression giving the rate of change in GW frequency (recall that f g w = πω) for a binary system with orbital period ω. In the context of PTA, we are sensitive to frequencies around f g w ∼ 10 -8 Hz and binaries with a chirp mass M c ∼ 10 9 M . For such systems, the coalescence time τ c is in the order of 10 4-5 years so the change in frequency that occurs 9 For curiosity, many pulse profile templates can be found at : https://www.jb.man. within the period of observation T (in the order of 10 years) is expected to be small. Using the previous equation, we can evaluate this change of frequency ∆ f g w to first order :

∆ f g w d f g w d t × T 96 5
G 5/3 M 5/3 c c 5 f 11/3 g w × T. (1.77)
For the considered values of the parameters M c , f g w and T we have a ∆ f g w that is less than 10 -9 Hz, hence below PTA frequency resolution. To that extent, the CGW sources can be considered stable and monochromatic (i.e. non evolving).

This approximation is no longer true for signals with higher frequencies. The very strong dependence of ∆ f g w on f g w can produce differences that are resolvable by the array. On figure 1.11, we plot ∆ f g w (not approximated as in equation 1.77) as a function of the initial signal frequency f g w , for 10 years of observation and three different values of chirp mass M c , to compare it with a typical PTA resolution of 1 nHz. It is clear that discrepancies between true and measured signal frequencies will appear around 100 nHz (or less for large M c ) where ∆ f g w exceeds the PTA resolution. We even start to see very strong frequency evolution for the largest values of chirp mass. It means that the non-evolving source approximation no longer stands. To search for signal at frequencies higher than ∼100 nHz, we must take into account the evolution of the SMBHB. If we do not, 100 nHz should be a reasonable limit up to which the ∆ f g w can be negligible for M c ≤ 10 9 M .

Millisecond pulsars

Pulsars appear as the last state of evolution of stars with masses around 9M . When nuclear fusion cycles are over within its core, the produced iron core collapses under the gravitational pull producing a gigantic explosion known as a supernova of type II. During this event, the whole star shrinks down into a small volume, leaving behind an extremely compact and highly magnetized neutron star with a diameter of about 10km. The angular momentum of the star before the supernova is conserved and transferred to the tiny residual neutron star, giving it a high rotational speed. The combination of both high magnetization and high spin velocity produces beams of electromagnetic emissions that we perceive as pulses, hence the name, Pulsar (Pulsating staR). They were originally detected in the radio frequency band. Some of them are extremely stable in their rotation, which is why they are interesting objects to study, opening a whole spectrum of very precise astrophysical measurements.

Due to their emissions, pulsars loose energy with time. As a consequence, they loose angular momentum and their rotational speed decreases. This "spin-down" can cause the pulsar to fade away and not be detectable by radiotelescopes. If the pulsar is in a binary system, it can get re-accelerated by stealing angular momentum from its companion through gas accretion. This process is called the "recycling" of pulsars [START_REF] Alpar | A new class of radio pulsars[END_REF]).

Once equilibrium is reached, such pulsars can have periods of rotation in the order of the millisecond (i.e. millisecond pulsars). Moreover, due to their history and age, they are much more stable than young pulsars that just went through a supernova (more prone to glitches/starquakes).

MSPs in our galaxy are monitored by several collaborations throughout the world using their radio observatories to collect timing data. Lately, efforts were made to put in common all these datasets which gave birth to the international collaboration.

The International Pulsar Timing Array

The The combination of all data is superior to the datasets of each collaboration. We have already observed the improvement in the detection of the common red noise process in Antoniadis 2022 by using IPTA DR2. This upgrade in sensitivity can be explained by :

• a better sky coverage providing better localization of GW signals,

• a better decoupling and identification of noise components due to the increased number of backends giving better constraints on noise parameters,

• the reduction of the number of gaps in the data due to absence of observations.

Today, the IPTA also counts the Indian PTA (InPTA) collaboration [START_REF] Tarafdar | The Indian Pulsar Timing Array: First data release[END_REF] and the South-African MeerKAT radio-telescope (Bailes et al. 2020) within its collaborators. For now, the only available combination is the DR2 of which we make extensive use in this manuscript. The IPTA collaboration is already working on the third data release that will combine latest datasets to give the best possible sensitivity to date.

PTA response to GW signal

We have presented in section 1.2 a solution to the linearized Einstein field equations in the TT gauge that gave GWs. The first order perturbation to the flat-space metric for a GW propagating along the ẑ axis was :

h T T µν = A=+,× A µν h A =       0 0 0 0 0 h + h × 0 0 h × -h + 0 0 0 0 0       , ( 1.78)
where we see the two independent polarization states h + and h × .

For a pulsar in the sky, let n be the unit vector pointing from the Earth to the pulsar and σ µ the four-momentum of a light-like signal of frequency ν travelling from the pulsar to the Earth. We can write the geodesic equation for σ µ as :

d σ µ d λ + Γ µ αβ σ α σ β = 0, (1.79)
with

σ µ = (ν, -ν n) = d X µ d λ = d d λ (t , x). (1.80)
In theory, the GW should influence the position of the pulsar n but this effect is considered small and neglected here. Expressing the Christoffel symbols for g µν = η µν + h µν and keeping only the first order terms, we get :

Γ µ αβ = 1 2 g µδ ∂ α g βδ + ∂ β g δα -∂ δ g αβ 1 2 η µδ ∂ α h βδ + ∂ β h δα -∂ δ h αβ + O(h 2 ) (1.81)
We plug this expression of Γ µ αβ into the geodesic equation and focus on the time component of σ µ :

d σ 0 d λ + Γ 0 αβ σ α σ β = 0 ⇔ d ν d λ + 1 2 η 00 -∂ 0 h αβ σ α σ β = 0 ⇔ d ν d λ + 1 2 ∂h i j ∂t ν 2 ni n j = 0.
(1.82)

We only have spatial components in h µν so we denote the metric perturbation by h i j , omitting the h 0ν and h µ0 terms that are null. To account for the fact that GWs propagate at the speed of light, we have to define the delayed wave solution :

h i j = h i j (t -k • x), (1.83)
with k the direction of propagation of the GW. This gives for the total derivative of h i j with respect to λ:

d h i j d λ = ∂h i j ∂t d t d λ + ∂h i j ∂( k • x) d ( k • x) d λ = ∂h i j ∂t ν - ∂h i j ∂t k • d x d λ = ∂h i j ∂t ν - ∂h i j ∂t k • (-ν n) = ∂h i j ∂t ν 1 + k • n , (1.84)
and for the time derivative :

∂h i j ∂t = d h i j d λ 1 1 + k • n 1 ν . (1.85)
Now, combining (1.82) and (1.85), we get :

1 ν d ν d λ = - 1 2 ni n j 1 + k • n d h i j d λ . ( 1.86)
This equation shows how the frequency of a signal is affected by the presence of a GW.

In fact, if we consider a pulsar that is emitting pulses of light at rate ν 0 , we can integrate this last expression between the times of emission of the pulse at pulsar t p and reception t on Earth. This yields the redshift (or blueshift) induced in the TOAs of a pulsar :

ln ν(t ) ν 0 = - 1 2 A=+,× ni n j 1 + k • n A i j (h A (t ) -h A (t p )) (1.87)
with ν(t ) the frequency at reception that we can define as ν(t ) = ν 0 + δν for small variation of ν giving ln ν 0 + δν /ν 0 δν/ν 0 .

We can relate the time at the pulsar t p with the time on Earth t using the relation :

t p = t - L c (1 + k • n), (1.88)
where L is the distance Earth-pulsar. The right hand side term represents the delay due to the propagation of the pulse, taking into account the geometry of the system (see figure 1.14). We see that for k• n = -1 we have t p = t . In that particular case, the pulses and GW propagate on the same axis at same velocity c (the emitting SMBHB is right behind the pulsar), canceling the delay between the measured time at Earth and pulsar, hence the frequency shift induced in the TOAs. This is called the "surfing" effect. We can also calculate the scalar product between the geometrical prefactor function of n and the A i j using the definition of the polarization basis ( ê1 , ê2 ) and the expression of vectors in spherical coordinates :

k = -sin θ cos φ x -sin θ sin φ ŷ -cos θ ẑ ê1 = -sin φ x + cos φ ŷ ê2 = -cos θ cos φ x -cos θ sin φ ŷ + sin θ ẑ (1.89)
This gives the antenna pattern functions F A for both polarizations, expressing the response of the PTA as a function of the sky position n of the pulsar and the direction of propagation of the GW k. We can write them as :

F + ( k) = 1 2 ( ê1 • n) 2 -( ê2 • n) 2 1 + k • n , (1.90
)

F × ( k) = ( ê1 • n)( ê2 • n) 1 + k • n , (1.91)
allowing us to rewrite equation (1.87) :

δν ν 0 = - A=+,× F A ( k) h A (t ) -h A (t p ) (1.92)
To get the timing residuals δ t GW induced in the TOAs by the passing GW signal, one needs to integrate again δν/ν 0 to calculate the phase shift produced by this frequency modulation of the pulsar signal :

t 0 d t δν ν 0 = - A=+,× F A ( k) t 0 d t h A (t ) -h A (t p ) . (1.93)
This requires to integrate the h i j terms, which cannot be done analytically. However, it is possible to approximate this integral by neglecting the evolution of the orbital frequency ω(t ) of the SMBHB (see section 1.4.1) producing the GWs over the course of PTA observation (of about 20 years). Given equations in (1.69), for non polarized signal (ψ = 0)

we have :

h + ∝ ω(t ) 2/3 cos(2Φ(t )) h × ∝ ω(t ) 2/3 sin(2Φ(t )), (1.94)
and for the integrals, neglecting the evolution of ω(t ) and using the chain rule to relate Φ and t :

d t h + ∝ ω 2/3 d Φ d t -1 d Φ cos(2Φ(t )) d t h × ∝ ω 2/3 d Φ d t -1 d Φ sin(2Φ(t )), (1.95)
with d Φ/d t = ω, we finally get :

d t h + ∝ 1 2 ω -1/3 sin(2Φ(t )) d t h × ∝ - 1 2 ω -1/3 cos(2Φ(t )).
(1.96)

After re-rotating the polarization basis by an angle 2ψ, we can write the timing residuals δ t GW induced in the TOAs of a pulsar by an individual SMBHB emitting CWs :

δ t GW = A=+,× F A ( k) s A (t ) -s A (t p ) s + (t ) = M 5/3 c G 5/3 d L c 4 ω(t ) -1/3 (1 + cos 2 ι) cos 2ψ sin(2Φ(t )) + 2 cos ι sin 2ψ cos(2Φ(t )) s × (t ) = M 5/3 c G 5/3 d L c 4 ω(t ) -1/3 (1 + cos 2 ι) sin 2ψ sin(2Φ(t )) -2 cos ι cos 2ψ cos(2Φ(t )) t p = t - L c (1 + k • n).
(1.97)

Note that we have two terms acting in δ t GW , the Earth term s A (t ) and the pulsar term s A (t p ) (for both polarizations s + and s × expressed above). We see that the amplitude of the timing residuals depends on ω(t ) -1/3 . It means that the sensitivity of the array decreases for signals with higher frequencies where the array experiences mainly white noise. We can expect an improved sensitivity for the low frequency signals but it is also mitigated by the presence of time correlated noises (red noise, dipsersion measure noise) as well as our limited time of observation (see section 1.4.1). The negative effect of the ω(t

) -1/3
factor is reflected in the behaviour of the PTA sensitivity curve on figure 1.9 where we see its degradation as we go higher in frequency.

Even though we have neglected the evolution of ω(t ) to approximate the integral in (1.95), it is not possible to neglect its evolution over the course of the whole travel of a pulse from pulsar to Earth. In section 1.4.1, we have explained why it is a reasonable approximation to consider the CGW sources to be monochromatic and non evolving during the time of observation (∼10 years) for low frequency signals. However, as demonstrated here, PTA is sensitive to the difference in GW amplitude h felt by the pulsar and our radioobservatories on Earth. Pulsars are on average at a distance of about 1kpc, meaning that it takes ∼ 1000 years for their radio emissions to reach us. For such timescale, it is no longer possible to ignore the frequency evolution of the SMBHB system. Therefore, the pulsar term in δ t GW can have different frequency and phase than the Earth term. For pulsar α, the frequency ω(t p,α ) and phase Φ(t p,α ) of the pulsar term are obtained using equations • ∆ f α = ω(t p,α ) -ω(t ) /2π < 1/T where the pulsar term and the Earth term have nearly the same frequency and their difference is below the resolution 1/T of the array. However, the phase Φ(t p,α ) will be different for each pulsar because of their disparate distances and sky locations, destroying the phase coherency of the pulsar terms, while the Earth phase Φ(t ) will be the same for all pulsar.

• ∆ f α = ω(t p,α ) -ω(t ) /2π > 1/T where the pulsar term and the Earth term have different frequencies that the array can resolve. In that case, they can be treated as separate components. We can then include the pulsar term in our analysis or consider it to be an extra source of noise because of its random phase and frequency (different for all pulsars), focusing on the Earth terms that are actually in phase.

In Sesana and Vecchio 2010, they estimate the distribution of ∆ f α for 100 isotropically distributed pulsars across the sky at 1kpc distance using a realistic simulation of SMBHB population (figure 1.15). They found that the frequency difference between pulsar and Earth term is mainly above the resolution of the array. To that extent, they can be treated separately.

Multi-band gravitational wave astronomy

As we have seen in the previous section, PTA aims to detect very low frequency GW signal in the nanohertz band. But the efforts of the scientific community gave birth to other detectors allowing to probe the universe for the presence of GW at higher frequencies (between 1-10 3 Hz). Some of these detectors are already functioning, some others are still being designed. We will give a brief overview of the main projects that are currently being developed (see figure 1.16).

Ground based detectors : LIGO/Virgo

The GW detectors of the LIGO (Laser Interferometer Gravitational-Wave Observatory)

and Virgo collaborations follow the same concept. The Galactic Binaries (GBs) present in the Milky way (white dwarf binaries) will be detectable by LISA. They share similarities with the PTA CGWs in the sense that they will be seen as deterministic components resolved on top of a stochastic signal generated by the large population of GBs (the stochastic "foreground"). The equations governing their evolution are the same as the CGWs, but here for higher frequencies (∼ 10 -3 Hz) and lighter systems (M c ∼ 1M ). In chapter 4, we develop and test a data analysis method using simulated LISA data for the detection of GBs (see section 4.4).

the samples in δ t are independent and Gaussian distributed with a standard deviation σ i that is determined by parameters θ A1 , the likelihood writes as :

L( δt | θ A , M A ) = i 1 2πσ i ( θ A ) exp - 1 2 δt 2 i σ 2 i ( θ A ) = 1 |2πN| exp - 1 2 δt N -1 δt , (2.1)
with |.| denoting the determinant and N = diag{σ i ( θ A )} the diagonal matrix containing the σ i ( θ A ).

Intuitively, we want the values of parameters θ A that maximize the likelihood to match the data δ t that was actually observed. Finding the maximum likelihood thus gives estimators of the parameters θ A . However, working with fixed values does not provide information about the probability distribution of the θ A .

Likelihood ratio test

We can compare two models A and B by taking the ratio of their respective likelihood functions. We often work with the logarithm of this ratio that we note Λ :

ln Λ = ln L( δt | θ A , M A ) L( δt | θ B , M B ) . ( 2.2)
This quantity can be used to measure the statistical significance of a model with respect to another. In the case of signal detection, we compare a model with noise + the anticipated signal to the model with only noise. We refer to the latter as the "null hypothesis" often noted H 0 . In the particular case where the two models are nested (i.e. when θ A = θ B θ s with θ s the parameters associated with the additional signal present in A)

and for large number of samples, this ratio follows a chi squared distribution (see Willk's theorem in [START_REF] Wilks | The Large-Sample Distribution of the Likelihood Ratio for Testing Composite Hypotheses[END_REF]).

Posterior distribution

In the Bayesian framework, we set up a model and define a prior probability distribution for each of the parameters. This prior probability is a key point to the approach because model parameters are now treated as random variables with probability distributions.

We want to update our prior knowledge π( θ A |M A ) of parameters θ A for a model M A , by expressing the probability of observing them given some data δ t . Using the rules of conditional probability we get the posterior probability :

p( θ A | δt , M A ) = Pr( δt | θ A , M A )Pr( θ A |M A ) Pr(δ t |M A ) = L( δt | θ A , M A )π( θ A |M A ) Z A .
(2.3) with θ A the vector of model parameters, L( δt | θ A , M A ) the likelihood, π( θ A |M A ) the prior probability and Z A = p(δ t |M A ) the evidence, expressing the probability of measuring data δ t given model M A , acting as a normalizing constant.

In order to assess the statistical significance of a model M A with respect to another model M B , we need to compare the odds Pr(M A |δ t ) and Pr(M B |δ t ) of the considered models given data δ t . Using conditional probabilities (Bayes theorem) we get the odds ratio :

Pr(M A |δ t ) Pr(M B |δ t ) = Pr(δ t |M A )Pr(M A ) Pr(δ t ) Pr(δ t |M B )Pr(M B ) Pr(δ t ) = Z A × p(M A ) Z B × p(M B )
.

(2.4)

The evidence Z corresponds to the integral of the posterior probability distribution over whole parameter space (fully marginalized posterior) :

Z A = d θ A p( θ A | δt , M A ) = d θ A L( δt | θ A , M A )π( θ A |M A ). (2.5)
Assuming we are agnostic about models M A and M B , we can take the two probabilities p(M A ) and p(M B ) to be equal. In that particular case, the odds ratio reduces to the Bayes factor :

B AB = Z A Z B , (2.6)
where B AB is the Bayes factor, the quantity that we use for model selection. The bigger the B AB , the preferred model A is over model B . In GW data analysis, we use it to search for the presence of GWs in the data by comparing a model with a signal against another model without that signal. For large values of B AB , we can say that we have detected that signal. To decide which model is best, we can refer to Jeffrey's scale [START_REF] Jeffrey | The Theory of Probability[END_REF]):

log 10 B AB Strength of evidence 0 to 1/2 Not decisive 1/2 to 1 Substantial 1 to 2 Strong >2 Decisive
It is important to note that other scales exist (Kass and Raftery 1995), so deciding whether the Bayes factor is significant may seem rather arbitrary. To that extent, it is possible to claim detection only when its value is high enough (>10 2-3 ). We must remain cau-tious when dealing with intermediate values of B AB that are sensitive to our own "bias" factors.

This quantity is often difficult to compute since we can only approximate the integral in (2.5) that yields the evidence Z, because of our poor knowledge of the posterior distribution. To estimate the latter, we resort to sampling algorithms allowing us to explore the parameter space. • The spin frequency and spin-down rate measuring the deceleration of the pulsar due to its emissions,

Modelling

• The sky position and proper motion of the pulsar in the sky,

• The dispersion measure due to the propagation of the pulses through the interstellar medium with which they can interact during their journey to Earth,

• The binary orbital effects when the pulsar orbits a massive companion that can induce GR effects like Shapiro delay or Einstein delay (or simply Doppler shift due to the orbital motion).

The geomoetry of the problem is expressed in the Solar System Barycenter (SSB) that we estimate using Solar System Epheremides (SSE) providing the position of the planets at a given time. This translation calls for an additional delay, the Roemer delay given by :

∆ R = - r • R B B c , (2.7)
with r the position of the radio observatory in the SSB frame and R B B the unit vector pointing to the sky position of the pulsar in the SSB frame. If this effect is not accurately quantified (i.e. errors in the SSE estimations hence in the position of the SSB), it can leave some artifacts in the timing residuals that will appear as a signal correlated between pulsars [START_REF] Tiburzi | A study of spatial correlations in pulsar timing array data[END_REF]).

For our analyses, we used a pre-fitted dataset including the timing residuals δ t with the associated values of timing model parameters for all pulsars so we did not have to perform the timing model fit ourselves. However, there are still unmodelled noises in the δ t that we will present in the next section.

Noises in PTA

To use PTA as a detector we need to identify all the sources of noise that will interfere with the desired signal (Verbiest and Shaifullah 2018). Each noise will be represented by a set of parameters in the model. One of the drawbacks of the PTA being a pulsar based detector is that it is not possible to characterize noises beforehand. They should be set as free parameters while running detection scripts, often leading to high dimensional parameter space and heavy computational cost.

We We assume that the noises are gaussian and the residuals should behave as a normal distribution if the model is matching well with the data. Including these white noise parameters allows us to rescale the residuals so they fit this assumption.

σ2 k = E 2 f ,k × σ 2 T O A,k + E 2 q,k (2.8) Each E f ,k (EFAC), E q,k (EQUAD)
Another white noise component is ECORR. It takes into account the jitter-like effect (pulse to pulse variation) that is correlated between TOAs for a given observation epoch.

Note that this noise is only modelled in the NANOGRAV data and is not to be considered in other datasets2 [START_REF] Lentati | From spin noise to systematics: stochastic processes in the first International Pulsar Timing Array data release[END_REF]).

In general, including white noise leads to a lot of additional parameters because of the many different backend-receiver systems in radio-telescopes. For that reason, it is common practice to evaluate the maximum likelihood values of EFAC, EQUAD and ECORR for the TOAs of each pulsar separately, and use these values as fixed for GW analysis that employs the whole PTA array. This allows to considerably reduce the size of the parameter space and the computational cost (Aggarwal 2019).

Red noise

S RN ( f ) = A 2 RN 12π 2 yr 3 f yr -1 -γ R N (2.9)
with A R N the amplitude and γ RN the spectral index.

In

Dispersion measure variation

The photons emitted from the pulsars get scattered by the interstellar medium on their way to Earth. The variability of the electron density on the axis Earth-pulsar introduces additional time-correlated noise in the data which resembles a low frequency red noise. However, photons with different frequencies will interact differently with the interstellar medium. Thus, observing incoming photons at different frequency bands allows us to distinguish this dispersion noise from a regular red noise. Its power spectral density writes as follows :

S D M ( f ) = A 2 D M 12π 2 yr 3 f yr -1 -γ D M 1400MHz f ν 2 (2.10)
with A D M the amplitude, γ D M the spectral index and f ν the photon frequency 3 .

The dispersion measure (DM) noise is individual to each pulsar because it depends on the density of the interstellar medium along the axis Earth-pulsar.

Bayesephem

The physics of PTA is expressed in the solar system barycenter (SSB) frame.

To translate from Earth's reference frame to SSB, we need a precise knowledge of Earth's position with respect to the SSB [START_REF] Hobbs | Manchester "tempo2, a new pulsar-timing package -I. An overview[END_REF][START_REF] Vallisneri | Modeling the Uncertainties of Solar System Ephemerides for Robust Gravitationalwave Searches with Pulsar-timing Arrays[END_REF]). We get this information via solar system ephemeris (SSE). If the translation is not performed correctly, one can observe timing residuals in the TOAs due to errors in the determination of planetary motion (thus, an error in the determination of the SSB). Small deviations can lead to noticeable delay and interfere with GW detection. Bayesephem is a method that proposes to mitigate the systematic errors in the SSE by including deterministic corrections in the timing residuals (Roemer delay). Gaussian weights are given to small deviations in the orbital parameters of planets, allowing to absorb uncertainties on planetary position and correct errors on the determination of the SSB.

Single pulsar likelihood

To get the timing residuals δt α of one pulsar, we need to fit a timing model to the observed TOAs, t obs,α . This timing model depends on many parameters (pulsar position, pulsar rotation frequency, spin down, parallax, ...) and gives a predicted value of the TOAs, t pr ed ,α [START_REF] Hobbs | Manchester "tempo2, a new pulsar-timing package -I. An overview[END_REF]). The difference between the observed and predicted values of the TOAs gives the timing residuals. Small imperfections are present in the timing model and assumed to be linear so first order corrections are included for better precision (Lentati et al.

2013

), with M the design matrix (containing the first derivatives of δ t pr ed ,α with respect to the timing model parameter errors) and the vector of small offset in timing parameters.

δt α = t obs,αt pr ed ,α = M + n.

(2.11)

The noises n are assumed to be Gaussian and stationary. It follows that the timing residuals δ t α behave as a normal distribution as developped in Haasteren and Vallisneri 2014. Given the noises that we presented in the previous section, we can write for n :

n = n W N + n RN + n D M , (2.12)
where we have the correlated RN and DM noises 4 and the uncorrelated WN. We model the noises n as Gaussian processes (Haasteren and Vallisneri 2014) and we can express them in time domain as a Fourier decomposition to the N f th harmonic :

n( t obs,α ) = N f l X l cos(2π f l t obs,α ) + Y l sin(2π f l t obs,α ) = Fa, (2.13)
with

a =             X 1 Y 1 X 2 . . . X N f Y N f             and F =        cos(2πt 0 f 1 ) sin(2πt 0 f 1 ) . . . cos(2πt 0 f N f ) sin(2πt 0 f N f ) cos(2πt 1 f 1 ) sin(2πt 1 f 1 ) . . . cos(2πt 1 f N f ) sin(2πt 1 f N f ) . . . . . . . . . . . . . . . cos(2πt N t -1 f 1 ) sin(2πt N t -1 f 1 ) . . . cos(2πt N t -1 f N f ) sin(2πt N t -1 f N f )        (2.14)
where N t is the number of samples in t obs,α , X l and Y l are the Gaussian distributed weights, the cos(2πt f l ) and sin(2πt f l ) are the basis functions and the f l are the N f linearly spaced frequencies starting from lowest frequency f 0 = 1/T with T the total timespan of observation of the pulsar (Haasteren and Vallisneri 2014) T = max t obs,αmin t obs,α .

In theory, a Fourier basis is orthogonal when the t obs,α are evenly spaced. This is not the case for a real PTA dataset because of stochasticity in the sampled values of TOAs. We still consider that the basis is approximately orthogonal so we can write the PSD Σ as a function of the Fourier coefficients a :

Σ =< a a >= diag{S n ( f L , θ H P )}/T, (2.15
π(a, Σ) = exp -1 2 a (Σ) -1 a |2πΣ| . ( 2.16)
We can also write the covariance matrix for noise n as :

< n n >= F < a a > F = F ΣF. (2.17)
This result will be useful to express the total covariance matrix of the likelihood. The final expression for the residuals δt is :

δt α = δt α -M -F RN a RN -F D M a D M (2.18)
Assuming Gaussian likelihood for the δ t α with uncorrelated white noise (therefore assuming a likelihood of the form 2.1), we have for the posterior distribution :

p(δ t α | θ α ) = exp -1 2 δt α N -1 δt α |2πN| exp -1 2 a (Σ) -1 a |2πΣ| π( )π( θ H P ). (2.19)
with N the diagonal matrix of the σi presented in equation (2.8), the Gaussian prior on a, π( ) the prior on the timing model parameter offset and π( θ H P ) the prior on the hyper-parameters.

The posterior distribution can be analytically marginalized with respect to the a thanks to their Gaussian behaviour. However, π( ) is a uniform distribution on the whole real axis and requires a different approach for marginalization. It is possible to mimick a uniform π( ) using a Gaussian prior in the limit :

lim λ→∞ exp -1 2 (λI) -1 |2πλI| = π( ). ( 2
p(δ t α | θ α ) = L( δt | θ α , C α )π( θ H P ) = 1 |2πC α | exp - 1 2 δt α C -1 α δt α π( θ H P ), (2.21)
with 5 :

C α = N + C T M ,α + C RN ,α + C D M ,α = N + M λIM + F RN Σ RN F RN + F D M Σ D M F D M . (2.22)
Bare in mind that this timing model marginalisation "trick" is an approximation as λ cannot be numerically set to ∞ and we use instead very large values (10 40 ). It is more convenient to work with this expression of the posterior for Bayesian analyses because the covariance matrix can be written as :

C α = N + T ΦT. (2.23) T =     M F RN F D M     Φ =     λI 0 0 0 Σ RN 0 0 0 Σ D M     .
(2.24)

In this form, the covariance matrix can be inverted using the Woodbury lemma :

C -1 α = (N + T ΦT) -1 = N -1 -N -1 T(Φ -1 + T N -1 T) -1 T N -1 , (2.25)
which is used to optimize numerical computation. In ENTERPRISE, the argument of the exponential δ t α C -1 α δ t α is computed by first solving the linear system for x and taking the dot product with δ t α N -1 T from the left :

(Φ -1 + T N -1 T) x = T N -1 δ t α ⇔δ t α N -1 T • x = δ t α N -1 T(Φ -1 + T N -1 T) -1 T N -1 δ t α , (2.26) giving, δ t α C -1 α δ t α = δ t α N -1 δ t α -δ t α N -1 T • x (2.27)
Finally, if we want to add the contribution of any deterministic signal to our model (CGW signal, BAYESEPHEM), we need to transform the residuals as follows :

δ t α → δ t α -s( t obs,α , θ) (2.28)
where s( t obs,α , θ) is a deterministic signal, function of time.

Multi pulsar likelihood

In the case of deterministic GGW signal, the likelihood for an array of pulsars is simply the product of the single pulsar likelihoods. For N pulsars indexed by α, we have the total timing residual vector :

δ t = α δt α . (2.29)
When we are using the whole array, the Fourier basis that is used for GP noises is defined from the total observation timespan of the array (T = max α δ t obs,α -min α δ t obs,α )

and the covariance matrix is : In this thesis, we are interested only in CGW signal so the covariance matrix we use is always block diagonal and the total likelihood factorises as the product of single pulsar likelihoods :

C ab =        C 0 0 0 0 0 C 1 0 0 0 0 . . . 0 0 0 0 C N        . ( 2
L P T A ( δt | θ, C ab ) = α L( δt | θ α , C α ), (2.31)
where the C α appear in the diagonal of the C ab . The GWB has a correlation matrix C GW B similar to that of a common red noise but with additional correlation terms χ ab that relate two pulsar a and b given their angular separation ζ ab (see figure 2.1) :

Correlated signals

C GW B = χ ab C RN , (2.32)
with

χ ab = 1 if a = b.
Dipolar correlation Dipolar correlations can originate from errors in the SSE [START_REF] Tiburzi | A study of spatial correlations in pulsar timing array data[END_REF]). The position of the observatories on Earth with respect to the SSB must be known to high precision, otherwise measurable (and correlated between pulsars) effects will appear in the timing residuals when the translation to the SSB is made (see section 2.2.1). In that case, the correlation function is of the form :

χ ab ∝ cos(ζ ab ). ( 2

.33)

A dipolar correlation is somewhat similar to a HD curve. This resemblance can cause problems when searching for GWB because the two signals may interact in the model. A correct modelisation of the SSE is then crucial to PTA efficiency [START_REF] Vallisneri | Modeling the Uncertainties of Solar System Ephemerides for Robust Gravitationalwave Searches with Pulsar-timing Arrays[END_REF].

Monopolar correlation A signal with monopolar correlation is fully correlated between

pulsars. There is no dependence in the angular separation ζ, hence the off-diagonal terms of the covariance matrix are equivalent. We have :

χ ab = 1. ( 2

Markov Chain Monte Carlo

When dealing with high dimensional models and complex datasets, we often work with difficult posterior probability distributions that are not analytically invertible or for which finding maxima is not trivial. Furthermore, it is almost impossible to visualize the proba-

Markov process

A Markov process (or Markov chain) is a stochastic process defining a series of events in which each event only depends on the previous. Such a process is said to be "without memory" or simply "Markovian". Let X = {A, B,C , ...} a random variable following the rules of a Markov process and p(t ) a vector expressing the probabilities of finding X at states {A, B,C , ...} at discrete time t . For a set of samples {X T }, the previous no-memory condition is satisfied by :

Pr(X T |X T -∆t , X T -2∆t , X T -3∆t , ...) = Pr(X T |X T -∆t ). (2.36)
Because of this condition, we can define a transition matrix Q that is independent of time, relating two successive events at times t and t + ∆t in the chain so that the components of probability vector p(t + ∆t ) change as :

p j (t + ∆t ) = i Q i j p i (t ). (2.37)
To preserve probabilities p j (t + ∆t ) to be between 0 and 1 and p(t + ∆t ) to be normalized, Q must be positive and respect the condition :

j Q i j = 1 (2.38)
When vector p(t + ∆t ) is different from p(t ), we are in transient regime and probabilities are function of time. But when vector p(t + ∆t ) is equal to p(t ), we have a stationary process for which we will denote probability vector as P so that :

P j = i Q i j P i .
(2.39)

The time independence of stationary processes is a key element that is utilized in the MCMC algorithms. Nonetheless, stationary does not necessarily mean equilibrium and it is custom to define a stronger condition to force stationarity.

Detailed balance

The detailed balance is a way of defining equilibrium by reversibility of a process. Indeed, if a process is reversible it is stationary, but not always the other way around. Detailed balance is expressed as :

Q i j P i = Q j i P j (2.40)
If this equality is respected, we have :

Q i j P i = Q j i P j ⇔ i Q i j P i = i Q j i P j ⇔ i Q i j P i = P j , (2.41)
where we have used (2.38) to go from second to last line. In other words, this ensures that going forwards or backwards in time has no effect on the system. A reversible Markov process must fulfill this condition and the core sampling algorithm of MCMC is entirely based on it.

The sampling algorithm

Until now, we have been working with a finite number of components for probability vector p i and transition matrix Q i j . We can easily generalize to continuous functions by saying :

p i → p( θ) Q i j → q( θ| θ * ), (2.42)
where θ is a vector in a d-dimensional parameter space following probability distribution p( θ). In that case, the detailed balance condition is :

q( θ| θ * )p( θ) = q( θ * | θ)p( θ * ) (2.43)
The idea behind the sampling algorithm is to use the properties of a stationary Markov process to draw samples from a target distribution p( θ). In other words, starting from an initial sample θ 0 , we want to build a Markov chain (i.e. a set of N samples { θ N }) that performs a random walk in the parameter space following the probability of p( θ). By ensuring that the target p( θ) is the stationary distribution of the Markov process, we ensure that the θ are sampled from p( θ). The resulting set of samples { θ N } can later be post-processed to study the probability distributions of the model parameters.

For any probability distribution p( θ), we want to find a transition matrix q( θ * | θ) that always respects the detailed balance condition. We first introduce a proposal distribution g ( θ * | θ) from which it is possible to draw new samples θ * given θ, the current sample (or current state of the Markov chain) 6 . Then we define an acceptance probability α( θ * , θ)

that decides if the new sample is accepted or rejected. We can write the transition kernel as the product of the two :

q( θ * | θ) = g ( θ * | θ)α( θ * , θ), ( 2.44)
and plugging this expression in (2.43), we get :

α( θ * , θ) α( θ, θ * ) = p( θ * )g ( θ| θ * ) p( θ)g ( θ * | θ) . ( 2

.45)

We can say that if a new jump going from θ to θ * is not accepted, then the same jump but reversed would always be accepted, that is to say, α( θ, θ * ) = 1. Thus, we have :

α( θ * , θ) = min 1, p( θ * )g ( θ| θ * ) p( θ)g ( θ * | θ) , (2.46)
where the min ensures that 0 ≤ α( θ * , θ) ≤ 1. This is called the Metropolis-Hastings choice [START_REF] Hastings | Monte Carlo Sampling Methods Using Markov Chains and Their Applications[END_REF]).

The acceptance probability will entirely determine our random walk in the parameter space, which by construction, respects the detailed balance condition for target distribution p( θ). Starting from sample θ, we jump to sample θ * with a probability α( θ * , θ) or stay at sample θ with probability 1 -α( θ * , θ). The process is repeated N times and obtained samples at each step are concatenated into a chain (Markov chain) that forms our final set of samples { θ N } following the target probability distribution p( θ).

Let us summarize the sampling algorithm in steps :

• For current (or initial) sample θ 0 compute p( θ 0 )

• Draw new sample θ * from proposal distribution g ( θ * | θ 0 ) and compute p( θ * )

• Compute the acceptance probability α( θ * , θ 0 ) = min 1,

p( θ * )g ( θ 0 | θ * ) p( θ 0 )g ( θ * | θ 0 )
• Accept new sample with probability α( θ * , θ) and set θ 0 = θ * OR reject new sample with probability 1 -α( θ * , θ) and set θ 0 = θ

• Append θ 0 to chain of samples and repeat

Hypermodels for model selection

To compute the Bayes factor, we can use a method presented in Hee et al. 2015 called hypermodels. A hypermodel is the union of two or several submodels. Its vector of parameters θ is :

θ = N -1 i =0 θ i (2.47)
with θ i the parameters for submodel i and N is the total number of submodels.

We introduce an additional integer parameter j with uniform probability distribution n( j ) = {0, 1, 2, ..., N -1}. Throughout sampling, we draw an integer k from n( j ) that determines which submodel is going to be active for the likelihood computation so that :

L( θ, j = k|δ t ) = L k ( θ k |δ t ) (2.48)
where L k ( θ k |δ t ) the likelihood of submodel k. In that case, the posterior distribution of the hypermodel is :

p( θ, j |δ t ) = p( θ 0 ∪ θ 1 ∪ ... ∪ θ i ∪ ... ∪ θ N -1 |δ t )n( j ) = p N -1 i =0 ∪ θ i |δ t n( j ) = L N -1 i =0 ∪ θ i |δ t N -1 i =0 π( θ i )n( j )/Z, (2.49)
where Z is the evidence of the hypermodel acting as a normalizing constant and π( θ i ) the prior probabilities for parameters θ i of submodel i . Fixing j = k, we get :

p( θ, j = k| δt ) = L k ( θ k | δt ) N -1 i =0 π( θ i )n( j = k)/Z, ( 2.50)
which the active posterior distribution when an integer k is drawn for n( j ). If we marginalize this expression for all θ, we obtain :

p( j = k|δ t ) = d θp( θ, j = k| δt ) = d θL k ( θ k | δt ) N -1 i =0 π( θ i )n( j = k)/Z = d θ k L k ( θ k | δt )π( θ k ) N -1 i =0 i =k d θ i π( θ i )n( j = k)/Z. (2.51)
Using the definition of the evidence Z k and the fact that prior probabily distributions π( θ i ) are normalized :

d θ k L k ( θ k | δt )π( θ k ) = Z k N -1 i =0 i =k d θ i π( θ i ) = 1, (2.52)
equation (2.51) greatly simplifies to :

p( j = k|δ t ) = Z k n( j = k)/Z ⇔ Z k = p( j = k|δ t ) n( j = k)
Z.

(2.53)

The p( j = k|δ t ) corresponds to the odds ratio for parameter j to be at value k. So for a chain of N samples, if n k is the number of samples corresponding to j = k, we have :

p( j = k|δ t ) = n k N , (2.54)
and for two different models k and k , we directly see that :

Z k Z k = n k n k , ( 2.55)
where we have used n( j = k) = n( j = k ) because n( j ) is a uniform distribution. Recall from section 2.1.3 that the ratio of evidences corresponds to the Bayes factor which is the quantity that we used to do model selection. Therefore, using hypermodels, the evaluation of the Bayes factor between models A and B is :

B AB = n A n B (2.56)
with n A and n B respectively the number of final samples for which submodel A was active and submodel B was active. They represent the time spent by the sampler on one model. Comparing them shows which model is preferred according to the MCMC rules.

The big advantage of this technique is that it requires no direct evidence estimates which are usually complicated multidimensional integrals 7 . However, it can sometimes be a bit long to run because we need the sampler to efficiently explore every submodel.

This requires many samples for which we need to run our MCMC for longer.

CHAPTER 3 . We achieved the best upper limit to date at low and high frequencies of the PTA band thanks to improved effective cadence of observations. In our analysis, we have taken into account the recently discovered common red noise process, which has an impact at low frequencies. We also find that the peculiar noise features present in some pulsars data must be taken into account to reduce the false alarm. We show that using custom noise models is essential in searching for continuous gravitational wave signals and setting the upper limit.

Introduction

The goal of the Pulsar Timing Array (PTA) collaborations is to detect gravitational wave We do not yet know the nature of this process, and its interpretation as GW background is inconclusive: the data is not informative enough to resolve the Helllings-Downs spatial correlations [START_REF] Hellings | Upper limits on the isotropic gravitational radiation background from pulsar timing analysis[END_REF], which should be present in the case of the GW signal.

In this work, we search for continuous GWs which could be present in the data in addition to the stochastic GWB. Following the steps of previous studies [START_REF] Babak | European Pulsar Timing Array limits on continuous gravitational waves from individual supermassive black hole binaries[END_REF] signal produced by a SMBHB binary in a circular orbit. In our study, we neglect the pulsar terms during the search and setting an upper limit on GW amplitude. However, we do an in-depth analysis of the (weak) candidate events identified as plausible GW signals. In the followup analysis on the restricted parameter space (frequency and sky position), we extend our model to include (i) pulsar term, (ii) eccentricity, (iii) extend the model beyond the assumption of a single source. For the first time, we have included in the analysis the common red component as a part of the total noise model.

The main results of the paper can be summarized as follows. We did not detect any CGW signal and set an upper limit on GW amplitude. We have found that the noise model plays a crucial role in interpreting PTA observations. The detailed analysis performed on the most promising candidate event revealed that it could be explained by a timecorrelated high-frequency noise in one of the pulsars.

This chapter is organized as follows. In the next Section, we will briefly describe the IPTA DR2 dataset and the data model used in the analysis. Most of the material needed for this Section is available in the literature, and we heavily refer to it, keeping only parts which are necessary for further presentation. In Section 3.4 we describe the methodology which we have followed to get our results presented in Section 3.5. In Section 3.6 we give a detailed follow-up study of a most promising GW candidate event and demonstrate the importance of noise modelling at high frequencies. We conclude with Section 3.7.

Throughout this chapter, we adopt geometrical units G=c=1.

3.2

The noise and signal model

The noise model

In all that follows, the noise model that we use is based on what is presented in section 2.2.2. That is to say, we use a PTA likelihood with marginalized timing model, including red noise and dispersion measure noise. The covariance matrix is block diagonal because we do not include common correlated processes between pulsars. However, we do include a common red noise process (except for the ranking) because it was found in the IPTA DR2 with good statistical significance (Antoniadis 2022). To summarize, the covariance matrix C α for the noise model of one pulsar is :

C α = N α + C T M ,α + C D M ,α + C RN ,α + C C RN , ( 3

The signal model

For the signal model, we consider monochromatic CGW signal as presented in section 1.4.4. The timing residuals expected in the PTA data for a CGW signal coming from a circular SMBHB:

s a (t , Ω) = - A F A ( Ω)[s A (t ) -s A (t -τ a )] (3.2)
with τ α = tt p,α the delay between Earth time t and pulsar time t p,α and :

s + (t ) = M 5/3 d L ω(t ) 1/3 sin[2Φ(t )](1 + cos 2 ι) cos 2ψ +2 cos[2Φ(t )] cos ι sin 2ψ , (3.3) s × (t ) = M 5/3 d L ω(t ) 1/3 sin[2Φ(t )](1 + cos 2 ι) cos 2ψ -2 cos[2Φ(t )] cos ι sin 2ψ , (3.4)
where M is the chirp mass, d L the luminosity distance, ω(t ) the CGW orbital angular frequency, ι is the orbital inclination to the line of sight, ψ is a polarization angle and Φ(t

)
is the phase of CGW.

The

F + ( Ω) = 1 2 ( m • p) 2 -( n • p) 2 1 + Ω • p , (3.5
)

F × ( Ω) = ( m • p)( n • p) 2 1 + Ω • p . (3.6)
In this work 1 , we neglect the pulsar term considering it as a part of the noise, assuming that the source has evolved sufficiently over τ α to move the pulsar term off the earth-term frequency. Including pulsar term should improve the parameter estimation but comes with a huge price of the increase in the complexity of the likelihood surface and the dimensionality of parameter space (2 additional parameters per pulsar for phase and frequency of pulsar term, e.g. see Corbin and Cornish 2010). We foresee the possibility of following up the candidate events (identified using the earth term only) with the extended signal model (pulsar term, eccentric orbit) on the reduced parameter space. We also neglect the evolution of the GW frequency (ω o = 2π f g w ) over the observation time.

The frequency evolution becomes potentially measurable for the heavy sources emitting at frequency ≥ 10 -7 Hz, neglecting the frequency evolution does not prevent us from detecting the sources but introduces a bias in the measured GW frequency (overestimating 1 Based on the submitted IPTA collaboration paper where I am corresponding author.

CGW parameter

Range log 10 h [-18, -11]

f g w (Hz) [10 -9 , 10 -7 ] φ 0 [0, 2π] cos ι [-1, 1] ψ [0, π] θ [0, π] φ [0, 2π]
Φ(t ) = ω o t + φ 0 /2, (3.7)
where φ 0 is initial orbital phase. Finally, the CGW amplitude h is a function of M, d L and f g w given by

h = 2M 5/3 (π f g w ) -2/3 d L . ( 3.8)
We consider the model containing only one CGW signal. This model still detects multiple CGW if they are present in the data at the non-overlapping Fourier frequencies (see Babak and Sesana 2012 for discussion). If we find more than one candidate with sufficient statistical significance as potential GW sources, we will conduct additional investigations extending our model to include several CGWs. We start the analysis with 1 CGW source characterized by 7 parameters summarized in table 3.1 together with their prior range (we always assume a uniform prior2 ).

Ranking the pulsars

Within the pulsars that make up the PTA, some may be more sensitive to the presence of GW signal, depending on their noise properties and time of observation. To characterize the contribution of individual pulsars we can evaluate the signal to noise ratios (SNR)

separately for each pulsar after injecting fake GW signal in their TOAs. It is interesting to rank them according to their SNR to see which contribute the most to the sensitivity of the array. We eventually discard the worst pulsars to lighten the data and the computational cost without significant loss in detectability.

Ranking method

In L. Speri 20223 we present two methods (Chimera and SNR B maximization, see reference and appendix 6.5 for details about the methods) to rank pulsars with respect to GWB, which is a stochastic and correlated signal. CGW are deterministic signals and their analysis has been treated separately from the stochastic GWB. CGWs are included in the model as a periodic delay applied to the timing residuals δt while the effect of the GWB is included in the covariance matrix C of the likelihood. This fundamental difference between the two signals and their mathematical description calls for a different ranking method.

Here, we want to rank pulsars according to their response to a CGW signal. One way to proceed is to inject a large number of fake CGW signals with randomized parameters except for fixed frequency and amplitude [START_REF] Babak | European Pulsar Timing Array limits on continuous gravitational waves from individual supermassive black hole binaries[END_REF]. Then, for each pulsar, the CGW signal-to-noise ratio is computed for each injection and averaged numerically. In this way, we have the average response of each individual pulsar in the array at a given frequency of the CGW signal. This averaging can also be done analytically, as shown in the following paragraph. Note that we refer to the signal-to-noise ratio of CGWs using the acronym SNR. However, we use the symbol ρ to distinguish the SNR of CGWs from the previously defined SNRs.

In the likelihood of Eq.(2.21), the inclusion of a deterministic signal is performed by changing the timing residuals as δ t → δ ts( θ) where s( θ) is the signal template we aim to measure. In that case, the likelihood can be rewritten as:

ln L = - 1 2 [(δ t |δ t ) + ( s| s) -2(δ t | s) + Tr ln 2πC] , (3.9)
where we have introduced the noise weighted inner product ( x| y) = x T C -1 y.

We can now calculate this expression for the hypothesis of the presence of a CGW (H 1 )

versus its absence (H 2 ) (see section 2.1.2). The expectation value of the log-likelihood ratio becomes:

〈ln Λ〉 H 1 = ln( p(δ t | s) p(δ t | 0) ) H 1 = 〈(δ t | s) - 1 2 ( s| s)〉 H 1 = 1 2 ( s| s) , (3.10)
where ρ Opt = ( s| s) is the optimal SNR for the CGW source.

Since the source parameters are not known a priori, we average ρ 2 Opt over gravitational wave polarization ψ, initial phase φ 0 , inclination ι, and sky location (θ, φ). To do so, we analytically compute the integral over the defined bounds of the CGW parameters:

ρ 2 = π 0 d ψ π 2π 0 d φ 0 2π -1 1 d cos ι 2 -1 1 d cos θ 2 2π 0 d φ 2π ( s| s) . (3.11)
Using the formula for a CGW signal from a circular SMBHB, s(t , Ω), as presented in section 1.4.4, the SNR 2 averaged over CGW parameters takes this simple form (see appendix for calculus development):

ρ 2 (h, f ) = 4 15 h 2π f 2 (cos 2π f t | cos 2π f t ) + (sin 2π f t | sin 2π f t ) , (3.12) with h = 2M 5/3 (π f ) 2/3 d L , ( 3.13)
where f and h are the gravitational wave frequency and amplitude, M is the chirp mass and d L is the luminosity distance. For pulsar a, we evaluate ρ 2 a at the TOAs t a . We consider an Earth-term SNR because the inclusion of the pulsar term requires assuming a chirp mass and having a good knowledge of pulsar distance to estimate frequency evolution. Furthermore, the pulsar term will only promote pulsars at large distances with typically low SNRs (large-distance pulsars typically have a worse timing and their radio signal is weaker, leading to a degradation in the distance estimation). We are interested in the relative SNR contribution of pulsars at fixed signal frequency thus including pulsar term constitutes additional complexity that is beyond our purposes. If pulsar term and earth term are at same frequency (but different phase), their averaging yields twice the same average SNR 2 therefore leaving the relative contribution unchanged. The correlated noises (e.g. intrinsic and dispersion measure noises) are taken into account in the covariance matrix C of the noise-weighted inner product of the cosine and sine terms.

Common (correlated) processes were not included in our noise model, so the covariance matrix is block-diagonal. In this way, the likelihood can be factorized and SNR 2 s can be computed independently for each pulsar. In practice, we should incorporate the common processes in the noise model, but this adds another level of complexity that is irrelevant for the goal of the selection procedure 4 . The ultimate goal is identification of the best pulsars for CGW detection, and therefore, only the intrinsic properties of the pulsars were considered.

We estimate the relative contribution of one pulsar to the total SNR of the array using the normalized SNR 2 :

ρ2 a ( f ) = ρ 2 a (h, f ) b ρ 2 b (h, f ) , (3.14)
Note that the amplitude h cancels out in this expression and the CGW frequency f remains the only parameter. Therefore we can fix h to any value without affecting the ranking.

We construct the cumulative sum of the normalized SNR 2 s of the pulsars ranked from best to worst. We fix a threshold value for the SNR 2 cumulative sum above which pulsar

Ranking results

We now test the performance of the CGW ranking method using noise-parameter values previously extracted from individual pulsar noise analyses of the latest IPTA data release [START_REF] Perera | The International Pulsar Timing Array: second data release[END_REF]) and a "realistic" simulated EPTA dataset consisting in 40 pulsars with same RMS and sky locations as real EPTA dataset (Desvignes 2016) but with white noise only, ignoring intrinsic red noise and dispersion measure noise, and injected GWB (L. Speri

2022).

Because the ranking method is based on an exact noise-averaged formula, it is unnecessary to simulate noise realizations to test its performance. However, we still want to prove that the selected pulsars recover most of the total SNR in the presence of a true (i.e. non-averaged) signal. We test this by comparing the fraction of total SNR 2 obtained using the CGW ranked pulsars to that obtained from a random pulsar selection. For an array of N pulsars, the fraction of total SNR 2 , given a list of M < N pulsars, is defined as:

ρ 2 M = M a=1 ρ2 a , with 0 < ρ 2 M < 1 , (3.15)
where ρ2 a is the normalized SNR 2 defined in Eq (3.14). After extracting the list of best pulsars, we test the selection procedure as follows: The distributions are obtained with the list of best pulsars and random selection for a real IPTA dataset and the realistic EPTA dataset.

• We draw the CGW signal parameters θ from a uniform distribution with bounds defined as in the integral of Eq. (3.11), and with frequency between 1 and 100 nHz.

As pointed out in section 3.3.1 the strain amplitude has no influence on the ranking and therefore we fix it to h = 10 -14 .

• We compute the non-averaged optimal SNR ρ Opt = ( s| s) for each pulsar for a CGW signal s(t , θ) and we use this quantity to calculate the normalized ρ2 a defined in Eq. (3.14).

• We compute ρ 2 M -CGW for the list of best selected pulsars and ρ 2 M -rand for a random subset of pulsars of random size M .

• We repeat the previous steps one thousand times.

This gives us 1000 values of ρ 2

M -CGW and ρ 2 M -rand that we plot as histograms on Figure 3.3. For the IPTA dataset, the distribution of fractional ρ 2 M -CGW for the selected pulsars is narrowly peaked around a mean value 0.97. The random selection ρ 2 M -rand gives an almost uniform distribution with 0.50 mean value. The distribution is not uniform because ρ 2 a is not uniform and a few ρ 2 a values are much bigger while many others are very small. Similar results are obtained for the realistic EPTA dataset. We find that the number of pulsars which gives 95% of the SNR 2 is 22 for both datasets. Now we briefly discuss the comparison between the CGW and GWB selection methods. Focusing on the realistic EPTA dataset, we find an overlap between the identified best pulsars with the CGW method and GWB method as shown in Table 3.2. This time we run the Chimera and SNR B -maximization ranking for brief explanation of the methods) without fixing the six initial pulsars of the EPTA. We find that 17 pulsars are common to all three selection methods (highlighted in bold).

In summary, when true CGW signals are injected in the data, the CGW ranking method selects the pulsars which provides most of the SNR of the array, whereas a random selection is inefficient. This method extracts the few best pulsars to optimize the search for a CGW signal.

Sensitivity curve CGW strain amplitude for ρ = 1

It is possible to plot the PTA sensitivity to CGW signal using the SNR formula in equation 3.12. Indeed, the total SNR 2 of the array ρ 2 is the quadratic sum of the individual SNRs ρ 2 a of pulsars. We can invert the expression 3.12 to express the amplitude h of the CGW as a

function of ρ 2 ρ 2 (h, f ) = a ρ 2 a (h, f ), (3.16) h( f ) = 15π f ρ( f ) (cos 2π f t | cos 2π f t ) + (sin 2π f t | sin 2π f t ) , ( 3.17)
giving h as a function of the CGW frequency f and the noise parameters included in the covariance matrix of the noise weighted inner products.

By setting ρ to be equal to 1 we obtain the value of CGW amplitude corresponding to the threshold above which the signal starts to become louder than the noise 5 . This dependence on the noise becomes obvious when we plot the sensitivity curve : the threshold amplitude h as a function of frequency f .

Effect of noises on the sensitivity

Naturally, we expect the chosen noise model to have a strong influence on the sensitivity of the array. This is a major concern when dealing with PTA data because we cannot characterize the noises present in the array beforehand, we have to fit for them. A wrong modelization of the noises could lead to fake detections.

In Bayesian analysis, we treat noise parameters as random variables with probability distributions. In equation 3.17 we explicitly see the dependence of h on noise model parameters present in the noise weighted inner products (denominator of the expression), meaning that the threshold amplitude of detection should itself have a probability distribution. We propose to test two things :

• Study the influence of different values of noise parameters for the same model,

• Examine the influence of the size of the Fourier basis for the noises modelled as Gaussian processes (RN and DM, see section 2.2.2).

For the first point, we use previous MCMC noise analyses that were performed on individual pulsars to recover the posterior distributions of the noise parameters. We produce a 100 different sensitivity curves using 100 different samples randomly drawn from the posteriors of the noise parameters. We use the noise model presented in section 2.2.2 with a timing model, white noise, red noise and dispersion measure noise. For the second point, we plot two different sensitivity curves with same noise model but different number of frequency bins (30 and 100) for RN and DM with maximum likelihood values of parameters. We clearly see on figure 3.5 that using a short basis has an impact on the predicted sensitivity at high frequencies. In section 2.2.2 we explained that for one pulsar the Fourier basis is made of the N harmonic frequencies of 1/T , T being the total time of observation for that pulsar. Logically, choosing a N that is too small will brutally turn off the effect of the noise after the frequency N /T , producing a "step" in the sensitivity that we see on figure 3.5. This effect will be even more pronounced for red noises with low spectral index. For very steep RN, the white noise rapidly takes over and the step will not appear. All things considered, having a short Fourier basis will not have a strong impact on signal searches at low frequency but it is crucial at higher frequencies.

Method

We work within the Bayesian framework and start with running the search for the CGW signal. As mentioned above, we sample parameters of CGW together with the noise pa-rameters for RN and DM. We keep the white noise parameters fixed and marginalize over the timing model and BAYESEPHEM parameters. We made two runs: with and without CRN, to check how much it affects our result. We always use Markov chain Monte-Carlo (MCMC) sampler (Ellis and Haasteren 2017), and we use ENTERPRISE [START_REF] Ellis | Enhanced Numerical Toolbox Enabling a Robust PulsaR Inference SuitE Zenodo[END_REF] software to construct the models and compute the likelihood and prior probability.

We extensively use the single pulsar noise explorations runs performed before the main analysis. We have converted posteriors for the RN and DM into 2D histograms ({A R N ,D M , γ RN ,D M }) and use them as one of a proposal for those parameters. This empirical proposal improves the efficiency of MCMC and reduces the autocorrelation length of the chain (Aggarwal 2019).

During the search, we compute the Bayes factor (BF), comparing the null model (noise only) against the model where we have a CGW signal on the logarithmically spaced frequencies. In the absence of the detection, we proceed to setting an upper limit, building a 95% sensitivity curve. During the upper limit analysis, we used a uniform prior on the amplitude of the GW signal.

IPTA data contains 53-millisecond pulsars; however, not all of them are equally sensitive to the CGW. We have selected 21 pulsars which, on average, recover 95% of the array's total signal-to-noise ratio (SNR) to CGW. The selected pulsars are depicted on the projected sky map in Figure 3.6 as red stars; we have used large green stars and annotation for the four best timers. The ranking method is briefly outlined in Babak 2016 and in greater detail in L. Speri 2022. This significantly reduces the computational cost without much affecting the final results.

Historically we have performed the search with the noise model with the uniform settings across all pulsars in the array; namely, we have used 30 frequency bins for modelling both RN and DM processes (as it was done in preceding work [START_REF] Arzoumanian | Gravitational Waves from Individual Supermassive Black Hole Binaries in Circular Orbits: Limits from the North American Nanohertz Observatory for Gravitational Waves[END_REF][START_REF]The NANOGrav 11 yr Data Set: Limits on Gravitational Waves from Individual Supermassive Black Hole Binaries[END_REF]. Using this uniform setting we have obtained quite erroneous results, and after a long investigation, we realised that this noise model does not adequately describe the observational data (see section 3.6). We have switched to another noise model where we have used the custom-made noise model for the six best EPTA pulsars (see Chalumeau 2021 for details), and for other pulsars, we changed the number of used frequency bins: RN30DM100 -30 bins for the RN and 100 bins for DM modelling. We will continue using this short-hand notation for the noise model, showing explicitly the number of Fourier frequencies used by the base functions in the Gaussian process describing the corresponding noise. The results presented in the next section were obtained using the "custom" made noise model; we postpone the detailed discussion on the noise model selection and influence until Section 3.6.

Model selection

During the search, we consider two competing models: noise only and noise with 1 CGW signal. We compute the Bayes factor (BF) to measure which model is preferred by the IPTA DR2 dataset. In particular we employ hyper-model jumps to compute BF following the methods outlined in Hee et al. 2015 and implemented in the enterprise extensions [START_REF] Taylor | [END_REF]. In this method, we introduce a hyper-parameter which indexes the models; then, we perform sampling inside each model and in this parameter. The BF then is a fraction of the samples the chain spends in each model. For example, if we consider two models M A and M B with the hyper-parameter n, then the BF is the ratio (see section

) B A B = n A n B (3.18)
where n A , n B counts the number of samples in the chain corresponding to the models A, B . In our previous investigations, we have found that this method gives quite a reliable result, e.g. see Chalumeau 2021 where comparison is done against the evidence computation obtained with Dynesty nested sampling.

Later in this work we will compute BF between the models where we also vary the noise.

Upper limit

For obtaining the upper limit, we again use MCMC assuming a uniform prior distribution for log h within the bounds [-18, -11] on the fixed set of f g w . We use marginalized posterior probability distribution for the CGW amplitude p(h) to set 95% upper limit h 95 defined as :

0.95 = h 95 0 p(h)d h. (3.19)
We use a grid of 100 logarithmically spaced frequencies between 10 -9 and 10 -7 Hz. The lowest bound determined by the IPTA observational span 1/T obs that is ≈ 10 -9 Hz while the upper bound, 10 -7 Hz, is constrained by reduction in the sensitivity due to response (∝ f -1/3 ) and by our assumption that GW frequency does not evolve. We should take into account the frequency evolution of the GW signal above 10 The figure was obtained computing the 95% upper limit for h on 12 patches across the sky. Gaussian interpolation was used to smoothen it across boundaries.

Results

We compute BF between two models B cust om+CGW cust om

for a noise only model M cust om and noise+CGW signal model M cust om+CGW . The subscript "custom" corresponds to the noise model we have used for these main results and distinguishes it from other noise models considered in the next section. Similarly to the upper limit run, we have used a log-uniform prior on log h bounded by [-18, -11], and the BF was computed on the grid of 100 CGW frequencies ( f g w ) between 1 and 100 nHz.

The main result of this work is summarized in Figure 3.7. In the lower panel we plot the BF, and, as one can see, the noise model is usually preferred. There are few spikes where BF reaches 1, this is definitely not a detection, however PTA efforts probably should monitor carefully those frequencies in the future extended IPTA datasets. The excess in BF was also used as identification of the CGW candidates to follow. These results show that there is nothing to follow 6

: B cust om+CGW cust om ≤ 1.
The top panel of Figure 3.7 shows the upper bound of the 95% central credible region for CGW strain computed at the same set of frequencies. Note the "spiky" features at several frequencies (8.1, 13, 16) nHz corresponding to the outlier indicating potential candidates of CGW. However the corresponding BFs are (0.95, 0.33, 0.34) indicating no statistical significance in the analyzed data. We have computed the upper limit using uniform in GW strain prior with and without CRN on the fine frequency grid; the results are present in Figure 3.9. The upper limit slightly worsens at low frequencies when we add the CRN; this is understandable as we We have also overplotted the best CGW upper limit available to date based on the analysis of the NANOGRAV 11 years dataset (Aggarwal 2019) as a dashed (green) line. Note that only the nine-year NANOGRAV data set was included in IPTA DR2. As expected, our current results are better at very low frequencies thanks to the longer observation time.

Extended sky coverage, improved effective cadence of observation thanks to overlapping timing data (gaps coverage) and the addition of pulsars like J0437-4715 (only present in PPTA data Manchester 2013a), which is an excellent timer, is reflected in a much improved upper limit at 100 nHz, where we might expect the first detection of CGWs.

We want to point at the double-peak feature just above 10nHz (see Figure 3.7). This peak remains the same under the prior change (from the uniform in the amplitude to the uniform in the log-amplitude), which often corresponds to a signal present in the data.

The Bayes factors for those peaks are low (0.33, 0.34). Nonetheless, we should keep an eye on those frequencies in the next data releases.

Effect of noise modeling on the CGW search

In this subsection, we consider several noise models and compute BF between those models of noise with and without CGW. The main results presented in the previous section were obtained with the custom-made noise model for the best EPTA pulsars (Chalumeau 2021) and with RN30DM100 choice of Fourier frequencies for other pulsars. The custom noise model modifies the number of Fourier frequencies as given in table 3.3 and includes additional noise components (like system noise). However, the Fourier basis for CRN is always fixed at 30 frequency bins.

Here we present results with what was considered "standard" noise settings before this work, namely RN30DM30 model. We have started analysis using this model, and the search quickly converged to a particular sky position at 51 nHz. The first peculiarity of these results is that 51 nHz is very close to Venus orbital frequency, and the second is

B R N 30_D M 30+CGW R N 30_D M 30
= 18, which is not very high; however, it seriously triggered our attention by being relatively well constrained in the parameter space.

We have launched a set of investigations trying to understand this event. Using samples taken from the candidate's posterior, we have checked the contribution to the SNR from each pulsar. It came out that the main part comes from J1012+5307, but a few other pulsars (J1713+0747 and J0437-4715) also contributed not negligibly. We have checked that the zero contribution from a very good timer J1909-3744 is expected given the presumed sky position of the event.

We also conducted a set of injections of CGW signals with parameters taken from the candidate's posteriors using a simulated IPTA data (same TOAs as the real) with white noise only (RN and DM are supposed to be sub-dominant at the candidate frequency).

We could not reproduce the observed results with injections even when we increased the amplitude of the simulated signal.

BAYESEPHEM does not include a contribution from the inner planets to the phenomenological model, so we consulted a group from Observatoire de la Côte d'Azur (IN-POP group), inquiring if there could be an error in the Venus orbit picked up in the CGW analysis. We were reassured that the Venus parameters are known with very high preci-sion: this is a simple planet without any moons. However, looking into the future, we probably should extend BAYESEPHEM (or alternative mitigation models) by including the perturbation of orbital elements for the inner planets. We have performed the analysis with a narrow prior around this event with the model including the pulsar term. The Bayes factor has slightly increased B F = 21; however, the parameter estimation did not change appreciably. We have tried several runs with an extended model that includes the orbital eccentricity (again with a narrow prior). The eccentric runs show a very poor convergence, but all of them suggest a relatively high eccentricity. Results of both models (eccentric and with pulsar term) imply that the power is not localized at one particular frequency but spread over some finite-size frequency band.

We have turned to several noise models. We have started with varying number of Fourier frequencies used in the Gaussian process for RN and DM and tried M RN 30D M 100 and M R N 100D M 100 . The BF for CGW with those noise models has increased tremendously.

Finally, we have tried the custom model for the six best EPTA pulsars. Most notable is the peculiar noise model for J1012+5307 (see table 3.3) which, in combination with the sky position of the candidate event (being next to it), suggests that the explanation might be in the time-correlated high-frequency noise present in that pulsar (see Appendix of Chalumeau 2021). This would also be consistent with the results of eccentric runs suggesting that this could be extended in the frequency feature. We do not know the origin of that noise; it was found empirically. Including the custom noise model reduced the BF for this event to B cust om+CGW cust om = 0.95 and the posterior samples are not anymore constrained (see Figure 3.11 and, please, note that we have used different range for the parameters).

We have also computed the Bayes factor for the noise-only models (assuming that GW signal is weak) and found that B RN 30D M 30 cust om =< 10 -2 and similar result for RN30DM100 model suggesting that the data prefers by far the custom model. Custom made noise models for PPTA pulsars were studied in Goncharov et al. 2020. However, the peculiar behaviour that we have found in IPTA data was mainly originating from J1012+5307 which is not observed by PPTA. For that reason, we have chosen to focus on noise models from Chalumeau 2021. On the other hand, the noise model for J0437-4715 was based on Goncharov et al. 2020 where the spectral index for RN is of 3 and the optimal number of frequency bins for DM is (at least) 91. For this pulsar, we chose 30 frequency bins for RN as it was the recommended value for spectral index > 1.5 (see [START_REF] Goncharov | On the Evidence for a Common-spectrum Process in the Search for the Nanohertz Gravitational-wave Background with the Parkes Pulsar Timing Array[END_REF]). This was a useful exercise that triggered a set of investigations we would have to do in case of any CGW candidate. In addition, this section shows the importance of custom modelling noise for the best pulsars in the array, especially the noise at high frequencies, which is often partially neglected, assuming that it is dominated by the white (measurements) noise and it does not affect the search for the stochastic GW signal (which is most pronounced at low frequencies). The considered event shows how unmodelled log 10 h = 13.36 +0.12 high-frequency noise could conspire for CGW signal.

Conclusion

We have searched for a continuous GW signal in the IPTA DR2 dataset. We have used the Bayesian approach and based detection criteria on the Bayes factor. We have shown that using a custom noise model for the six best EPTA pulsars is essential for the correct interpretation of the data. This is especially true for J1012+5307, which exhibits timecorrelated noise at high frequencies. We found no CGW candidates using this noise model and proceeded to set the upper limit on GW strain. The addition of CRN in the noise model slightly affects the upper limit by lowering the sensitivity of the array at low frequencies. The most sensitive frequency appears to be around 10 nHz with a 95% sky averaged upper limit for CGW amplitude h 95 = 9.1 × 10 -15 . The IPTA DR2 shows a much better upper limit than previously set at higher frequencies, making it a promising dataset to detect CGW.

During the analysis, we demonstrated the CGW candidate follow-up investigations program, which was an important exercise that should be used in subsequent PTA CGW analysis. The expected CGW signal has low SNR, and its SNR will be slowly accumulated as we get more pulsars and a longer observational span. Modelling noise in pulsar data is essential, especially at high frequencies.

This analysis was limited to circular SMBHBs using only the Earth term. The use of eccentric CGW signal and including the pulsar term might potentially improve the search; however, it brings signal complexity which might make harder the interpretation of the results and increases the parameter space. We are entering the era of very high quality and high cadence radio observations with new instruments like FAST [START_REF] Hobbs | The role of FAST in pulsar timing arrays[END_REF] or SKA [START_REF] Stappers | The prospects of pulsar timing with new-generation radio telescopes and the Square Kilometre Array[END_REF] with sophisticated data analysis techniques. Additional investigations of the best approach to detecting CGW have to be re-investigated, probably using simulated data and/or an extended CGW signal injection campaign.

Introduction

We live in the era of large physics and astrophysics projects and often have to deal with large and complex datasets. The data analysis usually requires large computing facilities and a single computation could sometimes last for weeks. Optimizing the analysis tech- proposal there and its efficiency depends on the interplay of number (and distribution) of hot chains (more chains is better) and computational demands (increase with the number of chains). Understanding the properties of the signal and the likelihood surface often leads to a custom proposal suitable for a particular problem.

Here we suggest another generic proposal based on the Kernel Density Estimation (KDE). The idea of using KDE has already been explored [START_REF] Ashton | BILBY-MCMC: an MCMC sampler for gravitational-wave inference[END_REF][START_REF] Farr | Kernel-density-based parallel ensemble sampler Astrophysics Source Code Library[END_REF]). The main result of this paper is in a particular implementation of KDE itself and its embedding into a sampler. Even though the proposed method is very generic we will mainly discuss its implementations in the gravitational waves (GWs) data analysis.

Let us summarize the key points of the KDE-based proposal.

• KDE is used together with the set of other proposals in building a Markov chain.

We assume an adaptive approach where we use the data accumulated in a chain to regularly rebuild the KDE. We repeat adjustments until the convergence criteria based on the Kulback-Leibler divergence is satisfied.

• In order to build KDE we split all parameters into several groups, where parameters in each group show evidence of mutual correlation. The KDE-based proposal with parameter h specifying the bandwidth of the kernel. We use latin subscripts from the first half of the alphabet to enumerate the samples in the set. The main idea is to sum smooth kernel functions of x centered on each sample (input) data point X a . The overlaps between neighbouring kernels will add-up, shaping the PDF for the set of samples { X }.

The choice of the kernel is arbitrary and we choose to work with a gaussian kernel of the form:

K g (x -X a , h) = d i =1 exp -1 2 |x-X a | 2 i h 2 i 2πh i , (4.2)
where the h i is the local bandwidth corresponding to i -th parameter |x -X a | i , and d is the dimensionality of the parameter space. We use the latin letter from the second half of the alphabet to enumerate particular parameters, and the vector notation corresponds to a vector in the parameter space.

Optimal bandwidth

A KDE has one free parameter which we want to tune, the bandwidth h. Its value should be adapted to the dataset we are working with. There is no direct way of estimating it and we use an optimisation method [START_REF] Jin | An Improved Variable Kernel Density Estimator Based on L2 Regularization[END_REF]) based on the minimisation of the mean squared error (MSE) 2 with respect to h :

2 = d x( f (x, h) -f (x)) 2 , (4.3) ∂ 2 ∂ h = 0, (4.4)
where f (x) is the true pdf that we want to approximate with the KDE. The numerical way of evaluating this integral is given in the appendix 6.3.

Instead of using a global bandwidth h, we can define a local bandwidth h a for each kernel K (x -X a , h a ) [START_REF] Terrell | Variable Kernel Density Estimation[END_REF]. In that case, our KDE f (x, h) is :

f (x, h) = 1 N N -1 a=0 K (x -X a , h a ). (4.5)
Intuitively we expect the local bandwidth h a to be scaled according to the local density of points. Indeed, the bandwidth is chosen so that it is narrow in the regions of parameter space where the samples are most dense and it is broad where we have fewer samples.

This can ensure good interpolation and overlap between kernels, in particular, in high dimensional problem where the sample points are very sparse.

In practice we use solution of equation 4.4 to find optimal local bandwidth. The global bandwidth (if needed) could be defined as an average over all local values (by setting the input parameter global bw = True).

In Appendix 6.3, we have shown that each local bandwidth h a = [h a,1 h a,2 ...] of the kernel centered on point X a can be approximated by solving a linear system for the k nearest neighbours X b of the form :

A( X b )      1 h 2 a,1 1 h 2 a,2 . . .      = B (k) (4.6)
Matrix A and vector B are given in Appendix 6.3 eqn. 6.33, where we provide a detailed description of the method. They depend on the position of nearest neighbours which are the k near points contained in a hypercube centered on the point X a in the parameter space. The edge ∆X i of the hypercube for each parameter X i is defined as :

∆X i = (max { X } X i -min { X } X i )/s (4.7)
where maximisation and minimisation are performed over the set of input samples and s is a scaling parameter we call "adapt scale".

It might happen that a hypercube contains no point (besides the central). In that case, we cannot compute the local bandwidth and the point is discarded. Its bandwidth is later set to the global bandwidth as defined above. In case of high dimensionality and if the points in the dataset are very sparse, we change parameter s iteratively decreasing by a factor 2 until we find non-empty hypercubes. However, if this happens the evaluation of the bandwidth will probably be flawed and there is not much we can do about it, except use bigger datasets with more sample points. Often the amount of additional points that is needed to cover all "holes" could be very large incurring unmanageable computational cost. That is why good parameter grouping is essential : reducing dimensionality without loss of correlated features in the data. This will be the main subject of the next section.

Method

The main idea is to build a KDE for a given d-dimensional set of sample points { X }. However, for a high dimensional KDE we are strongly affected by the "curse of dimensionality" [START_REF] Scott | Multi-dimensional density estimation[END_REF] because the sample sets are often limited in size leaving undercovered regions of the parameter space. In addition the efficiency of KDE is degrading if there are too many points since we place the kernel on top of each sample. For that reason, we will split a d-dimensional parameter space into several low dimensional subspaces grouping the most correlated parameters together. We assume that sub-groups are not correlated and build the KDE for each of them fα (x α , h α ), so the total KDE is the product of the low-dimensional KDEs:

F (x, h) = α fα (x α , h α ), (4.8)
where the greek indices enumerate the subgroups of parameters. Forming these subgroups relies on the assessment of the correlation between parameters based on the provided set of samples { X } and this is the main subject of the next subsection.

Parameter grouping

For a d-dimensional dataset { X } we want to split the parameters in several sub-groups.

Each sub-group will contain correlated parameters while parameters from different subgroups will be uncorrelated. We could use a covariance matrix to identify correlations, however it implicitly assumes Gaussian distribution and cannot account for any complex 2D structures between pairs of parameters. Instead we use a method based on the Jensen-Shannon divergence (JSD) [START_REF] Nielsen | On the Jensen-Shannon Symmetrization of Distances Relying on Abstract Means[END_REF]. JSD, similarly to the Kullback-Leibler (KL) divergence, measures the similarity between two distributions, but with the advantage of being symmetric and bound 0 ≤ J SD ≤ ln 2, moreover it does not make any assumptions about the distributions. For each pair of parameter (X i , X j ) with the joined probability distribution p(X i , X j ) we compute 0 ≤ JSD p(X i , X j )||p(X i )p(X j ) ≤ ln 2, (

where the second distribution p(X i )p(X j ) is a product of one-dimensional (marginalized) distributions obtained by shuffling the parameters {X i , X j } that supposed to destroy any correlations between parameters X i , X j (see figure 4.1). If the JSD is low, it means that the shuffling did not affect the dataset and the parameters did not exhibit correlation. On the other hand, if the JSD is large, the shuffling did change something and the parameters should be grouped together.

We define a JSD threshold (we usually use 0.1) which we consider that two parameters are correlated. Starting from one parameter, we iterate the process to extract all correlated (chained) pairs. Once no additional correlated parameter is found, we take the union of all correlated pairs of parameters to form a sub-group. This process is illustrated in figure 4.2.

In case of multimodality of the probability distribution function that we try to reproduce with KDE, we implemented an additional (optional) feature: clustering samples before building KDE. This feature is especially welcome when the modes are separated by very low probability valleys. We cluster the samples (using k-means method MacQueen 1967) and apply KDE building approach described above to each mode. This does not change the fundamental structure of the KDE but it helps the bandwidth adaptation.

Turning KDE into proposal

It is desirable in several application to use posterior points inferred for some parameters into a prior for another investigation. Let us give several examples:

• The data inferred from electromagnetic observations in form of samples is used as a prior for the GW experiment. In this case we can either build a joined likelihood Figure 4.2: Illustration of the parameter grouping process using a JSD matrix. For this example, the JSD threshold is set to 0.25, hence the parameter sub-groups will be [0, 2, 3, 4] and [START_REF] Abbott | Prospects for observing and localizing gravitational-wave transients with Advanced LIGO, Advanced Virgo and KAGRA[END_REF]. Starting from parameter 0, we find that JSD for parameters 2 and 3 are above the threshold so they are both correlated with parameter 0. Then we check for parameters 2 and 3 and find that 4 is correlated with 3 while 2 sees no additional correlation. The last step would have been to check for 4 and find that there are no additional correlations.

Therefore, 0, 2, 3 and 4 will form a sub-group of correlated parameters. The parameter "1" is the last parameter that does not correlate with others (according to the adopted threshold) and it will be a sub-group on its own. or, alternatively, build KDE on the external posteriors and use it as a prior while analyzing the GW data.

• In Pulsar Timing Array (PTA) data analysis we often first investigate data acquired for each pulsar and trying to build a noise model. Later this pulsar and associated noise model is plugged into "Array" of pulsars for searching for a GW signal. It is proven to boost significantly efficiency of the GW search if we use posteriors for the noise model inferred in the first step as a proposal in the global fit later on.

• Often the data is taken continuously and we want to analyse it "on the fly"; that is true for GW data analysis. In this case you want to increment the data with a certain cadence while using information about the sources acquired from the analysis of the past data. One possibility is to turn again posterior built from the analysis of, say, first half a year of data into a prior in the analysis of the whole year of data.

• The method in Hee et al. 2016 gives a practical suggestion on how to compute the Bayes factor comparing several models without computing the evidence for each model. In this approach we introduce hyper-parameter indexing the models and jump in this parameter (say, within MCMC) which corresponds to jumping between the models. For this method to be robust the exploration within each model must be very efficient otherwise it will lead to very long poorly converging runs or to spurious results. If we have posteriors for each (or some) model available, we can turn them into proposal and use in the hyper-model exploration.

The main idea of grouping parameters in building KDE implies that we probably also want to make jumps within each subgroup (or in some subgroups) while keeping other parameters fixed (Gibbs-like sampling). The subgroups for a current jump are chosen randomly assuming equal probability attached to each subgroup. The randomness implies reversibility and equal probability is a warrant that we jump in all parameters evenly (on average) while performing low dimensional jumps. Let us denote the number of sub-KDEs that is used for each jump n kd e , then the proposal probability is

Fn kd e (x, h) =

n kd e α fα (x α , h α ), (4.10)
where the subscript in x α implies that we vary only parameters that belong to that (α) subgroup. This probability1 is used to balance the chain in the Metropolis-Hastings step of the MCMC algorithm [START_REF] Hastings | Monte Carlo Sampling Methods Using Markov Chains and Their Applications[END_REF]). Choosing a point from a given sub-KDE fα (x α , h α ) is done by drawing a point from the randomly chosen kernel K (x -X a , h a ) of fα (x α , h α ) centered on X a :

X * b → X a + N ( 0, h a ), (4.11)
where N ( 0, h a ) is a normally distributed random variable with 0 mean and covariance matrix diag(h a ) that is the bandwidth of kernel K (x -X a , h a ). For a random set of n kd e sub-KDEs, the newly proposed point X * is the union of parameters from each subgroup

X * b : X * = b X * b . ( 4
.12)

Adaptive proposal

In case we do not have samples from the previous investigations, we still can build KDEbased proposal using the points accepted by a running MCMC. There are several caveats which need to be considered: (i) during the burn-in and even some time after the distribution of the accepted points is quite unstable that will reflect on the KDE (ii) we are breaking the rules of MCMC, the chain is only asymptotically Markov, so that at some point we should fix the KDE-based proposal and dismiss all samples accumulated before. The rest of this subsection gives the details of the practical implementation of the adaptation.

To build a KDE on the currently sampled chain, we select only a subset of the total N samples. In particular, we take n s uniformly spaced points, from burn-in to the last point, as illustrated in Figure 4.3. The burn-in is taken to be a fraction, q, of the total chain length N and those points are dismisssed (hence, n s = (1q)N ). As chain evolves, the burn-in could also grow until we get to the stationary distribution. The post-burn-in length is also increasing, and, since we keep n s fixed, the space ∆n between the selected samples grows. We expect that the quality of the KDE improves with increase in ∆n because of the reduced correlation in the samples taken for building KDE, and that is true until ∆n reaches typical autocorrelation length of the chain (Hogg and Foreman-Mackey 2018).

We re-build a new KDE after each N ad apt iterations (jumps) of the chain. We want to track the evolution of KDE and stop adapting when it has reached stability (that could be another indicator of the burn-in phase). We compare the new (re-built) KDE F1 (x, h 1 )

with the old F0 (x, h 0 ) by computing the KL divergence (Kullback and Leibler 1951):

KL( F0 || F1 ) = d x F0 (x)(log F0 (x) -log F1 (x)) (4.13)
This integral could be approximated as

KL( F0 || F1 ) 1 n s a (log F0 ( X 0a , h 0 ) -log F1 ( X 0a , h 1 )), (4.14)
where { X 0 } is the set of n s samples used to build F0 (x, h 0 ). This gives a measure of change in KDE between successive updates: ∆KL. We stop updating KDE in order to preserve the ergodicity of the process (Atchadé and Rosenthal 2005) as soon as convergence criterion

| < ∆KL > | < KL 2 > < 5% (4.15)
is satisfied. The angular brackets denote the averaging over the last 5 updates and we demand that the average change in KL is small compared to the average KL values. ∆KL can be negative or positive, depending on the evolution of KL. If KL does not converge to a specific value, this condition ensures that it is at least oscillating around a mean value.

Results

We consider two datasets and perform search/parameter estimation using MCMC with KDE-based proposal.

In first application we consider a dataset from International PTA collaboration and perform the noise analysis for each pulsar in the array (Perera 2019). The likelihood is expected to be quite broad and unimodal but the dimensionality of the parameter space is large as well as its overall volume.

In second application we work with the simulated LISA data and search/characterize small bandwidth with several Galactic white-dwarf binaries. The likelihood in this case has more complex structure with quite strong correlation between parameters.

We quantify performance of the KDE-based proposal using the following criteria:

• Closeness between built KDE and the true distribution using KL divergence2 .

• The acceptance rate when the KDE is used as a proposal with MCMC sampler

• The autocorrelation length (Hogg and Foreman-Mackey 2018) of the MCMC chain when the KDE is used as proposal, to evaluate improvement in the mixing of the chain.

IPTA dataset

We build the KDE using n s = 10000 samples from the chains generated by previous MCMC runs. We adopt the following choice of parameters for generating KDE:

• js threshold = 0.1

• adapt scale = 10

• use kmeans = False Increasing n kd e leads to the higher dimensionality of the jumps and has a strong impact on the acceptance rate as shown in Figure 4.6. The best results are achieved if we perform jumps in one subgroup at the time, the acceptance rate decreases exponentially with n kd e .

A high acceptance rate is not necessarily a sign of a good proposal as it has to be paired with the low autocorrelation length. For each run, we compute the autocorrelation lengths of all parameters and compare the maximum, minimum and mean autocorrelation lengths. The low autocorrelation length implies that the samples drawn/accepted are independent. Results are presented in Table 4.4. KDE performs very well for IPTA data.

Minimum autocorrelation length does not seem to be affected much by the choice of proposal but the maximum is reduced by a factor 2 when using KDE. Reducing the maximum is the most important because thinning the chain by this factor ensures that all samples are independent for all parameters. The mean autocorrelation length is just an indicator of the average performance of the proposal.

Like for the acceptance rate, we check the influence of n kd e on the autocorrelation length. Results are given in Table 4 ceptance rate of 0.22 that is very close to the expected optimal acceptance rate of 0.234 [START_REF] Gelman | Weak convergence and optimal scaling of random walk Metropolis algorithms[END_REF]). Based on these results, we recommended to work with values of n kd e between 1 and 10, especially when the dimensionality of the parameter space is high like in this IPTA example.

LISA dataset

Now we turn our attention to the simulated LISA data, in particular, we use "Sangria" dataset which is 1 year long, it contains about a dozen of merging massive black hole binaries and about 30 millions of Galactic binaries https://lisa-ldc.lal.in2p3.fr/.

Here we are interested in Galactic binaries in the very narrow frequency range around 4 mHz and we have removed all merging black holes. We have detected 6 sources in that frequency interval and we use KDE-based proposal together SCAM and DE 3 . This time we use home-made sampler M 3C 2 (https://gitlab.in2p3.fr/lisa-apc/mc3), specifically parallel tempering version of it. We describe this sample in details in the Appendix 6.4. This sampler is using Metropolis-Hastings acceptance-rejection step as well as slice sampling [START_REF] Neal | Slice sampling[END_REF].

Each Galactic binary is characterized by 8 parameters, so we have in total 48 parameters [START_REF] Littenberg | A detection pipeline for galactic binaries in LISA data[END_REF]). We expect some parameters (like amplitude and the orbital inclination angle) to correlate for each source, and, in addition, some parameters could correlate between the sources. Here we try to build a proposal on-the-fly. The likelihood surface for this problem is rather complex having many well separated maxima (reason for using parallel tempering). We will build KDE as we accumulate samples: adapting KDE proposal with the rate every 5000 samples and using n s = 5000 samples for each chain. Besides KDE we also use SCAM proposal and slice sampling. Note that the combination of SCAM jump with Metropolis-Hastings sampling and slice (being very independent) already significantly reduce the autocorrelation of the accepted points.

First, we consider the convergence of KDE adaptation. During the burn-in stage the KDE is changing quite violently. The correlation between parameters is quite unstable which leads to fluctuation in how parameters are grouped and in the number of subgroups. As burn-in proceeds we keep track of the grouping and fix the splitting in subgroups as soon as it stabilizes (when the same grouping appears at least 5 times). Once we have fixed sub-groups we check the KL divergence between the subsequent updates of KDE (as described above). The results of KDE adaptation are presented in As a next step we want to check the acceptance rate of the KDE-based proposal and compare it to SCAM and slice. Note that slice is not based on the Metropolis-Hastings acceptance/rejection algorithm, however we still can introduce an effective acceptance rate as a ratio of total number of slice calls to the total number of likelihood evaluations used by slice.

We consider two cases for the subgroup jumps n kd e = 1 and n kd e = 5. Figure 4.9 compares acceptance rate for n kd e = 1. One can see that it stabilizes around 0.31 very fast, and, despite that it is lower than what we had for the PTA application, it is a very decent acceptance rate. SCAM has a similar acceptance (but usually longer autocorrelation length), while slice is worse by a factor 10 (though it usually has low autocorrelation length). We observe that the mean value is slightly (about 17%) lower when we include KDE-based proposal and the maximum length remains the same. As we have already mentioned, mixing SCAM with Metropolis-Hastings and slice steps does reduce the autocorrelation already (compared to PTA example where we did not use slice sampling). In addition we use parallel tempering algorithm, where the hot chains could be seen as yet another jump proposal. All in all, KDE does not add much to already reduced autocorrelation run in the current analysis.

Conclusion

Bayesian formalism is a usual approach in nowadays gravitational waves data analysis.

The inference of the parameters posterior distribution is often done using MCMC and the efficiency strongly depends on the proposal it uses. In this article we have presented KDEbased proposal which can be either built on-the-fly during the extended burn-in stage or constructed using posterior points from another run. We tested this proposal by running MCMC on IPTA data, using a KDE that we have built from previous MCMC runs (i.e. non adaptive case). The advantage of KDE-based proposal was clearly seen in high acceptance rate with low autocorrelation length. We have found that using rather low value n kd e = 1 -5 (number of subgroups used in the jump simultaneously) seems to be optimal.

Another application of the KDE-based proposal was in running PTMCMC on the simulated LISA data searching for Galactic white dwarf binaries in a narrow frequency band.

In this case we have built KDE adaptively during an extended burn-in stage. The addition of KDE-proposal to the sampling had only a moderate impact: it shows a decent acceptance rate (about 31%) with only small improvement in the autocorrelation length.

Moreover we have shown that low-dimensional jumps are strongly preferred. Few things could be improved, most notably in the adaptation. The threshold for grouping parameters was chosen somewhat ad hoc, and correlation of some parameters could be close to the threshold. We did observe the fluctuation in choosing the sub- groups during the adaptation. One possibility could be to choose not one but two plausible grouping, build KDE for each and use two KDE proposals in a probabilistic manner.

The criteria for stopping adaptation was also chosen somewhat arbitrary, and might benefit from the further tuning. Finally, we should implement an adaptive tuning for optimal n kd e based on the acceptance rate.

Averaged signal to noise ratio

In this section, we propose to develop the analytical averaging of the signal to noise ratio formula presented in section 3.3.1. This formula is developed fromg the expression of the timing residuals induced in PTA by a CGW that was demonstrated in section 1.4.4. This calculation is rather trivial. We use for the timing residuals (only Earth term) The signal to noise ratio is, using the definition of the noise weighted inner product : ρ 2 = (s|s) = s T Σ -1 s (6.6)

s = F + s + + F × s × , ( 6
We will denote the averaging of a quantity A with respect to some parameter p with brackets and index < A|A > p . For the signal to noise ratio ρ (F + ) 2 + (F × ) 2 I 2 + J 2 (sin ωt | sin ωt ) + (cos ωt | cos ωt ) . (6.19) We now compute the average of the I 2 + J 2 term with respect to the inclination ι This partially averaged, sky location dependent SNR formula can be used to rank the pulsars according to frequency of CGW and sky location. However, we want the fully averaged ρ 2 so we need to express the antenna pattern functions in terms of the CGW sky location parameters To average over sky location parameters, we place ourselves in the frame where the earth-pulsar axis unit vector p = (0, 0, 1), which makes the calculation easy and gives :

< I 2 + J 2 > ι = 16 5 , (6.20
F + ( Ω) = 1 2 ( m • p) 2 -( n • p) 2 1 + Ω • p (6.
< (F + ) 2 + (F × ) 2 > θ,φ = 1 3 . ( 6

Kullback-Leibler divergence

We test the quality of our KDE by computing the KL divergence for each subgroup fα (x α , h α ).

Because we cannot know the true pdf for a set of samples { X }, we have to resort to binning. From sub-KDE fα (x α , h α ) built on subgroup of parameter samples { X α }, we draw a new set of samples { X * α }. We estimate the corresponding normalized histogram distributions P α and P * α using same grid of N bins to compute the KL divergence as Kullback and Leibler 1951 :

KL(P α ||P * α ) = N i =0
P α,i ln P α,iln P * α,i . (

To avoid divergence of the logarithms, we set every P α,i and P * α,i that are equal to 0 to the minimum found value in P α P * α that is not 0.

Optimal bandwidth of KDE

.

We start this appendix with defining few useful expressions that will be used in our derivations later.

• The overlap between two neighbouring kernels of same bandwidth is given by

d xK (x -X a , h)K (x -X b , h) = d i =1 exp -1 4 | X a -X b | 2 i h 2 i π2h i . ( 6.28)
Let us remind you that d is the dimensionality.

• If a set of samples { X } of size N drawn from the probability density function f (x), then we can approximate the averaging integral as :

d x f (x)g (x) 1
N a g (X a ) , (6.29) where the function g is evaluated at the sample points X a .

The main objective of this Appendix is to derive the optimal local bandwidth which is defined through the minimization of the mean square error: 2 πh a,i , (

where h a,i is i-th component of bandwidth attached to Gaussian kernel at X a and ∆X ab,i = ( X a -X b ) i is i-th component of the a vector connecting two samples in the parameter space.

Using now second bullet equation and excluding the actual sample from the sum (for improving stability and removing the possible bias, see "leave one out estimator" Jin et al.

2021

) we obtain for the second term As an alternative approach we can assume that the bandwidth is comparable to the distance to the neighnbours and define h 2 c, j = ∆X 2 c, j (1+ε c, j ), where ∆X 2 c, j = 1/k near b ∆X 2 bc, j is the average square distance (i -th component) to the points in vicinity of X c and assume that ∆X 2 bc, j /∆X 2 c, j ∼ 1 and ε c, j 1 for all b, c, j . This yields Solving this system at each point X c for each direction in the parameter space ( j ) gives us the desired local bandwidth h c .

The M3C2 sampler

The DE could be used as a proposal (snooker) described in Jasper A. Vrugt 2008 using multiple chains running in parallel or using the accumulated samples for each chain to propose the jump. There is no much difference between those two ways in case of well converged chains. However, the behavior and efficiency of those two implementations is very different during the burn-in stage.

In slice sampling, we use slicing of the parameter space either randomly or along the eigen directions of the covariance matrix, the choice is made with a probability set by user. In case of mixture of slice and MH, the frequency of each method is defined by a user specified weight. In addition to preset proposals available in M3C2, user can add custom jump-proposals using a common interface. The weights and proposals can be set individually for each running chain.

Besides running parallel independent chains, M3C2 sampler has also parallel tempering implementation with an adaptive temperature ladder following Vousden et al. 2016.

The adaptation is aiming at increasing the acceptance rate between the chains.

The multi-chain scheme of M3C2 can be easily deployed on the multicore CPU infrastructure. Data exchange between chains, in case of parallel tempering, is restricted to its minimum level (pairwise communication between the chains), to ensure a good scalability.

Pulsar ranking for correlated signals

In L. Speri 2022, several methods are presented for the ranking of pulsars according to their contribution to the total SNR of the array, for a given correlated signal. In this manuscript, we mention SNR-B and Chimera when we cross-check the ranking obtained with the CGW method to support its validity. Here we will briefly explain what the methods consist in.

SNR-B

The SNR-B method computes the optimal statistics that maximize the SNR of the array in the presence of correlated signal [START_REF] Rosado | Expected properties of the first gravitational wave signal detected with pulsar timing arrays[END_REF]. It is given by :

SNR 2 B = 2 a<b Γ 2 ab S 2 ( f )T ab P a ( f ) + S( f ) P b ( f ) + S( f ) + S 2 ( f )Γ 2 ab , (6.35)
where Γ ab is the correlation coefficient between pulsars a and b, T ab is the overlapping time of observation between pulsars a and b, S( f ) is the PSD of the expected signal and P a ( f ), P b ( f ) are the PSD of the intrinsic noises of pulsars (red noise, dispersions measure, etc).

Coupling matrix

The coupling matrix formalism seeks to quantify the overlap between different spatial correlation patterns between pulsars using spherical harmonics Y l m . Indeed, l = 0, 1, 2 describe respectively a monopole, a dipole and a quadrupole, corresponding to the correlations expected for the three main types of correlated signals in PTA (see section 2.2.5).

We decompose the sky distribution of pulsars on a spherical harmonic basis [START_REF] Efstathiou | Myths and truths concerning estimation of power spectra: the case for a hybrid estimator[END_REF]) and evaluate the overlap between the coefficients as : The weight w a is chosen according to the RMS value of pulsar a residuals (Roebber 2019).

K (l
For a set of pulsars, we evaluate the leakage between components with the ratio δ λ = λ mi n /λ max , λ being the eigenvalues of the matrix M l l obtained by averaging of K (l m)(l m)

over the components m. Selecting the largest ratio δ λ ensures that we have the best possible decoupling between different l .

Chimera

The Chimera method is a combination of the two previous methods. We define an effective score SC C hi mer a of SNR-B and δ λ (coupling matrix) according to which we rank the pulsars :

SC C hi mer a = SNR B δ λ . (6.38) In L. Speri 2022, the 25 best pulsars acquired with this method are saved as best pulsars of the array.

  tections. Le dernier chapitre présente une nouvelle méthode basée sur les Estimations par Noyaux (Kernel Density Estimation ou KDE) qui a été développée durant mon doctorat pour optimiser les techniques d'analyse de donnée pour la recherche d'OGs et en particulier pour améliorer le temps effectif d'échantillonage des méthodes MCMC. Le KDE est une méthode d'estimation de densité de probabilité des paramètres d'un ensemble de donnée. Pour estimer ces densités de probabilité, il est possible d'utiliser des histogrammes. Contrairement à ces derniers, le KDE a l'avantage d'être une fonction continue et définie dans tout l'espace. Pour qu'il décrive fidèlement la donnée, il faut l'adapter à celle-ci à l'aide de méthodes d'optimisation parfois coûteuses en temps de calcul. Nous présentons ici une nouvelle méthode d'optimisation rapide par linéarisation d'un problème de minimisation d'erreur quadratique. Ensuite, nous présentons une nouvelle méthode d'estimation des corrélations entre paires de paramètres basée sur la divergence de Jensen-Shannon à partir de laquelle un algorithme de groupement des paramètres est développé. Les KDEs produits avec ces méthodes sont utilisés en MCMC pour proposer de nouveaux points dans les régions intéressantes de l'espace des paramètres pour aider l'échantillonage. D'abord, nous avons testé les KDEs avec la donnée IPTA DR2 pour la recherche d'OGCs. Dans ce cas, le temps d'échantillonage effectif était réduit d'un facteur 2 ou 3 ce qui indiquait un bon fonctionnement de la méthode. Par la suite, nous avons testé une méthode adaptative du KDE avec de la donnée simulée de LISA (futur observatoire spatial d'OG) pour la recherche de signaux gravitationnels de binaires galactiques (systèmes binaires de naines blanches). Elle s'est avérée moins efficace avec la donnée LISA pour laquelle le temps effectif d'échantillonage n'était que légèrement réduit (de 17%).

Mots clés:

  Onde gravitationnelle -Pulsars -Trou noir supermassif -Analyse de données Ce manuscript est centré sur la détection d'Ondes Gravitationnelles (OG) basse fréquence en utilisant les données de chronométrage haute précision d'un réseau de pulsars milliseconde obtenues à l'aide de radiotéléscopes. Le premier chapitre donne un apercu du bagage théorique nécessaire pour comprendre l'origine des OGs. Après avoir présenté les bases de la Relativité Générale (RG), nous argumentons qu'un réseau de pulsars chronométrés (Pulsar Timing Array ou PTA) pourrait détecter des OGs générées par des systèmes de Trous Noirs Super Massifs Binaires (TNSMB) et nous dérivons une expression mathématique pour les résidus temporels qui seraient induits dans la donnée chronométrée. Lorsqu'un tel signal est stable et monochromatique, on le qualifie d'Onde Gravitationnelle Continue (OGC). Dans le second chapitre, nous présentons les méthodes d'analyse de donnée pour la détection d'OGs dans PTA. Nous expliquons brièvement l'inférence Bayesienne et nous détaillons comment les bruits intrinsèques aux pulsars peuvent être modélisés. Nous montrons que pour estimer la significativité statistique de nos modèles de signal et de bruits, nous calculons le Facteur de Bayes (FB) en utilisant des hypermodèles et des méthodes de Monte-Carlo par chaînes de Markov (MCMC). Dans le chapite trois, nous discutons des résultats obtenus avec le second ensemble de données de la collaboration PTA internationale (IPTA DR2) qui combine les données de : l'observatoire d'OG nanohertz d'Amérique du Nord (NANOGrav), le PTA Européen (EPTA) et le Parkes PTA (PPTA) en Australie. Nous avons cherché des OGC mais n'avons trouvé aucune preuve de leur présence dans la donnée. Nous avons établi une limite haute de détection pour l'amplitude des OGCs et avons trouvé que IPTA DR2 donne la meilleure sensibilité à ce jour pour des signaux dont la fréquence est plus haute que 11 nHz. Il a été découvert il y a peu qu'un bruit en basse fréquence et commun à tous les pulsars était présent dans la donnée. L'inclusion de ce dernier dans le modèle détériore la sensibilité aux OGCs du réseau. De plus, nous montrons qu'une mauvaise modélisation des bruits intrisèques aux pulsars peut mener à une fausse détection. Le dernier chapitre présente une nouvelle méthode basée sur les Estimations par Noyaux (Kernel Density Estimation ou KDE) qui a été développée durant mon doctorat pour optimiser les techniques d'analyse de donnée pour la recherche d'OGs. Cette méthode est utilisée avec MCMC pour améliorer son temps effectif d'échantillonage. Elle a été testée avec IPTA DR2 pour la recherche d'OGCs puis avec Sangria, une donnée simulée du futur projet d'observatoire spatial d'OG (le Laser Interferometer Space Antenna ou LISA), pour la recherche de Binaires Galactiques (BG). Nous démontrons que la méthode est efficace pour IPTA DR2 mais est moins performante avec LISA. Mots clés: Onde gravitationnelle -Pulsars -Trou noir supermassif -Analyse de données ABSTRACT : ENGLISH This manuscript focuses on the detection of low frequency Gravitational Waves (GW) using high precision timing observations of millisecond pulsars collected with radio telescopes.

: 1 THEORETICAL BACKGROUND Contents 1 . 1 1 . 3

 11113 Gravitational wave -Pulsar -Super massive black hole -Data analysis CHAPTER General relativity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.1 The notion of curved space . . . . . . . . . . . . . . . . . . . . . . . . 1.1.2 Mathematical description of a curved space . . . . . . . . . . . . . . 1.2 Gravitational waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Linearized Einstein field equations . . . . . . . . . . . . . . . . . . . 1.2.2 Polarization tensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.3 Quadrupole formula . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.4 Circular binary system . . . . . . . . . . . . . . . . . . . . . . . . . . . Supermassive black holes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1 Gravitational wave background . . . . . . . . . . . . . . . . . . . . . 1.3.2 Continuous waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Pulsar Timing Array . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.1 Sensitivity of the array . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.2 Millisecond pulsars . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.3 The International Pulsar Timing Array . . . . . . . . . . . . . . . . . 1.4.4 PTA response to GW signal . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Multi-band gravitational wave astronomy . . . . . . . . . . . . . . . . . . 1.5.1 Ground based detectors : LIGO/Virgo . . . . . . . . . . . . . . . . . . 1.5.2 Space based detector : LISA . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 . 1 :

 11 Figure 1.1: Illustration of the Ehrenfest paradox with on the left, a rotating disk with angular velocity ω seen from above. When it rotates, it should see its circumference C decrease due to relativistic length contraction while its radius R is unaffected. On the right, we have the resolution of the paradox showing that a curved geometry allows smaller circumference without change in radius R.

Figure 1 . 2 :

 12 Figure 1.2: The spherical surface of the the Earth (positive curvature) is fundamentally different from a plane (null curvature). Every planisphere of the Earth that we have are just approximations where length measurements are distorted.

Figure 1 . 3 :

 13 Figure 1.3: We can associate a tangent circle of radius R to each point of s(t). The inverse of the radius defines the curvature. This is the simplest way to think about curvature in 1D.

Figure 1 . 4 :

 14 Figure 1.4: In curved space, basis vectors e µ are function of coordinates x µ and straight lines are parametrized geodesics X µ (λ).

Figure 1 . 5 :

 15 Figure 1.5: The Riemann tensor contains the components of the deviation vectors R for all possible transportation of basis vector e α on a closed loop in directions µ and ν.

ψFigure 1 . 6 :

 16 Figure 1.6: Illustration of the geometry of the SMBHB. Here, L represents the rotation axis of the binary, k the direction of propagation of the GW with polarization basis { ê1 , ê2 }, the polarization angle ψ and finally ι the inclination (i.e. the angle between L and k).

Figure 1 . 7 :

 17 Figure 1.7: First images of SMBH captured by the Event Horizon Telescope.

Figure 1 . 8 :

 18 Figure 1.8: Illustration of the possible scenarios for MBH formation. Z is the metallicity of the halo, i.e. the presence of elements heavier than helium (figure taken from Volonteri 2010).

Figure 1 . 9 :

 19 Figure 1.9: A simulation of SMBHB population in the universe where y-axis is GW amplitude and x-axis is GW frequency. Each blue triangle represents a single SMBHB source. The total sum of these signals gives the characteristic strain of the GWB (black line) which can be approximated by a powerlaw (orange line) of index -2/3. The black stars are individual binaries resolved on top of the background (CGW). The dotted line displays the sensitivity of IPTA, showing that both GWB and CGW should be detectable at low frequencies. (Acknowledgement : we have used simulated population of binaries from A.Sesana to produce this plot.)

Figure 1 . 11 :

 111 Figure 1.11: The evolution of frequency ∆ f g w for 10 years of observation and three chirp masses compared to PTA resolution.

Figure 1 . 12 :

 112 Figure 1.12: The population of known pulsars given their spin period (x-axis) and spin period first time derivative (y-axis) (plot taken from Tiengo et al. 2011). MSPs lie on the bottom left corner after getting recycled. The dashed blue line is the limit of the pulsar graveyard below which no pulsar (or very few) is expected to radiate because of electronpositron pair production (Chen and Ruderman 1993).

Figure 1 . 13 :

 113 Figure 1.13: The list of IPTA telescopes and their locations in the various countries of the collaboration (picture taken from the NANOGRAV website).

Figure 1 . 14 :

 114 Figure 1.14: The geometry of the PTA in the radio observatory (Earth) frame. The SMBHB emitting a GW signal in the direction k is going to influence the TOAs of the pulsar observed on Earth. The angle between k and the unit vector n pointing to the sky position of the pulsar seen from Earth determines how the GW is felt by the pulsar according to the antenna pattern fuction F +,× . The distance Earth-pulsar L affects the delay in the pulsar term.

( 1 .Figure 1 . 15 :

 1115 Figure 1.15: Distribution of ∆ f α as shown in Sesana and Vecchio 2010. The vertical dotted line at ∼3 nHz represents the frequency resolution of the array.

Figure 1 .

 1 Figure 1.16: GW detectors sensitivities from Thermal noise in BEC-phononic gravitational wave detectors https://www.researchgate.net/publication/275974267_Thermal_ noise_in_BEC-phononic_gravitational_wave_detectors

  pairs are different for each radiometer backend k and are applied as global rescaling parameters to every TOA uncertainty σ T O A,k coming from backend k.

  The red noise is a time-correlated noise intrinsic to each pulsar dominating at low frequencies (Shannon and Cordes 2010, Haasteren and Levin 2012). It takes into account possible stochasticity in pulsar rotation. We describe it as a Gaussian process decomposed on a Fourier basis with finite number of components. Its power spectral density is expressed as :

  the case of an individual red noise, different values of A RN and γ RN are associated with each pulsar. To describe common red noise processes, we set the same values of A C R N and γ C RN for all pulsars, meaning that they share a red noise with same PSD. Lately, such common processes have been discovered by all PTA collaborations S Chen et al. 2021; Goncharov et al. 2021; Antoniadis 2022; Arzoumanian et al. 2020.

  ) with < . > denoting the averaging over the noise realizations, S n ( f L , θ H P ) the PSD of the modeled noise evaluated at Fourier frequency f L for hyper-parameters θ H P corresponding to, in the case of RN and DM, the amplitudes A RN , A D M and spectral indices γ R N , γ D M . Then, in the "weight-space" view of Gaussian processes (Haasteren and Vallisneri 2014), the Fourier coefficients a follow a Gaussian prior probability of :

  .20) with I the identity matrix. In that approximation, the are Gaussian weights and the design matrix M plays the role of basis functions. This method was first shown in Haasteren and Vallisneri 2014, allowing marginalisation on all parameters at once and giving the marginalised posterior distribution where the noise hyper-parameters θ H P (the amplitude and spectral indices of the Gaussian processes) appear in the covariance matrix C of the new marginalised likelihood L( δt α | θ, C α ) as

. 30 )Figure 2 . 1 :

 3021 Figure 2.1: Hellings and Down correlation with χ a function of angular separation ζ between two pulsars a and b (plot taken from Jenet and Romano 2014)

  When including correlated signals to the model, the covariance matrix presented above is no longer block diagonal. The off diagonal terms can follow various correlation patterns χ ab that are function of the angular separation ζ between pulsars.Hellings and Down PTA collaborations mainly target quadrupolar correlation patternsbecause they agree with the presence of GW signal. We can recover such correlations by calculating the overlap between the antenna pattern functions of two pulsars with angular separation ζ. This yields the Hellings and Down (HD) curve that is expected to be the signature of the GWB (see[START_REF] Hellings | Upper limits on the isotropic gravitational radiation background from pulsar timing analysis[END_REF].

. 34 )

 34 One possible origin for monopolar features is clock correction errors. The atomic clock systems in the radio-telescopes allowing local time measurements can accumulate errors that must be corrected in order to keep them synchronized. If the correction is not done correctly, it can cause a low frequency noise in the timing residuals of the TOAs measured with the affected radio-telescope[START_REF] Tiburzi | A study of spatial correlations in pulsar timing array data[END_REF]).

  bility distributions of individual parameters because of the tremendous amount of computation time that it would represent. For illustration, if we want to plot a posterior probability distribution p( θ) of a d-dimensional parameter space on a grid of N points per parameter, we would need N d evaluations of the probability function p( θ). Assume that one computation of p( θ) is fast (∼ 10 -3 s) and we want a grid of N = 10 points (low resolution) for a 10-dimensional parameter space. This plot would take : 10 -3 × 10 10 = 10 7 s 3 months and 3 weeks,(2.35) which is ridiculously long for such a low resolution plot. To tackle the issue, numerical algorithms were proposed, one of the most popular being Markov Chain Monte Carlo (MCMC). It allows to sample points from a target distribution p( θ) while focusing on the interesting regions, in a limited number of evaluations of the function. Before presenting the sampling algorithm, we must explain what a Markov process is in order to justify the motivation behind it. For more details, see[START_REF] Brooks | Handbook of markov chain monte carlo CRC press[END_REF]

Figure 2 . 2 :

 22 Figure 2.2: Illustration of the hypermodel method. Here, the hypermodel is the union of two different models A and B (red and blue) with posterior distributions P A and P B functions of parameters θ A and θ B . As the sampling goes on, the sampler can jump between models by selecting which posterior is active. It will spend more time on a model if it is favoured by the Metropolis-Hasting choice (model A in this example). The time spent on models A and B corresponds respectively to the total number of samples n A and n B .

(

  GW) signals in the nanohertz band, where we expect to see a gravitational wave background (GWB) produced by the superposition of GW signals from the population of supermassive black hole binaries (SMBHBs) (Maiorano et al. 2021). Some individual SMB-HBs might be brighter than most and stand above the stochastic signal; those are individually resolvable sources. The binaries detectable in the PTA band are in the orbits with the period from a few months to a few years and emit almost monochromatic GWs continuously during decades; we refer to those signals as continuous GWs (CGWs) (Ellis et al. 2012; Babak 2016; Aggarwal 2019; Corbin and Cornish 2010). The GWs affect propagation of the radio emission from millisecond pulsars leaving an imprint in the time of arrival (TOA) of pulses observed with the radio telescopes. CGWs impact TOAs from all observed millisecond pulsars in a deterministic manner characterized by parameters of the SMBHBs. In this work, we consider the data combined by the International Pulsar Timing Array (IPTA). IPTA is a consortium of NANOGRAV (NANOGrav Collaboration et al. 2015), European Pulsar Timing Array (EPTA) (Desvignes et al. 2016), Australian (PPTA) (Manchester et al. 2013b) and Indian Pulsar Timing Array (InPTA) (Tarafdar et al. 2022) collaborations. In particular, we analyze the second data released by IPTA (IPTA DR2) described in details in Perera 2019. Recently, PTA collaborations have reported on the discovery of the common red noise signal, that is the stochastic signal with the spectral shape common to all pulsars in the array. Its high statistical significance was demonstrated independently by three collaborations (Arzoumanian et al. 2020; S Chen et al. 2021; Goncharov et al. 2021) and, with even higher statistical confidence, was assessed using the IPTA DR2 (Antoniadis 2022).

. 1 )

 1 with N α the uncorrelated white noise, C T M ,α the term associated with the marginalized timing model, C D M ,α the dispersion measure noise, C RN ,α the individual red noise and C C R N the common red noise.

  F A are the antenna pattern functions (Babak and Sesana 2012; Sesana and Vecchio 2010; Ellis et al. 2012; Taylor et al. 2016) given as

Figure 3 . 2 :

 32 Figure 3.2: Normailzed ρ2 a of the five best pulsars of the IPTA DR2, at different CGW frequencies. The glitches at the right of the plots are due to the one year and half-year peaks.

Figure 3 . 3 :

 33 Figure 3.3: Distribution of the SNR 2 coverage for 1000 different sets of CGW parameters.

Figure 3 . 4 shows

 34 that amplitude h( f ) now becomes a broad distribution. It appears less constrained at low frequencies due to the stronger presence of RN and DM noise. The sensitivity we obtain with fixed maximum likelihood value of noise parameters seems to give a better sensitivity than the average distributions of h around 10nHz. Therefore, in the case of PTA, taking into account the variability of the noise parameters has a direct influence on the detectability of a signal.

Figure 3 . 4 :Figure 3 . 5 :

 3435 Figure 3.4: Binned sensitivity curve for 100 different realization of noise parameters. The grey colour scale represents the binning of h at a given frequency. The dashed line represents the sensitivity curve for maximum likelihood values of parameters.

- 7

 7 Hz during the analysis, as discussed in Petiteau et al. 2013. As mentioned above, we compute the upper limit for two models of noise; with and without CRN. The inclusion of CRN adds two more parameters to the model (an amplitude parameter A C RN and a spectral index γ C RN ), which we sample together with parameters of CGW and other noise components.

log 10 h 95 Figure 3 . 6 :

 9536 Figure 3.6: 2D sky sensitivity of IPTA for CGW signal around most sensitive frequency.

BFFigure 3 . 7 :

 37 Figure 3.7: On the top panel, the upper bound of the 95% central credible region on log h obtained with a log-uniform prior (with CRN). Bottom panel is the associated BF at given CGW frequency. The black dashed line shows where BF is equal to 1. The two peaks around 35 and 70 nHz with BF>1 are the 1yr and 2yr and should not be taken into account.

Figure 3 . 8 :Table 3 . 3 :Figure 3 . 9 :

 383339 Figure 3.8: CRN parameters for IPTA DR2 without CGW signal obtained in (orange) and with CGW signal at 1 nHz (blue).

Figure 3 . 10 :

 310 Figure 3.10: Posterior distribution for log h, cos θ and φ using noise model M RN 30 D M 30 . The black star indicates the sky position of J1012+5307. The dashed lines represent the median values of the parameters.

Figure 3 . 11 :

 311 Figure 3.11: Posterior distribution for log h, cos θ and φ using noise model M cust om . The black star indicates the sky position of J1012+5307. The dashed lines represent the median values of the parameters.

  niques and pipelines is then a key challenge of the data science associated with all large (astro)physical experiments.Nowadays, it is quite common to use Bayesian framework for analysing the data. It is especially convenient if we have a parameterized data model (or several competing models) describing the data. In this approach we treat all parameters as random variables with some prior based either on some physical principles or informed from the previous independent experiments. We use the observations at hands to refine our prior knowledge and infer a posterior probability distribution function for parameters of a model, or even perform a selection among several models. Often we have to deal with multidimensional parameter space with a non-trivial likelihood function which can be evaluated only numerically. One of the most used tool to perform the numerical sampling from a target distribution is Markov Chain Monte Carlo (MCMC). Building a Markov chain that represents the desired posterior requires two key ingredients: (i) proposal suggesting how to choose point X i +1 given the last point in the chain X i ; (ii) the detailed balance which ensures the reversibility of the chain. One of the most successful and frequently use proposal approach is to use parallel tempering: running several chains with logarithmic distributed temperature ladder (see for example[START_REF] Earl | Parallel tempering: Theory, applications, and new perspectives[END_REF]. The hot chains play role of a

Figure 4 . 1 :

 41 Figure 4.1: We plot three examples of datasets where we have on the left no correlation, in the middle linear correlation and on the right more elaborate features. On the top panels we have the corresponding values of the correlation coefficient based on simple evaluation of the covariance matrix. We see that it excels at finding the linear correlation but completely fails with the right panel features. On the bottom panels we have the same three datasets in red with their corresponding shuffled version that destroys correlations in blue. While the left panel remains unchanged, the others are affected and it is captured by the JS divergence.

Figure 4 . 3 :

 43 Figure 4.3: For a chain of total length N , we get rid of the burn-in, then we extract n s linearly space samples from the remaining fraction of the chain. These n s points are used to build the KDE.

Figure 4 . 4 :

 44 Figure 4.4: Distribution of the KL divergence for all sub-KDEs. 1-d sub-KDEs are of best quality.

Figure 4 . 7 :

 47 Figure 4.7: Mean, maximum and minimum autocorrelation lengths for several values of n kd e .

Figure 4 . 8 .

 48 The consistent grouping of parameters was reached after 26 updates as indicated by a dashed red line. Then we compute KL after each update and stop adaptation once the condition 4.15 is met (see the right panel ofFigure 4.8). From then on we keep KDE fixed and perform the actual sampling.

Figure 4 .

 4 Figure 4.10 compares acceptance rate for n kd e = 5. Increase in the dimensionality of the jumps has a drastic effect on the acceptance rate, its value drops to about 0.015. It is also interesting to compare behaviour of SCAM and slice for two runs. Slice shows very stable/consistent results, while SCAM has significant fluctuations though preserving the trend. The two-dimensional KDE jumps seem to be the best option in this application. Next we check the autocorrelation length when we run with and without KDE-based proposal. We restrict ourselves with the case n kd e = 1 since n kd e = 5 has a very poor acceptance rate.Figure 4.11 compares maximum and mean autocorrelation of two runs.

Figure 4 .

 4 [START_REF] Babak | European Pulsar Timing Array limits on continuous gravitational waves from individual supermassive black hole binaries[END_REF] compares maximum and mean autocorrelation of two runs.

Figure 4 . 8 :

 48 Figure 4.8: Evolution of KL and ∆KL. Left panel, evolution of KL before fixing parameter groups. Right panel, evolution of KL after fixing parameter groups, we start computing ∆KL 5 updates after the grouping was fixed. The thick red line indicates the ∆KL threshold level of 5% below which we reach convergence.

Figure 4 . 9 :

 49 Figure 4.9: Acceptance rate of all proposals for n kd e = 1. KDE and SCAM are on top panel, SLICE on bottom panel. Red plot and black dashed line shows where KDE finally converged and stopped updating.

Figure 4 . 10 :

 410 Figure 4.10: Acceptance rate of all proposals for n kd e = 5. KDE and SCAM are on top panel, SLICE on bottom panel. Red plot and black dashed line shows where KDE finally converged and stopped updating.

  max AL no kde mean AL kde mean AL no kde

Figure 4 . 11 : 6 APPENDIX Contents 6 . 1 6 . 3

 41166163 Figure 4.11: Maximum and mean autocorrelation lengths for n kd e = 1, comparing a run with KDE against a run without KDE. Black dashed line shows where KDE finally converged and stopped updating.

  .1) with s + (t ) = h ω sin[2Φ(t )]I cos 2ψcos[2Φ(t )]J sin 2ψ (6.2) s × (t ) = h ω sin[2Φ(t )]I cos 2ψ + cos[2Φ(t )]J sin 2ψ ,(6.3)where we have definedh = 2M 5/3 (π f ) 2/3d L (6.4)I = 1 + cos 2 ι J = 2 cos ι.(6.5) 123

5 h ω 2 (F + ) 2 +

 522 ) and we have< s|s > ψ,φ 0 ,ι = 4 (F × ) 2 (sin ωt | sin ωt ) + (cos ωt | cos ωt ) . (6.21)

22 )F

 22 × ( Ω) = ( m • p)( n • p) 1 + Ω • p , (6.23) ω = -(sin θ cos φ, sin θ sin φ, cos θ) m = (-sin φ, cos φ, 0) n =cos θ cos φ,cos θ sin φ, sin θ), (6.24)

. 25 ) 15 h ω 2 (

 25152 Finally, we have for the averaged signal to noise ratio to CGW signal :< s|s > ψ,φ 0 ,ι,θ,φ = 4 cos ωt | cos ωt ) + (sin ωt | sin ωt ) (6.26)

2 = 2 =

 22 d x( f (x, h)f (x)) d x f 2 (x, h) -2 d x f (x, h) f (x) + d x f 2 (x) = d x f 2 (x, h) -2 d x f (x, h) f (x) + const ,(6.30) where f (x) is the true PDF, ( f (x, h) is its KDE approximation and const is the term independent of the bandwidth h. We introduce local bandwidth h a attached to each sample point X a . Next we assume that all points in the vicinity of each point have similar bandwidth, in other words, h b ≈ h a for k near local points X b . Using these assumptions we can approximate the first term:N 2 d x f (x, h a ) f x, h b =

2 d x f (x, h a ) f x ≈ 2 N 2 N 1 , 1 ( 2 d / 2

 221122 N (N -1) ≈ N 2 . Combining these terms together gives us N 2 ( 2const ) of this expression by differentiating with respect to h c, j and equating it to zero: points X b in vicinity of X c and all components j . This assumption overestimates the bandwidth and therefore conservative: this is what is used in this paper. Expanding in this small parameters and retain only the terms quadratic in this small ratio we obtain the system of linear equations for 1/h 2 c, j :k near (2 d /2+1 -1)at each point X c for each direction in the parameter space ( j ) gives us the desired local bandwidth h c .

  we arrive at the system of linear equations for 1/h 2 c, j :

M3C2 1 (

 1 Multiple parallel Markov Chain Monte Carlo) is a python implementation of MCMC sampler. The aim of this tool is to improve the sampling robustness of complex posterior distribution, by running multiple chains in parallel. The cross check of the chain performance informs us about convergence (using Gelman-Rubin ratio Gelman and Rubin 1992). We have implemented two mechanisms of building the chain (1) using slice sampling (slice ) and (2) Metropolis-Hastings algorithm (MH) which could be used separately or together improving the mixture of the chains and reducing the auto-correlation length. Even though the sampler is very generic, we primarily use it within the context of GW data analysis. For the Metropolis-Hastings method we have implemented a set of proposal jumps: • SCAM (Single Component Adaptive Metropolis), jumps along one randomly chosen direction given by the eigen vectors of the covariance matrix Roberts and Rosenthal 2009; Haario et al. 2001 • DE (Differential Evolution), jumps along the direction given by difference of two randomly chosen samples of the chains, or (as in the classic implementation Ter Braak 2006; Jasper A. Vrugt 2008) by using state of different chains running in parallel. • ReMHA (Regional Metropolis Hastings Algorithm), the proposal represented by a mixture of several Gaussians distributions Roberts and Rosenthal 2009; Craiu et al. 2009. For SCAM we build the covariance matrix adaptively based on the accumulated samples. The use of the accumulated samples breaks the Markov property of the chain, making it asymptotically Markovian. The stability of the covariance matrix is yet another sign of the converged chain. One can stop adaptation after burn-in run. This proposal is suggested in Ellis 2013.ReMHA is similarly used adaptively. We use accumulated samples during the burn-in to estimate the number of clusters using Variational Bayesian Gaussian Mixture (skikit-learn package) and use this Gaussian mixture probablity as a proposal. This proposal is somewhat similar to the one suggested in[START_REF] Craiu | Learn From Thy Neighbor: Parallel-Chain and Regional Adaptive MCMC[END_REF]

  m)(l m) = d ΩY l m (Ω)W (Ω)Y (l m), giving the positions pa of pulsar with an associated weight w a .

  GWs emitted by merging stellar-mass black hole (or neutron star) binary systems because they can reach frequencies of ∼ 10Hz (see figure1.16), where the detector is most sensitive and the GW amplitude is loudest. SMBHBs are invisible to LIGO/Virgo because they merge before even getting to such high frequencies due their massive size.

	The first direct detection of a GW merger signal was announced by the collaborations
	on 11 February 2016. The development of a new generation of ground-based detectors
	(like Einstein telescope, Maggiore et al. 2020) is currently an active field of research.
	1.5.2 Space based detector : LISA
	The Laser Interferometer Space Antenna (LISA) is an international project of space-based
	GW detector in heliocentric orbit. Like LIGO and Virgo, it uses laser interferometry to

They are ground-based Michelson interferometers (LIGO in the USA and Virgo in Italy) with arms that are 3 to 4 km long. A GW passing through the detector modulates the length of the arms resulting in an interference pattern that can be measured

11

[START_REF] Abbott | Prospects for observing and localizing gravitational-wave transients with Advanced LIGO, Advanced Virgo and KAGRA[END_REF]

). This technology is suited for the detection of detect GWs, but with a different "triangle-shaped" configuration with arms that are ∼2.5 million kilometers long. This massive size brings technical difficulties that are worked around using Time-Delay Interferometry to measure the frequency shift δν/ν of the lasers where the signature of GWs will be imprinted (for details, see Amaro-Seoane et al. 2017).

It is most sensitive for GW signals with frequencies of ∼ 10 -3 Hz (see figure

1.16

).

PTA data 2.2.1 The timing model

In order to get the timing residuals δ t of a pulsar, we need to fit a timing model that gives a predicted value t pr ed of the TOAs t obs . This predicted value needs to account for all physical effects responsible for the delays induced in the TOAs. Here, we only give a brief overview of the mechanisms at work. For a complete description, please see

Edwards et al. 2006

. A least square fit is performed to estimate parameter values for :

. . 3.6 Effect of noise modeling on the CGW search . . . . . . . . . . . . . . . . . 3.7 Conclusion

	SEARCHING FOR CONTINUOUS GRAVITATIONAL WAVE SIGNAL
	WITH PTA
	The International Pulsar Timing Array 2nd data release is the combination of
	datasets from worldwide collaborations. In this study, we search for continu-
	ous waves: gravitational wave signals produced by individual supermassive
	black hole binaries in the local universe. We consider binaries on circular
	orbits and neglect the evolution of orbital frequency over the observational
	span. We find no evidence for such signals and set sky averaged 95% upper
	limits on their amplitude h 95 . The most sensitive frequency is 10nHz with
	h 95 = 9.1 × 10 -15

Contents 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 The noise and signal model . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 The noise model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 The signal model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Ranking the pulsars . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Ranking method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.2 Ranking results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.3 Sensitivity curve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 Model selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.2 Upper limit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 3 .

 3 1: List of the CGW parameters as defined in our model with their respected ranges. it), for more details see conclusion in Petiteau et al. 2013. So the CGW phase takes a very simple form:

Table 4 .

 4 1: Acceptance for various proposals (KDE run)

		Binned SCAM DE
	Acceptance rate	0.47	0.41	0.43

Table 4 .

 4 2: Acceptance for various proposals (binned proposal run) Tables 4.1, 4.2, compare acceptance rate of different proposals in two independent runs. One can see that KDE has the highest acceptance rate due to smart parameter grouping and interpolation between the samples incorporated in the KDE-based proposal. For the jumps we could use simultaneous jump in one or several (n kd e ) sub-groups.

Table 4 .

 4 .5. For high values of n kd e , even though we have decreasing acceptance rate, the autocorrelation drops too and mixing improves. However, this result should be take with a caution, using high n kd e could lead to a very low acceptance point as shown in Figure4.6. In Figure4.7 we show the autocorrelation length as a function of n kd e and it indicates that the optimal number is around n kd e = 5 with ac-Figure 4.6: Acceptance rate for several values of n kd e . 4: The maximum, minimum and mean autocorrelation lengths for three runs.

	n kd e KDE SCAM DE
	1	0.47	0.41	0.43
	5	0.22	0.36	0.41
	10	0.08	0.35	0.40
	Table 4.3: Acceptance for different n kd e

Table 4 .

 4

5: The maximum, minimum and mean autocorrelation lengths for different n kd e .

  2 , function of the GW param-=< s|s > ψ,φ 0 ,ι,θ,φ , (6.7)(s|s) = (F + ) 2 (s + |s + ) + (F × ) 2 (s × |s × ) + 2F + F × (s + |s × ),(6.8)In the following, we omit the h ω 2 prefactor. We have for each term in 6.8(s + |s + ) ∝ I 2 cos 2 2ψ(sin 2Φ(t )| sin 2Φ(t )) + J 2 sin 2 2ψ(cos 2Φ(t )| cos 2Φ(t )) × |s × ) ∝ I 2 sin 2 2ψ(sin 2Φ(t )| sin 2Φ(t )) + J 2 cos 2 2ψ(cos 2Φ(t )| cos 2Φ(t )) + ) 2 + (F × ) 2 I 2 (sin 2Φ(t )| sin 2Φ(t )) + J 2 (cos 2Φ(t )| cos 2Φ(t )) . = sin ωt cos φ 0 + cos ωt sin φ 0 (6.14) cos 2Φ(t ) = cos ωt cos φ 0sin ωt sin φ 0 , (6.15) yielding (sin 2Φ(t )| sin 2Φ(t )) = cos 2 φ 0 (sin ωt | sin ωt ) + sin 2 φ 0 (cos ωt | cos ωt ) + 2 cos φ 0 sin φ 0 (sin ωt | cos ωt ), (6.16) (cos 2Φ(t )| cos 2Φ(t )) = cos 2 φ 0 (cos ωt | cos ωt ) + sin 2 φ 0 (sin ωt | sin ωt ) -2 cos φ 0 sin φ 0 (sin ωt | cos ωt ). After averaging with respect to the initial phase we get < sin 2Φ(t )|2 sin Φ(t ) > φ 0 =< cos Φ(t )| cos Φ(t ) > φ 0 =

	eters, we have (s -2IJ cos 2ψ sin 2ψ(sin 2Φ(t )| cos 2Φ(t )), d ψ π d φ 0 2π d cos ι 2 d cos θ 2 d φ 2π ρ 2 + 2IJ cos 2ψ sin 2ψ(sin 2Φ(t )|2 cos Φ(t )), (s + |s × ) ∝ I 2 cos 2ψ sin 2ψ(sin 2Φ(t )| sin 2Φ(t )) -J 2 cos 2ψ sin 2ψ(cos 2Φ(t )| cos 2Φ(t )) (6.9) (6.10) -2IJ (sin 2 2ψ -cos 2 2ψ)(sin 2Φ(t )| cos 2Φ(t )). (6.11) We average the expression with respect to the polarization angle ψ to get < s|s > ψ = 1 2 h ω (6.12) In the case of a monochromatic, non evolving CGW signal, we can approximate 2Φ(t ) by : 2Φ(t ) → ωt + φ 0 , (6.13) so the sin 2Φ(t ) and cos 2Φ(t ) terms become 1 2 (sin ωt | sin ωt ) + 1 2 (cos ωt | cos ωt ), (6.18) So ρ 2 becomes sin 2Φ(t ) (6.17) 4 ω < s|s > ψ,φ 0 = 1 h 2

2 (F

A straight line is the shortest path between two points. The geodesic equation extends this deifinition to curved space.

Contraction with second index gives the same result with opposite sign and contraction with last index is null. This is due to the antisymmetric properties of the Riemann tensor.

At that stage, the trace reversed metric is introduced to render the "wave equation"-like properties of the linearized Einstein field equations.

Or De Donder gauge, or harmonic gauge... it has got many names.

This is true considering the point-mass approximation. In reality, for two objects of radius L, we can trivially show using r (t c ) = L that merger happens at t c = τ c 1 -L/r 0 4 which is less than τ c . We deduce from this that bigger objects merge sooner.

Note that this is not a true singularity but rather a coordinate singularity. It is possible to make it disappear by chosing an apprioriate coordinate system, like for instance, Lemaître coordinates.

This is a big approximation, for the sake of quickly demonstrating the expected frequency dependence of the characteristic strain. For a realistic derivation we would have to take into account redshift, mass distributions, mass ratios and even eccentricity of the binaries (see Taylor 2021).

Note true for legacy data because of gaps.

To be more precise, the GW affects the proper distance between the beam splitter and the mirrors.

In section

2.2.2, we show that in PTA noise models the σ i are function of two parameters EFAC and EQUAD characterizing the white noise.

For example, in the EPTA dataset, the jitter-like effect is considered to be included in the EQUAD term.

Not to be confused with the spin velocity of the pulsar sending photon packets (pulses) with a rate ν.The frequency f ν represents the "colour" of the photons measured on that frequency band.

In reality, many other noises can be added here, but only these two will be used in this manuscript

Note that we omit the index α for matrices N, M, F and Σ to avoid overcrowded equations but it should actually be there because they describe properties of pulsar α.

The proposal distribution can be anything, from a Gaussian distribution centered on θ to a uniform distribution, as long as the output sample is random.

Other methods like parallel-tempering or nested sampling can estimate these difficult integrals with fairly good accuracy.

For setting an upper limit we use uniform prior on the amplitude of GW strain

I am co-author of the paper for which I worked on the ranking method for CGW signals.

Furthermore, detectable CGW signals must be louder than the GWB. Since the GWB is stronger at lower frequencies, CGW signals are more likely to be found at high frequencies.

Here, the SNR is for a single frequency signal.

Note that the highest spike corresponds to frequencies close to f yr and should be discarded.

The probability Fn kd e (x, h) corresponds to the proposal probability for the jump q( θ| θ * ) presented in section

2.3.3

Because we do not know the true distribution, we have to use histograms to evaluate KL. The formula for binned KL is given in Appendix 6.2.

The DE introduced in Ter Braak 2006 is using population MCMC, here we rather use it on a single chain in the spirit described in the snooker proposal in Jasper A.[START_REF] Jasper | Differential Evolution Markov Chain with snooker updater and fewer chains[END_REF]

https://gitlab.in2p3.fr/lisa-apc/mc3

ACKNOWLEDGEMENTS

Lehenik eskertu nahi ditut ene famili osoa eta lagunak, bereziki ene aitamak, milesker sustengaturik eta beti hor izanik. I thank Stas and Antoine for their sympathy and always good advice throughout my PhD, as well as all the people in the PTA collaboration whom I have been working with.

I would also like to thank my examiners and reporters for accepting to evaluate my thesis and Alberto for wearing the flamingo shirt matching his hair at my defence.

Finally, I thank all my lab mates at APC and the shelf in our office.

The pulsars above the red dashed line contribute less than 5% of the total SNR 2 . This means only 12 pulsars out of 65 contribute on average to 95% of the total SNR 2 of the array at 5nHz. Note that, while only the best 22 pulsars in the array are shown in the figure, the normalized total SNR has been evaluated using all 65 pulsars in the array.

contributions to the total SNR 2 are not considered significant. This value was chosen to be 0.95. The process is illustrated in Figure 3.1 for pulsars from the IPTA second data release [START_REF] Perera | The International Pulsar Timing Array: second data release[END_REF]).

Due to the strong dependence of ρ2 a ( f ) on f , the resultant CGW pulsar ranking is also frequency dependent. This can be clearly seen from Figure 3.2. In our analysis, we use 100 log-spaced frequency bins between 10 -9 and 10 -7 Hz. Ranking lists were obtained separately for each frequency bin. In order to construct the final ranking catalog of best pulsars at a given frequency range, the lists at each frequency are merged together. This procedure ensures that we will gain at least, no matter the CGW frequency, 95% of the total SNR 2 of the array. Markov Chain Monte Carlo approach is frequently used within Bayesian framework to sample the target posterior distribution. Its efficiency strongly depends on the proposal used to build the chain. The best jump proposal is the one that closely resembles the unknown target distribution, therefore we suggest an adaptive proposal based on Kernel Density Estimation (KDE). We group parameters of the model according to their correlation and build KDE based on the already accepted points for each group. We adapt the KDEbased proposal until it stabilizes. We argue that such a proposal could be helpful in applications where the data volume is increasing and in the hypermodel sampling. We tested it on several astrophysical datasets (IPTA and LISA) and have shown that in some cases KDE-based proposal also helps to reduce the autocorrelation length of the chains. The efficiency of this proposal is reduces in case of the strong correlations between a large group of parameters.

is most efficient if the full parameter space could be split into many small uncorrelated groups. The performance drops significantly if dimensionality of a group is larger than 5.

• We have build KDE with the self-optimizing bandwidth based on the distribution of a sample points provided at the input.

Note that the adaptation breaks "Markovian" properties of the chain. Either one should dismiss the parts of the chain during the adaptation or, in case of uninterrupted adaptation, assume that the chain is only asymptotically Markov. We give detailed description of implementation in the next two sections (4.2, 4.3).

We have implemented the KDE-based proposal in a particular sampler https://gitlab. in2p3.fr/lisa-apc/mc3. We give a detailed description of this sampler in Appendix 6.4.

The main feature of this sampler is that it runs several chains either completely independently or as parallel tempering. Multi-chain run is used to compute Gelman-Rubin ratio (Gelman and Rubin 1992) to monitor the convergence.

We assess the performance of the suggested proposal in two applications to GW data analysis. In first one we analyse the data combined by International Pulsar Timing Array (IPTA) collaboration searching for a continous GW signal in the nano-Hz band. As the second dataset we use simulated LISA data publically available through LISA Data Challenge (2a) portal. We use the KDE-based proposal to infer parameters of 6 Galactic white dwarf binaries. We present the performance of our proposal for those two data analysis problems in Section 4.4, in particular we show that the KDE-based proposal allows to reduce the autocorrelation length while keeping high acceptance rate.

We conclude the paper with discussion on the limitation and possible extension of our method in Section 4.5.

Kernel Density Estimation

In this rather short section we describe our particular way of building KDE. We start with a short introduction to KDE and then give details of the bandwidth optimization that we use. KDE is a non-parametric method used to estimate a probability density function (pdf ) based on a finite set of sample points [START_REF] Parzen | On Estimation of a Probability Density Function and Mode[END_REF][START_REF] Rosenblatt | Remarks on Some Nonparametric Estimates of a Density Function[END_REF]). It is a smooth alternative to a histogram. The advantage of KDE is that it uses no binning and gives a continuous function interpolating (and extrapolating) across the whole parameter space.

For a d-dimensional dataset { X } of size N and kernel K (x, h), we have our KDE f (x, h) : to KDE. We search for a continuous gravitational wave signal while fitting for pulsar noise parameters (Aggarwal 2019). The KDE for the noise parameters has been built using posterior samples obtained from the preceding single pulsar analysis. We will compare three different runs:

• using a KDE-based + default jump proposals SCAM and DE (labeled as "KDE");

• using binned empirical distributions + default proposals SCAM and DE (labeled as "Binned");

• using only default proposals SCAM and DE (labeled as "None").

The binned empirical distributions are essentially 2-dimensional histograms based on the same posterior samples as used in building KDE [START_REF] Taylor | [END_REF]). Those carry a similar spirit to KDE, being pair-wise approximation to marginalized posterior, but at the same time are fundamentally different from the KDE in the sense that the KDE is a continuous function in space that interpolates between the sample points using smoothening kernel.

In addition KDE-based proposal makes grouping based on the parameter correlation that could lead to more than two dimensional group (see figure 4.4).

CHAPTER 5 SUMMARY AND PERSPECTIVES

In this PhD thesis, we have explored the low frequency realm of gravitational waves using the stability of millisecond pulsars with Pulsar Timing Arrays. More specifically, we have used the second data release of the International Pulsar Timing Array that is the combination of all the pulsar timing data collected for the past decades. It gathers high precision measurements from the European, North-American and Australian collaborations. We developed robust data analysis pipelines to search for the presence of continuous gravitational wave signals, i.e. monochromatic and deterministic gravitational wave signals induced by supermassive black hole binary systems. The manuscript was separated into four main sections : (I) the presentation of the theoretical background necessary to understand the physics of the pulsar timing array, (II) the data analysis methods that are used to search for gravitational wave signal in the data, (III) the results of this search obtained with the IPTA DR2 and finally, (IV) the presentation of a new method using kernel density estimation to optimize our data analysis pipelines.

The search for deterministic gravitational wave signal in the IPTA DR2 was conclusive.

We found no evidence for such signals so we have set an upper limit on amplitude h.

It appears that the IPTA DR2 is most sensitive for signals with a frequency around 10nHz

where its upper limit is h 95 = 9.1×10 -15 and provides the best sensitivity to date of all PTAs at high frequencies thanks to the improved cadence of observation. Following the recent discoveries of low frequency processes common to all pulsars in the array, we need to include it into the model. It came out that the latter correlates with the common red noise, slightly deteriorating the upper limit at low frequencies. Even though the search was not fruitful, we have stumbled upon gravitational wave like signals at high frequencies with relatively high Bayes factor that triggered our attention. These features were absorbed after changing the noise model. This result is crucial to understand how the noise model directly affects the observation and must be chosen with great care. This study led to a publication for which I am corresponding author. It was shared with the international collaboration and approved by the steering committee. It was submitted to MNRAS under the name : Searching for continuous Gravitational Waves in the second data release of the International Pulsar Timing Array.

We proposed a new method using kernel density estimation in gravitational wave data analysis. It came from a desire to simplify and yet optimize the data analysis pipelines that we use. In this chapter, we present a new fast semi-analytical method to optimize the bandwidths of the kernels. This optimal kernel density estimation is then used with MCMC as a proposal distribution. We quantify the quality of our proposal by estimating the autocorrelation length (i.e. the effective sampling time) of the produced Markov chains for different proposals. It turned out that our proposal was giving the best performance with IPTA data. We wanted to push it a step further by suggesting an adaptive proposal that would build itself as the MCMC is sampling. We tested it on simulated LISA data to search for galactic binary signals. The proposal was still reducing the effective sampling time but its performance was less probing. This work led to the writing of an article that was submitted of Physical Review D and pre-published on Arxiv (Falxa et al.

2022).

The search for gravitational wave signal is still an ongoing project. The recent detection of the common red noise process by all collaborations really has spawned renewed interest for the pulsar timing world. Hopefully, in the upcoming days, we will finally detect the signals that have been hunted for many years. The latest data releases together with the newly available radio-observatories provide the best sensitivity to date, raising hopes of gravitational wave detection in the next few years.