{'fraction_non_alphanumeric': 0.03506895788774984, 'fraction_numerical': 0.011468397374437642, 'mean_word_length': 4.383958705578718, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A possible signature of cosmic neutrino decoupling in the nHz region of the spectrum of primordial gravitational waves

In this paper we study the effect of cosmic neutrino decoupling on the spectrum of cosmological gravitational waves (GWs). At temperatures T>>1 MeV, neutrinos constitute a perfect fluid and do not hinder GW propagation, while for T<<1 MeV they free-stream and have an effective viscosity that damps cosmological GWs by a constant amount. In the intermediate regime, corresponding to neutrino decoupling, the damping is frequency-dependent. GWs entering the horizon during neutrino decoupling have a frequency f ~ 1 nHz, corresponding to a frequency region that will be probed by Pulsar Timing Arrays (PTAs). In particular, we show how neutrino decoupling induces a spectral feature in the spectrum of cosmological GWs just below 1 nHz. We briefly discuss the conditions for a detection of this feature and conclude that it is unlikely to be observed by PTAs.


Introduction
The presence in the Universe today of a stochastic background of gravitational waves (GWs) is a quite general prediction of several early cosmology scenarios. In fact, the production of gravitational waves is the outcome of many processes that could have occurred in the early phases of the cosmological evolution, like the amplification of vacuum fluctuations in inflationary [1] and pre-big-bang cosmology scenarios [2], phase transitions [3], and finally the oscillation of cosmic strings loops [4]. In most of these cases, the predicted spectrum of gravitational waves extends over a very large range of frequencies.
The detection of such primordial gravitational waves, produced in the early Universe, would be a major breakthrough in cosmology and high energy physics. This is because gravitational waves decouple from the cosmological plasma at very early times, when the temperature of the Universe is of the order of the Planck energy. In this way, relic gravitational waves provide us a "snapshot" of the Universe near the Planck time, in a similar way as the cosmic microwave background (CMB) radiation images the Universe at the time of recombination.
The extremely low frequency region (f 10 −15 Hz) in the spectrum of primordial gravitational waves can be probed through the anisotropies of the CMB. In particular, gravitational waves leave a distinct imprint in the so-called magnetic or B-modes of its polarization field [5,6]. The frequency region between ∼ 10 Hz and few kHz is probed by operating large scale interferometric GW detectors like LIGO [7] and VIRGO [8], that, although designed with the aim to detect astrophysical signals, can possibly also detect signals of cosmological origin [9]. The LISA space interferometer [10], that will hopefully operate in the 2020s, will probe the frequency region between 10 −4 and 1 Hz. Finally, pulsar observations can be used to obtain information on the stochastic GW background, through the technique known as pulsar timing. The so-called "Pulsar Timing Arrays" will probe the region from ∼1 to ∼ 100 nHz, with a maximum sensitivity between ∼3 and ∼10 nHz [11].
In order to compare the theoretical predictions with the expected instrument sensitivities, one needs to evolve the GWs from the time of their production to the present. This is done by assuming that gravitational waves propagate freely across the Universe. In this case, the only effect is that the amplitude of the wave decreases due to the expansion of the Universe. However, GWs are sourced by the anisotropic stress part of the energy-momentum tensor of matter, so that the above assumption is well-motivated only when anisotropic stress can be neglected. It is already known that the anisotropic stress of free streaming neutrinos acts as an effective viscosity, absorbing gravitational waves in the low frequency region, thus resulting in a damping of the B-modes of CMB [12,13,14,15,16]. ‡ In the present work we aim to study a possible signature of neutrino decoupling on the spectrum of cosmological GWs. During decoupling, the effective viscosity of neutrinos is increasing from zero (its value at high temperatures, when the neutrinos are tightly coupled to the cosmological plasma) to a finite value. Neutrino decoupling happens when the temperature of the Universe is T ≃ 1 MeV, or redshift z ≃ 10 10 ; GWs entering the horizon at that time have a frequency f ≃ 1 nHz. Since the viscosity results in a damping of the wave, we expect to have a frequency-dependent absorption of GWs in this frequency range. A spectral feature in the nHz range would be potentially interesting for PTAs. PTAs combine the fact that pulsars are very stable clocks with the fact that the time of arrival of the pulse also depends on the GW background between the pulsar and the Earth [11]. The use of a pulsar array allows to correlate the electromagnetic signals from different pulsars, thus eliminating variations in the arrival times that do not depend on the GW background. The dominant contribution to the GW background in the PTAs sensitivity range is that of supermassive black hole binaries following galaxy mergers [21,22,23]. This has an astrophyisical origin and would obviously not be affected by the decoupling of cosmological neutrinos. The dominant cosmological constribution to the GW signal in the nHz range is expected to be that of cosmic strings, that are topological defects left after symmetry-breaking phase transitions occurring in the early Universe. ‡ In Ref. [17], the evolution of cosmological GWs in the presence of neutrino free-streaming was studied up to second order in perturbation theory. Moreover, GW propagation is also affected by entropy-generating events (like the e + e − annihilation) occurring during the cosmological evolution; the effect on the GW spectrum has been studied in Ref. [18]. The authors of Ref. [18] also consider the effect of neutrino free streaming but assume an instantaneous neutrino decoupling. A similar analysis, including an accurate computation of the primordial spectrum as generated during inflation, and extended to higher frequencies, can be found in Ref. [19]. Finally, in Ref. [20], the role of a non-vanishing lepton asymmetry in determining the amount of absorption has been considered. (Note added in the arXiv version.) Other, possibly important, contributions in this frequency range are those of cosmic superstrings (i.e., topological defects arising in string-theory inspired inflationary models) and of inflation-generated GWs (through the usual amplification of quantum zero-point fluctuations). For our purposes, one important difference between the cosmic (super)strings-generated GWs on one side, and the inflationary GWs on the other, is that the former are created by causal processes at horizon and subhorizon scales, while the latter are generated by the stretching of quantum fluctuations to super-horizon scales. In the following we will concentrate on GWs that were outside the horizon at some time before neutrino decoupling, so that our results rigorously applies to inflationary GWs but not necessarily to string-generated GWs. It should be noted that the inflationary signal is, in the nHz range, the weakest of the cosmological signals mentioned above, and is probably beyond the reach of the ongoing PTA projects. It could be however within the reach of future instruments [11,24] like the Square Kilometer Array (SKA) [25].
The paper is organized as follows. In Sec. 2, we recall the basic equations describing the evolution of GWs in the presence of an imperfect fluid. In Sec. 3, we discuss the interaction between GWs and the neutrino fluid in the three regimes T ≫ 1 MeV, T ≃ 1 MeV and T ≪ 1 MeV, corresponding to f ≫ 1 nHz, f ≃ 1 nHz and f ≪ 1 nHz respectively. In Sec. 4, we present the results of the numerical integration of the Einstein-Boltzmann system, show the modification in spectrum of cosmological GWs and briefly discuss the conditions that could lead to its detection.

Basic equations
We shall use, all throughout the paper, natural units in which c = = k B = 1.
Let us consider a gravitational wave, propagating on the background of a flat (k=0) Friedmann Universe. In synchronous gauge, the spatial components of the perturbed metric are written as [we use the (− + ++) signature for the metric] while the other components are left unperturbed: g 00 = −1 and g 0i = 0. We will consider only the transverse traceless part of h ij . With this restriction, h ij only has two degrees of freedom, corresponding to the two polarizations of a GW. Here a(t) is the cosmological scale factor, that evolves according to the background Friedmann equation: whereρ is the background density of the cosmological fluid (in general, we will use overbars to denote background quantities). The Einstein equation for the evolution of a transverse traceless metric perturbation h ij takes the form of a wave equation [14,26,27]: where Π ij is the anisotropic stress, i.e. the traceless part of the three dimensional energy-momentum tensor T i j of the cosmological fluid, representing dissipative effects that are not present in a perfect fluid. It is then defined through the relation where P is the total pressure of the fluid (including possibly a small perturbation with respect to the background). We assume that the anisotropic stress is relevant only for one component of the cosmological fluid (the neutrinos), providing a densityρ ν = f νρ , while for the other components (e.g., the photons) it can be safely neglected.
In the following, we will find more convenient to use the conformal time η, related to the synchronous time t by dt = a(η)dη as our time variable. Moreover, since the perturbation equations are linear, it is very convenient to Fourier transform the spatial dependence of all the relevant quantities. Then the wave Eq. (3) rewrites as: where dots denote derivatives with respect to η, and H is the conformal Hubble constant H =ȧ/a. For the sake of simplicity, we use the same symbol for a given quantity and for its Fourier transform. In a perfect fluid, anisotropic stress is absent by definition (in fact, for a perfect fluid T i j = Pδ i j ), so that Eq. (3) takes the homogeneous form: In the radiation-dominated era, a ∝ η and H = 1/η. Then, considering the quantity χ ij ≡ kηh ij we find thaẗ This is a simple, undamped oscillator equation, that admits the solutions χ In the limit of kη → 0, i.e., for a wave that is far outside the horizon, we have that Then, the most general solution outside the horizon (during the radiation dominated era) is the superposition of a constant mode and of a mode that decays like 1/η. This is actually true also in the presence of anisotropic stress, because dissipative effects cannot operate on scales larger than the horizon, so that Π i j has to vanish in the limit kη → 0.
The macroscopic properties of the cosmological fluid, including the anisotropic stress, can be derived by the phase space distribution of its particles. The phase space is described by three positions x i and by their three conjugate momenta P i ≡ mdx i /ds, and evolves according to the Boltzmann equation: whereL ≡ d/dη is called Liouville operator, andĈ is the collision operator accounting for changes in the distribution function due to collisions between particles. Using the geodesic equation, the total derivative of f can be rewritten as: The link between the Einstein Eq. (3) describing the evolution of the metric perturbation h ij and the Boltzmann Eq. (9) is given by the following expression for the energy-momentum tensor in terms of the distribution function: where g is the determinant of the metric. Then, once the collision term is also specified, Eqs. (2), (3), (9) and (11) are all that is needed, at least in principle, to follow the propagation of a GW.
In general, the Einstein-Boltzmann system for the coupled evolution of h ij and f is an integro-differential system, because the source term in Eq. (3) is given by an integral over f . The system can be however reduced to an ordinary differential system in the case of massless particles. This is made by integrating out the dependence of f on energy and expanding its dependence on the direction of the particle momentum in Legendre polynomials. This method is widely applied in the numerical treatment of scalar perturbations [29]. The analogous equations for the tensor case have been derived in Ref. [28], to which we refer the reader for further details; here we just quote the relevant equations. The coupled Einstein-Boltzmann system is equivalent to the following system of infinite differential equations in k-space: In the above equations, the F (ℓ) 's and the C (ℓ) 's are related to the coefficients of the expansion in Legendre polynomials of f andĈ[f ], respectively (see Ref. [28] for details). The anisotropic stress is given by:

Interaction of gravitational waves with neutrinos
In this section we discuss how cosmological neutrinos affect the propagation of cosmological GWs. We restrict our attention to waves entering the horizon well before the time of matter-radiation equality, corresponding to a redshift z ≃ 10 4 . This corresponds to waves with a present frequency ν ≫ 10 −16 Hz. We recall that in the radiation-dominated era, a ∝ η and H = 1/η.
The thermal evolution of neutrinos can be divided into two distinct regimes. In the early Universe, when the temperature was sufficiently high, neutrinos were kept in thermal equilibrium with the cosmological plasma by frequent reactions like νν ↔ e + e − , νe ↔ νe, etc. As the Universe expanded and cooled down. the interaction rate Γ ν for these reactions eventually became smaller than the expansion rate H and the neutrinos decoupled from the cosmological plasma. From this point on, neutrinos can be considered as freely streaming through the Universe. The temperature T dec that separates the collisional from the non-collisional regime, defined as the temperature when the interaction rate is equal to the expansion rate, i.e. Γ ν (T dec ) ≃ H(T dec ), is found to be T dec ≃ 1 MeV, corresponding to z ≃ 10 10 and t ∼ a few seconds. A wave that enters the horizon at that time has a frequency f ∼ 1 nHz. Higher frequency waves enter the horizon before neutrino decoupling, and viceversa. Then, the following three regimes for the evolution of a GW can be correspondingly identified: • When T ≫ 1 MeV, the neutrinos are tightly coupled to the cosmological plasma, that behaves like a single, perfect fluid. The anisotropic stress of neutrinos (as well as of the other components) is negligible and the GW evolves according to the homogeneous Eq. (5). This regime is relevant for waves with f ≫ 1 nHz. • When T ≃ 1 MeV, the neutrinos are decoupling from the plasma. They cannot be considered as a perfect fluid because the mean free path of particles is getting large, so that the anisotropic stress cannot be neglected. However, the collisions have still to be taken into account. The GW evolution has to be calculated from the full system Eqs. (12)-(16) with the C (ℓ) 's modeling the collision processes. This regime is relevant for waves with f ≃ 1 nHz.
• When T ≪ 1 MeV the neutrinos are collisionless and behave as free particles. The anisotropic stress cannot be neglected: the GW evolution has to be calculated from the full system Eqs. (12)-(16) with C (ℓ) = 0 for all ℓ's. This regime is relevant for waves with f ≪ 1 nHz.
In the first, high-temperature regime, the metric perturbation simply evolves as h ij =h ij ≡ λ ij sin(kη)/kη (assuming the initial velocity vanishes), where λ ij is the initial value of the perturbation. In other words, the GW oscillates with a timedependent amplitude ∝ 1/kη. This can be seen as the energy loss of the wave due to the expansion of the Universe.
The third, low-temperature regime was first studied in Ref. [14], where it was realized that the anisotropic stress of neutrinos would lead to a partial absorption of the GW. In particular, for the standard value of the neutrino density, it was found that the amplitude of the wave when it is well inside the horizon is reduced to 80% of its value in the absence of stress. The intensity is reduced to (0.8) 2 = 64% of its value in the absence of stress. In Ref. [14], the Boltzmann-Einstein system was treated in its integro-differential form. We re-analyzed the problem in Ref. [28] using the multipole formalism, and allowing for different values of the neutrino density. In general, one has that the behaviour of the wave well inside the horizon (kη ≫ 1) is of the kind: where D is a damping factor, φ is a phase and f ν ≡ρ ν /ρ is the neutrino fraction, i.e. the ratio of the neutrino densityρ ν to the total densityρ. This is constant in the radiation-dominated era since bothρ ν andρ scale like a −4 , and is equal to f ν = 0.4052 for the standard case of three neutrino families with a present temperature T 0 ν = (4/11) 1/3 T 0 γ = 1.9 K. In the following we will neglect the phase shift φ, since we are only interested in the behaviour of the wave when averaged over many periods. The damping factor D only depends on f ν . The result of Ref. [14] can thus be written as D(f ν = 0.4052) = 0.8. In Ref. [28] we have found the slightly larger value D(f ν = 0.4052) = 0.88. We have also calculated the maximum amount of damping, corresponding to D(f ν = 1) = 0.75. It is worth stressing that the cosmological neutrino background has not been directly observed yet. Although strong deviations from the standard scenario with f ν = 0.4052 are unlikely, the possibility that f ν has a different value (due for example to the presence of additional free-streaming particles in the early Universe) should be taken into account.
It is worth stressing, at this point, that the behaviour of the wave in either of the two extreme regimes f ≫ 1 nHz and f ≪ 1 nHz does not depend on the frequency of the wave nor on the values of the cosmological parameters (at least for waves entering the horizon before matter-radiation equality). This is due to the fact that, when C (ℓ) = 0 and a ∝ η the explicit k dependence in the Einstein-Bolztmann system can be eliminated by considering the time variable u = kη. Also, the convenience of using u is that u = kη ∼ 1 corresponds to the time of horizon crossing. Moreover, by writing h ij = λ ijh (η), it can be seen that h(η) satisfies the same differential equations as h ij [upon a redefinition of the F (ℓ) 's and C (ℓ) 's], so that we can always take without loss of generality λ ij = 1. Taking all this into account, we have that the GW evolution for kη ≫ 1, in terms of χ ij = kηh ij = uh ij is: The fact that the (rescaled) amplitude of the wave is going from 1 for f > 1 nHz to D(f ν ) < 1 for f < 1 nHz already points out to the fact that one should expect a frequency dependence in the region f ≃ 1 nHz. In order to quantify this, one should solve the evolution equations for waves entering the horizon in the intermediate regime where neutrinos are neither a perfect fluid nor can be considered as freely streaming. The effect of collisions has to be taken into account. The exact computation of the collision term C[f ] (and consequently of the C (ℓ) 's) depends on the details of the interaction. However, a useful although rough approximation consists in writinĝ is the deviation of the neutrino distribution function f from its thermal equilibrium value f (0) , and τ is a characteristic time of the interactions that mantain the equilibrium. This form, albeit very simple, captures the main features that we expect from the actual collision term. When the characteristic time of the interactions is very large, τ → ∞ andĈ[f ] → 0, i.e. the fluid is collisionless. On the contrary, when the characteristic time is very small τ → 0 andĈ[f ] is very large unless δf → 0 too, i.e. the frequent collisions tend to mantain the equilibrium. This form forĈ[f ] gives [28]: In our calculations, we take τ to be the mean time between collisions, given by τ = (n σv ) −1 , where n is the number density of particles and σv is the thermally averaged cross section times velocity. The number density of particles in thermal equilibrium is n ≃ T 3 , while an adequate value for the cross section for weak interaction processes is σ ≃ G 2 F T 2 , where G F is the Fermi constant. Then:

Results and discussion
Equipped with the form (19) for the collision terms, with τ given by Eq. (20), we have integrated the Boltzmann-Einstein system for 71 values of f = k/2π between 10 −12 Hz and 10 −8.5 Hz ≃ 3 nHz, with a constant logarithm spacing ∆ log f = 0.05. We have assumed the standard value f ν = 0.4052. For each value of f , we have defined a frequency-dependent damping factor D f as: where the brackets . . . denote the average over many periods, taken when the wave is well inside the horizon (kη ≫ 1). It is clear that when h ij = sin(kη)/kη we have D f = 1. The damping factor basically quantifies how much the amplitude of a GW produced in the early Universe would be reduced due to the neutrino anisotropic stress, with respect to propagation in a perfect fluid. We express our results in terms of the density parameter of GWs Ω GW (f ) = (dρ GW /d ln f )/ρ c , where ρ GW is the energy density of GWs and ρ c is the critical density of the Universe. The quantity Ω GW is used by theorists to quantify the intensity of the cosmological GW background predicted in a given scenario, and experimental sensitivities are also often quoted in terms of it. The present intensity is usually computed simply by rescaling the GW amplitude at the source by a factor 1/kη to take into account the redshift due to the expansion. This should be corrected to take into account also the anisotropic stress of neutrinos. Since ρ GW ∝ h 2 ij , the correction factor is simply D 2 f . Thus, if we callΩ GW the value in absence of stress, we have that Ω GW = D 2 fΩ GW . In Fig. 1 we show D 2 f as a function of frequency. Frequency Hertz D f 2 Figure 1. The correction factor D 2 f = Ω GW /Ω GW . The highest frequencies in the plot correspond to waves entering the horizon when the neutrinos are tightly coupled to the plasma, so that there is no dissipation. The lowest frequencies correspond to waves entering the horizon when the neutrinos are freely streaming, so that D 2 f is constant and roughly equal to 0.75.
The fact that the collisions introduce a frequency dependence on the amplitude of the GW background (in addition to any dependence that could be already present at the source) is a more interesting case with respect to a frequency independent suppression, since the latter could not be disentagled from our ignorance on the amplitude of the original spectrum (that usually depends on the free parameters of the theory and thus cannot be determined a priori without any additional experimental input). On the other hand, a frequency-dependent effect can be disentangled by using this same dependence. To better quantify this, let us suppose that the GW spectrum at the source, has a featureless power law behaviour ∝ f α . After correcting for the expansion, the shape is still the same. Then the spectrum, once the effect of anisotropic stress has been taken into account, will be: The logarithmic slope of the spectrum in a given point is given by α ′ = d log Ω GW /d log f , so that, defining the deviation ∆α = α ′ − α, we have: In Fig. 2, we show ∆α as a function of frequency. We find that it has a maximum value of ∆α = 0.15 at f ≃ 0.1 nHz. The nHz region of the spectrum of GWs can be probed by PTAs, combining the fact that pulsars are very stable clocks with the fact that the time of arrival of the pulse also depends on the GW background between the pulsar and the Earth [11]. Correlating the electromagnetic signals from different pulsars could, in principle, enable a positive detection of the GW background in the galaxy. PTAs are sensitive to GWs with frequency between roughly 1 and 100 nHz, with a maximum sensitivity in the region between 3 and 10 nHz. The lower bound comes from the fact that standard pulsar timing techniques absorb any low-frequency signal, and so the time span of the data (currently ∼ 30 years) gives a lower bound on the observable frequencies.
The principal source of GWs in the region of maximum sensitivity are believed to be coalescing supermassive binary black-hole systems in the centre of merging galaxies. However, a signal of cosmological origin could also be present in that same frequency range, for example GWs produced during the inflationary era or from the oscillation of cosmic string loops.
Can the effect of anisotropic stress leave an imprint that could be, in the future, detected using PTAs? This could be possible, at some conditions. First of all, obviously, a signal of cosmological origin has to be present in the nHz range and it has to be strong enough to be detectable by PTAs. As noted in the introduction, the oscillation of cosmic (super)strings loops could give a detectable signal in the nHz range. However, GWs generated in this way are produced at scales smaller or equal than the horizon, so that the analysis made here does not rigorously apply. Although it is sensible to expect that GWs with wavelengths comparable to the horizon would behave in a similar way to the super-horizon GWs considered here, we cannot assess the effect on the overall spectrum and we defer this generalization to a future work. This leaves the GWs produced during inflation by the amplification of quantum vacuum fluctuations, that are too weak to be detected by ongoing PTA projects but could be within the reach of future experiments like SKA. In any case, an independent observation of the signal at larger frequencies, like those probed by interferometers, would also be useful to normalize the spectrum at the source. Secondly, the cosmological signal should be clearly separated from the astrophysical signal, like that produced by black hole binaries. The third point is the more problematic. We have shown that the effects of the damping are more evident at frequencies between 0.1 and 1 nHz. As explained above, the lower limit to the detectable frequency is given by the time span of the observations, currenly 30 years. In 70 years from now, with a total of 100 years of observations, the lower limit will be f = 1/(100years) ≃ 0.3 nHz, where the intensity of the wave is reduced by only 5% and the deviation from a featureless power law is ∆α = 0.15. The maximum change of the slope is at f ≃ 0.1 nHz, corresponding to 300 years of observations. It should also be noted that, even if we were willing to wait such a long time, it is not certain that it would actually increase the sensitivity, because of intrinsic low-frequency instabilities in the timing data. In conclusion, we think that the possibility to detect this effect are quite small, and that they will depend crucially on the existence of a sizeable inflationary-like GW background at the frequencies of interest and on the capability of estimating the logarithimic slope of the spectrum in the region just below 1 nHz.