COSMOS2015 photometric redshifts probe the impact of filaments on galaxy properties

TL;DR

This study uses COSMOS2015 photometric redshifts and simulations to show that galaxy properties vary with proximity to cosmic filaments, independent of local density, highlighting the influence of large-scale structure on galaxy evolution.

Contribution

It demonstrates that galaxy-filament proximity affects galaxy properties beyond local density effects, using photometric redshifts and simulations to establish this link.

Findings

Massive galaxies are closer to filaments than less massive ones.

Passive galaxies are nearer to filaments compared to active star-forming galaxies.

The filament influence persists after removing local density effects.

Abstract

The variations of galaxy stellar masses and colour-types with the distance to projected cosmic filaments are quantified using the precise photometric redshifts of the COSMOS2015 catalogue extracted from COSMOS field (2 deg$^{2}$). Realistic mock catalogues are also extracted from the lightcone of the cosmological hydrodynamical simulation Horizon-AGN. They show that the photometric redshift accuracy of the observed catalogue ($\sigma_z<0.015$ at $M_*>10^{10}{\rm M}_{\odot}$ and $z<0.9$) is sufficient to provide 2D filaments that closely match their projected 3D counterparts. Transverse stellar mass gradients are measured in projected slices of thickness 75 Mpc between $0.5< z <0.9$, showing that the most massive galaxies are statistically closer to their neighbouring filament. At fixed stellar mass, passive galaxies are also found closer to their filament while active star-forming galaxies statistically lie further away. The contributions of nodes and local density are removed from these gradients to highlight the specific role played by the geometry of the filaments. We find that the measured signal does persist after this removal, clearly demonstrating that proximity to a filament is not equivalent to proximity to an over-density. These findings are in agreement with gradients measured both in 2D or 3D in the Horizon-AGN simulation and those observed in the spectroscopic VIPERS survey (which rely on the identification of 3D filaments). They are consistent with a picture in which the influence of the geometry of the large-scale environment drives anisotropic tides which impact the assembly history of galaxies, and hence their observed properties.

Tables

Table 1. Table 1: Limiting mass ( M lim subscript 𝑀 lim M_{\rm lim} ), estimation of the comoving transverse width of the slice ( D trans subscript 𝐷 trans D_{\rm trans} ) , redshift errors ( Δ ​ z ​ [ M ∗ ∼ 10 10 ] Δ 𝑧 delimited-[] similar-to subscript 𝑀 superscript 10 10 {\Delta z}\left[M_{*}\sim 10^{10}\right] ) and corresponding thickness in Mpc ( D Δ ​ z ​ [ M ∗ ∼ 10 10 ] subscript 𝐷 Δ 𝑧 delimited-[] similar-to subscript 𝑀 superscript 10 10 D_{\Delta z}\left[M_{*}\sim 10^{10}\right] ), and chosen thickness of the slices in redshift ( Δ ​ z ​ [ D = 75 ​ M ​ p ​ c ] Δ 𝑧 delimited-[] 𝐷 75 M p c \Delta z\left[D=75{\rm Mpc}\right] ) and in Mpc ( D Δ ​ z subscript 𝐷 Δ 𝑧 D_{\Delta z} ). The limiting mass is estimated as in Laigle et al. ( 2016 ) . The comoving transverse distance is estimated at the upper limit of each redshift bin and corresponds to 1.4 deg (side of the COSMOS field). M ∗ > 10 10 ​ M ⊙ subscript 𝑀 superscript 10 10 subscript M direct-product M_{*}>10^{10}{\rm M}_{\odot} is the conservative limiting mass that we choose in this work to build a mass-selected sample at 0.5 < z < 0.9 0.5 𝑧 0.9 0.5<z<0.9 . Δ ​ z ​ [ M ∗ ∼ 10 10 ​ M ⊙ ] Δ 𝑧 delimited-[] similar-to subscript 𝑀 superscript 10 10 subscript 𝑀 direct-product \Delta z\left[M_{*}\sim 10^{10}M_{\odot}\right] corresponds to twice the median 1 ​ σ 1 𝜎 1\sigma redshift error determined by LePhare for the galaxies with masses 10 10 < M ∗ / M ⊙ < 10 10.5 superscript 10 10 subscript 𝑀 subscript 𝑀 direct-product superscript 10 10.5 10^{10}<M_{*}/M_{\odot}<10^{10.5} (the faintest galaxies in the bin, for which the redshift uncertainties are the highest). D Δ ​ z ​ [ M ∗ ∼ 10 10 ​ M ⊙ ] subscript 𝐷 Δ 𝑧 delimited-[] similar-to subscript 𝑀 superscript 10 10 subscript 𝑀 direct-product D_{\Delta z}\left[M_{*}\sim 10^{10}M_{\odot}\right] is the corresponding thickness in comoving Mpc at the upper limit of the redshift bin. Finally, Δ ​ z ​ [ D = 75 ​ M ​ p ​ c ] Δ 𝑧 delimited-[] 𝐷 75 M p c {\Delta z}\left[D=75{\rm Mpc}\right] is the maximum in each redshift bin of the redshift thicknesses chosen in this work corresponding to D Δ ​ z = 75 subscript 𝐷 Δ 𝑧 75 D_{\Delta z}=75 Mpc. It has been chosen such as to encompass two times the median redshift error in the lowest mass bin at z ∼ 0.9 similar-to 𝑧 0.9 z\sim 0.9 .

bin	$M_{\lim }$	$D_{\mathrm{trans}}$	$\Delta z\left[M_{\ast }\sim {10}^{10}\right]$	$D_{\Delta z}\left[M_{\ast }\sim {10}^{10}\right]$	$\Delta z\left[D=75Mpc\right]$	$D_{\Delta z}$
		(Mpc)	($2\times 1\sigma$)	(Mpc)		(Mpc)
0.50.6	8.72	54	0.019	58.3	0.024	75
0.60.7	8.92	62	0.021	61.6	0.025	75
0.70.8	8.95	70	0.024	68.0	0.027	75
0.80.9	9.17	75	0.029	75.2	0.029	75
0.91.0	9.19	82	0.032	80.4	–	–
1.01.1	9.23	87	0.045	105.1	–	–
1.11.2	9.37	93	0.059	131.4	–	–
1.21.3	9.44	98	0.059	128.1	–	–

Table 2. Table 2: The skeletons used in this work to assess the quality of the reconstruction or for gradients measurements. M lim subscript 𝑀 lim M_{\rm lim} refers to the mass limit of the galaxy sample used for the computation, σ 𝜎 \sigma to the persistence threshold and D Δ ​ z subscript 𝐷 Δ 𝑧 D_{\Delta z} is the thickness of the slice (for two-dimensional reconstructions). For skl DM 3 ​ D superscript subscript absent 3 D DM {}_{\rm 3D}^{\rm DM} the persistence threshold is not expressed in terms of σ 𝜎 \sigma . The quoted values are those used in the main text, but Appendices investigate how the reconstruction changes when these values vary.

Name	catalogue	Dim		$\sigma$ ∗	$M_{\lim }$ ∗	$D_{\Delta z}$ ∗	comments
skl${}_{}{}^{\mathrm{obs}}{}_{2D}{}^{}$	COSMOS2015	2D	galaxies with photo-$z$	2	1010M⊙	75 Mpc	density weighted by the photo-$z$ errors
skl${}_{}{}^{\mathrm{phot}}{}_{2D}{}^{}$	Horizon-AGN	2D	galaxies with photo-$z$	2	1010M⊙	75 Mpc	density weighted by the photo-$z$ errors
skl${}_{}{}^{\mathrm{spec}}{}_{2D}{}^{}$	Horizon-AGN	2D	galaxies with exact $z$	2	1010M⊙	75 Mpc	-
skl${}_{}{}^{\mathrm{DM}}{}_{3D}{}^{}$	Horizon-AGN	3D	dark matter particles	-	$8\times {10}^7$M⊙	-	density computed on a grid
skl${}_{}{}^{\mathrm{gal}}{}_{3D}{}^{}$	Horizon-AGN	3D	galaxies	5	1010M⊙	-	-

Table 3. Table 3: Summary of the results from the mass gradients towards filaments. The given values are the medians of the distributions and the standard errors from bootstrap resampling with 100 realisations expressed in comoving Mpc. Note that the median values in three- and in two-dimensions are not directly comparable because of projection effects. The bracketed values indicate the percentage of galaxies in each mass bin. The percentages are the same for the original and reshuffled samples. The reshuffled sample corresponds to a random signal preserving the galaxy mass-density relation.

bin	data	original	reshuffled
	Hz-AGN 3D	4.75$\pm$ 0.01 [54]	4.53$\pm$ 0.02
9.5log10	Hz-AGN 2D	1.23$\pm$ 0.01 [54]	1.01$\pm$ 0.01
	COSMOS2015	1.68$\pm$0.02 [47]	1.37$\pm$ 0.02
	Hz-AGN 3D	4.46$\pm$ 0.01 [27]	4.48$\pm$ 0.03
10log10.4	Hz-AGN 2D	0.70$\pm$ 0.01 [27]	0.92$\pm$ 0.01
	COSMOS2015	0.97$\pm$ 0.02 [23]	1.21$\pm$ 0.02
	Hz-AGN 3D	3.96$\pm$ 0.03 [14]	4.34$\pm$ 0.03
10.4log10.8	Hz-AGN 2D	0.60$\pm$ 0.02 [14]	0.89$\pm$ 0.02
	COSMOS2015	0.84$\pm$ 0.02 [18]	1.18$\pm$ 0.03
	Hz-AGN 3D	3.07$\pm$ 0.08 [05]	4.18$\pm$ 0.08
10.8log$M_{\ast }/M_{\odot }$	Hz-AGN 2D	0.42$\pm$ 0.03 [05]	0.86$\pm$ 0.03
	COSMOS2015	0.79$\pm$ 0.05 [12]	1.16$\pm$ 0.04

Table 4. Table 4: Summary of the results from the colour-type gradient signal towards filaments in COSMOS2015. The given values are the medians of the distributions and the standard errors from bootstrap resampling with 100 realisations expressed in projected comoving Mpc. The bracketed values indicate the percentage of galaxies within a mass bin in each population. The percentages are the same for the original and reshuffled samples. The reshuffled sample corresponds to a random signal preserving the galaxy type-density relation.

	original sample			reshuffled samples
bin	All	SF	passive	SF	passive
10log10.4	0.97$\pm$ 0.02	1.01$\pm$ 0.04 [72]	0.82$\pm$ 0.05 [28]	0.96$\pm$ 0.03	0.96$\pm$ 0.05
10.4log10.8	0.84$\pm$ 0.02	0.96$\pm$ 0.04 [58]	0.69$\pm$ 0.03 [42]	0.86$\pm$ 0.04	0.82$\pm$ 0.04
10.8log$M_{\ast }/M_{\odot }$	0.79$\pm$ 0.05	0.95$\pm$ 0.07 [34]	0.70$\pm$ 0.04 [66]	0.81$\pm$ 0.08	0.76$\pm$ 0.05
10log$M_{\ast }/M_{\odot }$ [All]	–	0.99$\pm$ 0.03 [59]	0.73$\pm$ 0.03 [41]	0.91$\pm$ 0.02	0.86$\pm$ 0.03

Table 5. Table 5: Summary of the results from the sSFR gradient signal towards filaments in Horizon-AGN using skl DM 3 ​ D superscript subscript absent 3 D DM {}_{\rm 3D}^{\rm DM} . The given values are the medians of the distributions and the standard errors from bootstrap resampling with 100 realisations expressed in comoving Mpc. Note that the median values found in three dimensions are not directly comparable to the two-dimensional ones because of the projection. The bracketed values indicate the percentage of galaxies within a mass bin in each population. The percentages are the same for the original and reshuffled samples. The reshuffled sample corresponds to a random signal preserving the galaxy sSFR-density relation in each mass bin. Galaxies with log sSFR/yr -1 < − 10.9 absent 10.9 <-10.9 are considered as passive, while those with log sSFR/yr -1 > − 9.8 absent 9.8 >-9.8 are star-forming.

	original sample				reshuffled samples
bin (log sSFR/yr-1)	All		$[-10.9,-9.8]$			$[-10.9,-9.8]$
10log10.4	4.46$\pm$ 0.01	4.38$\pm$ 0.12 [03]	4.45$\pm$ 0.03 [75]	4.70$\pm$ 0.05 [22]	4.54$\pm$ 0.14	4.50$\pm$ 0.03	4.52$\pm$ 0.05
10.4log10.8	3.96$\pm$ 0.03	3.87$\pm$ 0.26 [07]	4.06$\pm$ 0.04 [79]	3.94$\pm$ 0.18 [14]	3.99$\pm$ 0.11	4.02$\pm$ 0.04	4.02 $\pm$0.16
10.8log$M_{\ast }/M_{\odot }$	3.07$\pm$ 0.08	2.89$\pm$ 0.10 [21]	3.23$\pm$ 0.13 [74]	3.24$\pm$ 0.13 [06]	3.09$\pm$ 0.10	3.16$\pm$ 0.07	3.10 $\pm$ 0.22
10log$M_{\ast }/M_{\odot }$ [All]	–	3.69 $\pm$ 0.05 [06]	4.22$\pm$ 0.03 [76]	4.47$\pm$ 0.05 [18]	3.91$\pm$ 0.01	4.23$\pm$ 0.02	4.39 $\pm$ 0.07

Equations

$$p_{\mathrm{gal},i}={\int }_{z_{1,i}}^{z_{2,i}}{P}_{\mathrm{gal}}(z)dz/\int P_{\mathrm{gal}}(z)dz,$$

Full text

1 Introduction

A challenge for current galaxy formation models is to understand the role of the environment, its anisotropy and its dynamics in shaping the properties of galaxies across cosmic time. Galaxies form in the gravitational potential wells of dark matter haloes. Their mass assembly can thus be approached at first order by the gravitational history of halo assembly (Blumenthal et al., 1984; Navarro et al., 1995), while at a second order, the evolution of baryonic matter (gas accretion and cooling, star formation and feedback) induces additional levels of complexity (e.g. Hopkins et al., 2014, and reference therein). Starting from the dark matter only, the theory of cosmological structure formation predicts a strong dependency of halo mass assembly on its environment, driven by gravitational physics.

On the one hand, halo mass assembly history is encoded in the initial conditions of the matter that will end up in the halo (its so-called initial Lagrangian patch). At the simplest level, it is a function of its mass which determines its collapse time encoded in the spherical collapse model (Press & Schechter, 1974). An environmental dependency will occur when also considering the impact of the long wavelength density mode corresponding to the large-scale cosmic web. As encompassed by the peak-background split theory (Sheth & Tormen, 1999a), a large-scale overdensity induces a density boost, allowing the proto-halo to pass earlier the critical threshold of collapse and biasing the mass function in the vicinity of the large-scale structure (shifted Press-Schechter, Kaiser, 1984; Efstathiou et al., 1988; Cole & Kaiser, 1989; Bond et al., 1991; Mo & White, 1996; Sheth & Tormen, 1999b): the abundance of massive halos will thus be enhanced in overdense regions (Alonso et al., 2015), and these objects will be more strongly clustered. An observational consequence of this prediction is the mass-density relation (Oemler, 1974; Davis & Geller, 1976; Balogh et al., 1997), where the local density at a given redshift remains a tracer of the initial density boost.

On the other hand, in complement to the local density, the initial geometric distribution of neighbouring structures will also explicitly impacts the accretion and merger history of the proto-halo via the tides it applies onto the patch. It has been found that at a given mass, denser environments are on average populated by older haloes, where the age of a halo is defined as the epoch by which the halo has assembled one half of its current mass. This effect is called “assembly bias” (e.g. Sheth & Tormen, 2004; Gao et al., 2005; Wechsler et al., 2006; Croton et al., 2007) and can be interpreted as follows. At a given mass, if the halo is isolated, it can continue to grow by accreting matter from its neighbourhood, via the accretion of larger shells in the initial density field. If this halo lies in the vicinity of a massive structure, it will stop to grow earlier because the matter will recede towards that massive structure. Hence the large-scale structure around the proto-halo will affect the halo history via external tides (e.g. Dalal et al., 2008; Hahn et al., 2009; Wang et al., 2011) captured by the ellipsoidal collapse model (Bond & Myers, 1996). Those tides span a large range of scales. At small scale they may induce mass losses for haloes in the vicinity of a more massive halo, while on larger scale they imply statistically distinct mass assembly histories in different large-scale environment (see also Castorina et al., 2016), e.g. by delaying or quenching the mass inflow. An important ingredient at large scale is the anisotropic structure of the matter distribution, which is organised as the cosmic web (Bond et al., 1996), made of walls, filaments and nodes. The shape of the tidal tensor111Specifically, the traceless part of the tidal tensor, which orientation via its eigenvectors points along the main axes of the embedding cosmic web. Recall that gradients of the potential displace halos, while the tides distort and rotate them along the directions set by those vectors. will be precisely determined by this anisotropy. It will affect not only the accretion history of the proto-halo as described above, but also its dynamics via tidal torquing.

This dynamical impact of the anisotropic cosmic web has been investigated in dark matter simulations in relation to spin acquisition: Aubert et al. (2004b); Aragón-Calvo et al. (2007); Hahn et al. (2007); Codis et al. (2012); Aragon-Calvo & Yang (2014) report angular momentum alignment of low-mass haloes with filaments and a perpendicular orientation for more massive haloes. This alignment can be understood when measuring the vorticity content in walls and filaments inherited from flow-crossing (Libeskind et al., 2013; Laigle et al., 2015). Structures which are forming in these large scale flows will end up with a spin preferentially aligned with that vorticity and consequently with the filament. They will accrete matter through secondary infall with a coherent rotational motion, up to a specific transition mass corresponding to the Eulerian size of the quadrant, beyond which inflow will start to advect misaligned angular momentum. This dynamical connection can also be predicted from the initial conditions and its variation with the geometry of the environment is now well established in the context of the constrained tidal torque theory (Codis et al., 2015b). Interestingly, this dependency still holds for galaxies, despite the scale difference and the more complex baryonic physics in the circumgalactic-medium. Tempel & Libeskind (2013) have observed the alignment of galaxy’s spin with their closest filament in the SDSS survey (see also e.g. Hirv et al., 2016) and Dubois et al. (2014) have first confirmed this dependency for galaxies in hydrodynamical simulations (see e.g. also Codis et al., 2015a; Chisari et al., 2015, 2016; González et al., 2017). To summarise, it is expected that the large scale environment impacts both the mass and spin assembly history of the host haloes. Observations and simulations have shown that the dynamical dependency on the large-scale environment remains true also for galaxy spin. The key question now is: does it also affect scalar quantities such as the mass, colour and sSFR of galaxies?

If the effect of the tidal field on accretion is not completely erased by baryonic physics and feedback, the physical properties which depend on the accretion history should therefore be impacted by the corresponding anisotropic infall, and a specific dependency of galaxy masses, colours and star-formation rates on the geometry of the environment should be observable. However very few observational evidences sustain the prediction that the geometry of the environment impacts significantly galaxy accretion history.

The “galactic conformity” (Weinmann et al., 2006) is a first measured evidence of this effect. It quantifies a correlation between the quenching of the central and the quenching of its satellite galaxies: the fraction of quiescent satellites is higher around quiescent centrals than around star forming centrals. This effect has been detected up to redshift 2.5 for both low and high mass satellite galaxies (Kawinwanichakij et al., 2016). It suggests that galaxy properties depend not only on halo mass but also on the larger scales, and is likely connected to the halo assembly history. Hearin et al. (2015) observed this galactic conformity up to very large separation (4 Mpc): this effect is difficult to explain in semi-analytical models where galaxy properties depend essentially on halo mass. Attempting to make sense of it, they use the term of “halo accretion conformity” and connect this effect to the mutual evolution of haloes in the same large-scale tidal field. The effect of these large scale tides are not captured by measurements of the local density.

More recently, results from spectroscopic surveys have reported the dependencies of galaxy masses and types on the geometry of the environment using SDSS (Yan et al., 2013; Martínez et al., 2016; Poudel et al., 2017; Chen et al., 2017), GAMA (Alpaslan et al., 2016) and, at higher redshift, VIPERS (Malavasi et al., 2016), based on the three-dimensional reconstruction of the cosmic web. In particular, Malavasi et al. (2016) find a significant trend for galaxies with different stellar masses and type to segregate near the filaments: the most massive and passive galaxies are found closer to the filament center than the star-forming, low mass ones. Can photometric redshift surveys unravel such signal at higher redshift?

The recently achieved accuracy of current photometric redshifts opens the prospect of working now with projected two-dimensional slices which are sufficiently thin to unambiguously study the projected filaments. One of the assets of photometric surveys is their ability to probe different epochs of cosmic evolution to leverage their relative importance in building up galaxies, with complete sample at much lower cost than spectroscopic surveys. The potential of photometric surveys may have been initially underestimated because of their insufficient precision along the line-of-sight, making it seemingly too difficult to reconstruct the three-dimensional field at the required scale to trace the cosmic web. With the large and deep multi-wavelength surveys such as the Cosmological Evolution Survey222http://cosmos.astro.caltech.edu/ (COSMOS; Scoville et al., 2007), the upcoming Hyper Suprime-Cam survey (HSC; Miyazaki et al., 2012b) and in the long term the Large Synoptic Survey Telescope (LSST; Ivezic et al., 2008), it is timely and interesting to estimate the requested accuracy to achieve this goal via simulations; it is also timely to exploit the latest COSMOS2015 catalogue (Laigle et al., 2016) as a prototype of these upcoming surveys to try and detect the dependencies of galactic properties such as stellar mass and galaxy colours with respect to their position within the cosmic web. Previous studies on COSMOS have investigated the effect of local environment on galaxy properties, either relying on isotropic density estimators (e.g. Scoville et al., 2013; Darvish et al., 2016a), or dividing the field in morphological components (field, filament or cluster; Darvish et al., 2016b). However, we stress that our aim in this work is to investigate the specific impact of the anisotropic tides on galaxy properties, therefore disentangling explicitly the effect of the geometry of the large-scale environment on driving these tides (i.e the traceless part of the tidal tensor) from effects purely driven by the local density (i.e. the trace of the tidal tensor). In order to achieve this goal, we will focus on mass and colour gradients specifically towards filaments.

The paper is organised as follows. Section 2 describes the observed and simulated datasets, and the tool to extract the density and the skeleton. Section 3 assesses the reliability of the skeleton reconstruction using photometric redshifts, while relying on the realistic virtual photometric mocks built from the Horizon-AGN  hydrodynamical cosmological simulation (Dubois et al., 2014). Section 4 presents results from the data and our interpretation. Section 5 wraps up and outlines future works. Appendix A, B, C investigate how the reconstruction and the signal vary as a function of the thickness of the slice, the photometric redshift accuracy and the completeness of the sample, hence providing guidelines for computing a similar reconstruction in surveys with different redshift uncertainties. Appendix D assesses the robustness of the gradients towards filaments when removing the contribution of the nodes. We use a standard $\Lambda$CDM cosmology with Hubble constant $H_{0}=70.4$ km$\cdot$s${}^{-1}\cdot{\rm Mpc}^{-1}$, total matter density $\Omega_{\rm m}=0.272$ and dark energy density $\Omega_{\Lambda}=0.728$. All magnitudes are expressed in the AB system. Throughout this work, uncertainties on the distributions are computed from bootstrap resampling. The shaded area around the distributions correspond to 1$\sigma$ deviation to the mean for a bootstrap resampling with 100 realisations.

2 Datasets and methods

Before analyzing the relative distribution of galaxies within the projected two-dimensional cosmic web, let us first introduce the COSMOS2015 photometric catalogue, the mock catalogue extracted from the Horizon-AGN  simulation, and the ridge tracing algorithm used to trace the cosmic web.

2.1 Photometric catalogue from COSMOS

Let us first summarise the main properties of the COSMOS2015 photometric catalogue that we built in Laigle et al. (2016). This catalogue includes apparent magnitudes in 30 bands from ultra-violet (UV) to infra-red (IR): 0.25–8 $\mu$m. In optical, it contains the same broadband images as previous releases (Capak et al., 2007; Ilbert et al., 2009). The COSMOS-20 survey provides most of the optical coverage: six broad bands: $B,V,r,i,z^{++}$, twelve medium bands and two narrow bands taken with Subaru Suprime-Cam (Taniguchi et al., 2007, 2015). The $u^{*}$-band data is obtained from the Canada-Hawaii-France Telescope (CFHT/MegaCam). Intermediate bands are very useful to increase the spectral resolution of the measured spectral energy distribution (SED) in the apparent optical wavelength range, contributing in this way to the accuracy of photometric redshifts especially at low redshift (e.g. Cardamone et al., 2010). In the near-infrared (NIR), we rely on the new $Y,J,H,K_{s}$ images from the UltraVISTA survey (DR2, McCracken et al., 2012) and the $Y$ band of the Hyper Suprime-Cam at Subaru telescope (Miyazaki et al., 2012a). The photometry is extracted using SExtractor (Bertin & Arnouts, 1996) in dual image mode. Following a similar reduction procedure as described in McCracken et al. (2012), the detection image is a $\chi^{2}$ sum of the four NIR images of UltraVISTA DR2 and the $z^{++}$ band taken with Subaru Suprime-Cam. This combination ensures to increase the completeness of the catalogue. Mid-IR data in the four IRAC channels (i.e. in a wavelength range between $\sim 3$ and $8\,\mu$m) come from the SPLASH program (PI: Capak).

Photometric redshifts ($z_{\mathrm{phot}}$) are computed using LePhare (Arnouts et al., 2002; Ilbert et al., 2006) with a configuration similar to Ilbert et al. (2013). In particular, the principal nebular emission lines are implemented using an empirical relation between the UV light and the emission line fluxes as described in Ilbert et al. (2009). From a comparison with the spectroscopic samples zCOSMOS bright (Lilly et al., 2007) and faint (Lilly et al., in prep.) the uncertainty of our photometric redshifts turns out to be $\sigma_{\Delta z}/(1+z)=0.007$ at $i^{+}<22.5$ and $\sigma_{\Delta z}/(1+z)=0.03$ at $1.5<z<4$ . Stellar masses, absolute magnitudes, and star formation rates (SFR) are also derived with the LePhare code.

We remove also from the catalogue all the objects which are flagged as belonging to a bad area or for which the photometry is possibly contaminated by the light of saturated stars, meaning that we keep only objects in $\cal{A}^{\rm UVISTA}$&$\cal{A}^{\rm!OPT}$&$\cal{A}^{\rm COSMOS}$, according to the notations in Table 7 of Laigle et al. (2016).

Passive galaxies are identified using their locations in the colour-colour plane $NUV$-$r$/$r$-$J$ (Williams et al., 2009; Ilbert et al., 2013). Passive objects are identified as those with $M_{NUV}-M_{r}>3(M_{r}-M_{J})+1$ and $M_{NUV}-M_{r}>3.1$ (Figure 16 in Laigle et al., 2016). In particular, this technique avoids mixing the red dusty galaxies and passive ones. This criterion is expected to be more robust than a segregation based on the SFR derived from the SED-fitting. In the following, we call galaxy colour-type the galaxy type (passive or star-forming) determined based on the colour diagram.

The stellar mass completeness of the sample is estimated from the $K_{\rm s}$ magnitude following Pozzetti et al. (2010). The stellar mass limits up to redshift 1.3 are shown in Table 1 for the entire population of galaxies. In the following, our results are based on a galaxy sample with $M_{\rm lim}>10^{10}{\rm M}_{\odot}$ and $0.5<z<0.9$. This choice is justified in Section 2.4. The final galaxy catalogue contains 11 284 objects.

2.2 Virtual catalogue from Horizon-AGN

The Horizon-AGN​ simulation333http://www.horizon-simulation.org/ (see Dubois et al., 2014, for more details) is run with a $\Lambda$CDM cosmology with total matter density $\Omega_{\rm m}=0.272$, dark energy density $\Omega_{\Lambda}=0.728$, amplitude of the matter power spectrum $\sigma_{8}=0.81$, baryon density $\Omega_{\rm b}=0.045$, Hubble constant $H_{0}=70.4\,\rm km\,s^{-1}\,Mpc^{-1}$, and $n_{s}=0.967$ compatible with the WMAP-7 data (Komatsu, 2011). The simulation is run with the ramses code (Teyssier, 2002), and the initially coarse $1024^{3}$ grid is adaptively refined down to $\Delta x=1$ proper kpc, with refinement triggered in a quasi-Lagrangian manner: if the number of dark matter particles becomes greater than 8, or the total baryonic mass reaches 8 times the initial dark matter mass resolution in a cell. The size of the simulation box is $L_{\rm box}=100\,h^{-1}\rm\,Mpc$ on a side, and the volume contains $1024^{3}$ dark matter particles, corresponding to a dark matter mass resolution of $M_{\rm DM,res}=8\times 10^{7}\,\rm M_{\odot}$.

Heating of the gas from a uniform UV background takes place after redshift $z_{\rm reion}=10$ following Haardt & Madau (1996). Gas can cool down to $10^{4}\,\rm K$ through H and He collisions with a contribution from metals using rates tabulated by Sutherland & Dopita (1993). Star formation occurs in regions of gas number density above $n_{0}=0.1\,\rm H\,cm^{-3}$ following a Schmidt law: $\dot{\rho}_{*}=\epsilon_{*}{\rho_{\rm g}/t_{\rm ff}}$, where $\dot{\rho}_{*}$ is the star formation rate mass density, $\rho_{\rm g}$ the gas mass density, $\epsilon_{*}=0.02$ the constant star formation efficiency, and $t_{\rm ff}$ the local free-fall time of the gas. The stellar mass resolution is $M_{{\rm*,res}}\simeq 2\times 10^{6}{\rm M}_{\odot}$. The simulation follows feedback from stellar winds, supernovae type Ia and type II, as well as the formation of black holes, which can grow by gas accretion at a Bondi-capped-at-Eddington rate and coalesce when they form a tight enough binary. Black holes release energy either in quasar or radio mode when the accretion rate is respectively above and below one per cent of Eddington, with efficiencies tuned to match the black hole-galaxy scaling relations at $z=0$ (see Dubois et al., 2012, for details).

To better mimic observed surveys, a lightcone is produced on-the-fly at every coarse time step of the simulation. The lightcone generation is described in Pichon et al. (2010). For the lightcone extraction, we have replaced gas cells by gas particles, and treated them as for the stars and dark matter particles. The area of the lightcone is 5 deg2 below $z=1$, and 1 deg2 above.

We identify galaxies using the AdaptaHOP structure finder (Aubert et al., 2004a; Tweed et al., 2009), applied to the distribution of star particles in the simulated lightcone. Because AdaptaHOP works at fixed redshift, the extraction is done in thin slices of the lightcone, each slice beeing kept at fixed redshift. To avoid cutting galaxies, slices are overlapping, then duplicated galaxies are removed. Structures are selected using a local threshold of 178 times the average matter density, the local density of individual particles being computed from the 20 nearest neighbours. Only structures that have more than 50 particles are considered. We compute galaxy fluxes and magnitudes using the single stellar population (SSP) models from Bruzual & Charlot (2003), using a Salpeter (1955) initial mass function. We assume that each star particle behaves as an SSP, and compute its contribution to the total SED by logarithmically interpolating the models in metallicity and age and multiplying by the initial mass of the particle. In this work attenuation by dust is not taken into account. The computation of the photometry includes redshift distortions due to the galaxy peculiar velocities.

A critical assessment of the Horizon-AGN  simulation in term of the evolution of its galaxy populations is presented in Kaviraj et al. (2016) through comparison of colours, luminosity functions and mass functions. When taking into account uncertainties and systematics, they found that the simulation captures reasonably well the evolutionary trend of galaxy mass assembly as seen through the statistical quantities. The remaining point of tensions concern the high redshift Universe (underestimation of the mass function at the high mass end) and the low mass galaxies (overproduction of low mass galaxies in the simulation). These small discrepancies in these specific ranges of mass and redshift will not affect the conclusions drawn in our work. The separation between star-forming and passive galaxies is also less straightforward in Horizon-AGN than in data: a large fraction of passive galaxies are indeed moved towards the star-forming region in the colour diagram due to residual star formation, which is not completely extinguished in passive systems. For this reason relying on the colour diagram itself is not sufficient to properly separate passive and star-forming galaxies, and an additional dynamical or morphological criterium could be required (e.g. the $V/\sigma$ parameter, see Dubois et al., 2016). In this work we will rely on the specific star-formation-rate (sSFR) to characterise galaxy types. The sSFR in the simulation is computed from the stellar mass formed over the last 100 Myr. The sSFR-mass relation is in broad agreement with the observations (Dubois et al., 2016), although the low-mass galaxies at low redshift are generally too passive, and the high-mass galaxies contain residual star formation.

We extract from the Horizon-AGN  simulated lightcone (see Figure 2 below) a mock catalogue which has the same characteristics in terms of photometry and photometric redshift accuracy as COSMOS2015​. In particular, we divide the simulated galaxies in bins of redshift and magnitude, and we perturb their redshift using a random gaussian error corresponding to the observed one in the same redshift and magnitude bins. This mock catalogue is used to test the two-dimensional skeleton reconstruction and our ability to recover mass gradients in two dimensions.

2.3 Density and cosmic network identification

The filaments of the cosmic web are identified from the density field, which refers either to the number density of dark matter particles (in the case of the three-dimensional reference skeleton in the simulation) or to the number density of galaxies (for the observed and simulated two-dimensional skeletons). When computed from the galaxy distribution, the density is estimated from the Delaunay tessellation of the particles. Note that the density is not weighted by galactic mass.

To identify the cosmic network from the density, we use the persistence based filament tracing algorithm (DisPerSE, Sousbie, 2011; Sousbie et al., 2011). This method identifies ridges from the density field as the special lines connecting topologically robust saddle points to peaks. The identified filament network is by construction multi-scale, while the extraction is robust to noise. The set of all segments defining these ridges is called the skeleton (Pogosyan et al., 2009). DisPerSE has been already successfully used in two dimensions to characterise stellar filaments in the Milky Way (e.g. Arzoumanian et al., 2011; Panopoulou et al., 2016). Each filament is defined to be a set of connected small segments linking extrema together. Persistence is defined as the ratio of the density value at the two critical points in a topologically significant pair of critical points: maximum-saddle, saddle-saddle (in three dimensions only) and saddle-minimum. This ratio quantifies the robustness of the underlying topological feature characterised by this pair (the fact that the connected critical points are responsible for the appearance of a new topological feature such as connected component, tube, ball, in the corresponding excursion set). Expressed in terms of numbers of $\sigma$, persistence quantifies the significance of the critical pairs in the Delaunay tessellation of a random discrete Poisson distribution. In the following, a skeleton with a $2\sigma$ persistence level is a skeleton for which all the persistence pairs which have their probability to be found in a random discrete Poisson distribution below $2\sigma$ of the mean have been removed. Removing low-persistence pairs is a multi-scale non-local method to filter noise/low significance filaments. This method is particularly well adapted to noisy datasets such as redshift catalogues (Sousbie, 2011). A complementary notion is robustness (Weinkauf & Gunther, 2009), which is a local measure of how contrasted the critical points and skeleton segments are with respect to their local background.

Table 2 summarises the extracted skeletons both on observations and simulations and their main features that we detail also below.

To identify filaments in the two-dimensional density field from the COSMOS2015 catalogue, we proceed as follows. We first divide the galaxy sample in slices of constant thickness, the choice of which is justified in Section 2.4. To estimate the density in each slice, we compute the two-dimensional Delaunay tessellation from the galaxy distribution. To take into account the errors on photometric redshifts in this calculation, we associate to each galaxy in each slice a weight $p_{{\rm gal},i}$ which is the probability for this galaxy to be in the slice $i$ given the probability function distribution of its redshift $P_{\rm gal}(z)$ (as computed by LePhare):

[TABLE]

where $z_{1,i}$ and $z_{2,i}$ are the lower and upper redshifts of the slice. This weighting allows to reduce the pollution effects of foreground and background galaxies for which the probability to be in the slices is low. To correctly estimate the density and the topology close to boundaries, a surface of “guard” particles is added outside the bounding box and new particles are added by interpolating the estimated density computed on the boundary of the distribution. Note that the Delaunay tessellation is by construction a robust method to reconstruct reliably the density in a field which contains masks. The triangulation automatically connects galaxies from across the masked regions and thus interpolates the density (e.g. Aragon-Calvo et al., 2015). Filaments are then identified with DisPerSE from the two-dimensional Delaunay tessellation with a 2$\sigma$ persistence threshold. The properties of the noise are not necessarily the same in two dimensions as in three dimensions. This implies that although a persistence threshold above 3$\sigma$ would be more conservative in three dimensions (which guarantees that less than 0.3% of the critical pairs are spurious), the two-dimensional skeleton extracted with a 3$\sigma$ persistence would miss a large number of robust three-dimensional filaments. To compensate this effect, choosing a two-dimensional persistence threshold of 2$\sigma$ enables to drastically decrease the number of unrecovered three-dimensional filaments (by a factor $\sim$1.7), without sensibly increasing the number of spurious two-dimensional filaments (see also Malavasi et al., 2016, for a motivation behind this mapping; this choice is further justified in Appendix A using the simulation). Figure 1 shows the density of the COSMOS field in four slices at redshifts $z\sim 0.52$, $z\sim 0.60$, $z\sim 0.67$ and $z\sim 0.73$ estimated from the Delaunay tessellation and their corresponding skeletons in red. The intensity of the red colormap reflects the robustness of the segments, the more robust segments being darker.

We extract the two-dimensional skeletons in the simulation with mock photometric redshifts following the exact method used for the COSMOS field. In addition, a two-dimensional skeleton is also extracted in each slice with true (“spectroscopic” in the following) redshifts (hence without including the weighting with the photometric uncertainties). Finally, a three-dimensional skeleton (skl${}_{\rm 3D}^{\rm DM}$) is extracted from the dark matter distribution. To compute this skeleton, in order to save computational time, dark matter particles are first projected on a three-dimensional grid with a resolution of 250 comoving kpc using a cloud-in-cell algorithm. We perform a first extraction with a very low persistence threshold. Then the adopted persistence threshold is estimated to be $\sim 5$ times the rms level of density fluctuations for the low density filaments in the previous extraction. Note that our results are robust against the exact value chosen for this persistence threshold. This skeleton becomes the reference skeleton. It is indeed the most suitable, as the galaxy distribution itself is already a biased tracer of the underlying density field. In addition, dark matter particles are more numerous in the simulation, allowing us to extract more accurate filament positions, which will be all the more important when computing distances between two- and three-dimensional skeletons in Section 3. A three-dimensional view of the lightcone and the corresponding skeleton is shown in Figure 2.

2.4 Choice of slices

The thickness of the slices on which the projected density is computed and the skeleton extracted will impact our analysis. A volume-limited sample with a constant stellar mass cut and a constant slice thickness is required, in order to avoid possible systematics in the measurements due to different quality of reconstruction (because of increasing slice thickness with increasing redshifts). In two dimensions, there can be some confusion between nodes and projected three-dimensional filaments, or filaments and projected three-dimensional walls. Furthermore, projected filaments may link some overdensities which are in fact not truly linked in three dimensions. Hence if the slices are too thick, the skeleton derived from the projected density will have little correspondence to its three-dimensional counterpart, with a large number of fake two-dimensional filaments and non-recovered three-dimensional filaments. However, if the chosen thickness of the slices is too small, the two-dimensional skeleton will trace only fragments of the three-dimensional filaments and will be in turn a poor representation of the fully connected three-dimensional skeleton. Following Appendix C, we choose to work up to redshift 0.9 with slices of thickness 75 comoving Mpc, and by considering all galaxies more massive than $10^{10}$ M⊙. This final choice is guided both by the required thickness for an optimal reconstruction when exact redshifts are available and by the current redshift accuracy of our sample. As it is shown in Table 1, with a thickness of 75 Mpc, all galaxies of the final sample have their redshift uncertainties smaller than the thickness of the slice, where uncertainties are estimated from the $2\times 1\sigma$ errors444The $1\sigma$ error is defined as the value enclosing half of 68% of the probability distribution function of the photometric redshift of the galaxy. given by LePhare for each galaxy.

Tests on the virtual photometric catalogue detailed below confirm that this choice is good enough to reconstruct the skeleton. We do not consider galaxies below redshift 0.5 as the comoving transverse width of the field starts to be too small to extract reliably the large-scale structure. The adjacent redshift slices are overlapping and spaced by half of the thickness of the slices. There are 30 slices on the redshift range $0.5<z<0.9$.

3 Extraction assessment in mocks

3.1 Two-dimensional versus three-dimensional skeleton

Let us now use the Horizon-AGN simulated catalogue to estimate the quality of the skeleton reconstructed with photometric redshifts. For this assessment, we compare the two-dimensional photometric skeleton with the projected three-dimensional reference one. Figure 3 displays slices of thickness 75 Mpc around redshifts $z\sim 0.52$, $z\sim 0.60$, $z\sim 0.67$ and $z\sim 0.73$ in the Horizon-AGN simulation. The orange skeleton corresponds to the projection of our reference skeleton skl${}_{\rm 3D}^{\rm DM}$ computed from the dark matter distribution. The segments are blue when they do not have a two-dimensional counterpart within 0.5$\times D_{\rm int}$ where $D_{\rm int}$ is the mean inter-galaxy distance in the slice. The black skeleton corresponds to skl${}_{\rm 2D}^{\rm spec}$ while the red skeleton is skl${}_{\rm 2D}^{\rm phot}$. The two-dimensional skeletons are visually in good agreement, hence confirming by eye that photometric redshifts are sufficient to reconstruct the two-dimension skeleton. In more details, we see that most filaments of skl${}_{\rm 2D}^{\rm phot}$ have a reliable counterpart in skl${}_{\rm 3D}^{\rm DM}$. However some filaments of the former are not recovered in the latter and some filaments of the latter have no counterpart in the former. Hence the choice of the persistence threshold is important to mitigate the number of filaments in skl${}_{\rm 2D}^{\rm phot}$ which have no counterpart in skl${}_{\rm 3D}^{\rm DM}$ and which are thus potentially noisy detection. In addition, filaments detected in skl${}_{\rm 2D}^{\rm phot}$ but not in skl${}_{\rm 3D}^{\rm DM}$ may potentially be walls seen in projection. We note that the two-dimensional segments which have generally no counterpart in skl${}_{\rm 3D}^{\rm DM}$ have a low robustness.

Let us quantify this statistically by measuring the distribution of distances, $d_{\rm 2D\rightarrow 3D}$, between the segments of skl${}_{\rm 2D}^{\rm phot}$ and skl${}_{\rm 2D}^{\rm spec}$ skeletons and their nearest segments in the projected skl${}_{\rm 3D}^{\rm DM}$ skeleton (following the method introduced by Sousbie, 2011). These distributions give us an estimation of the reliability of the skeleton. These distributions are shown on Figure 4. Green and blue lines correspond to the distributions of the distances $d_{\rm 2D\rightarrow 3D}$ using spectroscopic and photometric redshift respectively for the computation of skl2D. Although the photometric skeleton has a slightly larger tail, the curves computed with photometric and spectroscopic redshifts are quite similar, which suggests that the photometric redshift uncertainties have little impact on our global ability to recover the three-dimensional structures, at least under the chosen parameters. However, the tail of the distribution implies that a certain number of two-dimensional filaments have no three-dimensional counterpart at all, i.e. they correspond to spurious two-dimensional filaments connecting points in two dimensions that are spatially disconnected in three dimensions. We find that, on average, 21.3% and 15.2% of segments have no counterpart closer than 1.5 projected Mpc in the photometric and spectroscopic skeleton respectively. If we keep only the 60% more robust segments in the two-dimensional skeleton, these numbers fall by a factor $\sim$2. Weighting the distribution of $d_{\rm 2D\rightarrow 3D}$ by the robustness allows to considerably reduce this tail, as shown by the red and yellow curves for the photometric and spectroscopic skeletons respectively, implying that the fake filament detections correspond generally to low-robustness filaments.

Note that skl${}_{\rm 2D}^{\rm phot}$ is degraded or biased compared to the projected skl${}_{\rm 3D}^{\rm DM}$ skeleton for three reasons. First, the photometric skeleton skl${}_{\rm 2D}^{\rm phot}$ will be noisier than the spectroscopic one skl${}_{\rm 2D}^{\rm spec}$ due to the galaxy redshift uncertainties. Secondly, even with exact redshifts, skl${}_{\rm 2D}^{\rm spec}$ is a biased reconstruction of the projected three-dimensional skeleton computed from galaxies due to the projection effect. Finally, it should be acknowledged that the three-dimensional galaxy skeleton skl${}_{\rm 3D}^{\rm gal}$ is biased compared to the three-dimensional dark matter skeleton skl${}_{\rm 3D}^{\rm DM}$. This bias, although relatively negligible, is dependent on the limiting mass of the galaxy catalogue and is investigated in Appendix B.

3.2 3D measurement from 2D at COSMOS2015 accuracy

Although the three-dimensional filaments can be reliably traced from two-dimensional slices at the COSMOS2015 redshift accuracy, it is not obvious that a distance-to-filament trend will be statistically recovered in two dimensions. The projection will inevitably blur the measurement: galaxies far from the filament but on the same line-of-sight will be counted as very close by projection. We will therefore test in the Horizon-AGN simulation if a signal existing in three dimensions may be recovered in two dimensions. Following the recent results from Malavasi et al. (2016), we expect to find mass gradients towards three-dimensional filaments. We first measure this trend in three dimensions. In a second step, we verify that the signal is also measurable in two dimensions.

Figure 5 presents the distributions of the three-dimensional distances to skl${}_{\rm 3D}^{\rm DM}$ for simulated galaxies in four different bins of masses in the Horizon-AGN lightcone. The uppest mass bin is chosen to contain galaxies above the knee of the observed galaxy mass function, which is found to be $\sim 10^{10.8}$M*⊙* in $0.5<z<0.8$ from the fit of the mass function with a double Schechter profile (Davidzon et al., 2017). The most massive galaxies are found closer to the filament center than the low mass ones, in agreement with Malavasi et al. (2016).

The same measurements are then carried in two dimensions with skl${}_{\rm 2D}^{\rm phot}$ on the Horizon-AGN  lightcone, in all slices between $0.5\!<z<\!0.9$. The projected transverse distances to the filament, measured in degrees, are then converted in comoving Mpc. The cumulative distributions in all slices are stacked and the mean distribution is presented in the left panel of Figure 6. The galaxy mass segregation is also found in two dimensions with photometric redshifts demonstrating our ability to recover the signal in observed data.

Let us now stress that we are looking for environmental effects driven specifically by filaments. As galaxies are more massive and passive in clusters than in the field, we expect mass gradients towards nodes driven by the effect of cluster alone. These gradients could contaminate the mass gradient-trend towards filaments. To minimise node contributions, we remove from the analysis galaxies closer to a node than 3.5 and 0.8 projected Mpc (see Appendix D) in three and two dimensions respectively and we reshuffle galaxy masses with respect to their position along the filaments. We show in Appendix D that this allows to correctly remove their contribution to the gradients towards filaments. We still recover mass gradients towards filaments both in three and two dimensions: more massive galaxies are more confined to filaments than less massive ones. These gradients are driven by filaments alone.

However, it is well known that denser environments contain more massive, redder and much less star-forming galaxies (e.g. Dressler, 1980; Balogh et al., 1997; Blanton et al., 2003; Kauffmann et al., 2004; Davidzon et al., 2016; Cucciati et al., 2016), even if we exclude cluster regions. Since the density is anisotropically distributed in filaments, such local mass-density relation could naturally lead to mass gradients towards filaments. To disentangle the contribution of the anisotropic large-scale environment from the effect of the local mass-density relation, we proceed following Malavasi et al. (2016). The galaxy catalogue is split in local density bins (40 density bins logarithmically distributed), where the density is computed from the three-dimensional galaxy Delaunay tessellation smoothed with a Gaussian kernel of 3 Mpc and 1 projected Mpc in three and two dimensions respectively (we checked that varying this smoothing length does not significantly change the conclusions). Then stellar masses are reshuffled between the galaxies in each bin. The distributions after reshuffling are shown as dashed lines in Figures 5 and 6. Because it preserves the mass-density relation, this reshuffling partially preserves the mass gradients towards filaments, however the reshuffled and original distributions are not identical. The median values of both original and reshuffled distributions are presented in Table 3. The gradient signal is much larger with the original distribution than the reshuffled one. In two dimensions, the median of the lowest mass bin distribution is higher than that of the other bins, but there is otherwise no significant gradient.

Let us summarise these results. First, mass gradients are found towards filaments in the Horizon-AGN simulation. These gradients are specifically induced by filaments and not driven by nodes. In addition, they cannot entirely be explained by the mass-density relation alone. The particular geometry of the large-scale environment itself is what drives the gradient. Secondly, we demonstrate that at the accuracy of the COSMOS2015 catalogue, mass gradients found in the three-dimensional mocks can also be identified in the two-dimensional ones. The method can therefore robustly be applied to the observed COSMOS2015 catalogue.

4 Galaxy properties in filaments

Having confirmed the feasibility to trace the filaments in two-dimensional slices with a thickness of $75$ Mpc and to recover the corresponding three-dimensional signal, let us now study in more details the distribution of galaxy masses and colour-types as a function of their distance to filaments using both the observed COSMOS2015 and simulated Horizon-AGN catalogues.

4.1 Galaxy mass gradients towards filaments

The right panel of Figure 6 presents the distributions of the distance to the filaments for the four bins of mass in the Horizon-AGN  lightcone (after including photometric redshift errors) and in COSMOS2015 respectively. The signal measured in all slices in the redshift range $0.5<z<0.9$ has been stacked. We remove the contribution of the nodes following the above-described prescription. Note that the lowest mass bin $9.5<\log M_{*}/{\rm M}_{\odot}<10.0$ should be considered with caution. These galaxies are not used to compute the skeleton, and their redshift uncertainties (2$\times$ 1$\sigma$) are generally higher than the thickness of the slices. In COSMOS, we also find that that massive galaxies are closer to the center of filaments compared to less massive ones, in agreement with the measurement in Horizon-AGN. We note that the galaxy mass gradients signal are stronger in the simulation in two dimensions than in COSMOS2015. A possible explanation is that uncertainties on galaxy masses are not modelled in the simulation, and not taken into account in the measurement from COSMOS2015. These uncertainties are likely to dilute the gradients. As low mass galaxies are more numerous and with larger uncertainties than high mass galaxies, they are likely to pollute the higher mass bins. This bias should dilute the signal especially in the highest mass bins, as we indeed observe in Figure 6. Node contribution is removed and we proceed with the above described reshuffling to disentangle the contribution of the anisotropic large-scale environment from the effect of the local mass-density relation. The distributions after reshuffling are shown as dashed lines in Figure 6. Consistently to the Horizon-AGN result, we find that the reshuffling almost completely dilute the gradient signal: the values of the median of each bin are consistent within the errobars, except the one of the lowest mass bin which is slightly higher (see Table 3). Varying the number of density bins used for the reshuffling does not change the trend. We conclude that, after removing the significant contribution of the nodes, most of the gradient signal towards filament is driven by the geometry of the large-scale environment itself.

4.2 Colour-type segregation towards filaments

Let us now investigate the effect of the anisotropic environment on galaxy star formation. Our COSMOS2015 galaxy sample has been divided in passive and star-forming populations based on the colour diagram, as explained in Section 2.1. Figure 7 presents the distributions of the distances to the filaments for star-forming and passive galaxies for different mass bins in COSMOS2015. All galaxies with $M_{*}>10^{10}~{}{\rm M}_{\odot}$ are considered. Once again, we remove the contribution of the nodes as explained previously. At fixed mass in all mass bins we find that passive galaxies are statistically closer to the filament center than star-forming ones. As we minimised the contribution of galaxies in nodes, this effect is specifically related to filaments. However, in accordance with the above mentioned mass-density relation, it has long been established that at low redshift galaxies are found to be much less star-forming in high-density regions relative to low-density regions. As the density is anisotropically distributed in filaments, the density itself could explain these colour-driven gradients towards filaments. To disentangle the effect of the geometry of the large-scale environment from the galaxy type-density relation, we proceed as previously. In each slice, the galaxy catalogue is split in local density bins (15 density bins logarithmically spaced) and stellar mass bins (4 stellar mass bins logarithmically spaced). Then galaxy types are randomised within these bins preserving the fraction of passive and star-forming galaxies in each bin. After reshuffling, no significant gradients are found whatever the considered mass bins. To better quantify the significance of the result, we look at the residual between the passive and the star-forming distributions ($D_{1}$ and $D_{2}$), expressed in number of $\sigma$: $\Delta_{1-2}=\Delta(D_{1}-D_{2})/\sqrt{(}\sigma_{1}^{2}+\sigma_{2}^{2})$, where the subscripts 1 and 2 refer to the passive and star-forming distributions respectively. We find that, for the original sample, the deviation generally exceeds $1\sigma$, and $2\sigma$ for the intermediate mass bin ($10^{10.4}\!<\!{M_{*}/\rm M}_{\odot}\!<\!10^{10.8}$) and the whole galaxy population. However the reshuffling generally completely cancels the signal. We conclude that the colour segregation within the filaments is not driven by the local density but by the geometry of the environment itself. Finally, we want to assess that these gradients are not dominated by the different mass distributions of the star-forming and passive populations even within a given mass bin. Therefore we reproduce the same measurements by selecting one galaxy sub-sample in each population with the same mass distribution. The results remain unchanged for the original sample while the residual gradients for the reshuffled one disappear.

The same measurement cannot be reproduced exactly in Horizon-AGN, since the separation between star-forming and passive galaxies is less straightforward in the simulation than in data. Nevertheless, we can compute the distribution of the distance to the closest filaments for galaxies in different sSFR bins. The result for all galaxies more massive than $10^{10}$ M*⊙* are presented in Figure 8. The measurement is computed in three dimensions, and the skeleton used is the reference dark matter skeleton skl${}_{\rm 3D}^{\rm DM}$. Three sSFR bins are chosen, such as to isolate an undoubtedly passive fraction (log sSFR/yr${}^{-1}<-$10.8) and a very star-forming ones ($-9.8<{\rm log\,sSFR/yr}^{-1}$). The intermediate bin corresponds to the remaining galaxies ($-10.8<{\rm log\,sSFR/yr}^{-1}<-$9.8). When considering together all the galaxies more massive than $10^{10}$ M*⊙*, we recover in the simulation the same trend as in COSMOS2015: passive galaxies are closer to the filament center than star-forming ones. The gradient signal is not completely cancelled by the reshuffling, but it is stronger for the original sample than for the reshuffled one. The results for individual bins of mass are presented in Table 5. Interestingly, depending on the mass bin, gradients are not only found between passive and star-forming populations, but also within the star-forming population. Thus there are significant gradients between the intermediate and star-forming populations for $10.0<\log M_{*}/{\rm M}_{\odot}<10.4$, while the local density alone does not drive such gradients. We also note that we do not find significant gradients for $10.4<\log M_{*}/{\rm M}_{\odot}<10.8$, although that mass bin displays the strongest signal in COSMOS2015. Finally, for the highest mass bin $10.8<\log M_{*}/{\rm M}_{\odot}$ a strong signal is also observed in the original sample, but the reshuffled one is consistent with no gradient. These measures are comforting the results from COSMOS2015 and highlight once again the role played by the geometry of the anisotropic environment in driving galaxy properties. Although not shown here, we find that, in the simulation, the strength of the signal is dependent on the choice of the sSFR bins as well as the robustness of the filaments included in the analysis (which indicates that the signal is probably multi-scale). A more detailed analysis of the galaxy property gradients in the simulation taking into account these variations will be the topic of future work.

4.3 Interpretation and discussion

Our results are quantitatively summarised in Tables 3, 4 and 5 which report the median values of each distribution. The normalisation of these values are quite sensitive to cosmic variance and to the chosen thickness of the slices for the two-dimensional reconstruction. Hence they have to be taken cautiously when comparing to other datasets. However the significant differences between the medians of the different populations, either divided in mass bins, sSFR bins or in terms of colour-types, is a robust result which allows us to report mass gradients towards filaments in COSMOS2015 and Horizon-AGN, and colour-type gradients towards filaments in COSMOS2015. In Horizon-AGN, we also find sSFR gradients in agreement with COSMOS2015. They corroborate results found in other fields by independent methods of analysis (Alpaslan et al., 2016; Malavasi et al., 2016; Poudel et al., 2017; Chen et al., 2017).

Theoretical considerations, together with measurements both in observations and simulations give us clues to how galaxy properties depend on the geometry of their environment. Galaxy harassment (Moore et al., 1996), ram-pressure stripping of the gas (Gunn & Gott, 1972; Quilis et al., 2000), strangulation of low-mass satellite galaxies (Larson et al., 1980; Balogh et al., 2000) or cumulative galaxy-galaxy hydrodynamic/gravitational interactions (Park & Hwang, 2009) may explain the depletion or heating of the cold gas in quenched galaxies in group or cluster environment, i.e. mainly at the node of the cosmic web. At intermediate density ranges, galaxy mergers, while temporarily modestly increasing star formation (e.g. Jogee et al., 2009; Kaviraj et al., 2015), may also halt star-formation and turn star-forming galaxies into passive ones (e.g. Dubois et al., 2016). But none of these processes can be responsible for galaxy segregation depending only on distance to filaments taking place in regions of low to intermediate densities, where galaxy interactions are limited.

Since the signal found in this work is sensitive to the large-scale geometry of the matter distribution itself, it must be driven by the large scale (traceless part of the) tidal tensor, which quantifies this anisotropy. The physical properties of dark matter halos (Codis et al., 2012) and their host galaxies (Dubois et al., 2014; Welker et al., 2014) indeed correlate with the dynamics of their embedding anisotropic large scale structure. This has identified kinematic signatures: angular momentum advection (Codis et al., 2015b) and internal kinetic anisotropy (Faltenbacher & White, 2009) are also set by the larger scale tides of the environment. The induced tides not only impact the merger tree and accretion history of the host, but also the filamentary flow of cold gas connecting to the host, hence its coherent gas supply. The past and present efficiency of star formation, as traced by the mass, colours and specific star formation rate, depends critically on the infalling rate and impact parameter of cold gas. It does not only depend on the timing of infall (in connection with the accretion history of the host), but also on its geometry and physical content, in connection with the dynamics of the matter flows (e.g. in plane gas-rich and co-rotating, or dry mergers with random orbital parameters). As discussed in the introduction, the vorticity-rich large-scale filaments are indeed the locus where low-mass galaxies steadily grow in mass via quasi-polar cold gas accretion (Welker et al., 2015), with their angular momentum aligned with the host filament. Galaxies are expected to accrete more efficiently cold gas when their angular momentum is aligned with the preferential direction of the gas infall, i.e. aligned with the filament (Pichon et al., 2011). Hence one expects star formation efficiency to be strongest wherever the alignment is tightest. Maximal alignment occurs in the highest vorticity regions (Laigle et al., 2015), i.e at the edge of filaments (whereas haloes terminate their mass assembly in the filament core while drifting towards nodes following the large-scale flow). The locus of this induced excess of star formation should therefore have measurable signatures in observations when quantified in the metric of the filament, as discussed e.g. in Codis et al. (2015b). This is fully consistent with what is observed in Figure 7 and 8.

In the present study we have attempted to disentangle the specific signature of this anisotropy while relying on the physical properties of galaxies at fixed mass. The various tests presented in this paper confirm that the anisotropy of the environment traced by transverse gradients with respect to filaments is detected in COSMOS2015 catalogue, as anticipated and quantified from mocks derived from the light-cone of Horizon-AGN. We provide observational evidence that what has been a pillar of galaxy formation, i.e. the assumption that the average effects of the environment can be approximated by the effect of the (in particular azimuthally) average environment is only valid at first order. This assumption is obviously strictly true for linear dynamic and linear observables, which does not capture the full complexity of what builds up the properties of galaxies. This past simplification was partially driven by the community’s strong focus on two-point functions, which cannot capture the anisotropy of the cosmic web. Yet it had long been known that for instance angular momentum – which is undoubtedly built from anisotropic tides – is a key underlying property driving morphology, which correlates strongly with colour and star formation efficiency. Hence one could attempt to extend, e.g. assembly bias theory (Sheth & Tormen, 2004) to try and explain the observed diversity at a given mass, so long as it is understood that mass assembly is in fact also driven by anisotropic large scale tides, which will impact gas inflow towards galaxies, hence their properties (Yan et al., 2013; Tramonte et al., 2017).

More elaborate – but still restricted in scope – bias theories such as excursion set peak theories involving a moving barrier account for the other invariants of the tidal tensor (see Desjacques et al., 2016, and reference therein), but fail to incorporate the effect of the orientation of its eigen-vectors. Semi-analytical models (e.g. Lacey & Silk, 1991) do account for the timings encoded in merging trees measured in simulations– which reflect in part the cosmic web – but are restricted to following dark halo only and typically leave out the configuration of their mergers.

An improved model for galaxy properties should therefore explicitly integrate the diversity of the topology of the large environment on multiple scales (following, e.g. Hanami, 2001) and quantify the impact of its anisotropy on galactic mass assembly history, and more generally on the kinematic history of galaxies (e.g. extending Codis et al., 2015b). What also remains to be understood is what process dominates at what redshift: does the offset of merger and accretion rate imposed by the large scale turbulent flow explain most of the environment dependence in observed physical properties (Aragon-Calvo et al., 2016), or is the geometry of gas inflow within filaments prevalent in feeding galactic discs coherently (Pichon et al., 2011; Stewart et al., 2011)? The kinematic of the large scale flow is neither strictly coherent nor fully turbulent. Both processes will in turn modify the efficiency of star formation and feedback, which is critically non-linear, currently poorly understood and fine tuned. In order to disentangle these effects, the next step will involve studying the connection between galaxy properties, angular momentum and distance to filaments in the Horizon-AGN simulation and in connection with the COSMOS2015 catalogue.

5 Conclusion

Using the photometric catalogue COSMOS2015, we investigated how galaxy properties depend on their anisotropic environment. A realistic photometric catalogue extracted from the Horizon-AGN hydrodynamical simulation was used to quantify our ability to recover the filamentary structures from two-dimensional slices at the same photometric redshift precision, and to predict the expected gradients. The main findings of this work are the following:

• •

Photometric surveys with a photometric precision similar to COSMOS2015 offer the prospect to reliably study the 3D properties of cosmic web from projected two-dimensional slices;

• •

We observe both in Horizon-AGN and in COSMOS2015 a significant tendency for massive galaxies to be closer to the filament center than for less massive galaxies;

• •

We also observe segregation of passive and star-forming galaxies near the filaments both in COSMOS2015 and Horizon-AGN. At fixed mass, passive galaxies are more confined in the core of the filament than star-forming ones. These two signals persist when minimising the contribution of nodes and cannot be explained by the local mass-density relation alone.

These findings are qualitatively consistent between simulation and data, and in agreement with studies relying on spectroscopic redshifts (e.g. Malavasi et al., 2016). They underline the specific role played by the anisotropy of the large scale cosmic web in shaping galaxies properties. The geometry of the large scale environment drives tides, which impacts the merger and accretion history of galaxies, their spin and in turn their observed properties. Mass and density alone are not sufficient to characterise the observed properties of galaxies aways from the cosmic web.

More observations are needed to disentangle the effects of all these competing processes, while relying on a wider fields and (photometric) redshift baseline, towards the peak of the cosmic star formation, where cold streams are expected to be more efficient. This will be challenging as one will need to carefully disentangle the evolutionary effects and selection biases. An alternative observational strategy is to rely on ongoing (SAMI and MANGA Croom et al., 2012; Bundy et al., 2015) or upcoming (Hector, Bland-Hawthorn, 2015) integral field spectroscopy surveys and explore spin properties of galaxies (e.g. Harrison et al., 2017) such as their physical and kinematic axis orientation as a function of the large scale environment. Simulated end-to-end catalogues such as those produced here will remain key to qualify these investigations.

Acknowledgments

CL was supported by the ILP LABEX (under reference ANR-10-LABX-63 and ANR-11-IDEX-0004-02) during the first part of this work and is now supported by a Beecroft Fellowship. JD and AS acknowledges funding support from Adrian Beecroft, the Oxford Martin School and the STFC. OI acknowledges the funding of the French Agence Nationale de la Recherche for the project “SAGACE”. This work relied on the HPC resources of CINES (Jade) under the allocation 2013047012 and c2014047012 made by GENCI and on the Horizon Cluster hosted by Institut d’Astrophysique de Paris. We warmly thank S. Rouberol for running the cluster on which the simulation was post-processed. This research is part of Spin(e) (ANR-13-BS05-0005, http://cosmicorigin.org) and horizon-UK. This research is also partly supported by the Centre National d’Etudes Spatiales (CNES). This work is based on data products from observations made with ESO Telescopes at the La Silla Paranal Observatory under ESO programme ID 179.A-2005 and on data products produced by TERAPIX and the Cambridge Astronomy Survey Unit on behalf of the UltraVISTA consortium. We thank D. Munro for freely distributing his Yorick programming language and opengl interface (available at http://yorick.sourceforge.net/).

Appendix A Choice of persistence in 2D

As stated in Section 2.3, the extraction of skeletons with the same persistence thresholds $\sigma$ in two and three dimensions will typically not lead to the detection of the same filaments. Indeed, one must significantly decrease $\sigma$ with respect to its three-dimensional value to recover in two dimensions most of the three-dimensional filaments. Let us estimate the optimal two-dimensional $\sigma$ which recovers most of the filaments without introducing spurious detections. This optimal $\sigma$ is estimated from the simulated catalogue of Horizon-AGN, by comparing the number of unmatched segments between skl${}_{\rm 2D}^{\rm spec}$ and skl${}_{\rm 3D}^{\rm DM}$. We compute first the distances between the two skeletons. In order to compare equitably the percentage of unmatched segments for different slice thicknesses, the distances between the two skeletons are divided by the mean two-dimensional inter-galaxy distance $D_{\rm int}$ (which increases with the slice thickness), taken as the square root of the area divided by the number of galaxies in the slice. Then we define the percentage of unmatched segments as the fraction of segments such as their closest counterpart is further than $D_{\rm int}$ (changing this threshold would move Figure 9 vertically). The percentage of unmatched segments is shown in Figure 9. In order to reduce as much as possible the number of skl${}_{\rm 3D}^{\rm DM}$ segments which do not have a counterpart in skl${}_{\rm 2D}^{\rm spec}$, without sensitively increasing spurious filament extraction in skl${}_{\rm 2D}^{\rm spec}$, we find that a 2$\sigma$ persistent threshold for skl${}_{\rm 2D}^{\rm spec}$ offers a reasonable balance between the percentage of unmatched segments in skl${}_{\rm 2D}^{\rm spec}$ and skl${}_{\rm 3D}^{\rm DM}$, whatever the slice thickness. Although skl${}_{\rm 2D}^{\rm spec}$ is computed for Figure 9 with a mass limit of $10^{10}{\rm M}_{\odot}$, this result does not significantly change with a different mass limit.

Appendix B Galaxy versus DM skeletons

The skeleton computed from dark matter particles is the most suitable to stand as the reference skeleton, since the galaxy distribution is a biased tracer of the underlying density field. In this Appendix, we compare the three-dimensional galaxy and dark matter skeletons and we highlight how their relative difference depends on the mass limit of the galaxy sample. To get values directly comparable to Figure 4, we measure the distances between the galaxy and dark matter skeletons in projection, although both skeletons are computed in three dimensions.

The distribution of the distance between the dark matter and galaxy skeletons skl${}_{\rm 3D}^{\rm gal}$ and skl${}_{\rm 3D}^{\rm DM}$ are plotted in Figure 10 for different mass limits. The solid lines refer to the distance from the galaxy skeleton towards the dark matter one, while the dashed lines refer to converse. All solid lines are almost overlapping. These distribution are peaked at very close separation, but still display a long tail. This suggests that whatever the chosen mass limits for the galaxy skeleton reconstruction, the galaxy filaments have generally a very close counterpart in the dark matter skeleton, although a small fraction has no counterpart at all. We note a slight improvement in recovering the DM skeleton when working with increasing masses. However, when decreasing the mass limit of the galaxy sample, the number of unmatched filaments in the dark matter skeleton decreases strongly (the tail of the dashed distribution decreases). To summarise, working only with massive galaxies allows us to recover the main dark matter filaments quite accurately, but at the price of missing a non negligible fraction of dark matter filaments. A reasonable balance involves choosing a mass limit between 109.5 and 10${}^{10}{\rm M}_{\odot}$. Finally we note that these variations are relatively small; overall the width of these distributions are smaller than those displayed in Figure 4.

Appendix C Choice of slice thickness

We aim now to determine if there is an optimal slice thickness to extract the cosmic web and estimate transverse gradients. We proceed in two steps. First, given an ideal survey (without photometric redshift uncertainties), we estimate how the accuracy of the reconstruction varies with the thickness of the slices. Secondly, we estimate how the quality of the reconstruction varies as a function of mass limits given a photometric redshift sample. Given the redshift accuracy of the COSMOS2015 catalogue,we choose then the optimal slice thickness and mass limit.

C.1 Optimisation of slice thickness with exact redshifts

The optimal slice thickness should lead to the best match between the two-dimensional skeleton and its projected three-dimensional counterpart. The agreement is quantified from the distribution of the distances between the projected skl${}_{\rm 3D}^{\rm DM}$ skeleton and skl2D. From the distribution, we measure the median and the width which encompasses 68% of the distribution around the median. Good agreement between the skeletons requires both the median and the width to be small. Let us first test the single effect of the projection for different thicknesses, independently of the redshift accuracy.

The distances are normalised by the mean inter-particle distance $D_{\rm int}$. Figure 12 shows three slices of skl${}_{\rm 2D}^{\rm spec}$ and the projected skl${}_{\rm 3D}^{\rm DM}$, corresponding to three slice thickness: 30 Mpc, 75 Mpc and 170 Mpc around redshift $z\sim 0.59$. In overall the more robust two-dimensional filaments trace existing three-dimensional filaments. For the lowest thickness slice, most of the three-dimensional filaments appear to have a counterpart in the two-dimensional skeleton, but there are a lot of fake filaments, generally with a low robustness because the skeleton traces only fragments of filaments. For the highest thickness slice, this fraction of fake two-dimensional filaments decreases strongly, but much more three-dimensional filaments have no counterpart. We statistically quantified it in Figure 11, which shows the distribution of the distances between skl${}_{\rm 2D}^{\rm spec}$ (M${}_{\rm lim}=$1010M*⊙) and skl${}_{\rm 3D}^{\rm DM}$, for slices of different thicknesses ranging from 30 Mpc to 200 Mpc. It therefore appears that for exact redshifts, there is no optimal slice thickness within the considered range. Working with thin slices ensures accurate filaments identification, but a non negligible fraction of two dimensional filaments have no counterpart at all in the three dimensional skeleton. This fraction diminish significantly when increasing the slice thickness, but at the price of less accurate filament identifications. In addition, as displayed by the black line in Figure 13, the fraction of unmatched three-dimensional segments increases significantly with the thickness. Hence a reasonable choice could lie between 60 and 100 Mpc. But our final choice of the slice thickness will be in any case limited by the redshift accuracy of the sample. Finally, Figure 12 indicates that filtering low robustness segments would be a way to remove the filaments which have probably no counterpart in three dimensions. This could be of interest when working with thicker slices. However, as the robustness is an estimation of the density of the filament with respect to the local background, keeping only the more robust filament would bias the measure towards the denser filaments. In addition, keeping only segments above a certain robustness threshold would destroy the connectivity of the cosmic web by removing entire sections of filaments. It would bias statistical measurements on the cosmic web, hence prevent accurate cosmological measurements based on such statistics. Note that although only the results for a mass limit of 1010M⊙* are shown in this Section, the trend is preserved when the mass limit of the sample varies.

C.2 Robustness w.r.t. photometric redshift accuracy

We found (Figure 10) that decreasing the mass limit of the sample is desirable to bring a better agreement between the galaxy and dark matter skeletons. At fixed mass limits, the optimal choice of the slice thickness broadly stands between 60 and 100 Mpc (Figure 13), although increasing the slice thickness is possible at price of missing a larger number of three-dimensional filaments. Let us now consider how the accuracy of the reconstruction is affected by the redshift uncertainties of the sample. This will drive our final choice both for the mass limit and the slice thickness, knowing that lower mass galaxies have higher redshift uncertainties. Figure 13 shows the percentage of unmatched segments within $D_{\rm int}$ between skl${}_{\rm 3D}^{\rm DM}$ and skl${}_{\rm 2D}^{\rm phot}$ as a function of the slice thicknesses for different mass limits. For comparison, the black line displays the distribution for skl${}_{\rm 2D}^{\rm spec}$ computed with a mass limit of 10${}^{10}{\rm M}_{\odot}$. When the mass limit decreases, more galaxies with high redshift uncertainties are included in the computation. It appears therefore that it is not reasonable to use a mass limit of 109.5M*⊙, as the fraction of unmatched two-dimensional filaments is very consequent (always greater than 20% for thickness lower than 100 Mpc). Increasing the mass limit to 1010M⊙* brings this fraction below 10%, while the fraction of unmatched three-dimensional segments does not increase significantly. Given the redshift uncertainties of COSMOS2015, we finally choose to work with galaxies more massive than 1010M*⊙*, which guarantees both that the galaxy skeleton is not strongly biased compared to the dark matter one, and that the fraction of unmatched two-dimensional segments is not too high. Within the allowing thickness range, we pick the value of 75 Mpc: it ensures that all the galaxies in the sample have their redshift uncertainties (2$\times$1$\sigma$) lower than the thickness of the slice.

Appendix D Filaments or nodes gradients?

Let us finally assess the robustness of the gradients towards filaments after removing the nodes. Nodes do indeed also drive mass gradients. These gradients could be either isotropic around the overdensity, or anisotropic along the filament, or a combination of both. In the first case, gradients towards nodes could mimic gradients towards filaments, even if they play no role in driving them.

To probe independently gradients towards nodes and gradients towards filaments, we will show in the following that we are able to cancel gradients towards nodes without canceling gradients towards filaments, and vice-versa. For each galaxy, the distance to the closest node is defined along the filament (i.e. the distance to the closest node is taken as that of its closest segment in the skeleton555The distance to the node is chosen to be that along the filament instead of the (euclidian) shortest distance to the node, because pure gradients towards filaments could otherwise also mimic gradients towards nodes: galaxies closer to filament would be also on average geometrically closer to the nodes. This cannot be the case if the distance to the node is taken as the one computed along the filament.).

As stated in the main text, in order to minimise the effect of the nodes in driving gradients, we first exclude from the sample galaxies which are closer to the nodes than a given distance $d_{\rm min}$. This distance has to be large enough to effectively mimimise the effect of the nodes, but small enough in order to keep a significant number of galaxies in the sample. We find that $d_{\rm min}=3.5$ Mpc in three dimensions and $d_{\rm min}=0.8$ (projected) Mpc in two dimensions offer a good compromise. The left panel of Figure 14 shows the mass gradients towards nodes in three dimensions after excluding the regions around nodes. Even after excluding a large region around nodes, significant gradients are found in the direction of the nodes.

We first try to cancel these gradients by preserving the gradients towards filaments. We therefore reshuffle galaxy masses with respect to their position along the filament (in given bins of distance to the filament). Hence the gradients towards filaments are preserved by construction (Figure 5). However, we find that, after reshuffling, no significant gradients towards nodes are measured (middle panel of Figure 14): we are able to cancel gradient towards nodes without cancelling gradients towards filaments. In a second step we try to cancel filament gradients by preserving the gradients towards nodes. We therefore reshuffle galaxy masses with respect to their distance to the filaments (in given bins of distance to the nodes). Hence the gradients towards nodes are preserved by construction. However, we find that, after reshuffling, no significant gradients towards filaments are measured (right panel of Figure 14): we are able to cancel gradient towards filaments without cancelling gradients towards nodes. This would not have been possible without excluding galaxies within $d_{\rm min}$ around the nodes. From these tests we conclude that gradients towards nodes and gradients towards filaments outside these exclusion regions are independent.

In two dimensions, it is not straightforward a priori that the nodes extracted from the two-dimensional skeleton are reliable tracers of their three dimensional counterparts. In the observed dataset we have only access to the two-dimensional nodes, but in the simulation, we also have access to the projected three-dimensional nodes. We use first the mock catalogue from the simulation to exclude regions around two-dimensional nodes and to reshuffle the position of the galaxies along the filaments as it would be possible to do with the observed catalogue. We compute then the distribution of galaxy distances by bin of mass around the projected three-dimensional nodes, in order to check that this method effectively cancel the signal. Figure 15 shows the mass gradients towards projected three-dimensional nodes, after excluding the regions around two-dimensional nodes, before and after the reshuffling along the two-dimensional filaments. After reshuffling, no significant gradients towards three-dimensional nodes are found, while the gradients towards filaments are preserved by construction (Figure 6). Similarly to the three-dimensional case, we also check that we can suppress gradient towards filaments but preserving those towards nodes. We therefore reshuffle galaxy masses with respect to their distance to the filaments in given bins of distance to the nodes. This operation partially preserves gradients towards three-dimensional nodes (as two-dimensional and projected three-dimensional nodes are not completely identical). However, it completely cancels gradients towards filaments. We conclude that gradients towards projected three-dimensional nodes and gradients towards two-dimensional filaments are independent.