Imprint of baryons and massive neutrinos on velocity statistics

We explore the impact of baryonic effects (namely stellar and AGN feedback) on the moments of pairwise velocity using the Illustris-TNG, EAGLE, cosmo-OWLS, and BAHAMAS suites of cosmological hydrodynamical simulations. The assumption that the mean pairwise velocity of the gas component follows that of the dark matter is studied here at small separations, and we find that even at pair separations of 10-20 $h^{-1}\mathrm{Mpc}$ there is a 4-5% velocity bias. At smaller separations, it gets larger with strength varying depending on the subgrid prescription. By isolating different physical processes, our findings suggest that the large scale velocity bias is mainly driven by stellar rather than AGN feedback. If unaccounted for, this velocity offset could possibly bias cosmological constraints from the kinetic Sunyaev-Zel'dovich effect in future cosmic microwave background (CMB) surveys. Furthermore, we examine how the first and the second moment of the pairwise velocity are affected by both the baryonic and the neutrino free-streaming effects for both the matter and gas components. For both moments, we were able to disentangle the effects of baryonic processes from those of massive neutrinos; and below pair separations of 20 $h^{-1}\mathrm{Mpc}$, we find that these moments of the pairwise velocity decrease with increasing neutrino mass. Our work thus paves a way in which the pairwise velocity statistics can be utilised to constrain the summed mass of neutrinos from future CMB surveys and peculiar velocity surveys.


Introduction
Over the last decade or so, cosmology has evolved to a state where we are able to precisely constrain the cosmological parameters with the help of galaxy redshift surveys (e.g. eBOSS Collaboration 2020), gravitational lensing surveys (e.g. Heymans et al. 2020) and cosmic microwave background (CMB) experiments (e.g. Planck Collaboration 2020). Some of the outstanding questions which remain are regarding the dark sector, including determining the nature of dark energy and the summed mass of neutrinos. In order to answer these questions, peculiar velocity surveys provide a complementary avenue to further our understanding. Forthcoming peculiar velocity surveys, such as the Taipan galaxy survey 1 (da Cunha et al. 2017), the Widefield ASKAP L-band Legacy All-sky Blind Survey 2 (WAL-LABY, Koribalski et al. 2020), and the Westerbork Northern Sky HI Survey (WNSHS), promise to be competitive as cosmological probes for very low redshifts with the respect to current galaxy clustering surveys (Koda et al. 2014;Howlett et al. 2017).
The current lower limit of sum of neutrino masses, M ν = m ν 0.06 eV, comes from the neutrino oscillation experiments (e.g. Forero et al. 2014;Gonzalez-Garcia et al. 2016;Capozzi et al. 2017;de Salas et al. 2017). Massive neutrinos, unlike dark matter or baryons, have large thermal velocities which imprint distinct signatures on the cosmological observables. Leveraging this by combining different cosmological probes serves as an upper limit on the neutrino mass constraints. Depending on which datasets are combined and how the analysis is done, the current upper limit on the summed Send offprint requests to: Joseph Kuruvilla, e-mail: joseph.kuruvilla@universite-paris-saclay.fr 1 https://www.taipan-survey.org 2 https://www.atnf.csiro.au/research/WALLABY/ neutrino mass ranges from 0.12 eV up to ≈0.4 eV (e.g Di Valentino et al. 2016;Vagnozzi et al. 2017;McCarthy et al. 2018;Poulin et al. 2018;Palanque-Delabrouille et al. 2020;Ivanov et al. 2020;Planck Collaboration 2020). The impact of massive neutrinos on clustering statistics in real and redshift space has been studied (e.g. Saito et al. 2008;Wong 2008;Castorina et al. 2015;Villaescusa-Navarro et al. 2018;García-Farieta et al. 2019). Forthcoming redshift surveys will provide tighter constrains on M ν using two-point and three-point galaxy clustering statistics (Chudaykin & Ivanov 2019). Furthermore, the bispectrum should help breaking neutrino mass and σ 8 degeneracy (Hahn et al. 2020). In addition to clustering statistics, the one-point probability distribution function of the total matter has been shown sensitive to neutrino mass and could provide strong constraints (Uhlemann et al. 2020).
In this era of precision cosmology, it is important to consider the effects of baryons and processes associated with galaxy formation (e.g. cooling and feedback) on cosmological observables, particularly as we push the analyses to smaller, 'non-linear' scales. It has been shown, for example, that dark energy constraints can be biased by baryonic effects if they are unaccounted for (e.g. Semboloni et al. 2011;Copeland et al. 2018). So far, much of the attention has been focused on the impact of baryons on the clustering statistics, for example in the case of two-point statistics in Fourier space (e.g. van Chisari et al. 2018Chisari et al. , 2019Schneider et al. 2019;van Daalen et al. 2020) and in configuration space (van Daalen et al. 2014). Mummery et al. (2017) has shown that the effects of baryon physics (feedback) and neutrino free-streaming are separable (multiplicative), to typically a few percent accuracy, in their effects on the clustering statistics, even on deeply non-linear scales. The effect of baryons on the matter bispectrum have also recently been exam-Article number, page 1 of 10 arXiv:2010.05911v1 [astro-ph.CO] 12 Oct 2020 A&A proofs: manuscript no. main ined (e.g. Foreman et al. 2020). These studies were done with the aid of cosmological hydrodynamical simulations. Note that an alternative approach is to use the so-called "baryonic correction model", where the baryonic effects are parameterised based on physically-motivated parameters and used to modify the outputs of cosmological N-body simulations (e.g. Schneider & Teyssier 2015;Schneider et al. 2019;Aricò et al. 2020b,a).
The main aim of the present paper is to understand the effects of baryonic processes and massive neutrino effect on velocity statistics, namely on the first two moments of the pairwise velocity statistics, at pair separations below 20 h −1 Mpc. Relatively little attention has been devoted to the impact of baryons and neutrinos on the velocity statistics to date, particularly in comparison to the numerous studies on the spatial distribution of matter/haloes. As we will describe in the following paragraphs, the pairwise velocity has applications mainly in three areas of cosmology: (i) Galaxy clustering: the observed positions of galaxies are perturbed from their true positions due to their peculiar velocities, an effect which is known as 'redshift-space distortions' (RSD). These distortions can be leveraged to accurately constrain the growth rate of structure, and hence cosmological parameters, by measuring correlation functions in redshift space (Percival & White 2009). In configuration space clustering, the state-of-the-art modelling is based on the 'streaming model' (Peebles 1980;Fisher 1995;Scoccimarro 2004;Kuruvilla & Porciani 2018;Vlah & White 2019), recently generalised to npoint correlation function in redshift space (Kuruvilla & Porciani 2020). In two-point clustering, it provides a framework to map the two-point correlation function in redshift space, which is obtained as the integral of the real-space isotropic correlation function with the line-of-sight (los) pairwise velocity distribution. The key element in this streaming model framework is the pairwise los velocity distribution. Thus, understanding how the pairwise velocity statistics are affected by baryons and neutrinos will further help in modelling small-scale redshift-space clustering statistics. Within the streaming model framework, Aviles & Banerjee (2020) have recently studied the effects of neutrinos on pairwise velocity statistics and redshift-space correlation function using Lagrangian perturbation theory, above scales of 20 h −1 Mpc.
(ii) Peculiar velocity surveys: direct measurements of the peculiar velocity can be achieved through redshifts and distances determined through scaling relations, such as the Tully-Fisher (Tully & Fisher 1977) or the Fundamental Plane relations (Djorgovski & Davis 1987;Dressler et al. 1987). These direct peculiar velocity surveys are shallow and thus offers an opportunity to probe the peculiar velocities in the nearby Universe. In Dupuy et al. (2019), the mean pairwise velocity estimator was used to constrain the growth rate of the structure using the Cosmicflows-3 dataset (Tully et al. 2016).
(iii) Kinetic Sunyaev-Zeldovich (kSZ) effect: a secondary anisotropy where CMB photons are scattered off free electrons which are in motion. This results in a Doppler shift, thus preserving the blackbody spectrum of the CMB (Sunyaev & Zeldovich 1972, 1980. The fluctuation in the CMB temperature can be written as where σ T is the Thomson scattering cross section, n e is the electron number density, u e is the velocity of the free electrons,n is the unit vector along the line-of-sight, and τ = σ T dl n e is the optical depth. It is also one of the techniques through which we can measure the peculiar velocities of objects at cosmological distances. However, the signal from kSZ effect is very weak, hence detections for individual objects have proven to be difficult so far. Currently, detections of kSZ are mainly limited to the mean pairwise velocity, as it can be measuring through stacking techniques to boost the signal. The first detection of kSZ effect through the pairwise mean velocity was by Hand et al. (2012) using the pairwise velocity estimator developed by Ferreira et al. (1999). Further evidence for kSZ using pairwise velocities were presented in Planck Collaboration 2016; Soergel et al. (2016); De Bernardis et al. (2017); Li et al. (2018). It has been shown that the mean radial pairwise velocity measured from the kSZ effect is capable of constraining alternative theories of gravity and dark energy (Bhattacharya & Kosowsky 2007, 2008Kosowsky & Bhattacharya 2009;Mueller et al. 2015a), in addition placing constraints on the summed mass of neutrinos (Mueller et al. 2015b). Alternatively, kSZ effect have been detected by correlating CMB maps with reconstructed velocity field (e.g. Schaan et al. 2016;Tanimura et al. 2020;Nguyen et al. 2020 (Abazajian et al. 2016) and CMB-HD (Sehgal et al. 2019a,b), will be able to measure the kSZ effect, and in turn the pairwise velocity statistics, much more precisely. As already noted, the aim of this paper is to disentangle the effects baryonic processes and massive neutrino on the first two moments of the pairwise velocity. We also examine the typical assumption that the pairwise velocity of gas follows that of the dark matter for the mean pairwise velocity. This assumption been tested for pairwise kSZ signal in Flender et al. (2016) using halos from N-body simulations and adding a gas profile following a model introduced in Shaw et al. (2010). However in this paper we follow the gas particles directly from hydrodynamical simulations.
The paper is structured as follows. In Sect. 2, we briefly summarise the various hydrodynamical simulations employed in this work. In Sect. 3, we introduce the radial pairwise velocity. We introduce its first moment, the mean radial pairwise velocity, and how it is impacted by different baryonic processes in Sect. 4. We focus on how massive neutrinos affect the first moment in Sect. 4.2, and the second moment in Sect. 5. Finally, we summarise our findings in Sect. 6.

Simulations
In this work, we make use of four suites of hydrodynamical simulations to measure the pairwise velocity statistics: Illustris-TNG, EAGLE, cosmo-OWLS and BAHAMAS. We briefly describe these simulations below.
1. Illustris-TNG: 'The Next Generation Illustris Simulations' Marinacci et al. 2018;Naiman et al. 2018;Pillepich et al. 2018a;Nelson et al. 2018) is a suite of cosmological simulations run using the moving mesh code arepo. It is a successor to the Illustris simulation (Vogelsberger et al. 2014a,b;Genel et al. 2014;Sijacki et al. 2015). The subgrid physics has been updated from the original Illustris with changes in AGN feedback, galactic winds and inclusion of magnetic fields, which are described in detail in Weinberger et al. (2017) and Pillepich et al. (2018b). The feedback processes were calibrated to Table 1. Characterisation of the various simulations used in this work. BAHAMAS (0) refers to the reference simulation with zero neutrino mass. The length of the simulation is denoted by L box . While m DM and m b denotes the mass of the dark matter and baryon species respectively.

Simulation
Hydrodynamical code Gadget 67.8 6.7 × 10 6 1.8 × 10 6 Planck 2013 Illustris-TNG100 arepo 75.0 5.1 × 10 6 0.9 × 10 6 Planck 2016 Illustris-TNG300 arepo 205.0 3.9 × 10 7 7.4 × 10 6 Planck 2016 cosmo-OWLS Gadget 400.0 3.7 × 10 9 7.5 × 10 8 WMAP7 BAHAMAS (0) Gadget 400.0 3.8 × 10 9 7.6 × 10 8 WMAP9 roughly reproduce several observed properties, such as the galaxy stellar mass function and the stellar-to-halo mass relation (see Pillepich et al. 2018b for details). This suite has simulations with three different volumes 50 3 , 100 3 and 300 3 Mpc 3 . In this work, we make use of the simulation boxes with side lengths of 100 and 300 Mpc, which have 1820  Table 1. Unlike the other suites used here, no attempt was made to calibrate the feedback to match particular observations with cosmo-OWLS. It was aimed at exploring the impact of large variations in the subgrid physics, including turning on or off physics such as radiative cooling and AGN feedback. The simulation suite adopts a WMAP7 cosmology, which is given by {Ω m , Ω b , Ω Λ , h, n s , σ 8 } = {0.2720, 0.0455, 0.7280, 0.7040, 0.9670, 0.8100}. 3. Evolution and Assembly of GaLaxies and their Environments (EAGLE, Schaye et al. 2015;Crain et al. 2015) is a set of cosmological hydrodynamical simulations evolved using gadget-3. The implemented subgrid physics is descended from OWLS but with several improvements as detailed in Schaye et al. (2015). The stellar and AGN feedback was calibrated to reproduce the present-day galaxy stellar mass function and the size-mass relation of galaxies. The hydro solver scheme was also modified from classic SPH to the pressureentropy 'Anarchy' scheme, also described in the above ref- erences 2017,2018). This was also run using gadget-3. It follows the evolution of 1024 3 DM and gas particles. Hence the mass resolution is lower than EAGLE or Illustris-TNG but is approximately the same as cosmo-OWLS. And like cosmo-OWLS, it follows significantly larger volumes than EAGLE or Illustris. The subgrid physics is based on the OWLS and cosmo-OWLS projects. However, unlike OWLS and cosmo-OWLS, the feedback was explicitly calibrated to reproduce the observed present-day galaxy stellar mass function and the amplitude of the hot gas-halo mass relation of groups and clusters. As BAHAMAS has the most realistic representation of baryons on large scales (including the gas fractions of massive groups and clusters), we expect the impact on large-scale structure to be more realistic for BAHAMAS. The reference simulation we use adopts the WMAP9 cosmology, which is given by .2793, 0.0463, 0.7207, 0.7000, 0.9720, 0.8211} with massless neutrinos.
We also use an extension of BAHAMAS that includes massive neutrinos (see McCarthy et al. 2018 for details). It consists of four simulations ranging from the lowest summed neutrino masses (M ν ) of 0.06 eV up to 0.48 eV in factors of 2. The massive neutrinos were implemented keeping all the cosmological parameters fixed apart from σ 8 (note that A s , the amplitude of the primordial power spectrum was kept fixed at the CMB value, consequently the inclusion of massive neutrinos lowers the σ 8 ) and the cold matter density (Ω cdm ), which was was decreased to ensure that the Universe is flat, where Ω Λ + Ω m = 1 and Ω m = Ω cdm + Ω ν + Ω b . The neutrino density (Ω ν ) is related to M ν by the relation Ω ν = M ν /(93.14 h 2 eV). Thus, the BAHAMAS simulation explores Ω ν for a range of 0.0013 to 0.0105. This suite allows us to study the degeneracy between baryonic physics and massive neutrino effects on the pairwise velocity statistics in a systematic way.
As noted above, the different sets of simulations have been calibrated employing various strategies. The characterisation of the simulations in terms of the box size, number of particles and mass resolution also varies between suites, as shown in Table 1. Despite the different underlying cosmological models in these simulations, we neglect the impact of cosmology on the baryonic effects. Previous work based on extensions of the BAHAMAS suite (e.g., Mummery et al. 2017;Stafford et al. 2020;Pfeifer et al. 2020) have shown that the effects of baryon physics are separable from changes in cosmology at the few percent level for most statistics. We have also verified (below) that the impact of fixed baryon physics on the pairwise velocity statistics is unaffected as the cosmology is changed to increase the summed mass of neutrinos.
It should be noted that there are corresponding collisionless simulations for each of the hydrodynamical runs above, including all of the massive neutrino cases.

Radial pairwise velocity
The observed galaxy velocities provide a biased view of the unbiased (and unobserved) total matter velocity field, u m . This unbiased velocity field can be defined as the fractional sum of its basic components where are the cold dark matter, the baryon and the neutrino fraction respectively. The velocity of the cold dark matter and the neutrino are denoted by u cdm and u ν respectively. Whereas the velocity of the baryons, u b , is further obtained as a fractional sum of the velocities of gas, stars and black holes (BH) where f i represents the fraction of gas, stars and BH. The radial component of the pairwise velocity is This can be measured directly from the simulations. In order to build the radial pairwise velocity distribution function (RPVDF) from the simulations, we randomly sample 192 3 tracer particles. In Fig. 1, we plot the RPVDF of the various components from the BAHAMAS simulation for pairs with a separation of (1, 2) h −1 Mpc and (40, 41) h −1 Mpc on top and bottom panels respectively. The solid lines denote the PDF for the matter from the corresponding collisionless simulation. The dashed, dotted and dash-dotted lines are for the gas, stars and dark matter species from the hydrodynamical simulation. It is evident from the PDF that the pairwise velocity information, i.e. all the moments, derived from the gas and DM particles are different at the scales shown. This is important as many studies of the kSZ effect normally assume that the gas perfectly traces the dark matter on large scales. It should be noted that w r < 0 denotes the pairs which are infalling towards each other, while w r > 0 implies that they are moving away from each other. RPVDF of both components are visibly skewed left, while the tails are much heavier for the dark matter component when compared to the gas component.

First moment of radial pairwise velocity
In this section, we compute the first moment of the radial pairwise velocity from the simulations. In the single stream regime, the mean radial velocity can be defined as: where δ represents the mass density contrast and r gives the pair separation vector. Using standard perturbation theory at leading order, it can be shown that (Fisher 1995) where f is the growth rate, j 1 (x) = sin(x)/x 2 − cos(x)/x and P(k) is the linear matter power spectrum. The particles in a pair tend to approach each other on average, i.e. w r (r) < 0, due to gravitational attraction. In Fig. 2, we explore the mean radial pairwise velocity for matter, DM and gas components from all the simulations mentioned before. Similar to building the RPVDF, we randomly sample 192 3 particles for each tracer (except in the case of BH particles) to compute the moments. To quantify the uncertainty of our measurements, we create three such catalogs and use the standard error of the mean. The top panel shows the effect of different subgrid physics on the gas radial pairwise velocity. The curves show that all four simulations follow a similar qualitative trend, in that the pairwise velocity of the gas is suppressed relative to the collisionless dark matter, particularly on small scales. However, the magnitude of this effect varies strongly from simulation to simulation. At intermediate scales of 1-10 h −1 Mpc, BAHAMAS (dotted lines) shows the maximal deviation of about 30% from the assumption that the gas follows the mean velocity of the dark matter. While both the Illustris-TNG runs show a maximal effect of 10-18% at the same scales. EAGLE shows an effect of about 10% (at most) on these intermediate scales. However on the smallest scales considered, EAGLE shows the largest effect, with the gas pairwise velocity deviating by up to 42% from the collisionless dark matter. It should be noted that on all scales considered here the ratio does not go to one, which implies that there is a velocity bias between the dark matter and gas component even on the largest scales that we measure. On below scales of 10 h −1 Mpc, the linear velocity bias approximation of mean pairwise velocity clearly breaks down. Intriguingly, this holds true for all the simulations we have considered with varying simulation volumes, thus suggest that this is robust to changes in the box size. The middle panel of Fig. 2 displays the ratio of the radial pairwise velocity of the dark matter component from the full physics run to the matter component from the collisionless simulation. This comparison tells us how dark matter responds to baryons in the full hydro runs. On pair separations of about 1-3 h −1 Mpc, the BAHAMAS shows a clear back-reaction effect whereby the dark matter component in the full physics run is infalling towards each other at a greater pace than its counterpart in the collisionless simulation. This trend is also seen in the Illustris-TNG100 simulation.
The bottom panel shows the effect of baryons on the total matter pairwise velocity. The Illustris-TNG runs are within 1% at scales above 2 h −1 Mpc. The matter mean radial pairwise velocity in EAGLE simulation behaves similar to Illustris-TNG at those scales and is affected by ≈ 1% at most. This is also in tandem with the findings in Hellwing et al. (2016), who find the effect of baryons on redshift-space clustering to be minimal in EAGLE. However, at small pair separations (≤ 1 h −1 Mpc), matter seems to infall towards each other faster around 0.5-1 h −1 Mpc and this trend reverses quickly at smaller scales. At the intermediate scales, BAHAMAS deviates at about the 2-3% level. This hints towards the possibility that the redshift-space clustering in BAHAMAS will be affected by baryonic effects to a larger degree than in EAGLE (as shown in Kwan et al., in prep). The fact that BAHAMAS produces a larger effect relative to EAGLE and Illustris-TNG is perhaps not that surprising, as the AGN is more effective at removing baryons from galaxy groups and clusters in BAHAMAS. This is a result of explicit calibration of the AGN feedback to reproduce the observed baryon fractions of massive systems, whereas neither EAGLE nor Illustris-TNG were calibrated on these data and predict baryon fractions in excess of that observed on mass scales of ∼ 10 14 M . By considering the matter pairwise velocity, we study the unbiased velocity field. To directly translate these effects to RSD measurements from redshift surveys, we would need to study the galaxy pairwise velocity statistics which we do not consider in this work. So far we have seen how the different baryonic models in the simulations affect the velocity statistics. However, we want to isolate the effect of different physical processes, such as AGN feedback. For this purpose, we use two different feedback runs from BAHAMAS with varying AGN subgrid heating temperatures, in addition to the reference run. The 'high-AGN' run has ∆T AGN = 10 8.0 K, while the 'low-AGN' run was run with ∆T AGN = 10 7.6 K. These values were chosen so that the simulations roughly bracket the upper and lower bounds of the observed hot gas fraction-halo mass relation inferred from X-ray observations (McCarthy et al. 2018). They therefore represent a kind of estimate of the allowed range of behaviours for models with AGN feedback. In Fig. 3, the solid and the dashed lines represent the ratio of matter and gas pairwise velocity with respect to the collisionless matter counterpart, respectively. The gas elements are pushed away from each other more strongly as the AGN temperature increases. This causes a stronger decrease in the gas radial pairwise velocity for the high-AGN model as can be seen. Similarly, the matter is also affected in the same manner. The high-AGN feedback causes the matter from the hydrodynamical simulation to deviate further from its counterpart in the collisionless simulation. It should, however, be noted that despite the fact that EAGLE has a higher AGN temperature, the effect of AGN heating is more prominent in BAHAMAS than in EAGLE. This can be attributed to differences in the mass resolutions of the two simulations, whereby each heating event in BAHAMAS deposits significantly more energy and thus results in a stronger expulsion. This has also been seen in the case of galaxy clustering information (Foreman et al. 2020). It would be interesting to run a high-resolution simulation such as EAGLE but to heat a similar volume/mass as that in BAHAMAS to see whether the effects are similar when the feedback is forced to operate in a similar way.
To further explore the impact of different physical processes, we also make use of the cosmo-OWLS simulations in Fig. 4. The dashed double-dotted line refers to the 'NoCool' simulation in cosmo-OWLS where there is no radiative cooling, star formation, stellar feedback or AGN feedback (there is only net photoheating from a UV/X-ray background). We see that in this case, the bias is nearly one on scales larger than 10 h −1 Mpc. This implies that it is indeed the physics of galaxy formation that is responsible for the velocity bias on large scales in the previously explored simulations. Turning on the cooling, star formation and stellar feedback, while keeping AGN feedback turned off ('NoAGN'), we see that this introduces a bias even at scales larger than 10 h −1 Mpc. This shows that physical processes like stellar feedback, prevents the gas from infalling. The fact that the bias on large scales is similar to that of runs that also include AGN feedback strongly suggests that it is stellar feedback, rather than AGN feedback, that is mainly responsible for the large-scale bias.
The dash-dotted line shows the effect of AGN feedback in cosmo-OWLS, two models of which have a much higher heating temperature than considered in the case of the BAHAMAS simulation. As a result, these models clearly expel gas away from each other to a much larger degree. We note, however, that the two highest heating temperature runs from cosmo-OWLS yield gas fractions significantly lower than observed on the scale of groups and clusters, implying that the feedback is somewhat too aggressive in those runs.

Redshift evolution
We also explore the effect of redshift evolution on the pairwise statistics. For this exercise we use the reference simulation from BAHAMAS with massless neutrinos and measure the mean pairwise velocities at redshifts 0.0, 0.5, 1.0 and 2.0. In Fig. 5, we show the effect of baryonic physics on matter mean radial pairwise velocity. The feedback is most efficient at higher redshifts at smaller scales (r < 1 h −1 Mpc), reaching a deviation of up to 9% for the matter fluid when compared to its collisionless matter counterpart. At z = 0 (denoted by solid line), we see the back reaction of DM having an effect on the matter. At scales above 1 h −1 Mpc, the ratio reaches a maximal deviation of ∼ 3%. Thus, these baryonic effects will be important to understand if we are to push the modelling of mean pairwise velocity to nonlinear scales and earlier times for forthcoming redshift surveys like Euclid. The gas elements show an even more pronounced effect when compared to the matter from the collisionless simulation, as shown in Fig. 6. At the highest redshift considered here (z = 2), the baryonic effects on the gas elements, denoted by dotted lines, strongly affect scales below 3 h −1 Mpc. Moving towards lower redshift, this effect is reduced in amplitude but more extended in scale, being seen on scales as large as ∼ 10 h −1 Mpc. It is again worth highlighting the fact that the velocity bias between the gas and the collisionless matter is below one at all scales considered here and at all times. This is a clear indication that one needs to be careful about the assumption that mean radial velocity of gas follows that of the dark matter at scales about 20 h −1 Mpc and below, especially for precise measurements in the future.
For comparison, the dashed double dotted line denotes the trend at z = 2 from the Illustris-TNG300 simulation. The effect of AGN feedback in this simulation is strongly reduced compared to BAHAMAS, as was also deduced from the z = 0 comparison previously.

The effects of massive neutrinos
Constraining neutrino mass is one of the primary objectives of forthcoming galaxy and CMB surveys. One of the main effects of neutrinos on the two-point clustering statistics in Fourier space (i.e. the power spectrum), is the damping of power on scales smaller than the free-streaming scale. Neutrinos will also affect velocity statistics (e.g. Mueller et al. 2015b). We focus on the mean radial velocity to exhibit the effects of neutrinos. Specifically, we show how they affect the matter mean pairwise velocity in Fig. 7. We see that the main effect of neutrinos is to reduce the mean pairwise velocity when compared with a massless neutrino simulation, implying that as the sum of neutrino mass increases, the infall of matter towards each other decreases. Physically, this is due to the fact that the neutrino component does not significantly cluster on scales below the free-streaming scale which in turn slows the collapse of the dark matter and baryons. Considering pair separation scales above 3 h −1 Mpc, we can see that the effect reaches approximately 20% on the matter component for the M ν = 0.48 eV. This will have important consequences for the RSD signal and hence on the redshift space clustering. Even for the most stringent of current constraints on the neutrino mass (M ν < 0.12 eV), the radial pairwise velocity of matter will be affected at the 3-5% level. This is also encouraging for the future peculiar velocity surveys, which might be able to provide independent constraints on the sum of neutrino masses.
However to decouple the effects of baryons and neutrinos using a single simulation Fig. 7 is non-trivial at small scales as the effects are intertwined with each other. Since we have a series of massive neutrino simulations from BAHAMAS (both hydro and collisionless for each neutrino mass), it is possible to disentangle the effects of baryonic physics and massive neutrinos. For this, we introduce the ratio statistics as follows.
where the velocity biases B (1) i (r) and N (1) (r) capture the effects of baryons and neutrinos, respectively, and the subscript i represents either the gas or the matter component. In Fig. 8, we show the velocity biases due to baryonic effect and massive neutrinos for the gas component, in the middle panel and the bottom panel respectively. The advantage of this approach is that we can treat these effects separately. In the future, one can build emulators for B (1) i (r) and N (1) (r) separately and combine them. The top panel shows the LHS of the equation (8). Similar to the matter, massive neutrinos reduce the mean pairwise velocity of gas component at scales above 3 h −1 Mpc and hence the gas velocity bias also decreases as the neutrino mass increases. The function B (1) gas (r) is roughly constant above 10 h −1 Mpc, below which the baryonic physics starts to have an effect. It can also be seen that the effect from the baryonic processes remain largely unchanged for different neutrino mass cosmology. We have also numerically verified that equation (8), holds true for the gas component. At the smallest pair separation considered the relative difference between the LHS and RHS is 0.1%, and at the largest separation it is roughly 10 −5 %.

Second moment of radial pairwise velocity
In this section, we focus on the second moment of the pairwise velocity and check how it is affected by massive neutrinos and the effects of baryons. In the single stream regime, we can define the second moment of the pairwise velocity as We are interested in the radial component which can be defined using standard perturbation theory at leading order as where ψ r (r 12 ) = f 2 2π 2 ∞ 0 j 0 (k r 12 ) − 2 j 1 (k r 12 ) k r 12 P(k) dk (11) is the radial velocity correlation function (Gorski 1988) with j 0 (x) = sin(x)/x, and is the one-dimensional velocity dispersion.
In Fig. 9, we show the direct measurement of the second moment directly from the BAHAMAS suite of simulations. In the case of the massless neutrinos, the matter in the hydrodynamical simulation has a smaller dispersion compared to the matter in the collisionless simulation. At the largest separation considered here, the second moment is reduced by 8-9%. As noted already in the first moment, the increasing neutrino mass decreases the velocity dispersion of the radial pairwise velocity. For the most massive neutrino case considered here, the pairwise dispersion is reduced by 25-40% when compared to the matter from the massless neutrino collisionless simulation. Understanding and accounting for this effect will be important in modelling RSD using the streaming model framework if we want to use clustering analysis at non-linear scales. To disentangle the effects of baryons and neutrino on the gas dispersion, we can write where B (2) gas (r) and N (2) (r) are the velocity biases due to baryons and neutrinos, in the context of the pairwise velocity dispersion. In Fig. 10, we show the LHS (top panel) and RHS (middle and bottom panels) terms of equation (13). The top panel shows that the pairwise velocity dispersion of the gas component in the massless neutrino cosmology is significantly less than that of the dark matter, by more than 40% at all scales considered. The dispersion decreases further as the neutrino mass is increased, as expected. This is encouraging as we can leverage the dispersion measure of the pairwise velocity from kSZ or peculiar velocity measurements to further constrain the summed mass of neutrinos.
In the middle panel, the effect of baryons is nearly invariant for the different neutrino cosmologies, although for the most massive neutrino case the baryonic effects differ by 1-3% at pair separations of around 10-20 h −1 Mpc from the massless neutrino case. The bottom panel effectively shows the impact of neutrinos on the velocity dispersion for dark matter species in the collisionless simulation. The most massive neutrino case causes a decreases of about 30% even at the largest separations, while a summed neutrino mass of 0.12 eV (dash-dotted line) and 0.24 eV (dotted line) show decreases in the pairwise velocity dispersion of approximately 4% and 8%, respectively.

Conclusions
In this study, we have focused on the imprint of baryons and neutrinos (and their interplay) on the first two moments of the radial pairwise velocity distribution. Understanding these effects will help us to alleviate any potential biases in constraining cosmological parameters, in particular the neutrino mass, from the future surveys.
The assumption that the mean pairwise velocity of gas component follows that of the dark matter is a crucial one undertaken in kSZ analyses. In Fig. 1, we demonstrated that even on large pair separations, r ∈ (40, 41) h −1 Mpc, the radial pairwise velocity distribution of the gas component differs from that of the dark matter.
Focusing on its first moment, we demonstrated that different subgrid models lead to different effects on the mean radial pairwise velocity statistics, especially on the very small scales below 1 h −1 Mpc, as can been seen in Fig. 2. We also see that even at pair separations of 15-20 h −1 Mpc, there is a pairwise velocity bias between the gas and dark matter. This indicates that the assumption that the mean pairwise velocity of gas follows that of dark matter breaks down at these scales.
We further studied the effect of AGN feedback in particular on the mean pairwise velocity in Fig. 3 using the BAHAMAS simulations, finding that more energetic AGN heating pushes the matter away leading to a decrease in the mean infall of material. In Fig. 4, we studied the effect of different baryonic processes using cosmo-OWLS suite of simulation. The assumption that the gas follows the dark matter (above scales of 10 h −1 Mpc) is valid only in the case when all non-gravitational physical processes like radiative cooling, star formation and stellar and AGN feedback are switched off. In the cases when those physical processes were switched on, the assumption breaks down for the pair separations we have considered. The source of the large-scale velocity bias appears to be driven by the stellar feedback rather the AGN feedback as suggested in Fig. 4. Turning AGN feedback on does not significantly alter this, but it does greatly affect intermediate scales. Thus the strength of the variation changes according to the subgrid physics considered.
The impact of baryonic process at different redshifts is studied using the BAHAMAS reference simulation in Figs. 5 and 6. We see that even at the highest redshift considered in our study z = 2, the baryonic processes introduce one percent level impact on matter mean pairwise velocity at scales above 10 h −1 Mpc. In the case of gas component, the impact is more prominent and introduces 4-5% change in the mean velocity with respect to the matter in a gravity-only calculation.
We studied the effect of massive neutrinos on mean radial pairwise velocity using the BAHAMAS suite of simulations. We showed that the matter mean pairwise velocity decreases as the summed neutrino mass increases, in Fig. 7. Though we studied the (unbiased) matter velocity field, these results suggest that the radial pairwise velocity could be used to potentially constrain neutrino mass from peculiar velocity surveys in the future, and in addition these effects could be important in modelling RSD using the streaming model framework in the presence of massive neutrinos. In Fig. 8, we disentangled the baryonic and massive neutrino effects on the mean radial pairwise velocity of gas component as introduced in equation (8). For the most massive neutrino considered in this work (M ν = 0.48 eV), we found that the mean radial pairwise velocity of gas decreases by roughly 20% when compared with dark matter in the massless neutrino simulation. The baryonic effect is nearly invariant when considering different neutrino mass simulations. Finally, we demonstrated the effect of neutrinos on the second moment of the radial pairwise velocity for both matter and and gas components, as shown in Figs. 9 and 10 respectively. Similar to the mean radial pairwise velocity, the second moment also decreases with increasing neutrino mass. The matter pairwise velocity is reduced by ∼15 % for M ν = 0.12 eV when compared to the massless neutrino case at pair separations of 10-20 h −1 Mpc. At the same separation, for the highest neutrino mass considered the impact is reduced by ∼35 %. This points towards the possibility of utilising the pairwise dispersion (as a function of pair separation) to constrain neutrino mass from either future peculiar velocity surveys or future CMB surveys using kSZ effect. The second moment would also be beneficial for breaking degeneracies between cosmological parameters. For example, the mean pairwise velocity scales as f σ 2 8 , while the second moment scales as ( f σ 8 ) 2 . Direct application of our results to either peculiar velocity surveys or kSZ would require us to study the effect separately on galaxies/haloes, which we reserve for a future study.
Thus, we have seen how different feedback models affect the moments of the pairwise velocity to varying degree. With the forthcoming peculiar velocity and CMB surveys, understanding these systematic effects from baryons and neutrinos will be essential for constraining the cosmological parameters using pairwise velocity accurately and precisely.