The role of the elaphrocentre in void galaxy formation

Voids may affect galaxy formation via weakening mass infall or increasing disk sizes, which could potentially play a role in the formation of giant low surface brightness galaxies (LSBGs). If a dark matter halo forms at the potential hill corresponding to a void of the cosmic web, which we denote the 'elaphrocentre' in contrast to a barycentre, then the elaphrocentre should weaken the infall rate to the halo when compared to infall rates towards barycentres. We investigate this hypothesis numerically. We present a complete software pipeline to simulate galaxy formation, starting from a power spectrum of initial perturbations and an N-body simulation through to merger-history-tree based mass infall histories. The pipeline is built from well-established, free-licensed cosmological software packages, and aims at highly portable long-term reproducibility. We find that the elaphrocentric accelerations tending to oppose mass infall are modest. We do not find evidence of location in a void or elaphrocentric position weakening mass infall towards a galaxy. However, we find indirect evidence of voids influencing galaxy formation: while void galaxies are of lower mass compared to galaxies in high density environments, their spin parameters are typically higher. For a fixed mass, the implied disk scale length would be greater. Tangential accelerations in voids are found to be high and might significantly contribute to the higher spin parameters. We find significantly later formation epochs for void galaxies; this should give lower matter densities and may imply lower surface densities of disk galaxies. Thus, void galaxies have higher spin parameters and later formation epochs; both are factors that may increase the probability of forming LSBGs in voids.


INTRODUCTION
The role of the void environment in galaxy formation is worth exploring. A particular case of interest is that of giant low surface brightness galaxies (LSBGs). The formation mechanisms of LSBGs (Sandage & Binggeli 1984;Bothun, Impey, Malin & Mould 1987) with high masses and low star formation rates remain unclear. Hoffman, Silk & Wyse (1992) presented a peaks-in-peaks structureformation calculation arguing that voids are likely to play a major role in the formation of giant low surface brightness galaxies. Here, we numerically investigate the role of voids in galaxy formation, with a particular emphasis of how voids may provide some of the characteristics leading to LSBG formation, by developing a reproducible pipeline that combines and builds on existing tools, which have developed considerably over the intervening three decades.
It has long been known that a high volume fraction of the Universe consists of large underdense regions (e.g. Gregory & Thompson 1978;de Lapparent, Geller & Huchra 1986). More recent measurements have increased our knowledge about these underdense regions, now known as cosmic voids or voids in the cosmic web. Sloan Digital Sky Survey (SDSS) measurements show that cosmic voids dominate the volume of our Universe, constituting 60% of the volume in Pan et al. (2012)'s analysis. Voids vary by about an order of magnitude in size, spanning scales from 10 to 200 Mpc/ℎ −1 , especially as found in various SDSS-related analyses, with higher numbers of smaller voids, and quantitative population estimates sensitive to the methods of defining and detecting voids (Hoyle & Vogeley 2002Pan et al. 2012;Nadathur & Hotchkiss 2014;Sutter et al. 2014;Pisani et al. 2014Pisani et al. , 2015bMao et al. 2017). Depending on the chosen catalogue and void detection algorithm, a modest fraction of galaxies can be considered as being located in voids. For example, Pan et al. (2012) estimate that about 7% of galaxies are located in voids. Since the fraction of volume occupied by voids may be as high as 80% (for theoretical estimates, see e.g. Colberg et al. 2008;Cautun et al. 2015), they are suspected to play a key role in relation to dark energy (e.g. Buchert et al. 2016).
In the cosmological galaxy formation context, we denote the gravitational potential hill corresponding to a cosmic void on megaparsec scales as an 'elaphrocentre' in order to emphasise its gravitational role in opposition to that of a barycentre 1 . The impact of voids on galaxy formation remains largely unknown. Whereas galaxies in clusters and the walls of the cosmic web typically undergo a gravitationally very active evolution with many mergers, void galaxies tend to have few mergers. Galaxy mergers usually lead to bursts of star formation, making galaxies briefly much brighter. It is typical in modelling the formation of a galaxy in a gravitational barycentre -a knot of the cosmic web -to approximate the surrounding universe via the Newtonian iron-spheres theorem. Relativistically, the conditions of Birkhoff's theorem, for spherical symmetry and an asymptotically 4-Ricci flat universe (Birkhoff & Langer 1923), are not satisfied for a flat Friedmann-Lemaître-Robertson-Walker (FLRW) model, but nevertheless provide heuristic motivation for applying the iron-spheres theorem. However, the situation in a void is different: the density contrast is much weaker, and there is a mass deficit below the mean density, rather than a mass excess above the mean density. Approximating a strong overdensity as being embedded in a surrounding empty universe that is Newtonian is likely to be less inaccurate than approximating a modest (in terms of linear mean density) underdensity in the same way, since gravity is attractive. What is effectively antigravity in voids -in comparison to the surroundings -is unlikely to be well modelled by a Newtonian approximation. Indeed, relativistically, to reach turnaround, an overdensity has to pass through a strongly positive spatial curvature phase (Roukema & Ostrowski 2019;Ostrowski 2020;Vigneron & Buchert 2019), after which it virialises at an overdensity of a few hundred times the mean density (e.g. Lacey & Cole 1993). Since a void tends to have negative spatial curvature (for a flow-orthogonal spacetime foliation), an overdensity inside a void will have difficulty forming. If it forms nevertheless, the negative spatial curvature environment will tend to weaken matter infall, weakening the star formation rate. A void can also be thought of as approximated by a spatially compact domain in a relativistic Milne model -empty of matter and spatially hyperbolic (often called "open") -in which structure forms more slowly than in the idealised background FLRW model.
A pseudo-Newtonian way of thinking about this is that compared to a background FLRW model, a galaxy in the emptiest parts of a void -the 'elaphrocentre' -feels a weak antigravitational environmental force around it, since the void is underdense compared to idealised average regions of space. Another toy model way of thinking about the elaphrocentric effect is as follows. Let the core of a dark matter halo be modelled as forming via linear theory followed by the standard pseudo-Newtonian spherical collapse approximation in a background FLRW model, and then by slow uniform mass infall induced by the elaphrocentric environment. This would appear as a history of mass infall that is spread out in time and gradual, rather than fast and sudden. If slow gradual infall is interpreted as the late collapse of the outer parts of the halo, then this corresponds to collapse at late epochs, when the FLRW critical density has dropped, implying a greater virial radius for a fixed total halo mass.
Based on these heuristic arguments we assume that if a dark matter halo forms at or near an elaphrocentre, then the merger rate of small haloes into that halo and the overall infall rate of dark matter and gas into the halo should be weaker than the usual infall rates towards a halo of similar mass at a barycentre. This effect should tend to create lower mass dark matter haloes at elaphrocentres compared to barycentres, which is generally the case: the most massive haloes form at the knots of the cosmic web. For a high-mass dark matter halo at the present epoch of a fixed mass, this elaphocentric effect should tend to increase the probability of the halo having grown in mass after its initial collapse by weak infall and a weak merger rate over a long period, rather than by an initial short burst of mass accumulation. These effects on the host haloes could lead, for a fixed mass and a fixed efficiency factor for converting available baryons into stars, to galaxies closer to the elaphrocentre preferentially forming with weak star formation rates over long time scales.
While Hoffman et al. (1992) proposed that voids play a major role in LSBG formation, here we briefly present a broader overview of LSBG formation scenarios. LSBGs were discovered more than three decades ago (Sandage & Binggeli 1984;Bothun et al. 1987), leading to the question of their formation mechanisms (Schombert, Maciel & McGaugh 2011;Schombert, McGaugh & Maciel 2013;Schombert & McGaugh 2014a,b). Interest in LSBGs was recently reignited by van Dokkum et al. (2015), who found a high abundance of large LSBGs in the Coma cluster, that they called ultra diffuse galaxies (UDGs), to distinguish them from traditional LSBGs. It is not yet clear whether LSBGs and UDGs share the same formation scenario, especially since this is dependent on their definitions. A common scenario is that UDGs form in high-spin haloes. Rong et al. (2017) find from simulations that UDGs form naturally in high-spin haloes within the ΛCDM model; Kim (2015) has started investigating this observationally. Chan et al. (2018) show that they can reproduce the observed quantities of red UDGs by imposing quenching, without assuming high spin haloes. Di Cintio et al. (2017) showed that UDGs can be produced in isolated dwarf galaxy haloes with stellar feedback and episodes of gas outflow. Jiang et al. (2019) extended this work by investigating field UDGs. Compared to other galaxies, LSBGs tend to be more spatially isolated, i.e. they tend to be somewhat elaphrocentric. Rosenbaum et al. (2009) quantified this on a 2-5 Mpc scale, finding that LSBGs are located in regions with lower galaxy number densities than those in which high surface brightness galaxies are located. LSBGs typically have low H surface densities, below around 5 ⊙ pc −2 , yielding weak star formation, with star formation rates that are approximately constant with cosmological time, rather than the exponentially declining star formation rates typically associated with high surface brightness galaxies (Di Paolo & Salucci 2020, Sect. 7.1, 7.2).
Here, we focus on the degree to which the elaphrocentric location may contribute to low surface brightness of void galaxies for a given host halo mass, via (i) a total-matter infall rate closer to being constant rather than being exponentially declining, and (ii) an enlarged disk size of a galaxy due to high spin (Rong et al. 2017) and/or an increase in the typical virial radius. (iii) We estimate the magnitude of elaphrocentric acceleration as a basis for more detailed studies.
To study this hypothesis, we present a highly reproducible (Akhlaghi et al. 2021) galaxy formation simulation and analysis pipeline (Peper, Roukema & Bolejko 2019). The packages in the pipeline are free-licensed packages, and should only require a POSIX-compatible operating system with sufficient memory and disk space for reproducing the full calculations, tables and figures, generating values that are statistically equivalent to those published here. We present the software packages and the pipeline in Sect. 2.1.
In Sect. 2.2 we propose a Voronoi-cell based definition of the 'elaphrocentre' and discuss alternative definitions of the 'centre' of a void from the literature. We present two different parameters characterising void membership and propose a criterion for use in global population comparisons between void and non-void galaxies in Sect. 2.3.1. We describe how we study the dependence of the infall of matter into galaxies (Sect. 2.3.2) and the dependence of galaxy sizes (Sect. 2.3.3) on these parameters and on the global void membership criterion. The elaphrocentres themselves are the places that are the most difficult to study via particle distributions, so in Sect. 2.3.4 we describe how we investigate accelerations near the elaphrocentres in preparation for future studies of elaphrocentric effects that might help form LSBGs.
We present our results in Sect. 3, discussion in Sect. 4 and conclude in Sect. 5. The reproducibility package for this paper is available at zenodo.4699702 and in live 2 and archived 3 repositories. The commit hash of the version of the source package used to produce this paper is 97bd6dc. The package was configured, compiled and run on a Little Endian x86_64 architecture.

Software pipeline
We provide a highly reproducible software pipeline for generating a realisation of galaxies with merger-history-tree based galaxy disk formation histories (and star formation histories, though we do not analyse these in this work) starting from early universe initial conditions. This approach not only combines existing community tools, but can also help in improving those existing tools by embedding them in a controlled software environment. Our pipeline is intended to be modular, so that the well established cosmological software tools currently chosen, can, in principle, be replaced in a modular way, provided that the user manages the input and output formats correctly. Our results are intended to be statistically reproducible. Parallelisation in several steps of the computational pipeline currently prevents byte-for-byte reproducibility.
Dark matter haloes are identified in an -body simulation, a merger history tree is created, semi-analytical recipes are used to generate galaxy disks with mass infall histories, and voids in the dark matter distribution are detected in the -body simulation. In the following we give a brief description of the codes and the key parameters used in these models. URLs and SHA512 checksums for the upstream versions of software used to produce this paper are listed in the reproducibility package. (An SHA512 checksum is a 512-bit integer computed from the bytes in a file using the SHA512 algorithm, aiming to provide a data integrity check on the file contents that is sensitive to small changes in the file.) The reproducibility structure is based on the Maneage template that aims for a high level of reproducibility (Akhlaghi et al. 2021). We follow Rougier et al. (2017) for the definitions of the 'reproducibility' of a research paper -in which independent authors attempt to use the same input data and the same source code and analysis pipeline to obtain the paper's claimed results -versus the paper's 'replicability' -in which independent authors attempt to use different but similar input data and/or a different but equivalent analysis to obtain the claimed results. Using these definitions, we believe that it should be straightforward for the reader to verify the reproducibility of our results. We expect that our results will also be replicable. Version identities of the software packages in the text below include commit hashes. Modifications that we have made to upstream versions of codes are included as patch files in the reproducibility source package (zenodo.4699702).

Initial conditions
We use -0.3.19-4b78328 (Prunet et al. 2008) to generate a set of standard initial conditions for a flat-space -body simulation with the standard 3-torus topology (often called 'periodic boundary conditions'). The package is a well-tested, parallelised package that generates peculiar velocity offsets against an FLRW background model for a standard cosmological power spectrum with Gaussian random fluctuations using the Zel'dovich approximation (Zel'dovich 1970a,b). Checks are made that the amplitude of the numerically generated power spectrum matches that of the input power spectrum. We generate a simulation with = 128 3 particles. The comoving fundamental domain size, often called the 'box size', is box = 80 Mpc/ℎ. These parameters give a dark matter particle mass of 2.03 × 10 10 ⊙ , which is a reasonable mass resolution for modest RAM and CPU resources. Together with the minimum number of particles per halo (Sect. 2.1.3), this sets a minimum halo mass and indirectly a minimum galaxy mass in the simulation. We do not expect to detect dwarf galaxies in the results presented here. We use the ΛCDM model as a proxy model that fits many observations. The FLRW cosmological parameter settings include the current values of the matter density parameter Ω m0 = 0.3, the dark energy parameter Ω Λ0 = 0.7 and the Hubble-Lemaître constant 0 = 70.0 km/s/Mpc.

Simulations
For our -body simulation, we chose --0.0-482f90f, a fork of the widely used adaptive mesh code -3.0 (Teyssier 2002). The fork has modifications to comply with the MPI 3.0 recommended standards for inclusion of the MPI header file 4 and optional extensions related to scalar averaging (Roukema 2018). The adaptive mesh structure of is designed to allow fast calculations that can resolve detailed gravitational behaviour in high density regions. The maximum resolution, which effectively corresponds to a softening length, is set at levelmax = 12. We produce snapshots starting at i = 10 Myr with an equally spaced time step of Δ = 100 Myr, and convert to scale factor values using -0.3.8.2. Newtonian gravitational potentials against the FLRW background are calculated in using the Poisson equation and a cloud-in-cell algorithm to estimate the matter density. These potentials are used to decide how to accelerate and shift particles. We output these values for later calculation of accelerations in relation to the elaphrocentre (Sect. 2.2).

Halo detection
For detecting dark matter haloes, we use -0.99.9-RC3+-6d16969 (Behroozi et al. 2013a), which uses a 7-dimensional friends-of-friends (FoF) algorithm. In contrast to other FoF halo finders, not only uses particles' spatial locations and peculiar velocity information to decide on physically meaningful dark matter haloes, but also uses temporal information to identify groups of particles that are persistent in time, rather than only using the instantaneous characteristics of a virialised object. Knebe et al. (2011) found that performs excellently in recovering haloes from an -body simulation. We run using a linking length of 0.28 and a minimum of 5 particles per halo. We set the virial radius criterion for detection to "200c", i.e., 200 times the critical density.

Merger history trees
We construct halo merger trees from the simulations Roukema 1993;Roukema et al. 1997;Roukema & Yoshii 1993, and references thereof). In this paper, instead of using the original Fortran77 routines from 1992, we use a more modern package, -1.01-e49cbf0 (Behroozi et al. 2013b), which was designed to perform on outputs from . In order to use these merger history trees for simulating galaxy evolution using , which was developed for following up simulations such as the Millenium simulation, we need to convert the trees to the LH T format. We do the conversion with -0.0-522dac5.

Galaxy formation and matter infall
To form galaxies within our dark matter haloes, we use semi-analytical galaxy formation recipes Kauffmann et al. 1993;Roukema et al. 1997;Kauffmann et al. 1999). Again, rather than using the original code from 1992, we use (Croton et al. 2016). For an introductory review, see Baugh (2006). We use the built-in functions of for estimating the 'size' and the infall rate history of each galaxy at the final output time. All galaxies are assumed to form as disk galaxies in , with a disk radius of where vir is the virial radius (set in at "200c") and is the dimensionless spin parameter computed using 's halo properties. The introduction of the parameter is generally attributed to Peebles (1969, eqs (35), (37)). The parameter is now widely used (e.g., Mo et al. 1998, eq. (10)), and is often normalised by a √ 2 factor (Bullock et al. 2001, eq. (5)), which is the convention chosen in , and adopted in this work: where is the total angular momentum (calculated as the vector sum of individual particles' angular momenta by functions including is the total energy in the non-relativistic sense (calculated as the sum of individual particles' Newtonian kinetic and potential energies by functions including _ _ in . in ), is the Newtonian gravitational constant, and is the halo mass.
The infalling mass is defined at a given time step by where reion is the reionisation factor, which estimates the effect of ionisation of the intergalactic medium from early stars; the baryon fraction b = 0.20 determines what fraction of total matter is baryonic (assumed to be the same for any dark matter halo/galaxy pair); vir is the virial (total matter) mass of the halo; and tot is the sum of all reservoirs of baryonic matter from the previous timestep (except at the initial timestep, when it is zero). Cases where infall < 0 are interpreted to mean that baryonic mass is ejected from the galaxy. We extended in order to estimate infall rate, star formation rate (SFR) and outflow rate histories. The history of any of these parameters for a given galaxy at a given time, traced backwards in cosmological time, is assumed to be the sum of the histories of all the separate pre-merger progenitors of the galaxy, appropriately matched by cosmological time. At each merger event, this physically corresponds to the components (dark matter, hot gas, cold gas, stars) of the progenitors being conserved in the merger. The summed star formation rate traced back for a given galaxy is what was originally used together with evolutionary stellar population synthesis to calculate galaxy spectral energy distributions and absolute magnitudes (Roukema et al. , 1997; we do not carry out evolutionary stellar population synthesis in this work. To evaluate these sums, we identify all galaxies present at the present epoch, ( ) = 1, and trace their progenitors' history back in time along the merger tree. For example, if a progenitor is itself the result of a merger at an earlier timestep, then its own history is the sum of its own progenitors. This procedure is continued recursively back in time in the merger tree for a given present-epoch galaxy. As in the original implementation ), other effects from mergers than conservation of mass, such as merger-induced starbursts, are also assumed in , but with a power law dependence on the satellite mass rather than direct proportionality (Croton et al. 2016, eq. (27)).

Voids
We identify the void environment of galaxies using the watershed void finder based on (Neyrinck 2008;Nadathur et al. 2019), which provides a nearly parameter-free void finder that does not require assumptions about void shapes. In contrast to other works, we use the full dark matter (DM) particle distribution as tracers. Nadathur & Hotchkiss (2015a) showed that the voids identified using the galaxy distribution as tracers differ from those traced using a randomly subsampled particle distribution. The authors recommend using the galaxy distribution as tracers in order to match observations. However, since our priority here is gravitational effects, we use the DM particle distribution rather than the galaxy distribution. It is likely that we will detect voids that are smaller and more numerous than those observed in the galaxy distribution (e.g. Mao et al. 2017), since the full DM particle distribution will show positive fluctuations in the density field that may be too weak to form DM haloes and galaxies, but will be detected by the watershed void finder and interpreted as boundaries of voids.
We introduce several small changes into . In addition to Nadathur & Hotchkiss (2015b)'s definition of the circumcentre, we calculate the position of the elaphrocentre, as defined below in Sect. 2.2. We add a routine to read in simulation data in G -2 format (Springel 2005), the default output format that we chose for the -body simulation. We output lists of particle identities of the DM particles that constitute each void. This information is needed to decide the extent to which a galaxy's host halo is located in a void.
We adopt the definition of the effective radius of a void, eff , which is not directly related to the choice of a definition of the void centre. The effective radius eff is defined by as the radius of a hypothetical sphere that has the same volume as the total volume of all the Voronoi cells that constitute the void, i.e. eff := 3 4 ( ) 1/3 where are the volumes of the Voronoi cells that determine the void.

Elaphrocentre and other definitions of void centres
To investigate if a galaxy's position in a void -its elaphrocentric location -has a significant effect on the formation and evolution of the galaxy, we first need to clarify earlier terminology regarding void centres from the literature, and we need to define the elaphrocentre.
Nadathur & Hotchkiss (2015b, Sect. 2.3) define the "circumcentre" for a given void using the Voronoi cell with the lowest density and the three lowest density adjacent Voronoi cells. The intersection of these four Voronoi cells determines the circumcentre. By construction, the circumcentre is the centre of the largest sphere that can be inscribed in the tetrahedron determined by the particles in these four (neighbouring) Voronoi cells, and the centre of the largest empty sphere that can be inscribed in the void. Nadathur et al. (2017, Sect. 2.3 (ii)) rename this the "void centre" and show that it correlates strongly with the local maxima of the gravitational potential with respect to the background FLRW model.
We define the elaphrocentre similarly. In a void identified by the watershed algorithm, we identify the particle at which the potential is highest, using the potentials estimated by (Sect. 2.1.2). We then identify the three adjacent Voronoi cells whose particles have the highest potentials. Together with the cell of the highest potential particle, we again form a tetrahedron between the four particles that respectively define the four cells. The centre of this tetrahedron is the elaphrocentre. While we expect a strong spatial correlation between elaphrocentres and circumcentres, they will differ in general, in particular for small, highly non-spherical voids. By definition, elaphrocentres are appropriate for studying elaphrocentric effects on galaxy formation. A group of test particles at an elaphrocentre will, in the Newtonian sense, be accelerated away from the elaphrocentre, and disperse rather than cluster together. Thus, the elaphrocentre would seem to be a good environment to form a large, diffuse galaxy, provided that the mass that forms the future galaxy is low compared to the mass deficit determining the gravitational properties of the void as a whole.
A third centre commonly defined in void studies is the 'macrocentre' or 'volume-weighted barycentre'. This is defined (Sutter et al. 2015, Sect. 3, eq. (4); Nadathur & Hotchkiss 2015b, Sect. 2.3, eq. (2)) as the volume-weighted mean of the position vectors ì of all DM particles identified as being in the void, i.e., ì vwb := ( ì ) / , where is the Voronoi cell volume associated with the -th particle. Since the position is not mass-weighted, it is unrelated to the normal Newtonian definition of a barycentre for a particle distribution, ì m := ( ì ) / , where is the mass of the -th particle. In the continuous limit to arbitrarily high particle resolution (assuming a continuous fluid), the volume-weighted barycentre ap- Clearly, in the continuous limit, ì vwb is the geometrical centroid of the domain D, which in geometry is often termed the 'barycentre'. In general, this corresponds to the astronomical barycentre only if the density distribution in the domain D is uniform. In other words, the volume-weighted barycentre contains no information about the density distribution within D apart from numerical noise.
By definition, apart from discretisation and numerical effects, ì vwb only determines the geometrical mean position (the centroid) of the overall shape of the void, as defined by the outermost particles of the void. Adding a few particles with big Voronoi cells adjacent to a single side of a void would significantly modify the global shape (union of Voronoi cells) of the void. This would shift the volume-weighted barycentre -the centroid -significantly. Thus, Nadathur & Hotchkiss (2015b) are correct that ì vwb depends on the presence rather than the absence of tracers. However, the fundamental problem with using ì vwb in the context of cosmological voids is that it indicates a void centre that (apart from numerical effects) has no dependence on the density variations in the interior of the polyhedron (union of all Voronoi cells) that defines the void; only the void boundary affects ì vwb . Since there is no reason for the centroid of the polyhedron bounding a void to have any tight relation with the position of minimal density if the void is even mildly asymmetrical, it is unsurprising that Nadathur & Hotchkiss (2015b) found the density at ì vwb to be higher than at the circumcentre. Thus, in the context of cosmological voids, the physical relevance of ì vwb is unclear, and if used, we recommend that it be described by the term 'geometrical centroid' (or 'boundary centroid') rather than 'macrocentre' or 'volume-weighted barycentre'.

Analysis
To study our hypothesis that the elaphrocentric location of galaxies in voids plays a significant role in their evolution, we analyse the simulated haloes, galaxies and voids produced by our pipeline in relation to the voids' elaphrocentres as follows.

Void membership and elaphrocentric distance
Identifying which galaxies are located in a void is non-trivial. For a given galaxy, we could find the void that gives the shortest elaphrocentric or circumcentric distance, and consider the galaxy to be a member of the void if the circumcentric or elaphrocentric distance is below a given fraction of the void effective radius. A distance between two positions in this work is calculated using the shortest of the multiple 3-torus (R 3 /Z × Z × Z ≡ S 1 × S 1 × S 1 ) spatial geodesic comoving distances. This is often described more loosely as 'the comoving distance with periodic boundary conditions'.
However, identifying galaxy membership in a void by the elaphrocentric or circumcentric distance would only be accurate for voids that are spherically symmetric. Although voids tend to evolve to become more spherical, as shown analytically by Icke (1984) by reverting a simple toy model for collapsing density perturbations and numerically by Sheth & Van de Weygaert (2004) in -body simulations, a void will in general be non-spherical. Moreover, the elaphrocentre will not, in general, coincide exactly with the circumcentre. Thus, a more accurate way of deciding on void membership should, in principle, be possible by using knowledge of the particle positions.
The void membership criterion proposed here, as with the > 50% merging identity criterion initially published in 1993 for merger history trees (Roukema & Yoshii 1993, Sect. 3;Roukema et al. 1997, Sect. 2.2.1), is a simple proposal that we expect to be improved upon later. Our voids are detected as a union of Voronoi cells -each containing a DM particle -by . Thus, for any given void, we have a list of particles that approximately define the void. Since we know from which particles are members of a halo in which a given galaxy forms, we can check which of these halo particles are present in the list of void member particles for any given void. We restrict the list of void particles to those that are within ≤ 2.0 eff from the elaphrocentre, where eff is the void effective radius calculated by (the radius of the hypothetical sphere having the same volume as the void's Voronoi cells; see Sect. 2.1.6). The ≤ 2.0 eff restriction should remove some of the sharpest regions adjacent to the knots of the cosmic web and exclude the outermost regions of voids of high ellipticity. In this sense, it will counteract the space-filling nature inherent to any watershed voidfinder to some degree.
Our void membership criterion is that we require that a fraction H∩V strictly greater than min H∩V = 0.50 of the particles in a halo H be members of a void V for a galaxy in that halo to be considered a member of the void. As in the case of the > 50% merging identity criterion, which prevents a halo from having multiple descendants, this void membership criterion is strong enough to prevent a galaxy in a halo that lies on a wall, filament or knot from being allocated to more than one void. This criterion could be strengthened to force a selection of galaxies that are placed further in the interior of the voids, at the cost of reducing the total number of galaxies recognised as being members of voids. A galaxy that does not satisfy this criterion is considered to be a non-void galaxy.

Infall dependence on environment
We wish to see if infall rates -of dark and baryonic matter in general -are affected by the host halo's location in a void. The infall history should affect the star formation rate, which requires baryonic matter to first collapse into the centre of its host dark matter halo's potential well. We consider the infall rate traced backwards in time for any halo at the final output time. This infall rate is the sum of the mass accumulation histories of the component haloes of the final halo's merger tree. We can write this backtraced history, over predecessor haloes destined to merge together, as the mass evolution assigned to the final halo, ( ), so that d /d ≈ Δ /Δ ( ) is the infall rateof small haloes and diffuse matter together.
The hypothesis that the elaphrocentre of a void (corresponding to a spatially compact part of a hyperbolic, super-Friedmannian expanding region) would weaken the infall rate can now be formalised. For a given mass , the average (mean) infall rate is, by definition, independent of location. To distinguish the archetypal case of a typical disk galaxy, with an initial burst of star formation followed by an exponential decay, from that of an archetypal LSBG, with an approximately flat star formation rate, we attempt to fit d /d by an exponential, of the form where is the infall amplitude and is a decay time scale. Typical disk galaxies should have low time scales , while we hypothesise that void galaxies should on average have higher , corresponding to approximately flat infall rates.
The exponential form of the fit for the matter infall rate is a heuristic choice inspired by the traditionally used declining exponential for modelling the star formation rate (Bruzual A. 1983). In reality, merger histories are more complex than these simplified extremes, so we expect a fair fraction of automated fits to fail, especially since we apply this to all timesteps defined in Sect. 2.1, starting from the first time step in which a galaxy is modelled to form in a dark matter halo. Nevertheless, this automated fit procedure should be able to distinguish whether the infall of matter is closer to a brief, quickly weakening series of early events or rather a more steady infall extended over a long period. We expect that a flat infall rate would correspond to a roughly constant star formation rate over time. The early, brief, burst scenario of the infall of matter would allow a high star formation rate immediately after the infall, whereas a constant rate of matter infall should yield an approximately constant star formation rate. Long cooling times would modify this relation, especially for galaxies in high-mass haloes in voids (Hoffman et al. 1992, Einstein-de Sitter case).
We first compare and for galaxies depending on their classification as void or non-void galaxies. The infall histories are calculated using our modification of . The prediction for creating LSBGs is that void galaxies should tend to have low amplitude and a high (slow) decay rate , and vice versa for non-void galaxies.

Galaxy size dependence on environment
The other parameter that may indicate an elaphrocentric contribution to a disk galaxy becoming an LSBG is the disk scale length disk (for a density profile ∝ exp(− / disk )). This can be converted to a disk half-mass radius using the relation 1/2 = disk , where solves ( + 1) exp(− ) = 1/2 (e.g. Kravtsov 2013, = 1.678). As stated above (Sect. 2.1.5; Eq. (1)), calculates disk . Since voids are underdensities, dark matter haloes forming in voids will tend to collapse somewhat later than in overdensities (e.g. Lacey & Cole 1993, App. A). In an expanding FLRW universe, the critical density crit decreases, so for a fixed virialisation overdensity threshold and fixed mass, vir increases with time. Thus, modelling galaxy disk scale lengths as being proportional to the halo virial radius (Eq. (1)), it would be reasonable to expect void galaxies to have greater disk than non-void galaxies, for a fixed value of the spin parameter . Elaphrocentric galaxies will typically undergo a more isolated evolution than barycentric galaxies, with fewer merger events. D'Onghia  6), using the full dark matter particle distribution, which leads to voids much smaller than typically observed or found in simulations when traced by galaxies. Histogram of elaphrocentric distance of galaxies identified as being located in voids. The distribution at / eff < 1 is typical of observational and simulated void profiles, with most galaxies located near the effective radius. The distribution at / eff ≥ 1 can be interpreted as showing galaxies in the outer parts of voids that are generally quite asymmetrical. A perfectly symmetrical void would have all its member galaxies at / eff < 1.
(2008) found, based on -body simulations, that the spin parameter of haloes in equilibrium is not influenced by merger events. If the role of the spin parameter is indeed weak, then void galaxies should be marginally larger than non-void galaxies. We consider both disk directly, and vir and individually.

Elaphro-acceleration
As a complement to the direct analyses of simulated galaxies via , we also investigate accelerations near the elaphrocentres, as preparation for future studies of elaphrocentric effects that might have an effect on galaxy formation in voids and might help form giant LSBGs. We estimate the acceleration (compared to the FLRW reference model) of test particles directed away from the elaphrocentre, in the direction of the boundaries of the void. This acceleration can be thought of as counteracting the self-gravity of an overdensity that is destined to collapse into a dark matter halo and allow the formation of a galaxy within the halo. This acceleration is the effective antigravity that we have assumed could enlarge the size of a galaxy and decrease its infall rate, provided that the simulation has a detectable galaxy near the elaphrocentre.
Here, we describe how we estimate this 'elaphro-acceleration' without requiring the presence of a halo or galaxy. We will compare our results with the Newtownian estimate for the gravitational pull towards the centre of a halo of a Malin-1-like galaxy (Bothun et al. 1987). In this simplified model, we assume a high-mass test halo at the elaphrocentre of a void identified in the simulation, without modifying the underlying DM distribution. In the real Universe and in simulations, it is rather unlikely for a galaxy to form exactly at an elaphrocentre. Nevertheless, we feel that this calculation will be a useful guide, since the elaphro-acceleration should be maximal in amplitude at the elaphrocentre. For our canonical high-mass test halo we adopt parameters that are motivated by observations of Malin 1 (Seigar 2008;Junais et al. 2020). For an order of magnitude estimate, we adopt test = 10 12 M ⊙ for the mass of the halo and test = 1.20 Mpc for the region from which dark matter originated.
For any given void, we interpolate the potential linearly, and calculate the acceleration as the gradient of the potential, ∝ −∇ . We use the full dark matter particle distribution, since gravitationally this should be more accurate than that of collapsed haloes (or galaxies) alone. We use the gravitational potential estimates calculated by (Sect. 2.1.2). Since only provides potential estimates at particle positions, the resolution limit implied by using these is determined by the particle number density. By definition, the number density is very low inside a void, so there are very few particle positions available for interpolating the potential. We sample the elaphro-acceleration at six positions which lie on the sphere with radius test , i.e. at ((± test , 0, 0), (0, ± test , 0), (0, 0, ± test ), where the elaphrocentre is the origin of the coordinate system. For a given void, we calculate the radial (signed) and tangential (amplitude) components of the velocity vector ì for each of the six positions, and find the mean values and ⊥ , respectively. As stated in Sect. 2.1.6, some of the six positions may fall outside of a void when the void is too small; we ignore the void in such cases.

Reproducibility versus cosmic variance
We present results below (Sect. 3) with a preference for reproducibility over cosmic variance. Our pipeline, using the Maneage template (Akhlaghi et al. 2021), includes a step for verification, in the sense of verifying that when the reader recalculates our complete pipeline, s/he should obtain statistically equivalent results to our original results, within some tolerances. Estimation of these tolerances effectively requires an approximate estimate of cosmic variance in the parameters of interest. In principle, we could use these repeated full runs of the pipeline to obtain mean or median estimates of our main parameters of interest, rather than presenting the values from a single simulation, and the random uncertainties derived from that simulation. However, that would reduce the reproducibility of our results, in the sense that readers wishing to run our pipeline would also have to perform multiple full runs.
Thus, here we favour reproducibility over cosmic variance, using the latter only for reproducibility purposes, not for obtaining results of the paper. We first run a fixed version of our full package 10 times, for different random seeds and with physical randomisation induced by parallel computation. We obtain a full set of the results presented in this paper from each run. For any given parameter, we calculate the standard deviation cv of the values of the parameter, where the subscript 'CV' refers to 'cosmic variance'. The particular version of the source package used for these verification runs has the string e93569c as its commit identity (commit hash). We use . Differential halo number counts 2 / d /d versus halo mass for the haloes that host galaxies in the simulation, where is the number density of haloes in a given mass interval and is the mean mass density of the simulation. The solid, black (top) curve represents all haloes; the dashed, red (middle) curve is for haloes that are not in a void; and the dash-dotted, blue (bottom) curve is for haloes that are identified as being located in a void. The gradual decline in the numbers of haloes towards the lower mass scales and the sharp cut at the lowest mass scale are artefacts of the limited resolution of our simulation. these repeat runs only for verification in the reproducibility sense, not for results. A given parameter in a fresh realisation that aims to reproduce our results can then be verified for consistency by requiring that it agree with our published value to within 5.0 √ 2 times the stated random error ran , where the √ 2 factor represents the assumption that both have independent Gaussian errors drawn from the same distribution. For parameters with high cv (i.e., those known to fail this verification), we require agreement between the fresh value and the mean from the ensemble of runs within 5.0 times cv , and we state cv as an additional uncertainty, using the notation '± cv '. We use this formalised verification procedure for the parameters that we judge to be the more physically relevant.
The simulation presented below was chosen randomly.

Simulation pipeline
In the final time step of our = 128 3 -particle simulation, we detected 5329 haloes, with 4817 galaxies evolved along the merger history trees for these haloes. Among these galaxies, 3848 have virial mass in the range 10 11 -10 13 M ⊙ /ℎ, which we study with the aim of seeking a factor in the formation of high-mass galaxies. Applying the H∩V > 0.50% definition in Sect. 2.3.1, we identify 1998 galaxies in voids, among which 1588 of these galaxies have a mass in the selected range. This fraction of galaxies identified as void galaxies is larger than the 7% estimated by Pan et al. (2012). A more detailed analysis of the galaxies that are identified to be in voids is given below. We investigate key quantities dependent on the fraction of their host haloes particles in a void H∩V and dependent on their relative distance to the voids centre / eff . We use Theil-Sen robust linear fits (Theil 1950;Sen 1968) on each key quantity to see if it has a statistically significant dependence on either of the two void location parameters. There are 2198 voids in the final time step. The void size distribution is shown in Fig. 1. Since we detect voids physically, following the dark matter particular distribution, as described in Sect. 2.1.6, it is unsurprising that the void population is dominated by voids a few Mpc/ℎ in size, with a correspondingly higher number density than that typically seen in void catalogues calculated using galaxies as tracers. The distribution of void galaxy elaphrocentric distances is shown in Fig. 2.

Infall rate
As described in Sect. 2.3.2, we first compared infall rates for galaxies as separate void and non-void populations. For each galaxy mass infall history d /d ( ), we first find a linear least-squares best fit to log 10 (d /d ) versus for time steps where d /d ( ) > 0. The optimal parameters of this fit are used to find a non-linear least-squares best fit of d /d ( ) to a decaying exponential (Eq. (4)  time steps with d /d = 0. As stated above, merger histories are complex, and many galaxies' infall histories are poorly fit by this procedure. This applies both to void and non-void galaxies, and we do not attempt to analyse these cases. The limitations of this simplified approach should similarly affect both populations and should not affect our comparison of the successfully fit subsets of the two populations. We find v = 1263 valid fits for the void galaxies, and − v = 325 fits that are rejected either as failed fits ( − vf = 1) or as having an unreasonably high amplitude, > 1000 ⊙ /yr, indicating a physically unrealistic fit ( − vu = 324 cases). For non-void galaxies we find v = 1877 valid fits and − v = 525 fits rejected either as failed fits ( − vf = 1) or as having unphysically high amplitudes ( − vu = 524 cases). For galaxies in host haloes with virial masses in the range 10 11 -10 13 M ⊙ /ℎ, we find the medians listed in Table 1, where the uncertainties are standard errors in the median. Throughout this work, uncertainties in the median are given as the standard error in the median, unless otherwise stated. The median host halo mass for the void galaxies is (4.1± 0.4) × 10 11 ⊙ , lower than that of the non-void galaxies, (5.7 ± 2.6) × 10 11 ⊙ . Figure 3 shows the differential mass function of the haloes in the form of the halo multiplicity function. We show the differential halo masses for all haloes, for haloes that are classified as being associated with a void and for haloes that not associated with any void. As stated in Sect. 2.1.1, we cannot resolve dwarf galaxies; this is seen in the sharp cut at the lowest masses. The gradual decline in low-mass haloes and the sharp cut are both consistent with the use of the halo finder (Sect. 2.1.3). The void and non-void halo mass functions clearly differ in the higher mass ranges. The most massive haloes form in voids very rarely, while among the lowest mass haloes, the probabilities of forming in a void or not are of similar magnitude.
We find (Table 1) no significant difference in either the amplitude or the time scale of infall between the void and non-void galaxy populations. The dispersion in infall patterns within each population is as great as both and themselves. We found that, to very high significance, void galaxies typically form later than non-void galaxies, as expected, since they form in underdensities. We find median collapse epochs (in standard FLRW cosmological time) of f v = 4.1 ± 0.1 Gyr and f nv = 3.3 ± 0.1 Gyr for the void and non-void galaxies, respectively. By the collapse epoch of a galaxy, we mean the first epoch at which the mass infall rate calculated by for the galaxy is non-zero.
We investigate the results for and more closely by checking if either or has a dependence on either the fraction of a host halo's particles that are identified as void particles, H∩V , or on the host halo's elaphrocentric distance / eff . Figures 4-7 show the dependence of and on H∩V and / eff for all void galaxies. The H∩V axis is shown with H∩V decreasing from left to right, so that the galaxies that are best qualified as void galaxies are shown towards the left in all four figures. There is no visually obvious dependence of the infall parameters on H∩V . However, Fig. 4 does show a modestly significant non-zero slope, i.e., the median of the infall amplitudes is somewhat higher for galaxies better identified as void galaxies (having a higher value of H∩V ). This would tend to oppose the hypothesis of a general tendency to form LSBGs in voids. The other slopes of the best fit linear relations, using robust

Galaxy Sizes
While we do not detect significant elaphrocentric effects on infall rates, effects on galaxy sizes could play an important role in forming large diffuse galaxies. As stated above (Sect. 2.3.3), to check the size of a galaxy at the final output time step, we use the disk scale length disk provided by . The results for the galaxies divided into void and non-void populations are shown in Table 2, where we list the disk exponential scale length disk , the spin parameter and the virial radius vir . Table 2 shows a significant difference for the overall scale length disk . Our results show that, as a population, void galaxies form significantly smaller disks, both for our selected mass interval and for the full sample. Although this might seem to support van de Weygaert et al. (2011, fig . 2, left)'s finding that, for a given absolute magnitude, void galaxies tend to be smaller than the general galaxy population, we do not (yet) model stellar populations and estimate absolute magnitudes, so this qualitative agreement is promising but not conclusive.
A likely explanation for the smaller sizes of void galaxies is shown by the other rows in Table 2: the void galaxy population has a much smaller median virial radius vir than the non-void population, but an insignificantly higher spin parameter . The slightly greater spin parameter appears insufficient to compensate or override the lower vir of the void galaxies. The values of the spin parameter are reasonable in relation to standard values in the literature. Our non-void host halo values of listed in Table 2 are consistent with the friendsof-friends halo estimate of Zjupa & Springel (2017, Sect. 2.4, para. 8), B,FOF = 0.0414, while our void host halo values of are slightly higher.
As stated above, the void galaxy host haloes are typically somewhat less massive than the non-void haloes and the collapse epochs of void galaxies are significantly later. These two parameters should have opposite effects on the halo virial radii. In Table 2, we see that vir is significantly larger for non-void galaxies, showing that the higher mass of non-void galaxies plays the dominant role.
To see if a general trend of disk also exists as a function of a galaxy's void location, Figs 8 and 9 examine the dependence of disk on H∩V and elaphrocentric distance for void galaxies. The slope of the best fit in Fig. 8, d disk /d H∩V = −1.44 ± 0.46 ± cv 0.42 kpc/ℎ, is not significantly non-zero when we take into account cosmic variance (the distribution of d disk /d H∩V over repeated runs includes a strong tail of values that are not significantly nonzero). The dependence on elaphrocentric position (Fig. 9) is not significant either. In Figs 10-13, we investigate whether the spin parameter or the virial radius vir is more responsible for the modest reduction in the disk scale length of void galaxies as indicated in Table 2. The sharp lower limit in Figs 12 and 13 is an artefact of the detection threshold of dark matter haloes in the -body simulation. The virial radius is calculated by from the halo mass. Out of the four figures (Figs 10-13), two show highly significant slopes: Figs 10 and 12. The slopes in Figs 11 and 13 are not significantly different from zero.
The slopes in Figs 10 and 12, d /d H∩V = 0.021 ± 0.004 ± cv 0.005 and d vir /d H∩V = −100 ± 9 kpc/ℎ, respectively, are both very strong, but opposed. Galaxies better identified as void galaxies have higher spins, but also lower virial radii and lower masses. The overall effect, as shown in Fig. 8, is that the two effects more or less cancel, in contrast to the full-population results shown in Table 2. Figure 9 and its fit show that overall, the elaphrocentric distance Table 3. Median radial and tangential ⊥ accelerations at test = 1.20 Mpc/ℎ from the elaphrocentre. median acceleration 0.09 ± 0.01 km/s/Gyr ⊥ 1.90 ± 0.04 ± cv 0.29 km/s/Gyr / eff has only a weak effect on galaxy disk scale lengths disk . The halo size and spin parameter both have weak, though apparently again opposite, dependences on disk , with d /d ( / eff ) = 0.007 ± 0.002 ± cv 0.002 in Fig. 11 and d vir /d ( / eff ) = −15.8 ± 5.9 ± cv 6.7 kpc/ℎ in Fig. 13.
In summary, the trends for disk , , and vir found in Table 2 are similar, but strengthened, when void location of a galaxy is quantified by H∩V , and insignificant when void location is quantified by / eff . The lack of a significant dependence of these parameters on the elaphrocentric distance, / eff , is somewhat surprising, since one might expect H∩V and / eff to be proxies for one another, equally valid for defining how high a galaxy is on the potential hill of a void. We discuss this counterintuitive result further in Sect. 4.2.
Since vir is obtained from vir by on the assumption of a detection threshold of 200 times the critical density, Fig. 12 can be qualitatively compared with observational estimates of the masses of void galaxies. Keeping in mind the fixed lower limit in mass resolution, the robust best fit relation can be used to describe the galaxies best located in a void ( H∩V = 1) as having host haloes with vir ∼ 100 kpc/ℎ, and those best located in walls ( H∩V = 0) as having higher mass host haloes, with vir ∼ 200 kpc/ℎ. Thus, the masses of galaxies' host haloes located in the walls should be typically eight times those of galaxies located in voids. Weistrop et al. (1995) found that 12 H -emitting galaxies in the Boötes void were mostly quite luminous, and Szomoru et al. (1996) found that a sample of 16 targetted Boötes void galaxies (along with 21 companion galaxies) appeared to be similar to corresponding late-type, gasrich field galaxies and of similar masses. These Boötes void surveys would appear to be inconsistent with the mass difference found here. However, the more recent and bigger survey of 60 void galaxies in the Sloan Digital Sky Survey (SDSS) by Kreckel et al. (2012) found that these have moderately low stellar masses, mostly around 10 9 -10 10 ⊙ . The SDSS void galaxy survey would appear more likely to be consistent with our results. Future work, in which the remaining steps in galaxy formation and evolution modelling are added to the pipeline presented here, should enable quantitative comparisons to see if the Boötes and SDSS void galaxies (see also Pan et al. 2012;Nadathur & Hotchkiss 2015b;Sutter et al. 2015) are consistent with those modelled here.

Elaphro-acceleration
Elaphro-accelerations were calculated as described in Sect. 2.3.4, at elaphrocentric distances of test = 1.20 Mpc. We found 1570 voids that allowed valid estimates. The median radial acceleration for a test particle at these positions is given in Table 3. The amplitudes of these two values are not directly comparable, because is the median of signed values, while ⊥ is non-negative by construction. The Newtonian estimate for the gravitational pull at test away from the centre of the canonical high-mass test halo of mass test , a barycentre, is an inward-pointing acceleration, i.e. test-halo = −3.05 km/s/Gyr.
The median estimate of given in Table 3 is an outward-pointing elaphro-acceleration to high statistical significance. This is consistent with what could be expected of the elaphrocentre, defined as the location of the global maximum in the potential of a void, with mass typically moving away from the elaphrocentre. The amplitude is more than an order of magnitude weaker than that of the barycentric acceleration towards our canonical high-mass halo. Figure 14 shows that the full spread of radial elaphro-accelerations is wide, including many negative values, and that dependence on the void effective radius eff is weak. Together, these properties imply that, at least with the numerical techniques and simulation parameters adopted in this work, a systematic antigravitational effect at the elaphrocentre helping to oppose infall is likely to be modest. This is consistent with our infall results above. The median tangential acceleration ⊥ is given in Table 3 and the individual estimates and fit are shown in Fig. 15. These values are about an order of magnitude higher in amplitude than those of , and similar to that for our canonical high-mass halo. This supports the argument that rotational properties of the fluid flow such as shear and vorticity are likely to be significant in understanding voids. Figures 14  and 15 do not show significant dependence of and ⊥ on the size of a void.

Infall rates
We generally found a lack of statistically significant trends in the two infall parameters on H∩V and / eff ( Table 1, Figs 4-7), for those galaxies whose host haloes' infall rates could be fit with an exponential best fit. A moderately significant non-zero dependence is that of on H∩V , shown in Figure 4, opposite to that expected: void galaxies have slightly higher amplitudes of their best fit infall rate histories. The simplest interpretation is that in a halo destined to collapse with a given final mass, the infall history of matter into that halo is nearly independent of the environment. The modest amplitude of the median acceleration outwards from the elaphrocentre -the elaphro-acceleration (Table 3) -is probably the main reason for this. This is something like a Newtonian numerical equivalent of Birkhoff-like (Birkhoff & Langer 1923) or 'finite infinity' (Ellis 1984;Wiltshire 2007) arguments for modelling galaxies in isolation from their environment. Our hypothesis that the void environment helps to form giant LSBGs by providing slow, weak infall is not supported by our numerical results.
Although the parameters of our simplified infall fit are not affected by the position of a galaxy, we found a significant difference in the median galaxy collapse epoch (the first time step with a non-zero, non-negative infall rate) between void galaxies and nonvoid galaxies. While we did not expect this to play a major role in galaxy formation, it should. The median collapse epochs (in standard FLRW cosmological time) that we found were f v = 4.1 ± 0.1 Gyr and f nv = 3.3 ± 0.1 Gyr for the void and non-void galaxies, respectively (Sect. 3.2). Thus, to very high significance, void galaxies typically form later than non-void galaxies. We interpret this as a result of their formation in underdensities. Haloes that collapse later should, according to the standard spherical collapse model, have a lower matter density. Galaxies should thus form with lower dark matter and baryonic matter densities (M ⊙ /kpc 3 ), which may lead to lower surface densities (M ⊙ /kpc 2 ) of galaxy disks and lower surface brightnesses (L ⊙ /kpc 2 ) induced by star formation.
This characteristic of void galaxies agrees with Rong et al. (2017), in the sense that these authors found that UDGs have a later formation time than typical dwarfs, and assuming that we associate UDGs as being located in voids. Rong et al. (2017) found that UDGs at the current epoch have a median age of 7.1 Gyr compared to typical dwarfs with a median age of 9.6 Gyr. For the current epoch estimated at 0 ∼ 13.8 Gyr, these correspond to median UDG and ordinary dwarf formation epochs of 6.7 Gyr and 4.2 Gyr, respectively. Since (i) we consider high-mass host haloes and correspondingly high-mass galaxies, rather than the more typical UDGs and dwarfs, and (ii) we separate populations by their location in voids rather than by continuing through to stellar population synthesis, more precise correspondence with Rong et al. (2017)'s results would be unlikely with our current pipeline. The qualitative agreement that void galaxies typically form later than non-void galaxies by about a Gigayear (our result) and that UDGs form about three Gigayears later than ordinary dwarf galaxies (Rong et al. 2017) is a promising sign of progress towards a cohesive theory of LSBG formation.

Galaxy sizes
The most massive galaxies and their host haloes (with the largest virial radii vir ) form in the tight knots of the cosmic web. Although we found that galaxies in voids tend to be smaller when comparing the overall void to non-void populations (Table 2) in the dependence of the galaxy disk scale length disk on H∩V , nor on on / eff . In contrast, we did find significant dependences of the two contributing parameters to disk , the spin parameter and the virial radius vir , on H∩V , but insignificant dependence on / eff . While the dependence on vir is clearly the dominant effect, we find that the spin parameter considered alone, which tends to form large galaxy disks, is higher for void galaxies. This is seen most significantly in Fig. 10, via the dependence of on H∩V . The higher spin parameter effect could be interpreted as the result of either fewer merger events weakening the spin parameter, or of gravitational effects inside the void. A likely candidate for the latter is the tangential acceleration ⊥ (Table 3). This is typically of the same order of magnitude as the gravitational pull of the source region of the highmass test halo that is, in its properties, inspired by Malin 1. However, in this work we have focussed on overall population properties and reproducibility of the pipeline. Continuation through to disk surface densities, for comparison with Di Paolo & Salucci (2020, Sect. 7.1, 7.2), and to stellar population evolution, remains a task for future work. Moreover, the rare high-mass galaxies that are well identified as void galaxies may require high numbers of simulations, if realised randomly, since, by definition, they are rare. Alternatively, a small number of big simulations may provide qualitative clues, as in the Malin 1 analogue found in the IllustrisTNG simulation by Zhu et al. (2018).
It may seem somewhat surprising that these significant dependences on H∩V do not translate into significant dependences on / eff . The explanation most likely lies in the fact that traces voids using Voronoi tessellation and the watershed algorithm, and voids are in general aspherical. Galaxies can lie in fairly empty parts of a void, with high H∩V , while lying, for example, in the far ends of a prolate void, where / eff > 1. It would be useful to quantify the relation between H∩V and / eff . Figure 16 shows visually that there is no obvious relation between H∩V and / eff . As indicated in the caption, the best robust fit indicates no statistically significant non-zero linear slope relating the two parameters. The fact that most of the void galaxies lie at / eff > ∼ 1 is consistent with the explanation suggested above. This can also be thought of as follows. Voids are defined by the absence of particles. Galaxies are generally not found in the interior of voids, because then the void shapes would be defined differently, shifting those galaxies' host haloes from void status to near-boundary status in the redefined voids. A halo located at approximately / eff > ∼ 1 in a highly aspherical void is not necessarily located in a locally low density region, so it is not constrained to contain a high number of particles identified as void particles.
Thus, H∩V and / eff appear to be statistically independent parameters. The significant trend of disk on H∩V , and the fact that H∩V has a more local physical meaning than / eff , suggest that H∩V is the more physically useful parameter to choose. Numerical and observational studies that measure the local number density around a given density will tend to correspond to the use of H∩V as a parameter for characterising the void nature of a galaxy.
To see if earlier definitions of void centres, as discussed in Sect. 2.2, could have more significant effects on galaxy sizes, we repeated our calculations for distances from the circumcentre and from the geometrical centroid instead of the elaphrocentre. Table 4 shows the robust best fits to the dependence of disk , vir and on void-centric distance for the three definitions. The full scatter plots (not shown) are visually indistinguishable from Figs 9, 11 and 13. The differences between the three cases are statistically negligible. Given that the geometrical centroid (macrocentre or volume-weighted barycentre) only encodes information about the void's periphery, with no information from its interior, it may seem surprising that this gives similar results to the other two centres. However, this is probably explained by the near-total absence of galaxies located within the central half-radius of the void; the detected galaxies' radial distances from the centre only vary mildly with different definitions of the centre.

Elaphro-acceleration
We found a significantly non-zero positive median acceleration towards the edges of a void. However, this median outwards acceleration, (Table 3), is of the order of only a few percent of the inward gravitational pull at test = 1.20 Mpc that the source mass excess for our canonical high-mass host halo would create, test = −3.05 km/s/Gyr. Moreover, Fig. 14 shows a wide scatter between outward and inward accelerations from the elaphrocentre. Thus, while a modest average effect in opposing infall could be expected for galaxies that are close to the elaphrocentre of a void, the effect would be sensitive to the wide distribution in values and to relations between the elaphrocentric acceleration and other dynamical parameters. Future work in placing a test halo near an elaphrocentre, with the assumption that the test halo has no dynamical effect on the underlying void properties, may use these results as a guide to judging the likely strength of the effect. For example, the probability that a test halo of a given mass placed at the elaphrocentre of a random void has its infall rate weakened sufficiently to make it a candidate LSBG could be estimated. This could be compared with the infall rate behaviour from haloes from the -body realisation itself, as we presented in Sect. 3.2, at elaphrocentric positions far from the elaphrocentre.
The median tangential acceleration ⊥ is much higher than the radial acceleration, and might be used to study the higher spin parameter found for void galaxies when identified by H∩V , as shown in Fig. 10. Since this is of the same order of magnitude as our canonical radial acceleration, test , it is likely that ⊥ could play an important role for galaxies forming in voids.

Future extensions
An obvious further development, not yet included in the present work, would be to analyse the star formation rate histories and to extend the pipeline with evolutionary stellar population synthesis methods. This would allow us to identify LSBGs in our galaxy population in a way more closely comparable to observational results, while continuing to benefit from the reproducibility and modularity of the pipeline presented in this work. The inclusion of metallicity evolution, in particular that of O/H, would allow comparison with the populations of extremely metal-poor gas-rich, dwarf galaxies that seem to characterise a difference between void and non-void galaxy formation (Pustilnik et al. 2020).
Another extension would be to extend or replace the gravitational simulation. Using a relativistic simulation, rather than a standard (Newtonian) simulation, would provide a major theoretical improvement towards more realistic results. The scalar averaging extensions provided by through the / front end to check background-independent dynamical properties (Roukema 2018), using the relativistic Zel'dovich approximation Buchert & Ostermann (2012); Buchert, Nayet & Wiegand (2013), could easily be added. Other options could include using either (Adamek et al. 2016) or the fully relativistic E T (Bentivegna & Bruni 2016;Macpherson, Lasky & Price 2017). Hydrodynamical simulations would also be useful for comparison. Given the aims of this project in providing a reproducible pipeline with modular, free-licensed components, it should, in principle, be straightforward to replace any of the pipeline steps or to start the pipeline at an intermediate step, such as analysing pre-calculated -body simulation outputs. The present form of the pipeline assumes 2 format for the -body simulation output snapshots. Alternatives in the statistical analysis of infall rate histories would also be useful to explore. Here, we chose to fit the infall rate history with decaying exponentials, which include nearly constant rates as a special case with very long time scales , but the reality of the mass infall rate history, and the corresponding star formation rate history, is in generally much more complex, depending especially on merger events. A more general quantitative way of characterising the global population of mass infall or star formation rate histories would bring this pipeline closer to physical reality.

CONCLUSION
We have presented a complete, ab initio, reproducible galaxy formation pipeline starting from a standard post-recombination-epoch spectrum of initial perturbations, aiming to identify key factors in void galaxy formation that might contribute to the formation of giant low surface brightness galaxies in voids (Sect. 2.1). We introduced the term elaphrocentre to clarify its opposite physical nature to the barycentre and we clarified the confusing use of the latter term in void studies (Sect. 2.2).
We did not find statistically significant numerical evidence that the elaphrocentre, or the void location of a galaxy more generally, plays a key role in forming major populations of large diffuse galaxies -LSBGs -via the parameters that we considered as the most likely to play a strong role -and (Figs 4-7). Since gravity is attractive, there is an asymmetry between the sharp nature of barycentres (potential wells) and the wide spread of elaphrocentres (potential hills), which could explain the lack of a strong effect.
We found that the fractional elaphrocentric distance of a void galaxy / eff is, statistically, a less useful independent variable than H∩V , the fraction of a galaxy's host halo particles that are identified as being in a single void. This is important for observational studies of void galaxies. The characterisation of galaxies as void galaxies by H∩V , which should roughly correspond to a low local dark matter density, or by relative elaphrocentric radius / eff , which would require identification of voids in a catalogue, will in general give uncorrelated results; the two parameters show no significant linear correlation (Fig. 16).
A serendipitous result is that void galaxies were found to be significantly smaller in virial radius (host halo mass) than non-void galaxies ( Table 2,  . This complicates the question of giant LSBG formation, because the disk scale length disk , as calculated in Eq. (2), is dominated by the virial radius. We did find that galaxies better identified in voids have a higher spin parameter. This finding of a higher spin parameter for high H∩V is qualitatively consistent with Rong et al. (2017)'s result that a higher spin is a key feature of UDGs and thus indirectly supports the hypothesis of void location constituting a significant factor in LSBG formation. Higher resolution simulations, extending to lower mass galaxies, would be needed to see if the higher spin of UDGs is quantitatively explained as a consequence of void location.
We also found that the median galaxy collapse epoch differs to very high statistical significance between void and non-void populations ( f v = 4.1 ± 0.1 Gyr and f nv = 3.3 ± 0.1 Gyr for void and non-void galaxies, respectively; Sect. 3.2). For a standard spherical collapse model, the later collapse of void galaxies should lead to these galaxies being less dense, quite likely resulting in lower surface densities and star formation rates.
In summary, despite not finding direct numerical evidence for LSBG formation in our overall populations, the higher spin parameter for the overall population of void galaxies, especially when characterised by H∩V , and the later formation epoch of void galaxies, are qualitatively consistent with Rong et al. (2017)'s findings for UDGs, assuming that the extension to lower masses remains valid. Together with these two key features that contribute to the formation of diffuse galaxies, the smaller size of void galaxies suggests that, in contrary to our hypothesis of giant LSBG formation in voids, the role of voids is to preferentially form diffuse, somewhat smaller galaxies. Moreover, we hope that by providing our complete software pipeline 5 using the Maneage template that aims at full reproducibility (Rougier et al. 2017;Akhlaghi et al. 2021), rather than only giving the names and URLs of cosmological software packages, our work will encourage the community to avoid unnecessary effort spent in guessing the precise details of the computational software used in this and other extragalactic research.

DATA AVAILABILITY STATEMENT
As stated in the introduction, the full reproducibility package for this paper, including all input data (parameters), is available at zenodo.4699702 and in live 6 and archived 7 repositories. The numerical output data for our main results is available as a DOIidentified record at zenodo.4699702/voidgals_infall.dat. The specific version of the source package used to produce this paper can be identified by its commit hash 97bd6dc.