Probing Modified Gravity with Pulsar Timing Array experiments and Gravitational Waves

A Stochastic Background of Gravitational Waves in the nHz frequency band has been recently detected by the Pulsar Timing Array (PTA) collaboration, with an evidence between 2 and 4 σ. The origin of such a signal is still unclear: both Super-Massive Black Hole Binaries and some cosmological sources can fit the data. Among the cosmological sources, Scalar-Induced Gravitational Waves (SIGWs) show the largest Bayes factor, so they look particularly promising, also because they are strictly connected to Primordial Black Holes, which represent possible candidates for Dark Matter. This ground-breaking result opens the door to the possibility to use pulsars to probe Modified Gravity (MG), in an indirect way. Hence, this thesis aims to study such SIGWs in a general theory of gravity. In order to do so, we have first derived the evolution equations for cosmological fluctuations by perturbing the spatially flat Friedmann-Lemaître-Robertson-Walker metric and the stress-energy tensor, working in the Poisson gauge. While at linear order scalar, vector and tensor perturbations evolve independently, at second order they mix, and scalar fluctuations become a source for (second-order) tensor perturbations. We have computed the second-order GW equation of motion in General Relativity (GR), keeping the anisotropic stress sourced by an imperfect fluid and including the contribution coming from scalar-tensor mode mixing. In doing so, we noticed that there is a coupling between the linear anisotropic stress and the scalar potential. Consequently, any mismatch between the two scalar modes will manifest in the source term of the SIGWs equations and can potentially influence the spectrum. Besides the anisotropic stress sourced by an imperfect fluid, another source of anisotropic stress can arise in MG theories due to a purely geometrical effect. Therefore, after reviewing SIGWs in GR, we have investigated second-order tensor modes sourced by terms quadratic in first-order scalar perturbations in f (R) gravity. Specifically, we have computed the perturbed evolution equations in a general f (R) theory and found how the source term gets modified in the induced GWs equation. Then, we applied our general set of equations to the specific model f (R) = R + αR^2. We used an analytical approach to study the first-order corrections in the field equations introduced by a given f (R) model with respect to GR, without relying on the so-called quasi-static approximation. We considered the radiation-dominated epoch and used a log-normal primordial power spectrum, assuming Gaussian curvature fluctuations. Within this framework, we provided a formula for calculating the SIGWs spectral density which accounts for the usual GR part and the first-order correction introduced by the MG model in terms of the α parameter. In our analytical approach, we noticed that to compute the post-GR corrections we shall account for the sub-dominant matter contribution that is present during radiation dominance. Finally, focusing on PTA scales, we compared the spectral density of SIGWs in standard GR and f (R) gravity.