Probing modified gravity in cosmic filaments

Multiple modifications of general relativity (GR) have been proposed in the literature in order to understand the nature of the accelerated expansion of the Universe. However, thus far all the predictions of GR have been confirmed with constantly increasing accuracy. In this work, we study the imprints of a particular class of models --"screened"modified gravity theories -- on cosmic filaments. We have utilized the $N$-body code ISIS/RAMSES to simulate the symmetron model and the Hu-Sawicky $f(R)$ model, and we post-process the output with DisPerSE to identify the filaments of the cosmic web. We investigated how the global properties of the filaments -- such as their lengths, masses, and thicknesses -- as well as their radial density and speed profiles change under different gravity theories. We find that filaments are, on average, shorter and denser in modified gravity models compared to in $\Lambda$CDM. We also find that the speed profiles of the filaments are enhanced, consistent with theoretical expectations. Overall, our results suggest that cosmic filaments can be an effective complementary probe of screened modified gravity theories on Mpc scales.


Introduction
The current standard model of cosmology is ΛCDM. By introducing two major unknowns in the theory, cold dark matter (CDM) and a cosmological constant (Λ), this model provides an excellent fit to many observations on large scales, such as the cosmic microwave background (Melchiorri et al. 2000;Netterfield et al. 2002;Bennett et al. 2013;Planck Collaboration et al. 2016) and the late-time acceleration of the Universe as measured from Type Ia supernovae (Riess et al. 1998;Perlmutter et al. 1999;Tonry et al. 2003). There are many alternatives to the ΛCDM model (see, e.g., Amendola & Tsujikawa 2010), and one of the simpler ones is to introduce a scalar field to the Einstein-Hilbert Lagrangian. When coupled to matter, this scalar field gives rise to an additional gravitational force componentoften called a fifth force (see, e.g., Mota & Shaw 2006;Clifton et al. 2012;Amendola et al. 2013). Gravitational tests on Earth and within the Solar System, however, constrain the strength of a potential fifth force to be orders of magnitude lower than its Newtonian counterpart (for a review of local constraints see Will 2006). This implies that if a scalar field coupled to matter exists in the Universe, the fifth force must be somehow suppressed in certain environments to satisfy the current constraints. There are many such screening mechanisms that hide the fifth force in high-density regions (see, e.g., Khoury 2010;Joyce et al. 2015, for reviews).
In this paper, we study two different types of screening mechanisms: the symmetron model (Hinterbichler et al. 2011), and the chameleon screening mechanism (Khoury & Weltman 2004). Both models suppress the fifth force in high-density regions but mediate a long range force in low-density regions. However, while the former relates the vacuum expectation value and the coupling strength between the scalar and matter, the latter alters the scalar field mass -and, thus, the range of the fifth force -based on the density of the region. In order to detect such a screening mechanism, it is necessary to probe scales beyond the solar system and regions which are not screened. It is therefore worthwhile to search for signatures of these models in the cosmic web of large scale structure.
The cosmic web is the complex, hierarchical distribution of matter on large scales (Bond et al. 1996;Libeskind et al. 2018). It is composed of over-dense dark matter halos, in which galaxies reside; long, filamentary structures that connect and feed the halos; flat walls or "pancakes"; and vast, under-dense voids. The web-like arrangement of dark matter is reflected in the galaxies, which are biased tracers of the underlying density field, as has been observed by, for example, the Sloan Digital Sky Survey (SDSS; Tegmark et al. 2004) and the 2MASS redshift survey (Huchra et al. 2012).
Different structures of the cosmic web have previously been explored as potential testbeds for modified gravity. In particular, dark matter halos and clusters have been studied extensively in the literature, including their mass functions (Schmidt 2010;Ferraro et al. 2011;Brax et al. 2012Brax et al. , 2013Davis et al. 2012;Hagstotz et al. 2018), density profiles (Lombriser et al. 2012), velocity statistics (Arnold et al. 2014;Hellwing et al. 2014;Gronke et al. 2015a) and even internal properties, such as their spin (L'Huillier et al. 2017). There has been some attention to the effect of screening on the general structures of the cosmic web (Falck et al. 2014(Falck et al. , 2015, and more recently, cosmic voids have been investigated a potentially powerful probe of modified gravity theories (Cai et al. 2015;Zivick et al. 2015;Voivodic et al. 2017;Falck et al. 2018;Cautun et al. 2018).
In this paper, we explore the effect of symmetron and chameleon screening mechanisms on the filaments of the cosmic web. The filaments are good candidates to probe modified gravity, as their relatively low mass-density means they are unlikely to be screened. Detecting signatures of the fifth force may be easier in filaments compared to halos, where the enhancement Article number, page 1 of 11 arXiv:1807.07287v2 [astro-ph.CO] 14 Sep 2018 is mostly seen near the edge of the halos (see, e.g., Falck et al. 2015). In spite of this, filaments were previously mostly ignored as a probe of (screened) modified gravity.
N-body simulations are essential tools for studying the evolution of the cosmic web due to its non-linear nature (see reviews, e.g., by Bagla 2005;Dehnen & Read 2011). Various gravitational N-body simulation codes have been developed which solve for the evolution of an additional gravitational force component stemming from different modified gravity theories (see Winther et al. 2015, for a comparison of modified gravity simulation codes). We will use the N-body code ISIS (Llinares et al. 2014), which takes into account a fifth force mediated by an additional scalar field and includes both the symmetron (Hinterbichler et al. 2011) and Hu-Sawicky f (R) models (Hu & Sawicki 2007).
The identification of cosmic filaments in N-body simulations is a non-trivial task as filaments are harder to define and identify than halos or voids. Filament identification techniques vary greatly and can be based on graph theory, stochastic point processes, or topological methods (Sousbie 2011;Alpaslan et al. 2014;Tempel et al. 2016;Libeskind et al. 2018). For this work, we will use DisPerSE (Sousbie 2011;Sousbie et al. 2011) which identifies persistent topological features in discretely-sampled distributions.
The paper proceeds as follows. In Sect. 2, we briefly introduce the modified gravity theories analyzed in this paper. In Sect. 3, we describe the N-body simulations used and our methods for identifying and defining cosmic filaments in the distribution of dark matter particles. We then examine the filament properties and their radial density and speed profiles in Sect. 4, and we investigate how these differ in modified gravity compared to ΛCDM. Finally, we conclude the paper in Sect. 5.

Modified gravity
We investigated extended theories of gravity which include an extra canonical scalar degree of freedom, φ, in the gravity sector (Amendola & Tsujikawa 2010). These can in general be written in the Einstein frame as: where M Pl and H are the the reduced Planck mass M −2 Pl ≡ 8πG, κ 2 = 8πG, R is the Ricci scalar, and S M is the matter action. We note that the matter fields couple to the Jordan frame metricg µν which is connected to to the Einstein frame through g µν = A 2 (φ)g µν , where A(φ) is the coupling function. Here and henceforth, we indicate quantities in the Jordan frame with˜.
When minimizing this action with respect to g µν , an additional force, referred to as the fifth force, arises. The strength of this force can be quantified as γ ≡ |F Fifth |/|F N |, where F N is the Newtonian gravitational force. In particular, local gravity experiments constrain γ 1 in the solar system (Will 2006;Bertotti et al. 2003). However, on larger scales the current constraints are weaker, allowing γ ∼ 1.
We consider two gravity models with two different screening mechanisms, both of which suppress the fifth force on solar system scales, therefore avoiding gravity bounds (Khoury 2010;Joyce et al. 2015), while still modifying gravity in lower density environments on cosmological scales.

Symmetron
In the symmetron model (Hinterbichler et al. 2011), the coupling strength between the scalar field and matter is proportional to the vacuum expectation value (VEV). In regions of low mass density, the VEV becomes large, increasing the coupling strength. On the other hand, in regions of high mass density, the VEV becomes small and decouples from matter, thus the fifth force is screened in these regions.
The action of the symmetron models is given by Eq.
(1) with the symmetric potential, and the coupling function, where µ and M are mass scales and λ is a dimensionless constant. These choices of V(φ), A(φ) imply that the effective potential has a parabolic shape if the local matter density is larger than a critical threshold and a "sombrero" shape otherwise. The resulting minima at φ = 0 and offset from it yield a screened and unscreened fifth force, respectively. As in Winther et al. (2012), we can rewrite the parameters (µ, M, λ) to more physical ones. The range of the fifth force in a vacuum can be written as the scale factor at the time of the symmetry breaking is given by and finally, the coupling strength can be written as In the symmetron model, the maximum enhancement of gravity is given by γ max = 2β 2 1 − (a ssb /a) 3 (e.g., Gronke et al. 2015b).

2.2.
Hu-Sawicky f (R) gravity f (R) gravity is perhaps the most well-studied modified gravity model. In this model, one adds a function, f (R), of the Ricci scalar to the Einstein-Hilbert action in the Jordan frame: In this paper, we consider the Hu-Sawicky model (Hu & Sawicki 2007), which makes use of chameleon screening. In this model, f (R) is given by The parameters c 1 , c 2 , and n are positive constants, and m 2 = H 2 0 /Ω m0 . By requiring that the model yields late-time accelerated expansion, in the form of an effective cosmological constant, the parameters c 1 and c 2 may be reduced down to a parameter f R0 (see, e.g. Llinares et al. 2014;Arnalte-Mur et al. 2017), given as The choice of f R0 and n fully specifies our models, and they relate to the range of the fifth force in the cosmological background today as The maximum enhancement of gravity for all the f (R) models is γ max = 4/3 (e.g., Clifton et al. 2012).

Simulations
We use the N-body simulation code ISIS (Llinares et al. 2014), a modification of RAMSES (Teyssier 2002). ISIS follows the evolution of a scalar field in the quasi-static limit and changes the gravitational force accordingly. Each simulation has 512 3 dark matter particles in a box with side length 256 Mpc h −1 . The parameters used for the background cosmology are (Ω m0 , Ω Λ0 , H 0 ) = (0.267, 0.733, 72 km s −1 Mpc −1 ), and the corresponding dark matter particle mass is m p = 9.3×10 9 M h −1 .
In addition to the ΛCDM model, we run simulations with three sets of parameters for the f (R) model and four for the symmetron model. The parameters used for these models are listed in Table 1. Here, we list also the range of the fifth force in vacuum, λ φ , and the maximum enhancement of the fifth force, γ max . All simulations are analyzed at redshift z = 0.

Identifying cosmic filaments using DisPerSE
DisPerSE 1 (Discrete Persistent Structures Extractor; Sousbie 2011) identifies persistent topological features in discretelysampled distributions such as N-body simulations and galaxy catalogs Laigle et al. 2018;Malavasi et al. 2017). It is based on discrete Morse theory (Forman 2002), which partitions space into a series of n-dimensional domains defined by the gradient of a function (e.g., the density field) and the connectivity of critical points (maxima, minima, and saddle points). Structures are identified as components of the Morse-Smale complex of the density field. DisPerSe defines filaments as one-dimensional lines connecting a maximum (critical point of third order) and a two-saddle (critical point of second order).
DisPerSE assigns a density to each discrete tracer using the Delaunay tessellation (Delaunay 1934;Schaap & van de Weygaert 2000), which divides the three-dimensional volume into tetrahedra with vertices at the particle positions. The density at each dark matter particle's location is then given by the volume of its surrounding tetrahedra. In DisPerSE, a random subsample of the N-body particles can be input to the Delaunay tessellation for computational convenience while retaining the main features of the filaments ; in this paper, we use a subsample of 188 3 particles.
The finite sampling of the density field by discrete tracers such as dark matter particles induces Poisson noise. In DisPerSE, noisy features are removed by filtering the Morse-Smale complex using persistence theory (Edelsbrunner et al. 2002), which measures the robustness of pairs of critical points. The relative importance of a given critical pair is assessed by its persistence ratio, defined as the relative density contrast of the critical pair, which can be given a significance value in units of "σ" (see Eq. 4 of Sousbie 2011). In this work, we choose a value of 4σ. Sousbie (2011) has shown that for a 2σ threshold, the probability of a filament to be noise is roughly 5%; increasing this to a 4σ threshold reduces the noise probability down to roughly 0.006%. To determine which σ value to use, we compared the filament distribution for different values of sigmas at different particle subsamples and found that the 4σ distribution closely resembled the case using a 256 3 particle subsample at 5σ. Figure 1 illustrates the filaments identified by DisPerSE, which trace the underlying density field of the N-body simulation. The endpoints of the filaments connect at clusters, and the endpoint which is defined as the maximum usually resides within a halo identified by AHF (Knollmann & Knebe 2009), while most of the two-saddles are found not to be within AHF halos. The output of the DisPerSE code is a list of the coordinates of critical points that define filament segments. In order to measure filament properties, such as mass and density, to study the effect of the fifth force on these filaments, we next need to determine which dark matter particles belong to each filament.

Assigning particles to filaments
The basic steps of our particle-collection procedure are as follows: firstly, for a subset of particles around each filament, calculate their distance to the filament and their distance along the filament; secondly, filter out particles near the end-points of the filaments that are likely to be gravitationally-bound to nearby halos; and finally remove particles further away than the "edge" of the filament, determined by its radial density profile. We note that for this procedure and for the analysis of the results, we use the full set of 512 3 particles in each simulation: the sub-sampling was only for the identification of the filament segments. Each filament in the slice is given a different color, and they are plotted over the background density determined from the dark matter particles.

Particle-filament distances
In order to reduce computing time, we first define a rectangular box around each filament that extends 6 Mpc h −1 beyond the maximum extent of the filament, and we only compute distances for particles within this box. Each filament consists of a set of segments, defined by their endpoints s 1 and s 2 . We want to compute the distance of a particle to the line defined by these endpoints, which can be parameterized as s(t) = s 1 +t(s 2 −s 1 ), where t ∈ [0, 1]. The minimum distance from this line to a particle located at p is then D = |s(t s ) − p| with Since many particles will have a minimum distance to a point on the line outside the bounds given by the endpoints, we enforce that t s is within the range [0, 1] by setting all t s < 0 to 0 and t s > 1 to 1. The final distance to the filament is then the minimum distance for all segments, which we characterize by t p = t s + N for the t s which minimizes D, where for the first segment N = 0, for the second N = 1, and so on.

Removing halo particles
As filaments identified by DisPerSE connect regions of high density, it is very likely that the filaments are contaminated with dark matter particles that are gravitationally bound to halos.
Since we are only interested in the particles defining the filament itself, and not the halos, we must find a way to filter out these particles. We first compute the number density profile along the filament by computing a histogram of the t p values of all particles in the sub-box around the filament. We then compute the average number of particles in all t p bins, and we remove the particles belonging to the bins closest to the endpoints of the filament (i.e., the halos) with a number density above this average. If the number density is higher than average in a bin that is not connected to the endpoints. For example, if there is a small density peak within the filament, it is not removed.

Defining the filament edge
Once we have removed particles near the endpoints of the filament that may belong to a halo, we next determine the "edge" of the filament in order to define its thickness (and also its mass). The edges of structures in N-body simulations are notoriously hard to define, even for halos (for a discussion, see Knebe et al. 2013). In the case of halos, the theory of spherical collapse (e.g. Lemaître 1958;Cooray & Sheth 2002) yields a theoretical foundation for a density threshold which divides the collapsed or virialized region from the background. Cosmic filaments, however, are -as opposed to halos -non-virialized structures, and a similarly motivated definition does not exist. We thus choose an approach resembling the the commonly-used definition of the virial radius of dark matter halos. The filaments are assumed to consist of segments with cylindrical shape. We then define the thickness as the radius at which the cumulative density profileρ(< r), that is, the average density within r, falls below 10ρ c0 , where ρ c0 denotes the critical density today. Filaments that never reach this density threshold are removed from the dataset.

Pruning the filament catalog
At this point we have a catalog of filaments (and their constituent particles), each with a defined length, thickness, and mass. In Section 4, we ignore filaments which fulfill one of the following conditions: 1. filaments which are smaller than 1Mpc h −1 , 2. filaments larger than 20 Mpc h −1 , 3. filaments with thickness larger than 10 Mpc h −1 , or 4. filaments with mass larger than 10 15 M h −1 .
The first point is meant to filter out filaments which are not well resolved due to the mass resolution of our simulation, and removes between 13 to 18% of all identified filaments in the simulation for the different gravity models. The latter three points remove filaments of which there are too few at these sizes (large mass, length, or thickness) to make statistically significant conclusions. We find that less than 1% of the filaments in our dataset are removed because of the latter three points. The fact that few large filaments are found in our simulations is supported by observations. Tempel et al. (2014) found that 0.01% of the filament volume, detected in the SDSS galaxy survey, have lengths less than 1 Mpc h −1 or larger than ∼30 Mpc h −1 .

Velocity components
Though most of our velocity analysis will be focused on the particle speed, that is, the magnitude of its velocity vector, we also decompose the velocity of each particle into radial, v r , and tangential, v t , components with respect to the segment of the filament that the particle is closest to. The radial component is orthogonal to the filament segment and is defined as positive when the particle is moving away from filament axis. The tangential component runs parallel to the filament segment: we define the positive direction as that which starts from the two-saddle and ends at the maximum point. We note that the third, "circular" component is ignored; results for radial and tangential speeds of particles within filaments are given in the Appendix.

Error computation
For the global filament properties, such as the number of filaments at a given length bin, we compute Poissonian errors: where N is the number of data points in the given bin. For stacked profiles, such as the speed profiles of different filaments, we compute the error on the mean, σ/ √ N, where σ is the standard deviation of the distribution and N is the number of filaments in the respective stack. Furthermore, the uncertainties for relative differences are computed by Gaussian error propagation.

Global properties
In this section, we present the distributions of filament length, mass, and thickness in ΛCDM and modified gravity simulations, as well as the correlations between these properties. The shaded regions of the global properties show the Poisson noise while the errors for the relative differences are computed by Gaussian error propagation.

Filament lengths
The distributions of the filament lengths in the symmetron and f (R) models are found in Fig. 2 and 3, respectively. Both figures show the length distributions in the top panel and the relative differences of the modified gravity models with respect to ΛCDM in the bottom panel. We note that the uncertainties are larger for longer filaments, due to the fact that there are fewer filaments with these lengths.
With the exception of fofr4, there are a larger number of short filaments in the modified gravity simulations than in ΛCDM, where "short" is L 10 Mpc h −1 for the symmetron models and L 7 Mpc h −1 for the f (R) models. This is especially true for the symmetron C case, which has more than 30% more filaments at length scales ∼1-3 Mpc h −1 . Despite the fact that fofr6 possesses the shortest fifth force range (cf. § 3.1), and is thus closer to ΛCDM than the other f (R) simulations, it has more short filaments than ΛCDM, while fofr4 is the closest to ΛCDM. While there are, overall, more filaments in each modified gravity model, long filaments appear to be more common in the ΛCDM case. The overall shape of the filament length distribution agrees well with the findings of Tempel et al. (2014) and Bond et al. (2010), which used different methods of filament-finding on the SDSS galaxy catalog. Figure 4 and 5 show the filament mass distributions of the symmetron and f (R) models, respectively. The figures also include, in the bottom panels, the relative differences of the modified gravity models with respect to ΛCDM. As above, due to the low number of filaments at large masses, the error becomes significantly larger.

Filament masses
The number of filaments, of all mass scales, are generally larger in the symmetron models. Symmetron C has the most filaments of small mass, while symmetron D has the most filaments of large mass. In the case of the f (R) models, there are more filaments of small mass in fofr6, and more of large mass fofr4, compared to in ΛCDM. A similar result was found for dark matter halos: more of small mass in fofr6, and more of large mass in fofr4, compared to in ΛCDM (Lombriser et al. 2013;Li & Hu 2011). We discuss this connection between halos and filaments in more detail in Sec. 5.  Fig. 2, but comparing the f (R) models instead of the symmetron models.

Filament thickness
In Fig. 6, we show the thickness distributions of the ΛCDM model and the symmetron models at the top, with the relative differences with respect to ΛCDM at the bottom. Similar results for the f (R) models are shown in Fig. 7. The thickness distribution can be compared to the findings of Bond et al. (2010), who define the filament width based on the shape of the density profile 2 , and did not use a density threshold as we did. As a consequence, their filaments are generally thicker, that is, they find more filaments with a larger thickness, between 5 Mpc h −1 to 10 Mpc h −1 . It is important to note that the thickness of our filaments depend on the density threshold parameter described in Sec. 3.3.3. The overall thickness will increase and decrease if the threshold is set lower and higher, respectively. While all modified gravity mod-    Fig. 4, but comparing the f (R) models instead of the symmetron models. els possess more filaments of thicknesses between ∼ 0.1 Mpc h −1 and 0.3 Mpc h −1 , the differences with respect to ΛCDM are negligible for filaments with thicknesses between 0.3 Mpc h −1 and 2 Mpc h −1 . The one exception is the symmetron C model, which has more than ∼20% filaments in all thickness ranges. The recovered thickness is naturally linked with the filaments' (orthogonal) density profile which we discuss in § 4.3.

Correlations between the filament properties
In this section, we show how the three global properties (mass, length, and thickness) correlate with each other. Intuitively, they should all be correlated. For instance, a long filament, or even a thick filament, should have an increased mass because it covers a larger amount of space and thus includes more particles in the filament. The errors are computed by the standard deviation,   while the error of the relative differences are still computed by Gaussian error propagation.
The average filament mass as a function of the filament thickness is shown in the top row of Fig. 8, and the relative differences with respect to the ΛCDM model are shown in the bottom row. As we can see, the average mass of the filament does indeed increase as the thickness of the filament becomes bigger. This is not surprising, because a thicker filament would contain more particles. This result also agrees with the findings of Cautun et al. (2014) (see, e.g., their figure 39), which compares the filament linear density with the filament diameter. The authors also conclude that thicker filaments have higher masses. This correlation is enhanced in all analyzed modified gravity models -with the exception of the symmetron C model.
The average mass as a function of the filament length is shown in the top row of Fig. 9, with the relative differences to ΛCDM in the bottom row. The average mass does indeed increase as the filament becomes longer, for the same reason as The average mass of the filaments with a given length. Bottom: The relative difference of the modifies gravity models with respect to the ΛCDM model an increased thickness: it contains more particles. However, the correlation between mass and length is much weaker than between mass and thickness. The difference between the average masses of the modified gravity theories with respect to ΛCDM varies in this case. Some models, such as fofr4 and symmetron D, have a larger average mass, while other models, such as fofr6 and symmetron C, have a smaller average mass. Finally, we compare the average thickness for filaments of different length in the top row of Fig. 10, and the relative differences with respect to ΛCDM are shown in the bottom row. Interestingly, the thickness decreases as the filament becomes longer. This anticorrelation can be explained due to the fact that, by definition, the filaments start and end at critical points in the density field, and even after potential halo particles are removed (as described in § 3.3.2), the endpoints of the filaments are still overdense compared to the interior regions. Shorter filaments have these overdense endpoints closer together, resulting in denser, and thus thicker, filaments. The filaments of almost all modi-  fied gravity models have, on average, a smaller thickness than in ΛCDM at all length scales, with fofr4 being the closest to ΛCDM.

Density profiles
We first study the radial density profiles of cosmic filaments and how they change under different gravity models. We only present the relative difference of the density profiles with respect to the ΛCDM because the profiles themselves are very similar. The average density profiles for all filaments in each model are shown in Fig. 11. This figure shows that, on average, the filaments have a larger density in the center than the outer edge. The displayed uncertainties of the density profiles are the error of the mean (of the corresponding stack). The uncertainties on the relative differences follow through Gaussian error propagation.
To determine whether screening affects filaments of different mass, we consider how the density profiles change at dif- The average densities are generally larger in the modified gravity models. The difference between symmetron C and ΛCDM becomes smaller as the filament mass becomes larger, and is practically zero for masses M ≥ 10 14 M h −1 . This is most likely due to the screening effect being active for these most massive -and most dense -filaments. In the case of symmetron D, the density profiles of the most massive filaments deviate more from ΛCDM in the outer regions of the filament. Interestingly, the models that are supposed to deviate the least from ΛCDM (e.g. fofr6 and symmetron A), have, on average, larger densities for filaments of small and intermediate masses. This implies that cosmic filaments may be a good probe of screened modified gravity models in which halos are already self-screened.

Particle velocities and speed profiles
In this section, we analyze the speeds of particles within filaments, as well as their radial speed profiles. Generally, we find that the radial component of the particle velocity accounts for most of the total speed, which means that the tangential speed will only give a small contribution. This also indicates that most particles will move toward the filament as opposed to along it. The results for both the radial and tangential speeds are found in Appendix A. The shaded regions in the plots show the uncertainties computed as for the density profiles.
We show the average speed of the particles within a filament as a function of the filament mass (such that each filament contributes one average particle speed) in the top two panels of Fig.  13, along with the relative differences with respect to ΛCDM in the bottom two panels. We find that more massive filaments contain particles with larger speeds, as expected. The models symmetron A and fofr6 deviate the least from ΛCDM, while symmetron D and fofr4 deviate the most, in their respective classes of models. In general, the difference with respect to ΛCDM decreases as the filament mass increases, which is likely a result of the screening beginning to be effective for filaments with larger mass. The total speed profiles of particles within filaments for filaments of similar mass, for the ΛCDM & symmetron and ΛCDM & f (R) models, are shown in Fig. 14 and 15, respectively. The top panels show the speed profiles while the bottom images shows the relative difference of the speed profiles with respect to the ΛCDM model. The mass ranges are M ∈ [10 12 , 10 13 ]M h −1 , M ∈ [10 13 , 10 14 ]M h −1 and M ∈ [10 14 , 10 15 ]M h −1 , as they were for the density profiles (see Fig. 12). Below, we refer to these mass ranges as filaments of small, intermediate and large masses.
For all mass ranges, the average speed is generally larger near the filament center, compared to the outer regions of the filament. Similar to speed profiles of clusters of galaxies, this is due to the fact that in the central region a larger fraction of the total energy is kinetic. This also agrees very well with the result of Fig. 11, which indicates that the density is larger, on average, near the center. As expected, by increasing the filament mass, the average speed will also increase because more massive filaments create deeper gravitational potentials.
The relative difference of the average speed profiles agrees mostly with our theoretical expectations. For instance, comparing the f (R) models with ΛCDM (see Fig. 15), fofr6 deviates the least from the ΛCDM, which is expected as the fifth force using this parameter has the smallest range. Fig. 15 also shows that the speed profile of the fofr6 model follows the typical behavior of screened modified gravity theories: the screening mechanism is active at small radii, leading to a speed profile similar to ΛCDM. Furthermore, this threshold between screened and unscreened radii shifts inward for more massive filaments, which also happens in dark matter halos (Falck et al. 2015). From the same figure, it also becomes clear that filaments are not screened at all in fofr4, which causes the deviation with respect to ΛCDM to become larger in all filament mass ranges. Figure 14 shows the average speed profiles for the symmetron model. Similarly to the f (R) results, the difference with respect to ΛCDM is smaller in large filaments and also smaller near the center of the filaments. The differences between the symmetron model and ΛCDM again coincide with the theoretical expectations. Unsurprisingly, the deviations become larger for larger coupling strengths (such as symmetron C compared to symmetron A) for the least massive filaments. Nevertheless, the deviation between symmetron D and ΛCDM is the largest for the symmetron models because in this case the screening mechanism becomes active at higher filament densities, compared to the background density.

Conclusions
In this paper, we have explored how filaments of the cosmic web change in different modified gravity models. First, we have analyzed the global properties of the filaments, such as their masses, lengths, and thicknesses, and studied how they correlate with each other. In addition to the global properties, we have also examined the average density and speed profiles.
We found that the length distribution is overall comparable to the findings of Tempel et al. (2014) and Bond et al. (2010). However, filaments are generally shorter (L 10 Mpc h −1 ) in all the modified gravity models studied. While the mass distributions were also comparable to observations (Bond et al. 2010;Tempel et al. 2014), the screened modified gravity models clearly altered this distribution at the 10% level. Interestingly, this deviation with respect to ΛCDM is different for each analyzed model, and can thus potentially be used to constrain the parameter space of screened modified gravity theories.
Generally, we found that the mass distribution of filaments resembles that previously found for halos (e.g., Ferraro et al. 2011;Lombriser et al. 2013). For instance, both the halo and filament mass distributions show that fofr4 and fofr6 have more filaments and halos of large and small masses, respectively. This suggests that these deviations in the respective mass functions have similar causes. In screened modified gravity models, halos become (self-)screened if their mass exceeds a certain threshold (µ 200 in the notation of Gronke et al. 2015b; also see Mitchell et al. 2018). Below this threshold, they grow at an accelerated rate (with the boost given by ∼ γ max ), and above it at the "normal" rate set by GR. This leads to a maximum increase in abundance right around this screening threshold where the change in growth rate is largest. For instance, amongst the f (R) models analyzed, the fifth force range -and, consequently, the screening threshold -is largest for fofr4 and smallest for fofr6 (cf. Table 1). Thus, within this simplistic picture we expect more small mass filaments, in agreement with our findings.
We found furthermore that the filament thickness distributions did not change dramatically amongst the models analyzed, and the greatest deviations with respect to ΛCDM were for filaments with thicknesses 0.1-0.3 Mpc h −1 . As our definition of the filaments' boundaries is a fixed density threshold, thicker filaments are also denser on average. Thus, the fact that there are no deviations for filaments with thickness R T 1 Mpc h −1 (with the exception of the "extreme" symmetron C model with γ max ∼ 7) is a signature of screening.
We also showed that the three global properties are correlated with each other. For instance, thicker filaments are also more massive (cf. § 4.2). However, these correlations change slightly between the different gravity scenarios -which has to be taken into account when interpreting the deviations of the observables.
Apart from these "global properties", we studied the radial density profiles of cosmic filaments, too. On average, the density near the center of the filaments is larger than at the edges. As expected, the filaments were denser in the modified gravity models due to the additional fifth force causing the particles to cluster more (with a ∼ 20% increase in matter density). However, some deviations between modified gravity and ΛCDM were contrary to the theoretical expectations. For instance, in the fofr6 model, that is, the f (R) model with the shortest fifth force range, the filaments were on average denser than in the fofr4 case. Similarly, for the symmetron models, the average density of symmetron D, in which the fifth force was active the earliest, was smaller than both symmetron A and symmetron B.
The velocity field is a direct probe of the underlying gravitational forces and thus a good proxy for the strength of the fifth force. We found that the filament velocity profiles, for the most part, coincided with our theoretical expectations. For instance, an increased filament mass correlates with a larger average speed, which is due to the fact that a larger filament mass would create a deeper gravitational potential and thus result in a greater kinetic energy. Furthermore, the differences in the speed profiles between the modified gravity models and the ΛCDM model were mostly consistent with the theoretical predictions. For instance, in the symmetron model, symmetron A deviates the least from ΛCDM, while symmetron D deviates the most. Interestingly, for the most massive filaments, the fifth force is screened in the symmetron C model near the center of the filament. This can be seen from both the density profiles and the speed profiles.
Overall, we found that the distribution, structure, and dynamics of cosmic filaments are (heavily) altered for some screened modified gravity models. Furthermore, our results show that the probes analyzed were not altered uniformly -for example, an overall increase for all models -but instead display a variety of signatures with properties being enhanced, constant, or even decreased with respect to ΛCDM. While this makes the deviations harder to interpret, it can be used as an advantage, for instance, to distinguish between different modified gravity scenarios. These non-linear changes of the filament properties are due to the fact that the filaments are less dense than halos and thus less likely to be screened. While these results suggest that cosmic filaments can be an excellent probe of modified gravity, the identification and analysis of filamentary structure in emission or absorption data is very challenging, and it remains to be seen whether the resulting data-sets will have much constraining power over the screened modified gravity parameter space. components We have in § 4.4 only shown the total speed of the particles within the filament. The total speed, however, does not describe the general direction of the particles themselves, such as whether the particles is moving away or toward the filament center or how the particles moves along the filament axis. By intuition, the particles themselves should mostly move toward the filament center. In addition to this, the average tangential component should be close to zero as the particles should not have any preferential direction along the filament axis. Although, this is not true for the particles near the end points, where the high density regions usually reside, but because they move in opposite directions, the average tangential speed should therefore cancel each other out.
In this section, we show the average radial and tangential speeds of different mass ranges. The average radial speed within a mass bin is shown in the top row of Fig. A.1, with the relative differences with respect to ΛCDM in the bottom row. The radial component is defined such that a positive value indicates that particles are moving outward the filaments. As expected, we see that the particles are, on average, moving toward the filament. The increase in radial speed due to an increased mass is caused by the fact that the gravitational potentials are deeper for the more massive filaments. The relative differences can also be compared to the ones in Fig. 13, and it is easy to see that the relative differences between the two figures are almost the same. This indicates that the radial speed will make up the majority of the total speed. Figure A.2 shows the average tangential speed at different mass bins. For the least massive filaments, the particles generally have no preferential direction along the filament axis. For the intermediate massed filaments, the particles tend to move toward the endpoint which is considered to be the maximum point. Meanwhile, for the most massive filament, most particles appears to move toward the endpoint which defines the two-saddle. We found that most maximum points resided within a halo defined by the AHF halo catalog, while most of the two-saddles did not. The maximum points are therefore more dense. However, it is puzzling that the particles in the most massive filaments tend to move toward the less dense endpoint.