\NeedsTeXFormat{LaTeX2e}
\documentclass[usenatbib]{mnras} % mnras.cls is available in e.g. Debian stable (9.10+)

%% This is a convenience variable if you are using PGFPlots to build plots
%% within LaTeX. If you want to import PDF files for figures directly, you
%% can use the standard `\includegraphics' command. See the definition of
%% `\includetikz' in `tex/preamble-pgfplots.tex' for where the files are
%% assumed to be if you use `\includetikz' when `\makepdf' is not defined.
\newcommand{\makepdf}{}

%% Import the necessary preambles.
\input{tex/src/preamble-mnras-free.tex}
\input{tex/src/preamble-style.tex}
\input{tex/build/macros/project.tex}
\input{tex/src/preamble-pgfplots.tex}
\input{tex/src/preamble-journals-astro.tex}
%\input{tex/src/preamble-biblatex.tex}
\input{tex/src/preamble-shortcuts-astro.tex}



\title[Elaphrocentre and LSBGs]{The
  role of the elaphrocentre in low surface brightness galaxy formation}

\author[Peper \& Roukema]{Marius Peper,$^{1}$
         Boudewijn F. Roukema$^{1,2}$
         \\
         $^1$ Institute of Astronomy, Faculty of Physics,
           Astronomy and Informatics, Nicolaus Copernicus
           University, Grudziadzka 5, 87-100 Toru\'n, Poland
         \\ %email: {\tt marius.peper@astro.umk.pl} % avoid this for antispam reasons if we can
         $^2$ Univ Lyon, Ens de Lyon, Univ Lyon1, CNRS, Centre de
           Recherche Astrophysique de Lyon UMR5574, F--69007, Lyon,
           France}

%Mohammad Akhlaghi$^{4,5}$
%$^4$ Instituto de Astrof\'isica de Canarias, Calle V\'ia L\'actea s/n, 38205 La Laguna, Tenerife, Spain.
%$^5$ Departamento de Astrof\'isica, Universidad de La Laguna, Avenida Astrof\'isico Francisco S\'anchez s/n, 38200 La Laguna, Tenerife, Spain.

\date{Accepted \ldots Received \ldots; in original form \ldots}
\pubyear{2020}

\begin{document}%\layout


\maketitle

%% Abstract
\begin{abstract}
%WCWC word count
%context
  {The formation mechanisms of giant low surface brightness galaxies (LSBGs), with high masses and low star formation rates, remain unclear.}
%aims
  {If a dark matter halo forms at the potential hill corresponding to a void of the cosmic web, which we denote the \enquote*{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.}
%method
  {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 reproducibility.}
%results
  {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 LSBG formation: while void galaxies are of lower mass, 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.}
%conclusions
  {Thus, void galaxies have higher spin parameters and later formation epochs; both are factors that may increase the probability of LSBGs forming in voids.}
%WCWC word count
\end{abstract}

\begin{keywords}
  methods: numerical, galaxies: evolution, cosmology: dark matter
\end{keywords}


\section{Introduction}

%context
The formation mechanisms of giant low surface brightness galaxies (LSBGs; \citealt*{Sandage84,Bothun87}) with high masses and low star formation rates remain unclear.
A promising hypothesis is that the void environment plays a significant role in LSBG formation.
Here, we numerically investigate some aspects of this hypothesis, by developing a reproducible pipeline that combines and builds on existing tools.
In the cosmological galaxy formation context, we denote the gravitational potential hill corresponding to a cosmic void on megaparsec scales as an \enquote*{elaphrocentre} in order to emphasise its gravitational role in opposition to that of a barycentre\footnote{From ancient Greek: \href{https://en.wiktionary.org/wiki/\%E1\%BC\%90\%CE\%BB\%CE\%B1\%CF\%86\%CF\%81\%CF\%8C\%CF\%82\#Ancient_Greek}{$\varepsilon\lambda\alpha\phi\rho o \varsigma$} (light) and \href{https://en.wiktionary.org/wiki/\%CE\%B2\%CE\%B1\%CF\%81\%CF\%8D\%CF\%82\#Ancient_Greek}{$\beta\alpha\rho\upsilon\varsigma$} (heavy), respectively.}.
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.
Overall, 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 generally 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.
A barycentric spatial location would tend to lead to a halo having early rapid mass accumulation and later weak infall is less likely.
These effects on the host haloes could lead, for a fixed mass, to galaxies closer to the elaphrocentre preferentially forming with weak star formation rates over long time scales.
Thus, elaphrocentric location might help to explain the formation of rare, giant low surface brightness galaxies.

The discovery of galaxies with low surface brightness is over three decades old \citep{Sandage84,Bothun87}.
The discovery led to many questions on how low surface brightness galaxies (LSBGs) can form \citep*{Schombert2011LSBG,
  Schombert2013LSBG,
  SchombertMcGaugh2014a,
  SchombertMcGaugh2014b}.
Observations by \cite{vanDokkum15} recently reignited interest in LSBGs.
These authors found a high abundance of large LSBGs that they called ultra diffuse galaxies (UDGs).
Van Dokkum et al.\/ (\citeyear{vanDokkum15}) proposed the new term and definition because the newly found galaxies do not seem to share the classic features expected for LSBGs.
The discovery of UDGs in the Coma Cluster again raises questions on formation scenarios for LSBGs in general.
It is not yet clear whether LSBGs and UDGs share the same formation scenario, especially since this is dependent on the choice of their definitions.
A common scenario for UDGs to form is to assume that these galaxies form in haloes with large spin components.
\citet{Rong2017} find from simulations that UDGs form naturally in high-spin haloes within the $\Lambda$CDM model; \citet{JiHoon15} has started investigating this observationally.
\citet{Chan2018} show that they can reproduce the observed quantities of red UDGs by imposing quenching, without assuming high spin haloes.
\citet{DiCintio2017} showed that UDGs can be produced in isolated dwarf galaxy haloes with stellar feedback and episodes of gas outflow.
\citet{Jiang2019} extended this work by investigating field UDGs.
Compared to other galaxies, LSBGs tend to be more isolated, located at the edges of or well inside voids of the cosmic web, i.e. they tend to be elaphrocentric.
\citet{Rosenbaum09} 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.
As summarised by \citet[][\SSS{}7.1, 7.2]{DiPaoloSalucci20}, LSBGs typically have low surface densities, below around $5 M_{\odot}\mathrm{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.

Thus, it is a fair assumption that the location of a galaxy in an underdense region benefits the formation of large, diffuse galaxies.
It has long been known that a high volume fraction of the Universe consists of large underdense regions \citep*[e.g.][]{GregoryThompson78,deLappGH86}.
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 \citet{Pan2012voids}'s analysis.
Voids vary by about an order of magnitude in size, spanning scales from $10$ to $200$~{\hMpc}, 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 \citep{HoyleVOg02voids, HoyleVog04voids, Pan2012voids,NadHot2013, SutterLWWWP14, PisaniLSW14stack, Pisani15voidcounts, PisaniSW15vpec, Mao2017}.
Depending on the chosen catalogue and void detection algorithm, a modest fraction of galaxies can be considered as being located in voids.
For example, \citet{Pan2012voids} estimate that about $7$\% of galaxies are located in voids.
Since the fraction of volume occupied by voids may be as high as 80\% \citep[for theoretical estimates, see e.g.][]{Colberg08voidcomparison,Cautun15voidfrac80pc}, they are suspected to play a key role in relation to dark energy \citep[e.g.][]{BCKRW16}.

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, the evolution of void galaxies will tend to involve many fewer mergers.
Galaxy mergers usually lead to bursts of star formation, making galaxies briefly much brighter.
Here, we hypothesise that a void will contribute to the creation of large galaxies of low surface brightness.
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 \citep{BirkhoffLanger23}, are not satisfied for a flat Friedmann--Lema{\^i}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 wrong than approximating a modest (in terms of linear mean density) underdensity in the same way, since gravity is attractive.
What is effectively a form of 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 \citep{RO19flatness,Ostrowski19}, after which it is expected to virialise at an overdensity of a few hundred times the mean density \citep[e.g.][]{LaceyCole93MN}.
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 the infall of matter, and thus weaken the star formation rate.
This can also be thought of as a void being approximated by a spatially compact domain in a relativistic Milne model -- almost 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 \enquote*{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.
Given that LSBGs and UDGs are diffuse, typically having low star formation rates, and tend to be located in low-density environments, it is thus possible that elaphrocentres play a significant role in their formation via these effects.

Here, we mainly focus on the degree to which the elaphrocentric location is responsible for these low, roughly constant star formation rates 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 \citep{Rong2017} and/or an increase in the typical virial radius.
We also (iii) estimate the magnitude of elaphrocentric acceleration as a basis for more detail studies.


To study this hypothesis, we present a highly reproducible \citep[][and references therein]{Akhlaghi2020maneage} galaxy formation simulation and analysis pipeline, which not only combines existing community tools for this particular purpose, but can also help in improving the existing tools by embedding them in a controlled software environment.
Our pipeline \citep*{PeperRoukemaBolejko2019} is constructed out of well established cosmological software tools, which can, in principle, be replaced in a modular way, provided that the user manages the input and output formats correctly.
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.
Parallelisation in several steps of the computational pipeline currently prevents byte-for-byte reproducibility.
We present the software packages and the pipeline in \SSS\ref{s-method-pipeline}.

In \SSS\ref{s-method-elaphro} we discuss definitions of the \enquote*{centre} of a void: we comment on the ambiguity between the mathematical and physical usage of the term \enquote*{barycentre} in relation to voids and propose some clarifications; we propose a Voronoi-cell based definition of the \enquote*{elaphrocentre}.
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 \SSS\ref{e-method-void-membership}.
We continue by describing how we study the dependence of the infall of matter into galaxies (\SSS\ref{s-method-infall}) and the dependence of galaxy sizes (\SSS\ref{s-method-size}) 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 \SSS\ref{s-method-elaphro-acc} 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 \SSS\ref{s-results}. We discuss these in \SSS\ref{s-discuss} and conclude in \SSS\ref{s-conclu}.
The full reproducibility package for this paper is available at {\projectzenodohref} and in a live {\sc git} repository\footnote{\projectgitrepository}.
%and archived\footnote{\projectgitrepositoryarchived} {\sc git} repositories. %% softwareheritage.org
The commit hash of the version of the source package used to produce this paper is \projectversion.


\section{Method}

\subsection{Software pipeline} \label{s-method-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.
Our pipeline is intended to be modular and our results are intended to be statistically reproducible.
Dark matter haloes are identified in an $N$-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 $N$-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.
The reproducibility structure is based on the Maneage template that aims for a high level of reproducibility \citep{Akhlaghi2020maneage}.
We follow \citet{RougierHinsen2017Repdefn} for the definitions of the \enquote*{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 \enquote*{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 {\sc git} commit hashes.
Modifications that we have made to upstream versions of codes are included as patch files in the reproducibility source package ({\projectzenodohref}).


\subsubsection{Initial conditions} \label{s-meth-initial-cond}
We use {\mpgraficname}-{\mpgraficversion} \citep{PrunetPichon08mpgrafic} to generate a set of standard initial conditions for a flat-space $N$-body simulation with the standard 3-torus topology (often called \enquote*{periodic boundary conditions}).
The {\mpgraficname} 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 \citep{ZA,ZA70b}.
Checks are made that the amplitude of the numerically generated power spectrum matches that of the input power spectrum.
We generate a simulation with $N={\Ncrootvalue}^3$ particles.
The comoving fundamental domain size, often called the \enquote*{box size}, is $L_{\mathrm{box}}={\Lboxvalue}$~Mpc/$h$.
These parameters give a dark matter particle mass of ${\sageparamPartMassvalue} \times 10^{10} ~ M_{\odot}$, which is a reasonable mass resolution for modest RAM and CPU resources.
We use the $\Lambda$CDM model as a proxy model that fits many observations.
The FLRW cosmological parameter settings include the current values of the matter density parameter $\Ommzero={\OmegaMvalue}$, the dark energy parameter $\OmLamzero={\OmegaLvalue}$ and the Hubble--Lema\^{\i}tre constant $H_0={\Hubblevalue}$~km/s/Mpc.

\subsubsection{Simulations} \label{s-meth-simulations}
For our $N$-body simulation, we chose \ramsesscalavname-{\ramsesscalavversion}, a fork of the widely used adaptive mesh code {\sc ramses}-3.0 \citep{Teyssier02}.
The {\ramsesscalavname} fork has modifications to comply with the MPI~3.0 recommended standards for inclusion of the MPI header file\footnote{\url{https://mpi-forum.org/docs/mpi-3.0/mpi30-report.pdf}} and optional extensions related to scalar averaging \citep{Roukema17silvir}.
The adaptive mesh structure of {\sc ramses} is designed to allow fast calculations that can resolve detailed gravitational behaviour in high density regions.
The maximum {\sc ramses} resolution, which effectively corresponds to a softening length, is set at {\ramsesparamLevelmaxname} $= \ramsesparamLevelmaxvalue$.
We produce snapshots starting at $t_{\mathrm{i}}={\TimeStartvalue}$~Myr with an equally spaced time step of $\Delta t={\TimeStepvalue}$~Myr, and convert to scale factor values using {\cosmdistname}-{\cosmdistversion}.
Newtonian gravitational potentials against the FLRW background are calculated in {\sc ramses} 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 (\SSS\ref{s-method-elaphro}).

\subsubsection{Halo detection} \label{s-meth-halo-det}

For detecting dark matter haloes, we use \rockstarname-{\rockstarversion} \citep{Behroozi13rockstar}, which uses a $7$-dimensional friends-of-friends (FoF) algorithm.
In contrast to other FoF halo finders, {\rockstarname} 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.
\citet{Knebe2011} found that {\rockstarname} performs excellently in recovering haloes from an $N$-body simulation.
We run {\rockstarname} using a linking length of $\rockstarparamFOFLINKINGLENGTHvalue$ and a minimum of $\rockstarparamMINHALOOUTPUTSIZEvalue$ particles per halo.
%The following human explanation of the virial radius definition will not change automatically, but the reader will be alerted if these are inconsistent.
We set the virial radius criterion for {\rockstarname} detection to {\tt \rockstarparamMASSDEFINITIONvalue}, i.e., $200$ times the critical density.


\subsubsection{Merger history trees}
We construct halo merger trees from the simulations \citep[][and references thereof]{1993ASPC...51..298R,RPQR97,1993ApJ...418L...1R}.
In this paper, instead of using the original Fortran77 routines from 1992, we use a more modern package, \ctreesname-{\ctreesversion} \citep{Behroozi13MHT}, which was designed to perform on outputs from {\rockstarname}.
In order to use these merger history trees for simulating galaxy evolution using {\sagename}, which was developed for following up simulations such as the Millenium simulation, we need to convert the trees to the {\sc LHaloTree} format.
We do the conversion with {\convertctreesname}-{\convertctreesversion}.

%\sometablesandfigures

\subsubsection{Star formation and galaxies} \label{s-meth-galaxies}

To form galaxies within our dark matter haloes, we use semi-analytical galaxy formation recipes \citep{1993ASPC...51..298R,KWG93,RPQR97,KauffCDW99}.
Again, rather than using the original code from 1992, we use {\sagename} \citep{Croton16SAGE}.
For an introductory review, see \citet{Baugh06review}.
We use the built-in functions of {\sagename} for estimating the \enquote*{size} and the infall rate history of each galaxy at the final output time.
All galaxies are assumed to form as disk galaxies in {\sagename}, with a disk radius of
\begin{equation}
  \rdiskscale = \frac{\lambda}{\sqrt{2}} R_{\mathrm{vir}},
\label{sage-disk-radius}
\end{equation}
% := \sqrt{\lambda_x^2 + \lambda_y^2 + \lambda_z^2}$ % putting this formula would require explaining what \lambda_x, \lambda_y, \lambda_z are
where $R_{\mathrm{vir}}$ is the virial radius (set in {\rockstarname} at {\tt \rockstarparamMASSDEFINITIONvalue}) and $\lambda$ is the dimensionless spin parameter computed using {\rockstarname}'s halo properties.
The introduction of the parameter $\lambda$ is generally attributed to \citet[][eqs~(35), (37)]{Peebles69spin}.
The parameter is now widely used \citep[e.g.,][eq.~(10)]{MoShudeWhite98spinparam}, and is often normalised by a $\sqrt{2}$ factor \citep[][eq.~(5)]{Bullock01spin}, which is the convention chosen in {\sagename}, and adopted in this work:
\begin{equation}
  \lambda := 2^{-1/2}\,J\,|E|^{1/2}\,G^{-1}\,M^{-5/2} \,,
  \label{e-defn-spin}
\end{equation}
where $J$ is the total angular momentum, $E$ is the total energy in the non-relativistic sense, $G$ is the Newtonian gravitational constant, and $M$ is the halo mass.

The infalling mass is defined at a given time step by
\begin{equation}
  \Delta M_{\mathrm{infall}}(t_i) := f_{\mathrm{reion}} \, f_{\mathrm{b}} \, M_{\mathrm{vir}} - M_{\mathrm{tot}},
\end{equation}
where $f_{\mathrm{reion}}$ is the reionisation factor, which estimates the effect of ionisation of the intergalactic medium from early stars;
the baryon fraction $f_{\mathrm{b}} = \sageparamBaryonFracvalue$ determines what fraction of total matter is baryonic (assumed to be the same for any dark matter halo/galaxy pair); %\sageparamBaryonFracname
$M_{\mathrm{vir}}$ is the virial (total matter) mass of the halo;
and $M_{\mathrm{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 $M_{\mathrm{infall}} < 0$ are interpreted to mean that baryonic mass is ejected from the galaxy.
We extended {\sagename} 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 \citep{1993ASPC...51..298R,RPQR97}; 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, $a(t)=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 \citep{1993ASPC...51..298R}, other effects from mergers than conservation of mass, such as merger-induced starbursts, are also assumed in {\sagename}, but with a power law dependence on the satellite mass rather than direct proportionality \citep[][eq.~(27)]{Croton16SAGE}.

\subsubsection{Voids} \label{s-method-voids}
We identify the void environment of galaxies using the {\revolvername} watershed void finder based on {\sc zobov} \citep{Neyrinck2008Zobov}, 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.
\citet{NadathurHotch15voidsII} 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 \citep[e.g.][]{Mao2017}, 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 {\revolvername}.
In addition to \citet{NadathurHotch15voidsI}'s definition of the circumcentre, we calculate the position of the elaphrocentre, as defined below in \SSS\ref{s-method-elaphro}.
We add a routine to read in simulation data in {\sc Gadget-2} format \citep{Springel05GADGET2}, the default output format that we chose for the {\sc ramses} $N$-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.

\subsection{Elaphrocentre and other definitions of void centres} \label{s-method-elaphro}
To investigate if a galaxy's position in a void -- its elaphrocentric location -- has a significant effect on LSBG formation, we first need to clarify earlier terminology regarding void centres from the literature, and we need to define the elaphrocentre.

\citet[][\SSS{}2.3]{NadathurHotch15voidsI} define the \enquote{\em{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.
\citet[][\SSS{}2.3 (ii)]{Nadathur2017VoidCentre} rename this the \enquote{\em{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 {\em elaphrocentre} similarly.
In a void identified by the watershed algorithm, we identify the particle at which the potential is highest.
We then identify the three adjacent Voronoi cells with the highest potentials, which together with the highest potential cell 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 \enquote*{macrocentre} or \enquote*{volume-weighted barycentre}.
% volume-weighted barycentre is a term that is used *in* revolver's code, see revolver/python_tools/zobov.py.
This is defined (\citealt{Sutter15VIDE}, \SSS{}3, eq.~(4); \citealt{NadathurHotch15voidsI}, \SSS{}2.3, eq.~(2)) as the volume-weighted mean of the position vectors $\vec{x}_i$ of all DM particles identified as being in the void, i.e., $ \vec{c}_{\mathrm{vwb}} := \left({\sum_i {{V_i}\vec{x}_i}}\right)/{\sum_i V_i}$, where $V_i$ is the Voronoi cell volume associated with the $i$-th particle.
Since the position is not mass-weighted, it is unrelated to the normal Newtonian definition of a barycentre for a particle distribution, $ \vec{c}_{\mathrm{m}} := \left({\sum_i {{m_i}\vec{x}_i}}\right)/{\sum_i m_i}$, where $m_i$ is the mass of the $i$-th particle.
% Chose a slightly different name to underline that the volume-weighted barycentre is different from the Newtonian definition of the barycentre. We maybe even remove the sentence about the Newtonian barycentre?
In the continuous limit to arbitrarily high particle resolution (assuming a continuous fluid), the volume-weighted barycentre approaches $\vec{c}_{\mathrm{vwb}} = {\int_{\cal D} \vec{x} \,\diffd V}/{\int_{\cal D} \diffd V}$ over a spatial domain ${\cal D}$.
Clearly, in the continuous limit, $\vec{c}_{\mathrm{vwb}}$ is the geometrical centroid of the domain $\cal D$, which in geometry is often termed the \enquote*{barycentre}.
In general, this corresponds to the astronomical barycentre only if the density distribution in the domain $\cal D$ is uniform.
In other words, the volume-weighted barycentre contains no information about the density distribution within $\cal D$ apart from numerical noise.

By definition, apart from discretisation and numerical effects, $\vec{c}_{\mathrm{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, \citet{NadathurHotch15voidsI} are correct that $\vec{c}_{\mathrm{vwb}}$ depends on the presence rather than the absence of tracers.
However, the fundamental problem with using $\vec{c}_{\mathrm{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 $\vec{c}_{\mathrm{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 \citet{NadathurHotch15voidsI} found the density at $\vec{c}_{\mathrm{vwb}}$ to be higher than at the circumcentre.
Thus, in the context of cosmological voids, the physical relevance of $\vec{c}_{\mathrm{vwb}}$ is unclear, and if used, we recommend that it be described by the term \enquote*{geometrical centroid} (or \enquote*{boundary centroid}) rather than \enquote*{macrocentre} or \enquote*{volume-weighted barycentre}.

\subsection{Analysis} \label{s-method-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.

\subsubsection{Void membership and elaphrocentric distance} \label{e-method-void-membership}
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 ($\mathbb{R}^3/\mathbb{Z}\times \mathbb{Z} \times \mathbb{Z} \equiv \mathbb{S}^1 \times \mathbb{S}^1 \times \mathbb{S}^1$) spatial geodesic comoving distances.
This is often described more loosely as \enquote*{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 \citet{Icke84} by reverting a simple toy model for collapsing density perturbations and numerically by \citet{Sheth04} in $N$-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 \mbox{$>50$\%} merging identity criterion initially published in 1993 for merger history trees (\citealt{1993ApJ...418L...1R}, \SSS3; \citealt{RPQR97}, \SSS2.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 {\revolvername}.
Thus, for any given void, we have a list of particles that approximately define the void.
Since we know from {\rockstarname} 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 $\relaphro \leq {\revolverparamiRadFracElaphrovalue}\, R_{\mathrm{eff}}$ from the elaphrocentre, where $R_{\mathrm{eff}}$ is the void effective radius calculated by {\revolvername}.
The effective radius $R_{\mathrm{eff}}$ is the radius of a hypothetical sphere that has the same volume as the total volume of all Voronoi cells constituting a void.
The $\relaphro \leq {\revolverparamiRadFracElaphrovalue}\, R_{\mathrm{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 $\haloinvoidfrac$ strictly greater than $\haloinvoidminfrac = {\HaloIsInVoidMinFractionvalue}$ of the particles in a halo $\cal H$ be members of a void $\cal V$ for a galaxy in that halo to be considered a member of the void.
As in the case of the \mbox{$>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.


\subsubsection{Infall dependence on environment} \label{s-method-infall}
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, $M(t)$, so that $\diffd M/\diffd t \approx \Delta M/\Delta t\, (t)$ is the infall rate -- of 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 $M$, 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 $\diffd M/\diffd t$ by an exponential, of the form
\begin{equation}
  \frac{\diffd M}{\diffd t} (t) = {\InfallAmplitude} \exp({-t/\InfallDecayRate})\,,
  \label{e-infall-model}
\end{equation}
where ${\InfallAmplitude}$ is the infall amplitude and ${\InfallDecayRate}$
is a decay time scale.
Typical disk galaxies should have low time scales ${\InfallDecayRate}$, while LSBGs should have high ${\InfallDecayRate}$, corresponding to approximately flat infall rates.

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 \SSS\ref{s-method-pipeline}, starting from the first time step in which a galaxy is modelled to form in a dark matter halo.

We first compare ${\InfallAmplitude}$ and ${\InfallDecayRate}$ for galaxies depending on their classification as void or non-void galaxies.
The infall histories are calculated using our modification of {\sagename}.
The prediction for creating LSBGs is that void galaxies should tend to have low amplitude ${\InfallAmplitude}$ and a high (slow) decay rate ${\InfallDecayRate}$, and vice versa for non-void galaxies.

\begin{figure}
  \includegraphics[width=0.9\columnwidth]{fig/void_population\revolverparamiRadFracElaphrovalue}
  \caption{Histogram of void effective radii $R_{\mathrm{eff}}$, using the full dark matter particle distribution, which leads to voids much smaller than typically observed or found in simulations when traced by galaxies.
\label{void_population}}
\end{figure}

\begin{figure}
\includegraphics[width=0.9\columnwidth]{fig/number_counts\revolverparamiRadFracElaphrovalue}
\caption{Histogram of {\CentreAdjectivevalue} distance of galaxies identified as being located in voids.
  The distribution at $r/R_{\mathrm{eff}} < 1$ is typical of observational and simulated void profiles, with most galaxies located near the effective radius.
  The distribution at $r/R_{\mathrm{eff}} \ge 1$ can be interpreted as showing galaxies in the outer parts of voids that are generally quite asymmetrical.
\label{results-galaxy-hist}}
\end{figure}

\subsubsection{Galaxy size dependence on environment} \label{s-method-size}

The other parameter that may indicate an elaphrocentric contribution to a disk galaxy becoming an LSBG is the disk scale length $\rdiskscale$ (for a density profile $\rho \propto \exp(-r/\rdiskscale)$).
This can be converted to a disk half-mass radius using the relation $r_{1/2} = \nu \rdiskscale$, where $\nu$ solves $(\nu+1) \exp(-\nu) = 1/2$ \citep[e.g.][$\nu=1.678$]{Kravtsov_2013}.
As stated above (\SSS\ref{s-meth-galaxies}; Eq.~\eqref{sage-disk-radius}), {\sagename} calculates $\rdiskscale$.
Since voids are underdensities, dark matter haloes forming in voids will tend to collapse somewhat later than in overdensities \citep[e.g.][App.~A]{LaceyCole93MN}.
In an expanding FLRW universe, the critical density $\rho_{\mathrm{crit}}$ decreases, so for a fixed virialisation overdensity threshold and fixed mass, $R_{\mathrm{vir}}$ increases with time.
Thus, modelling galaxy disk scale lengths as being proportional to the halo virial radius (Eq.~\eqref{sage-disk-radius}), it would be reasonable to expect void galaxies to have greater $\rdiskscale$ than non-void galaxies, for a fixed value of the spin parameter $\lambda$.
Elaphrocentric galaxies will typically undergo a more isolated evolution than barycentric galaxies, with fewer merger events.
\citet{DOnghia08} found, based on $N$-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 $\rdiskscale$ directly, and $R_{\mathrm{vir}}$ and $\lambda$ individually.

\subsubsection{Elaphro-acceleration} \label{s-method-elaphro-acc}

As a complement to the direct analyses of simulated galaxies via {\sagename}, we also investigate accelerations near the elaphrocentres, as preparation for future studies of elaphrocentric effects that might help form 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 \enquote*{elaphro-acceleration} without requiring the presence of a halo or galaxy.
In this simplified model, we assume that a Milky-Way--like galaxy forms at the elaphrocentre.
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.
We adopt a radius of the region from which dark matter originated, for a Milky-Way--like galaxy, as $\rMW=\ElaphrocentreHaloOriginalRadiusMpcvalue$~{Mpc}.
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 {\sc ramses} (\SSS\ref{s-meth-simulations}).
Since {\sc ramses} 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.
For any given void, we interpolate the potential linearly, and calculate the acceleration as the gradient of the potential.
We sample the elaphro-acceleration at six positions which lie on the 2-sphere with radius $\rMW$, i.e. at $((\pm \rMW,0,0),(0,\pm \rMW,0),(0,0,\pm \rMW)$, 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 $\vec{v}$ for each of the six positions, and find the mean values $\dot{v}_{\parallel}$ and $\dot{v}_\perp$, respectively.
As stated in \SSS\ref{s-method-voids}, some of the six positions may fall outside of a void when the void is too small; we ignore the void in such cases.

\subsection{Reproducibility versus cosmic variance} \label{s-method-reprod-cosvar}

We present results below (\SSS\ref{s-results}) with a preference for reproducibility over cosmic variance.
Our pipeline, using the Maneage template \citep{Akhlaghi2020maneage}, 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.
We formalise the cosmic variance of a given parameter by the standard deviation $\sigmacosvar$ (not the literal variance) over repeated runs of commit {\CosmicVarCommitID} of our source package.
We preferentially verify consistency of a given parameter in our published run by requiring that it agree with the mean from the ensemble of runs within $\NVerifySigmas$ times the stated random error $\sigmarandom$.
For parameters with high $\sigmacosvar$, we require agreement within  $\NVerifySigmas$ times $\sigmacosvar$, and we state $\sigmacosvar$ as an additional uncertainty, using the notation \enquote*{$\pmcosvar$}.
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.


\section{Results} \label{s-results}

\subsection{Simulation pipeline} \label{s-pipeline}

In the final time step of our $N={\Ncrootvalue}^3$-particle simulation, we detected {\NHaloesFinalTimevalue} haloes, with $\PlotsSizeiNGalsTotalvalue$ galaxies evolved along the merger history trees for these haloes.
Among these galaxies, $\PlotsSizeiNGalsInIntervalvalue$ have virial mass $M$ in the range ${\MassgapLowerLimitvalue}$--${\MassgapUpperLimitvalue} \mathrm{M}_\odot /h$, which we study with the aim of seeking a factor in the formation of massive LSBGs.
Applying the $\haloinvoidfrac > {\HaloIsInVoidMinFractionvalue}\%$ definition in \SSS\ref{e-method-void-membership}, we identify ${\PlotsSizeiNVoidGalsTotalvalue}$ galaxies in voids, among which ${\PlotsSizeiNVoidGalsvalue}$ 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 \citet{Pan2012voids}.
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 $\haloinvoidfrac$ and dependent on their relative distance to the voids centre $r/R_{\mathrm{eff}}$.
There are $\nvoidsiNvoidsvalue$ voids in the final time step.
The void size distribution is shown in Fig.~\ref{void_population}.
Since we detect voids physically, following the dark matter particular distribution, as described in \SSS\ref{s-method-voids}, it is unsurprising that the void population is dominated by voids a few Mpc/$h$ in size, with a correspondingly higher number density than that typically seen in void catalogues calculated using galaxies as tracers.

\begin{figure}
  \includegraphics[width=0.9\columnwidth]{fig/infall_frac_amp\revolverparamiRadFracElaphrovalue}
  \caption{Amplitude of infall rate ${\InfallAmplitude}$ versus fraction $\haloinvoidfrac$ of a galaxy's host halo composed of void particles.
    A Theil--Sen robust linear fit to the relation, ${\InfallAmplitude} = \left[(\InfallAmpFracZerovalue \pm \InfallAmpFracSigZerovalue) + (\InfallAmpFracSlopevalue \pm \InfallAmpFracSigSlopevalue \pmcosvar \InfallAmpFracSlopevalueCosmStdDev) \haloinvoidfrac\right]\,M_{\odot}/\mathrm{yr} $, is shown (the slope is sensitive to cosmic variance).
    Galaxies towards the left ($\haloinvoidfrac = 1$) are those best identified as being in voids.
    Plain text table for this figure through to Fig.~\ref{f-infall-rreff-tau}: \mbox{\href{\projectzenodofilesbase/voidgals_infall.dat}{\projectzenodoid/voidgals\_infall.dat}}.
    \label{f-infall-frac-amp}}
\end{figure}

\begin{figure}
  \includegraphics[width=0.9\columnwidth]{fig/infall_frac_tau\revolverparamiRadFracElaphrovalue}
  \caption{Infall decay rate ${\InfallDecayRate}$ versus fraction $\haloinvoidfrac$ of a galaxy's host halo composed of void particles, as in Fig.~\protect\ref{f-infall-frac-amp}, with a Theil--Sen robust linear fit ${\InfallDecayRate} = \left[(\InfallTauFracZerovalue \pm \InfallTauFracSigZerovalue) + (\InfallTauFracSlopevalue \pm \InfallTauFracSigSlopevalue \pmcosvar \InfallTauFracSlopevalueCosmStdDev) \haloinvoidfrac\right] \, \mathrm{Gyr}$, is shown.
    \label{f-infall-frac-tau}}
\end{figure}


\begin{figure}
  \includegraphics[width=0.9\columnwidth]{fig/infall_rreff_amp\revolverparamiRadFracElaphrovalue}
  \caption{Amplitude of infall rate ${\InfallAmplitude}$ versus elaphrocentric location $r/R_{\mathrm{eff}}$ of a galaxy's host halo.
    The fit displayed is ${\InfallAmplitude} = \left[(\InfallAmpReldistZerovalue \pm \InfallAmpReldistSigZerovalue) + (\InfallAmpReldistSlopevalue \pm \InfallAmpReldistSigSlopevalue) \,r/R_{\mathrm{eff}}\right] \,M_{\odot}/\mathrm{yr}$.
    As in Fig.~\protect\ref{f-infall-frac-amp}, galaxies towards the left are those best identified as being in voids, but voidness is characterised in this plot by a lower elaphrocentric distance $r/R_{\mathrm{eff}}$, instead of by a higher void fraction $\haloinvoidfrac$.
    \label{f-infall-rreff-amp}}
\end{figure}

\begin{figure}
  \includegraphics[width=0.9\columnwidth]{fig/infall_rreff_tau\revolverparamiRadFracElaphrovalue}
  \caption{Infall decay rate ${\InfallDecayRate}$ versus versus elaphrocentric location $r/R_{\mathrm{eff}}$ of a galaxy's host halo, as in Fig.~\protect\ref{f-infall-rreff-amp}, with robust fit ${\InfallDecayRate} = \left[ (\InfallTauReldistZerovalue \pm \InfallTauReldistSigZerovalue) + (\InfallTauReldistSlopevalue \pm \InfallTauReldistSigSlopevalue) \,r/R_{\mathrm{eff}}\right] \,\mathrm{Gyr}$.
    \label{f-infall-rreff-tau}}
\end{figure}


\subsection{Infall rate}  \label{s-results-infall}
As described in \SSS\ref{s-method-infall}, we first compared infall rates for galaxies as separate void and non-void populations.
For each galaxy mass infall history $\diffd M/\diffd t\,(t)$, we first find a linear least-squares best fit to $\log_{10}(\diffd M/\diffd t)$ versus $t$ for time steps where $\diffd M/\diffd t\,(t) > 0$.
The optimal parameters of this fit are used to find a non-linear least-squares best fit of $\diffd M/\diffd t (t)$ to a decaying exponential (Eq.~\eqref{e-infall-model}), starting from the first time step with $\diffd M/\diffd t\,(t) > 0$ and no longer excluding time steps with $\diffd M/\diffd t = 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.

\begin{table}
  \centering
  \caption{Exponential decaying fit parameters (medians and standard error in the median) for infall rates for void and non-void galaxies.
    A large standard error in the time scale, of the order of the age of the Universe, indicates that many fits represent nearly constant infall rates.
    \label{Infall-rate-table}}
  $\begin{array}{l c c}
    \hline
    & \mbox{in voids} & \mbox{not in voids} \\
    \hline
    \log_{10}(A) \,({M_{\odot}/\mathrm{yr}})
    & {{\InfallMedianValueAmplitudeVoidsvalue} \pm {\InfallRobustErrorAmplitudeVoidsvalue}}
    & {{\InfallMedianValueAmplitudeNonVoidsvalue} \pm {\InfallRobustErrorAmplitudeNonVoidsvalue}}\rule{0ex}{2.7ex}%strut
    \\
    \InfallDecayRate \,(\mathrm{Gyr})
    & {\InfallMedianValueDecayVoidsvalue} \pm {\InfallRobustErrorDecayVoidsvalue}
    & {\InfallMedianValueDecayNonVoidsvalue} \pm {\InfallRobustErrorDecayNonVoidsvalue} \\
    \hline
  \end{array}$
\end{table}

We find $\mvoidgood={\VoidFitiGoodParametervalue}$ valid fits for the void galaxies, and $\mvoidbadtot={\VoidFitiBadParameterTotalvalue}$ fits that are rejected either as failed fits ($\mvoidbadfit = {\VoidFitiBadParameterFailedFitvalue}$) or as having an unreasonably high amplitude, ${\InfallAmplitude} > {\VoidFitiParameterThresholdInputvalue}$~$M_{\odot}/$yr, indicating a physically unrealistic fit ($\mvoidbadthresh = {\VoidFitiBadParameterThresholdvalue}$ cases).
For non-void galaxies we find $\mnonvoidgood={\NonVoidFitiGoodParametervalue}$ valid fits and $\mnonvoidbadtot={\NonVoidFitiBadParameterTotalvalue}$ fits rejected either as failed fits ($\mnonvoidbadfit = {\NonVoidFitiBadParameterFailedFitvalue}$) or as having unphysically high amplitudes ($\mnonvoidbadthresh = {\NonVoidFitiBadParameterThresholdvalue}$ cases).
For galaxies in host haloes with virial masses in the range ${\MassgapLowerLimitvalue}$--${\MassgapUpperLimitvalue} \mathrm{M}_\odot /h$, we find the medians listed in Table~\ref{Infall-rate-table}, 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 $(\HaloMassVoidsMedvalue \pm \HaloMassVoidsStderrMedvalue) \times 10^{11} M_{\odot}$, lower than that of the non-void galaxies, $(\HaloMassNonvoidsMedvalue \pm \HaloMassNonvoidsStderrMedvalue) \times 10^{11} M_{\odot}$.

\begin{table}
  \centering
  \caption{Median disk galaxy scale length, spin parameter and virial radius for void and non-void galaxies in the mass interval ${\MassgapLowerLimitvalue}$--${\MassgapUpperLimitvalue} \mathrm{M}_\odot /h$ and on all mass scales.
    \label{Galaxy-size-table}}
  $\begin{array}{l c c}
    \hline
    &\mbox{in voids} & \mbox{not in voids} \\
    \hline
    %{\revolverparamiRadFracElaphrovalue} &
    \multicolumn{3}{c}{\mbox{restricted mass interval}} \\
    \hline
    \rdiskscale (\mbox{kpc}/h)
    & {\PlotsSizeiRAvgVoidsvalue} \pm {\PlotsSizeiErrorRAvgVoidsvalue}
    & {\PlotsSizeiRAvgNonvoidsvalue} \pm {\PlotsSizeiErrorRAvgNonvoidsvalue}\rule{0ex}{2.7ex}%strut
    \\
    \lambda
    & {\SpiniSpinAvgVoidsvalue} \pm {\SpiniErrorSpinAvgVoidsvalue}
    & {\SpiniSpinAvgNonvoidsvalue} \pm {\SpiniErrorSpinAvgNonvoidsvalue} \\
    R_{\mathrm{vir}} (\mbox{kpc}/h)
    & {\RviriRvirAvgVoidsvalue} \pm {\RviriErrorRvirAvgVoidsvalue}
    & {\RviriRvirAvgNonvoidsvalue} \pm {\RviriErrorRvirAvgNonvoidsvalue} \\
    \hline
    \multicolumn{3}{c}{\mbox{all mass scales}} \\
    \hline
    \rdiskscale (\mbox{kpc}/h)
    & {\PlotsSizeiRAvgVoidsTotalvalue} \pm {\PlotsSizeiErrorRAvgVoidsTotalvalue}
    & {\PlotsSizeiRAvgNonvoidsTotalvalue} \pm {\PlotsSizeiErrorRAvgNonvoidsTotalvalue}\rule{0ex}{2.7ex}%strut
    \\
    \lambda
    & {\SpiniSpinAvgVoidsTotalvalue} \pm {\SpiniErrorSpinAvgVoidsTotalvalue}
    & {\SpiniSpinAvgNonvoidsTotalvalue} \pm {\SpiniErrorSpinAvgNonvoidsTotalvalue} \\
    R_{\mathrm{vir}} (\mbox{kpc}/h)
    & {\RviriRvirAvgVoidsTotalvalue} \pm {\RviriErrorRvirAvgVoidsTotalvalue}
    & {\RviriRvirAvgNonvoidsTotalvalue} \pm {\RviriErrorRvirAvgNonvoidsTotalvalue} \\

    %  ${\revolverparamiiRadFracElaphrovalue}$ & ${\PlotsSizeiiRAvgNonvoidsvalue} \pm {\PlotsSizeiiErrorRAvgNonvoidsvalue}$
    %  & $ {\PlotsSizeiiRAvgVoidsvalue} \pm {\PlotsSizeiiErrorRAvgVoidsvalue}$ \\
    %  ${\revolverparamiiiRadFracElaphrovalue}$ & ${\PlotsSizeiiiRAvgNonvoidsvalue} \pm {\PlotsSizeiiiErrorRAvgNonvoidsvalue}$
    %& $ {\PlotsSizeiiiRAvgVoidsvalue} \pm {\PlotsSizeiiiErrorRAvgVoidsvalue}$
    \hline
  \end{array}$
\end{table}

We find (Table~\ref{Infall-rate-table}) no significant difference in either the amplitude ${\InfallAmplitude}$ or the time scale ${\InfallDecayRate}$ of infall between the void and non-void galaxy populations.
The dispersion in infall patterns within each population is as great as both ${\InfallAmplitude}$ and ${\InfallDecayRate}$ 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 $t^{\mathrm{f}}_{\mathrm{v}}=\FormationTimeVoidsMedvalue \pm \FormationTimeVoidsStderrMedvalue$~Gyr and $t^{\mathrm{f}}_{\mathrm{nv}}=\FormationTimeNonvoidsMedvalue \pm \FormationTimeVoidsStderrMedvalue $~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 {\sagename} for the galaxy is non-zero.

We investigate the results for ${\InfallAmplitude}$ and ${\InfallDecayRate}$ more closely by checking if either ${\InfallAmplitude}$ or ${\InfallDecayRate}$ has a dependence on either the fraction of a host halo's particles that are identified as void particles, $\haloinvoidfrac$, or on the host halo's elaphrocentric distance $r/R_{\mathrm{eff}}$.
Figures~\ref{f-infall-frac-amp}--\ref{f-infall-rreff-tau} show the dependence of ${\InfallAmplitude}$ and ${\InfallDecayRate}$ on $\haloinvoidfrac$ and $r/R_{\mathrm{eff}}$ for all void galaxies.
The $\haloinvoidfrac$ axis is shown with $\haloinvoidfrac$ 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 $\haloinvoidfrac$.
However, Fig.~\ref{f-infall-frac-amp} does show a modestly significant non-zero slope, with galaxies better identified as void galaxies (having a higher value of $\haloinvoidfrac$) tending to have somewhat higher amplitudes ${\InfallAmplitude}$.
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 statistics as above, indicated numerically in the figure captions, are not significantly non-zero.

%We find a modest but significant difference in amplitudes of the infall rates.
%The median amplitude of the infall rate is about ${\InfallMedianAmplitudeVoidsPercentLowervalue}$\% weaker for void galaxies than for non-void galaxies, significant at about the ${\InfallMedianAmplitudeSignifvalue}\,\sigma$ level (for a Gaussian interpretation of the standard errors).
%However, this difference is quite likely to be degenerate with the difference in typical host halo masses, which differ between the two populations.
%More massive haloes tend to form at the knots of the cosmic web.
%So there is both moderately lower infall in the void galaxies and a stronger factor in lower host halo mass in void galaxies.

\subsection{Galaxy Sizes} \label{s-results-size}

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 (\SSS\ref{s-method-size}), to check the size of a galaxy at the final output time step, we use the disk scale length $\rdiskscale$ provided by {\sagename}.
The results for the galaxies divided into void and non-void populations are shown in Table~\ref{Galaxy-size-table}, where we list the disk exponential scale length $\rdiskscale$, the spin parameter $\lambda$ and the virial radius $R_{\mathrm{vir}}$.
Table~\ref{Galaxy-size-table} shows a significant difference for the overall scale length $\rdiskscale$.
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 \citet[][fig.~2, left]{vdWeygaert11VGS}'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 the table: the void galaxy population has a much smaller median virial radius $R_{\mathrm{vir}}$ than the non-void population, but an insignificantly higher spin parameter $\lambda$.
The slightly greater spin parameter appears insufficient to compensate or override the lower $R_{\mathrm{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 $\lambda$ listed in Table~\ref{Galaxy-size-table} are consistent with the friends-of-friends halo estimate of \citet[][\SSS2.4, para.\/ 8]{ZjupaSpringel2017}, $\lambda_{\mathrm{B,FOF}} = 0.0414$, while our void host halo values of $\lambda$ are slightly higher.

\begin{figure}
  \includegraphics[width=0.9\columnwidth]{fig/galaxy_size_frac\revolverparamiRadFracElaphrovalue}
  \caption{Galaxy disk scale length $\rdiskscale$ versus fraction $\haloinvoidfrac$ of a galaxy's host halo composed of void particles.
    The fit is $\rdiskscale = \left[ (\GalsizeFracZerovalue \pm \GalsizeFracSigZerovalue) + (\GalsizeFracSlopevalue \pm \GalsizeFracSigSlopevalue \pmcosvar \GalsizeFracSlopevalueCosmStdDev) \haloinvoidfrac \right] \,\mathrm{kpc}/h$.
    Plain text table for this figure through to Fig.~\ref{f-galaxy-Rvir-rreff}: \mbox{\href{\projectzenodofilesbase/voidgals_infall.dat}{\projectzenodoid/voidgals\_infall.dat}}.
    \label{f-galaxy-size-frac}}
\end{figure}

\begin{figure}
  \includegraphics[width=0.9\columnwidth]{fig/galaxy_size_rreff\revolverparamiRadFracElaphrovalue}
  \caption{Disk scale length $\rdiskscale$ versus elaphrocentric location $r/R_{\mathrm{eff}}$, as in Fig.~\protect\ref{f-galaxy-size-frac}.
    The fit is $\rdiskscale = \left[ (\GalsizeReldistZerovalue \pm \GalsizeReldistSigZerovalue) + (\GalsizeReldistSlopevalue \pm \GalsizeReldistSigSlopevalue) \,r/R_{\mathrm{eff}} \right] \,\mathrm{kpc}/h$.
    \label{f-galaxy-size-rreff}}
\end{figure}

\begin{figure}
  \includegraphics[width=0.9\columnwidth]{fig/galaxy_spin_frac\revolverparamiRadFracElaphrovalue}
  \caption{Dimensionless spin parameter $\lambda$ (Eq.~\protect\eqref{e-defn-spin}) versus $\haloinvoidfrac$.
    The fit is $\lambda = (\GalspinFracZerovalue \pm \GalspinFracSigZerovalue) + (\GalspinFracSlopevalue \pm \GalspinFracSigSlopevalue \pmcosvar \GalspinFracSlopevalueCosmStdDev) \haloinvoidfrac$.
    \label{f-galaxy-spin-frac}}
\end{figure}

\begin{figure}
  \includegraphics[width=0.9\columnwidth]{fig/galaxy_spin_rreff\revolverparamiRadFracElaphrovalue}
  \caption{Dimensionless spin parameter $\lambda$ (Eq.~\protect\eqref{e-defn-spin}) versus elaphrocentric location $r/R_{\mathrm{eff}}$, as in Fig.~\protect\ref{f-galaxy-spin-frac}.
    The fit is $\lambda = (\GalspinReldistZerovalue \pm \GalspinReldistSigZerovalue) + (\GalspinReldistSlopevalue \pm \GalspinReldistSigSlopevalue \pmcosvar \GalspinReldistSlopevalueCosmStdDev) \,r/R_{\mathrm{eff}}$ (this is sensitive to cosmic variance).
    \label{f-galaxy-spin-rreff}}
\end{figure}

\begin{figure}
  \includegraphics[width=0.9\columnwidth]{fig/galaxy_Rvir_frac\revolverparamiRadFracElaphrovalue}
  \caption{Virial radius $R_{\mathrm{vir}}$ (Eq.~\protect\eqref{e-defn-spin}) versus $\haloinvoidfrac$.
    The fit is $R_{\mathrm{vir}} = \left[ (\GalRvirFracZerovalue \pm \GalRvirFracSigZerovalue) + (\GalRvirFracSlopevalue \pm \GalRvirFracSigSlopevalue) \haloinvoidfrac \right] \,\mathrm{kpc}/h$.
    \label{f-galaxy-Rvir-frac}}
\end{figure}

\begin{figure}
  \includegraphics[width=0.9\columnwidth]{fig/galaxy_Rvir_rreff\revolverparamiRadFracElaphrovalue}
  \caption{Virial radius $R_{\mathrm{vir}}$ versus elaphrocentric location $r/R_{\mathrm{eff}}$, as in Fig.~\protect\ref{f-galaxy-Rvir-frac}.
    The fit is $R_{\mathrm{vir}} = \left[ (\GalRvirReldistZerovalue \pm \GalRvirReldistSigZerovalue) + (\GalRvirReldistSlopevalue \pm \GalRvirReldistSigSlopevalue \pmcosvar \GalRvirReldistSlopevalueCosmStdDev) \,r/R_{\mathrm{eff}} \right] \,\mathrm{kpc}/h$.
    \label{f-galaxy-Rvir-rreff}}
\end{figure}

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~\ref{Galaxy-size-table}, we see that $R_{\mathrm{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 $\rdiskscale$ also exists as a function of a galaxy's void location, Figs~\ref{f-galaxy-size-frac} and \ref{f-galaxy-size-rreff} examine the dependence of $\rdiskscale$ on $\haloinvoidfrac$ and elaphrocentric distance for void galaxies.
The slope of the best fit in Fig.~\ref{f-galaxy-size-frac}, $\diffd \rdiskscale/\diffd \haloinvoidfrac = \GalsizeFracSlopevalue \pm \GalsizeFracSigSlopevalue \pmcosvar \GalsizeFracSlopevalueCosmStdDev$~kpc/$h$, is not significantly non-zero when we take into account cosmic variance (the distribution of $\diffd \rdiskscale/\diffd \haloinvoidfrac$ over repeated runs includes a strong tail of values that are not significantly non-zero).
The dependence on elaphrocentric position (Fig.~\ref{f-galaxy-size-rreff}) is not significant either.

In Figs~\ref{f-galaxy-spin-frac}--\ref{f-galaxy-Rvir-rreff}, we investigate whether the spin parameter $\lambda$ or the virial radius $R_{\mathrm{vir}}$ is more responsible for the modest reduction in the disk scale length of void galaxies as indicated in Table~\ref{Galaxy-size-table}.
Out of the four figures (Figs~\ref{f-galaxy-spin-frac}--\ref{f-galaxy-Rvir-rreff}), two show highly significant slopes: Figs~\ref{f-galaxy-spin-frac} and \ref{f-galaxy-Rvir-frac}.
The slopes in Figs~\ref{f-galaxy-spin-rreff} and \ref{f-galaxy-Rvir-rreff} are not significantly different from zero.

The slopes in Figs~\ref{f-galaxy-spin-frac} and \ref{f-galaxy-Rvir-frac}, $\diffd \lambda/\diffd \haloinvoidfrac = \GalspinFracSlopevalue \pm \GalspinFracSigSlopevalue \pmcosvar \GalspinFracSlopevalueCosmStdDev$ and $\diffd R_{\mathrm{vir}}/\diffd \haloinvoidfrac = \GalRvirFracSlopevalue \pm \GalRvirFracSigSlopevalue$~kpc/$h$, 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.~\ref{f-galaxy-size-frac}, is that the two effects more or less cancel, in contrast to the full-population results shown in Table~\ref{Galaxy-size-table}.

Figure~\ref{f-galaxy-size-rreff} and its fit show that overall, the elaphrocentric distance $r/R_{\mathrm{eff}}$ has only a weak effect on galaxy disk scale lengths $\rdiskscale$.
The halo size and spin parameter both have weak, though apparently again opposite, dependences on $\rdiskscale$, with $\diffd \lambda/\diffd \left(r/R_{\mathrm{eff}}\right) = \GalspinReldistSlopevalue \pm \GalspinReldistSigSlopevalue \pmcosvar \GalspinReldistSlopevalueCosmStdDev$ in Fig.~\ref{f-galaxy-spin-rreff} and $\diffd R_{\mathrm{vir}}/\diffd \left(r/R_{\mathrm{eff}}\right) = \GalRvirReldistSlopevalue \pm \GalRvirReldistSigSlopevalue \pmcosvar \GalRvirReldistSlopevalueCosmStdDev$~kpc/$h$ in Fig.~\ref{f-galaxy-Rvir-rreff}.

In summary, the trends for $\rdiskscale$, $\lambda$, and $R_{\mathrm{vir}}$ found in Table~\ref{Galaxy-size-table} are similar, but strengthened, when void location of a galaxy is quantified by $\haloinvoidfrac$, and insignificant when void location is quantified by $r/R_{\mathrm{eff}}$.
The lack of a significant dependence of these parameters on the elaphrocentric distance, $r/R_{\mathrm{eff}}$, is somewhat surprising, since one might expect $\haloinvoidfrac$ and $r/R_{\mathrm{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 \SSS\ref{s-discuss-sizes}.

%%\left(r/R_{\mathrm{eff}}\right)

\begin{table}
  \centering
  \caption{Median radial $\dot{v}_{\parallel}$ and tangential $\dot{v}_{\perp}$ accelerations at $\rMW = \ElaphrocentreHaloOriginalRadiusMpcvalue$~Mpc/$h$ from the elaphrocentre.
    \label{t-elaphro-accel}}
  $\begin{array}{l c}
    \hline
    & \mbox{median acceleration}\\
    \hline
    \dot{v}_{\parallel} &
    \ElaphroAcciAccRadMedvalue \pm \ElaphroAcciAccRadStderrMedvalue \,\mathrm{km/s/Gyr}\rule{0ex}{2.7ex} \\
    \dot{v}_{\perp}  &
    \ElaphroAcciAccTanMedvalue \pm \ElaphroAcciAccTanStderrMedvalue \pmcosvar \ElaphroAcciAccTanMedvalueCosmStdDev \,\mathrm{km/s/Gyr} \\
    \hline
  \end{array}$
\end{table}

\begin{figure}
  \includegraphics[width=0.9\columnwidth]{fig/accel_rad\revolverparamiRadFracElaphrovalue}
  \caption{Dependence of the radial elaphro-acceleration $\dot{v}_{\parallel}$ as a function of effective void radius $R_{\mathrm{eff}}$, together with a robust linear fit.
    The fit is $\dot{v}_{\parallel} = [(\AccRadZerovalue \pm \AccRadSigZerovalue) + (\AccRadSlopevalue \pm \AccRadSigSlopevalue \pmcosvar \AccRadSlopevalueCosmStdDev)\,h/\mathrm{Mpc}\, R_{\mathrm{eff}}]$~km/s/Gyr.
\label{f-acc-rad-rreff}}
\end{figure}

\begin{figure}
  \includegraphics[width=0.9\columnwidth]{fig/accel_tan\revolverparamiRadFracElaphrovalue}
  \caption{Dependence of the tangential elaphro-acceleration $\dot{v}_{\parallel}$ as a function of $R_{\mathrm{eff}}$, together with a robust linear fit.
    The fit is $\dot{v}_{\perp} = [(\AccTanZerovalue \pm \AccTanSigZerovalue) + (\AccTanSlopevalue \pm \AccTanSigSlopevalue \pmcosvar \AccTanSlopevalueCosmStdDev)\,h/\mathrm{Mpc}\,  R_{\mathrm{eff}}]$~km/s/Gyr.
\label{f-acc-tan-rreff}}
\end{figure}


\subsection{Elaphro-acceleration}  \label{s-results-elaphro}

Elaphro-accelerations were calculated as described in \SSS\ref{s-method-elaphro-acc}, at elaphrocentric distances of $\rMW = \ElaphrocentreHaloOriginalRadiusMpcvalue$~Mpc.
We found ${\ElaphroAcciCountedVoidsPotvalue}$ voids that allowed valid estimates.
The median radial acceleration for a test particle at these positions is given in Table~\ref{t-elaphro-accel}.
The amplitudes of these two values are not directly comparable, because $\dot{v}_{\parallel}$ is the median of signed values, while $\dot{v}_{\perp}$ is non-negative by construction.
The Newtonian estimate for the gravitational pull at $\rMW$ away from the centre of a Milky-Way--like halo of mass $M\tagMW = 10^{12} \mathrm{M}_\odot$, a barycentre, is an inward-pointing acceleration, i.e. $\dot{v}_{\parallel}^{\textrm{MW-halo}} = \MWaccel\, \mathrm{{km/s}/{Gyr}}$.

The median estimate of $\dot{v}_{\parallel}$ given in Table~\ref{t-elaphro-accel} is an outward-pointing elaphro-acceleration to high statistical significance.
The amplitude is more than an order of magnitude weaker than that of the barycentric acceleration towards our canonical Milky-Way--like halo.
Figure~\ref{f-acc-rad-rreff} shows that the full spread of radial elaphro-accelerations is wide, including many negative values, and that dependence on the void effective radius $R_{\mathrm{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 $\dot{v}_{\perp}$ is given in Table~\ref{t-elaphro-accel} and the individual estimates and fit are shown in Fig.~\ref{f-acc-tan-rreff}.
These values are about an order of magnitude higher in amplitude than those of $\dot{v}_{\parallel}$, and similar to that for our canonical Milky-Way--like 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~\ref{f-acc-rad-rreff} \& \ref{f-acc-tan-rreff} do not show significant dependence of $\dot{v}_{\parallel}$ and $\dot{v}_{\perp}$ on the size of a void.


\section{Discussion} \label{s-discuss}

\subsection{Infall rates} \label{s-discuss-infall}

We generally found a lack of statistically significant trends in the two infall parameters on $\haloinvoidfrac$ and $r/R_{\mathrm{eff}}$ (\SSS\ref{s-results-infall}), for those galaxies whose host haloes' infall rates could be fit with an exponential best fit, apart from a modestly significant trend of the amplitude $A$ to increase with greater $\haloinvoidfrac$.
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.
This is something like a Newtonian numerical equivalent of Birkhoff-like \citep{BirkhoffLanger23} or \enquote*{finite infinity} \citep{Ellis84relcosmo,Wiltshire07clocks} arguments for modelling galaxies in isolation from their environment.
Our hypothesis that the void environment helps to form LSBGs by providing slow, weak infall is not supported by our numerical results.


\subsection{Galaxy sizes} \label{s-discuss-sizes}

\begin{figure}
  \includegraphics[width=0.9\columnwidth]{fig/frac_rreff\revolverparamiRadFracElaphrovalue}
  \caption{Relation between fraction in void $\haloinvoidfrac$ and elaphrocentric position $r/R_{\mathrm{eff}}$, showing a robust linear (Theil--Sen) fit.
    The fit is $\lambda = (\FracRreffZerovalue \pm \FracRreffSigZerovalue) + (\FracRreffSlopevalue \pm \FracRreffSigSlopevalue \pmcosvar \FracRreffSlopevalueCosmStdDev) \,r/R_{\mathrm{eff}}$ (the slope is sensitive to cosmic variance).
    Plain text table for this figure: \mbox{\href{\projectzenodofilesbase/voidgals_infall.dat}{\projectzenodoid/voidgals\_infall.dat}}.
\label{f-frac-rreff}}
\end{figure}

Although we found that galaxies in voids tend to be smaller when comparing the overall void to non-void populations (Table~\ref{Galaxy-size-table}), this was not detected in the dependence of the galaxy disk scale length $\rdiskscale$ on $\haloinvoidfrac$, nor on on $r/R_{\mathrm{eff}}$.
We did find significant dependences of the two contributing parameters to $\rdiskscale$, the spin parameter $\lambda$ and the virial radius $R_{\mathrm{vir}}$, on $\haloinvoidfrac$, but insigificant dependence on $r/R_{\mathrm{eff}}$.
It may seem somewhat surprising that these significant dependences on $\haloinvoidfrac$ do not translate into significant dependences on $r/R_{\mathrm{eff}}$.
The explanation most likely lies in the fact that {\revolvername} 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 $\haloinvoidfrac$, while lying, for example, in the far ends of a prolate void, where $r/R_{\mathrm{eff}} > 1$.
It would be useful to quantify the relation between $\haloinvoidfrac$ and $r/R_{\mathrm{eff}}$.

Figure~\ref{f-frac-rreff} shows visually that there is no obvious relation between $\haloinvoidfrac$ and $r/R_{\mathrm{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 $r/R_{\mathrm{eff}} \gtapprox 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 $r/R_{\mathrm{eff}}  \gtapprox 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, $\haloinvoidfrac$ and $r/R_{\mathrm{eff}}$ appear to be statistically independent parameters.
The significant trend of $\rdiskscale$ on $\haloinvoidfrac$, and the fact that $\haloinvoidfrac$ has a more local physical meaning than $r/R_{\mathrm{eff}}$, suggest that $\haloinvoidfrac$ 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 $\haloinvoidfrac$ as a parameter for characterising the void nature of a galaxy.

\subsection{Elaphro-acceleration} \label{s-discuss-elaphro}

We found a significantly non-zero positive median acceleration towards the edges of a void.
However, this median outwards acceleration, $\dot{v}_{\parallel}$ (Table~\ref{t-elaphro-accel}), is of the order of only a few percent of the inward gravitational pull at $\rMW = \ElaphrocentreHaloOriginalRadiusMpcvalue$~Mpc that the source mass excess for our canonical Milky-Way--like host halo would create, $\dot{v}_{\parallel}\tagMW = \MWaccel\, \mathrm{{km/s}/{Gyr}}$.
Moreover, Fig.~\ref{f-acc-rad-rreff} 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 $N$-body realisation itself, as we presented in \SSS\ref{s-results-infall}, at elaphrocentric positions far from the elaphrocentre.

The median tangential acceleration $\dot{v}_{\perp}$ is much higher than the radial acceleration, and might be used to study the higher spin parameter $\lambda$ found for void galaxies when identified by $\haloinvoidfrac$, as shown in Fig.~\ref{f-galaxy-spin-frac}.
Since this is of the same order of magnitude as our canonical radial acceleration, $\dot{v}_{\parallel}\tagMW$, it is likely that $\dot{v}_{\perp}$ could play an important role for galaxies forming in voids.

\subsection{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.
%This would need, for example, $k$-corrections \citep{Wirtz1918} for calculating surface brightness estimates.

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 {\sc inhomog} through the {\sc ramses}/{\ramsesscalavname} front end to check background-independent dynamical properties \citep{Roukema17silvir}, using the relativistic Zel'dovich approximation \citet*{BuchRZA1,BuchRZA2}, could easily be added.
Other options could include using either {\sc gevolution} \citep{AdamekDDK16code} or the fully relativistic {\sc Einstein Toolkit} \citep*{BentivegnaBruni15,Macpherson17}.
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 $N$-body simulation outputs.
The present form of the pipeline assumes {\sc gadget2} format for the $N$-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 $\InfallDecayRate$, 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.

\section{Conclusion} \label{s-conclu}

In this work we have presented a complete, {\em ab initio}, reproducible galaxy formation pipeline starting from a standard post--recombination-epoch spectrum of initial perturbations, aiming to identify key factors in forming low surface brightness galaxies in voids (\SSS\ref{s-method-pipeline}).
We introduced the term {\em elaphrocentre} to clarify its opposite physical nature to the barycentre and we clarified the confusing use of the latter term in void studies (\SSS\ref{s-method-elaphro}).
Our most promising results were not those we expected.

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 -- ${\InfallAmplitude}$ and ${\InfallDecayRate}$.
More specifically, the effects on mass infall for those galaxies with acceptable exponential fits are mostly insignificant, whether characterised by the overall population, by $\haloinvoidfrac$, or by $r/R_{\mathrm{eff}}$ (Table~\ref{Infall-rate-table}, Figs~\ref{f-infall-frac-amp}--\ref{f-infall-rreff-tau}).
The only moderately significant slope, that in Figure~\ref{f-infall-frac-amp}, shows a dependence of ${\InfallAmplitude}$ on $\haloinvoidfrac$ opposite to that expected: void galaxies have slightly higher amplitudes of their best fit infall rate histories.
The modest amplitude of the median acceleration outwards from the elaphrocentre -- the {\em elaphro-acceleration} $\dot{v}_{\parallel}$ (Table~\ref{t-elaphro-accel}) -- is probably the main reason for this.
While the potential wells of barycentres are sharp, since gravitational collapse shrinks the spatial spread of the source material, the potential hills of elaphrocentres are wide, since gravity expands the spatial spread of the early-epoch source material.
With hindsight, our numerical results together with this qualitative interpretation suggest that elaphrocentric gravitational contributions to LSBG formation are likely to be modest on average.

Moreover, the fractional elaphrocentric distance of a void galaxy $r/R_{\mathrm{eff}}$ is, statistically, a less useful independent variable than $\haloinvoidfrac$, the fraction of a galaxy's host halo particles that are identified as being in a single void.
We found that the two parameters show no significant linear relation (Fig.~\ref{f-frac-rreff}).
This is important for observational studies of void galaxies.
The characterisation of galaxies as void galaxies by $\haloinvoidfrac$, which should roughly correspond to a low local dark matter density, or by relative elaphrocentric radius $r/R_{\mathrm{eff}}$, which would require identification of voids in a catalogue, will in general not be expected to give correlated results.

Our galaxy disk scale length results (Table~\ref{Galaxy-size-table}, Figs~\ref{f-galaxy-size-frac}--\ref{f-galaxy-Rvir-rreff}) do not directly show any support for LSBG formation.
This is because the disk scale length $\rdiskscale$, as calculated in Eq.~\eqref{e-defn-spin}, is dominated by the population difference between void and non-void galaxies.
The most massive galaxies and their host haloes (with the largest virial radii $R_{\mathrm{vir}}$) form in the tight knots of the cosmic web.
While this 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.~\ref{f-galaxy-spin-frac}, via the dependence of $\lambda$ on $\haloinvoidfrac$.

We thus have indirect support for a spatial location well inside a void constituting a significant factor in LSBG formation, by giving a higher spin parameter to a host halo than a location in a barycentric location.
This finding of a higher spin parameter for high $\haloinvoidfrac$ is consistent with \citet{Rong2017}'s result that a higher spin is a key feature of UDGs.
A followup step to confirm the consistency will be to add evolutionary population stellar synthesis to our pipeline in order to study the spin parameter of LSBGs identifed by their estimated surface brightness.
For the study of intermediate-mass LSBGs of a fixed host mass, we could reasonably expect that the spin parameter in galaxies whose host haloes are better identified as being located in voids would increase even more strongly with increasing $\haloinvoidfrac$.
Increased $\haloinvoidfrac$ would then appear to play an even stronger role in generating low surface brightness.

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 $\dot{v}_{\perp}$ (Table~\ref{t-elaphro-accel}).
This is typically of the same order of magnitude as the gravitational pull of the source region of a Milky-way--like halo.
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 \citet[][\SSS{}7.1, 7.2]{DiPaoloSalucci20}, 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 \citet{ZhuXuGaspari18Malin1}.

We also found another parameter whose median values are different to very high statistical significance between void and non-void populations: the galaxy collapse epoch.
While we did not expect this to play a major role in LSBG formation, it should.
For a standard spherical collapse model, with a fixed overdensity ratio of collapse compared to the critical density at that epoch, haloes that collapse later will have lower matter densities.
It is credible that the lower densities (by volume) will lead to lower surface densities of the disks that form, and correspondingly lower star formation rates, and, thus, LSBGs.
The median collapse epochs (in standard FLRW cosmological time) that we found were $t^{\mathrm{f}}_{\mathrm{v}}=\FormationTimeVoidsMedvalue \pm \FormationTimeVoidsStderrMedvalue$~Gyr and $t^{\mathrm{f}}_{\mathrm{nv}}=\FormationTimeNonvoidsMedvalue \pm \FormationTimeVoidsStderrMedvalue $~Gyr for the void and non-void galaxies, respectively (\SSS\ref{s-results-infall}).
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.
This should lead to haloes and galaxies forming with lower dark matter and baryonic matter densities (M$_{\odot}$/kpc$^3$) and may lead to lower surface densities (M$_{\odot}$/kpc$^2$) of galaxy disks and lower surface brightnesses (L$_{\odot}$/kpc$^2$) induced by star formation.

This is the second characteristic of void galaxies found in our results that agrees with \citet{Rong2017}, 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.
\citet{Rong2017} found that UDGs at the current epoch have a median age of $7.1$~Gyr compared to typical dwarfs with a median age of $t_0 - 4.2$~Gyr.
For the current epoch estimated at $t_0 \sim 13.8$~Gyr, these correspond to median UDG and ordinary dwarf formation epochs of $6.7$~Gyr and $9.6$~Gyr, respectively.
Since (i) we consider high-mass host haloes and correspondly 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 \citet{Rong2017}'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 \citep{Rong2017} is a promising sign of progress towards a cohesive theory of LSBG formation.

In summary, despite not finding direct numerical evidence for LSBG formation in our overall populations, the higher spin parameter $\lambda$ for the overall population of void galaxies, especially when characterised by $\haloinvoidfrac$, and the later formation epoch of void galaxies, are consistent with voids playing a significant role in LSBG formation.
Moreover, we hope that by providing our complete software pipeline\footnote{DOI-stamped record: {\projectzenodohref}; live {\sc git} repository: {\projectgitrepository}} using the Maneage template that aims at full reproducibility \citep{RougierHinsen2017Repdefn,Akhlaghi2020maneage}, rather than only giving the names and URLs of cosmological softwage 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.

%This work demonstrates how powerful our freshly build pipeline can be.
%We can analyse galaxies in our simulation in full detail and compare their properties based on position, mass or any other criterion.

%% Acknowledgements
\section*{Acknowledgments}
The authors wish to thank Krzysztof Bolejko, Mariana Jaber, Matteo Cinus, Mohammad Akhlaghi, Ra{\'u}l Infante--Sainz and M.\/ Ahsan Nazer for many very useful suggestions.
Work on this paper has been supported by the ``A next-generation worldwide quantum sensor network with optical atomic clocks'' project of the TEAM IV programme of the Foundation for Polish Science co-financed by the European Union under the European Regional Development Fund.
This paper has been supported by Polish MNiSW grant DIR/WK/2018/12.
This paper has been supported by Pozna\'n Supercomputing and Networking Center (PSNC) computational grant 314.
Part of this work was supported by Universitas Copernicana Thoruniensis in Futuro under NCBR grant POWR.03.05.00-00-Z302/17.

\input{tex/build/macros/dependencies.tex}

%% Tell BibLaTeX to put the bibliography list here.
\bibliographystyle{\bibtexbstfile}
\bibliography{\bibtexbibfiletmp}

%% Finish LaTeX
\end{document}

%% This file is part of a paper describing the Maneage workflow system
%%   https://gitlab.com/makhlaghi/maneage-paper
%
%% This file is free software: you can redistribute it and/or modify it
%% under the terms of the GNU General Public License as published by the
%% Free Software Foundation, either version 3 of the License, or (at your
%% option) any later version.
%
%% This file is distributed in the hope that it will be useful, but WITHOUT
%% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
%% FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
%% for more details.
%
%% You should have received a copy of the GNU General Public License along
%% with this file.  If not, see <https://www.gnu.org/licenses/>.