\section{Introduction}

In the $\Lambda$CDM paradigm, structure in the Universe arises from the initial density perturbations of an (almost) homogeneous dark matter distribution. Due to gravitational evolution, this leads to the appearance of collapsed structures, i.e. dark matter haloes. Some of the baryonic matter, following this process, cools down and settles at the centres of the gravitational potentials where it forms galaxies.

This mechanism has been studied through models of so-called spherical collapse \citep{Gunn1972, 1985ApJS...58...39B}, whose main prediction is the existence of a radius within which the material orbiting the halo is completely virialized. In general, this virial radius depends on cosmology and redshift, but both in numerical simulations and observations, fixed overdensity radii are widely used as proxies for this quantity. An example of this is $r_{200m}$, defined as the radius within which the average density is $200$ times the average matter density of the Universe, $\rho_\text{m}$. The corresponding enclosed mass is known as $M_{200m}$.

Halo mass functions constructed with these idealized definitions can capture the effects of cosmology \citep{1974ApJ...187..425P}, the nature of dark matter \citep{2013MNRAS.434.3337A}, and dark energy \citep{Mead_2016} on the growth of structure. In the real Universe, however, this picture is complicated by the triaxiality of haloes \citep{1991ApJ...378..496D, Monaco_1995} and the existence of clumpy (baryonic) substructure \citep{Bocquet_2015}.

Because the process of structure formation is hierarchical, massive haloes contain subhaloes, some of which host galaxies themselves. The resulting clusters of galaxies are the focus of this work. What makes these objects particularly unique is the fact that they are not fully virialized yet. To this day, they are still accreting both ambient material and subhaloes through filamentary structures surrounding them \citep{Bond_1996}. Because of their definition, however, traditional overdensity definitions of mass are not only affected by halo growth, but also by a pseudo-evolution due to  the redshift dependence of $\rho_\text{m}$ \citep{2013ApJ...766...25D}.

\cite{Diemer_2014} and \cite{More_2015} were the first to note that this growth process leads to the formation of a sharp feature in the density profile that separates collapsed and infalling material. This feature, therefore, defines a natural boundary of the halo. The location of this edge, i.e. the splashback radius $r_\text{sp}$, has an obvious primary dependence on halo mass, but also a secondary dependence on accretion rate. While this behaviour can be qualitatively explained using simple semi-analytical models of spherical collapse, none of the analytical models currently proposed \citep{Adhikari_2014, Shi2016} can fully describe its dependency on mass and accretion rate \citep{Diemer_2017b}. Despite this, the corresponding definition of halo mass is particularly suited to define a universal mass function valid for a wide range of cosmologies  \citep{diemer2020universal}.

In this paper, we try to bridge the gap between the theoretical understanding of the splashback feature and observational results, both past and future. The outer edge of clusters has already been extensively measured through different tracers: the radial distribution of galaxies from wide surveys \citep{2016ApJ...825...39M, Baxter_2017, Chang_2018}, but also their velocity distribution \citep{tomooka2020clusters, 2021MNRAS.503.4250F}, and in the weak-lensing signal of massive clusters \citep{Umetsu_2017, Chang_2018, Contigiani_2019b}. Furthermore, forecasts have already set expectations for what will be obtainable from near-future experiments \citep{Fong_2018, xhakaj2019accurately, 2020arXiv201011324W}. Despite the wealth of data and studies, however, not many \emph{splashback observables} have been proposed. The only robust application of this feature found in the literature is related to the study of quenching for newly accreted galaxies \citep{2020arXiv200811663A}.

To achieve our goal, we make use of hydrodynamical simulations of massive galaxy clusters, which we introduce in Section~\ref{sec:Hydrangea}. We focus mainly on $z=0$, but also make use of snapshots at redshifts $z=0.474$ and $z=1.017$. In Section \ref{sec:definition}, we start our discussion by introducing the physical interpretation of splashback and consider the connection between the galaxy and dark matter distributions. We then continue in Section \ref{sec:msrelation} and \ref{sec:redshift}, where we explain how galaxy profiles and weak-lensing mass measurements can be combined to construct a mass-size relationship for galaxy clusters. Finally, we summarize our conclusions and suggest future developments in Section \ref{sec:conclusions}.

\section{Hydrangea}
\label{sec:Hydrangea}

The Hydrangea simulations are a suite of $24$ zoom-in hydrodynamical simulations of massive galaxy clusters ($\log_{10}M_\text{200m}/\mathrm{M}_\odot$ between $14$ and $15.5$ at redshift $z = 0$) designed to study the relationship between galaxy formation and cluster environment \citep{Bahe2017}. They are part of the Cluster-EAGLE project \citep{Bahe2017,Barnes_2017} and have been run using the EAGLE galaxy formation model \citep{Schaye_2014}, which is known to reproduce galaxy observables such as colour distribution and star formation rates. To better reproduce the observed hot gas fractions in galaxy groups, the AGNdT9 variant of this model was used \citep{Schaye_2014}.

The zoom-in regions stretch to between 10 and 30 Mpc from the cluster centre, roughly corresponding to $\lesssim 10 r_{200m}$. For the definition of the cluster centre, in this work, we choose the minimum of the gravitational potential. We note, however, that this choice will not impact our conclusions since we will focus on locations around $r_\text{200m}$. The particle mass of $m \sim 10^{6}~\mathrm{M}_\odot$ for baryons and $m\sim10^{7}~\mathrm{M}_\odot$ for dark matter allows us to resolve galaxy positions down to stellar masses $M_\ast \geq 10^{8}~ \mathrm{M}_\odot$ and total masses $M_\text{sub} \geq 10^{9}~
\mathrm{M}_\odot$, respectively.

In Figure \ref{fig:comparison} we show the log-derivative of the stacked subhalo density $n_s(r)$ at large scales. This is the result of a fit obtained using the model of \cite{Diemer_2014}, and we refer the reader to the aforementioned paper and \cite{Contigiani_2019b} for a detailed explanation of the model and its components. \edit{The choice to employ this profile is based on its ability} to capture the sharp feature visible around $r_\text{200m}$, which is the focus of this work. We optimally sample its $8$-dimensional parameter space using an ensemble sampler \citep{Foreman-Mackey2013}.

In the same plot, we also include the stacked subhalo profile of the accompanying dark matter only (DMO) simulations, initialized with matching initial conditions. The two profiles match almost exactly, suggesting that baryonic effects do not alter this feature to a significant extent \citep[see also][]{2020arXiv201200025O}. While not shown, we report that the same conclusion can be reached by looking at the full matter distribution $\rho(r)$ in the two sets of simulations. Similarly, this feature is also visible in the number density of galaxies, $n_g(r)$. \edit{Due to our focus on all three of these profiles, we choose not to work with background subtracted quantities.}

For reference, we present a full list of the simulated clusters used in this paper and their relevant properties, some of them defined in the following sections, in Table~\ref{tab:clusterlist}.

\begin{table*}
	\centering
	\caption{The Hydrangea clusters used in this paper and their $z=0$ properties. $\Gamma_{0.3}$ is the accretion rate measured between $z=0$ and $z=0.297$. The three splashback radii $r_\text{sp}$, $r_\text{sp}^\text{g}$, $r_\text{sp}^\text{s}$ refer to the splashback radius measured, respectively, in the dark matter, galaxy, and subhalo distributions (see Section~\ref{sec:definition}). For two clusters, CE-28 and CE-18, the radius $r_\text{sp}$ is not used in this work because the dark matter distribution displays a featureless profile at large scales. All quantities are in physical units.}
	\label{tab:clusterlist}
	\begin{tabular}{lccccccc}
		\hline
		Name & $\Gamma_{0.3}$ & $M_\text{200m}$ & $r_\text{200m}$ & $r_\text{sp}$ & $r_\text{sp}^\text{g}$ & $r_\text{sp}^\text{sub}$\\
		 & & $[10^{14}~M_\odot]$ & [Mpc] & [Mpc] & [Mpc] & [Mpc]\\
		\hline
CE-0 & 0.8 & 1.74 & 1.74 & 2.98 & 2.72 & 2.60 \\
CE-1 & 2.0 & 1.41 & 1.63 & 1.71 & 1.56 & 1.79 \\
CE-2 & 0.5 & 1.41 & 1.63 & 2.36 & 3.27 & 2.36 \\
CE-3 & 0.8 & 2.04 & 1.84 & 2.60 & 2.72 & 2.72 \\
CE-4 & 2.8 & 2.19 & 1.89 & 1.63 & 1.87 & 1.79 \\
CE-5 & 2.0 & 2.24 & 1.90 & 2.36 & 2.60 & 2.48 \\
CE-6 & 1.1 & 3.31 & 2.16 & 2.60 & 2.48 & 2.60 \\
CE-7 & 1.2 & 3.39 & 2.17 & 3.13 & 2.60 & 2.85 \\
CE-8 & 1.8 & 3.09 & 2.12 & 2.26 & 2.48 & 2.06 \\
CE-9 & 1.1 & 4.27 & 2.36 & 3.76 & 3.76 & 3.27 \\
CE-10 & 0.8 & 3.55 & 2.21 & 3.13 & 3.13 & 2.98 \\
CE-11 & 1.4 & 4.27 & 2.34 & 3.13 & 2.85 & 2.72 \\
CE-12 & 0.1 & 5.13 & 2.49 & 3.43 & 3.76 & 4.13 \\
CE-13 & 1.5 & 5.25 & 2.52 & 2.26 & 3.13 & 2.72 \\
CE-14 & 2.1 & 6.17 & 2.66 & 2.60 & 2.72 & 2.48 \\
CE-15 & 4.2 & 6.76 & 2.73 & 1.96 & 2.26 & 2.48 \\
CE-16 & 2.7 & 7.59 & 2.84 & 1.42 & 4.13 & 3.43 \\
CE-18 & 1.1 & 9.12 & 3.03 & - & 3.76 & 3.76 \\
CE-21 & 3.7 & 12.30 & 3.34 & 2.36 & 2.85 & 2.60 \\
CE-22 & 1.5 & 16.98 & 3.72 & 4.53 & 4.53 & 4.33 \\
CE-24 & 1.5 & 15.49 & 3.61 & 3.27 & 3.27 & 4.33 \\
CE-25 & 3.4 & 19.05 & 3.87 & 3.43 & 3.43 & 3.43 \\
CE-28 & 1.9 & 21.88 & 4.06 & - & 3.94 & 3.27 \\
CE-29 & 3.5 & 32.36 & 4.61 & 3.94 & 4.13 & 3.94 \\
		\hline
	\end{tabular}
\end{table*}

\begin{figure}
    \centering
    \includegraphics[width=0.46\textwidth]{1comparison.pdf}
    \caption{The splashback feature visible in the average subhalo distribution of simulated high-mass clusters. We extract the logarithmic slope by fitting a smooth profile to the mean of the Hydrangea profiles rescaled by $r_\text{200m}$. We perform this operation both on the hydrodynamical simulations (Hydro) and their dark matter only counterparts (DMO). The minimum around $r_{200m}$ marks the halo boundary, and this figure highlights the lack of baryonic effects on the location or depth of this feature. The two logarithmic slope profiles are consistent with each other at the $1$ per cent level.}
    \label{fig:comparison}
\end{figure}

\section{Splashback}
\label{sec:definition}

\subsection{Definition}

For haloes that continuously amass matter, material close to its first apocentre piles up next to the edge of the multi-stream region, where collapsed and infalling material meets \citep{Adhikari_2014}. A sudden drop in density, i.e. the feature visible in the profiles of Figure~\ref{fig:comparison}, is associated with this process.

This intuitive picture leads to three characterizations of the splashback radius, depending on the approach used to measure or model it:
\begin{enumerate}
\item The location of the outermost phase-space caustic.
\item The point of steepest slope in the density profile.
\item The apocentre of recently accreted material.
\end{enumerate}
While these definitions have all been previously hinted at in the introduction, in this section we explicitly present them and discuss the connections existing between them. This also justifies our adopted definition, based on the density profile.

The first definition is clearly motivated in the spherical case but fails once it is applied to realistic haloes. The presence of angular momentum and tidal streams from disrupted subhaloes \citep[see, e.g.,][]{2011MNRAS.413.1419V}, smooth out this feature and make its description murky. The second definition was the first suggested in the literature. Introduced by \cite{Diemer_2014}, it is based on the study of dark matter profiles in N-body simulations and has been linked to the first, more dynamical, definition \citep{Adhikari_2014, Shi2016}. The third was first suggested by \cite{Diemer_2017}, who showed that this location can be calibrated to the second one \citep{Diemer_2017b} by choosing specific percentiles of the apocentre distribution.

To clarify the relationship between the outermost caustic and apocentre, it is educational to use a self-similar toy model based on \cite{Adhikari_2014} to show the phase-space distribution of a constantly accreting halo with an NFW-like mass profile \citep{Navarro_1997}.

In the absence of dark energy, we follow the radial motion of particles,
\begin{equation}
    \label{eq:motion}
    \ddot{r} = \frac{GM(<r, t)}{r^2},
\end{equation}
 between their first and second turnaround in the mass profile:
\begin{equation}
    \label{eq:motion2}
    M(r, t) = M(R, t) \frac{f_\text{NFW}(r/r_\text{s})}{f_\text{NFW}(R/r_s)}.
\end{equation}
We impose that the total mass evolves as $M(R, t) \propto t^{2\Gamma/3}$, $R\propto t^{2(1+\Gamma/3)/3}$, and the dimensionless NFW profile is defined as: $f(x) =\log(1+x)-x/(1+x)$. In this set of equations, $\Gamma$ is the dimensionless accretion rate, $R$ represents the turnaround radius, and the scale parameter $r_s$ is defined by the infall boundary condition
\begin{equation}
    \label{eq:boundary}
    \frac{d\log M}{d\log r}(R) = \frac{3\Gamma}{3+\Gamma}.
\end{equation}
This condition, combined with the turnaround dynamics, imposes that $M(R, t)\propto (1+z)^{-\Gamma}$ \citep{1984ApJ...281....1F}.

We point out that the dependence on the time-sensitive turnaround properties $M(R, t), R(t)$ can be factored out from the equations above, meaning that the entire phase-space at all times can be obtained with a single numerical integration.

In Figure~\ref{fig:gammavel} we show the result of this calculation, denoting the location of the outermost caustic as $r_\text{sp}^{c}$. The caustic is formed by the outermost radius at which shells at different velocities meet ($r/r_\text{sp}^c=1$ in the plot) and the location of shells at apocenter is defined by the intersection between the zero-velocity line and the phase-space distributions. From the figure, two things are noticeable: material at $r_\text{sp}^c$ has not reached its apocentre yet, and the ratio between these two locations depends on the accretion rate.

It is beyond the scope of this work to quantify this dependence since it depends heavily on the mass profile inside $r_\text{sp}^{c}$. Qualitatively, however, the difference between caustic and apocentre is easy to understand once the dynamical nature of this feature is considered: the halo is growing in size, and while some material is now reaching its apocentre, mass accreted more recently has the chance to overshoot it and form the actual caustic. In a static picture, this would not be the case.

In realistic haloes, this dependence on accretion rate is only one of many factors that biases and adds scatter to the relationship between the halo boundary and apocentres. Other factors include, e.g., non-spherical orbits and the presence of multiple accretion streams. Despite this, \cite{Diemer_2017} has shown that there is a clear link between the apocentre distribution and splashback. The percentile definition introduced there is particularly suited to theoretical investigations, but its usefulness in the very low-$\Gamma$ regime is still uncertain \citep{Mansfield_2017, xhakaj2019accurately}, and it has not been explored in the presence of modifications of gravity \citep{Adhikari_2018, Contigiani_2019}.

\begin{figure}
    \centering
    \includegraphics[width=0.46\textwidth]{2accretionratedep.pdf}
    \caption{The phase-space structure of accreting dark matter haloes depends on the accretion rate $\Gamma$. We employ a toy model of spherical collapse to describe the multi-stream region of NFW-like haloes. The figure shows that the material at the outermost caustic, $r_\text{sp}^c$, is not necessarily at apocentre (where $v=0$) and that the ratio of these two radii is a function of accretion rate. For ease of readability, we have rescaled the coordinates by $r_\text{sp}^c$, and the velocity of collapsed material at this point. }
    \label{fig:gammavel}
\end{figure}

For this work, we define the splashback radius as the location of the steepest slope as defined by a profile fit. In Table~\ref{tab:clusterlist} we report, for each cluster, this radius measured in the distribution of galaxies, subhaloes, and total matter ($r_\text{sp}^g, r_\text{sp}^s, r_\text{sp}$). The model is a modified Einasto profile \citep{Einasto1965} with the addition of a power-law to take into account infalling material \citep{Diemer_2014}. Regarding the goodness of fit, we find that up to and around $r_{200m}$ the standard deviation of the residuals is of order $10$ per cent. On the other hand, the presence of substructure superimposed on a shallow density profile results in normalized residuals of order $50$ per cent in the outer regions.

To further justify our approach, we show in Figure~\ref{fig:fit} how this simple definition of splashback radius is able to capture the phase-space boundary of different haloes, even when a sudden drop in density is absent. The main benefit of this definition is that it avoids the arbitrariness of the apocentre definition, or the bias induced by multiple caustics in the minimum slope definition \citep{Mansfield_2017}. Its main caveats, however, are that 1) it is computationally expensive since it requires high-resolution simulations and a multi-parameter fit procedure, and 2) it might not apply to low-mass clusters and galaxy groups. We leave this last question open for future investigations.

We wrap this subsection up by stressing that this definition of ``the'' splashback radius is, like any other, useful only to study its correlation with other properties, or quantify the impact of different physical processes. While the flexibility of the chosen model is not surprising given the number of free parameters, the clear connection between the phase-space and the log-derivative in individual haloes is a powerful and seemingly general result. Ultimately, however, the observational results focus on stacked projected density profiles, and so should the predictions.

\begin{figure}
    \centering
    \includegraphics[width=0.47\textwidth]{3fitexample.pdf}
    \caption{Fitting simulated subhalo profiles with a smooth model. In the top panels, we show the radial subhalo distributions of two clusters (CE-16, left and CE-9, right), together with the best-fit profiles used to reconstruct the log-derivative. In the bottom panels, we show how the inferred location of the log-derivative minimum (vertical line) identifies the phase-space edge of relaxed (left) and perturbed (right) galaxy clusters. In the phase-space plots, the cluster on the left is formed by collapsed particles, while the stream visible on the lower right is infalling material. The right panels demonstrate how our approach is effective even in the presence of an on-going merger when the splashback feature is not visible as a sharp transition in the density profile.}
    \label{fig:fit}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=0.46\textwidth]{4gammacorr.pdf}
    \caption{The splashback radius and its correlation with the accretion rate. The ratio between the splashback radius and the $200$m overdensity radius correlates with the accretion rate. We show that this correlation exists for the clusters studied in this work and compare it to the relations obtained in three other studies (see text for references). }
    \label{fig:gammarsp}
\end{figure}

\subsection{Accretion}

It is well established \citep{Diemer_2014, More_2015,Mansfield_2017,diemer2020part} that the location of the halo boundary correlates with the accretion rate
\begin{equation}
    \Gamma_{0.3} = \frac{\Delta \log M_\text{200m}}{\Delta \log (1+z)}.
\end{equation}
In this work, this ratio is calculated in the redshift range $z=0$ to $z=0.293$, since this time interval roughly corresponds to one crossing time for all clusters considered here, i.e. how long ago the material currently at splashback has been accreted \citep{Diemer_2017}. Although this choice is partially arbitrary, we have investigated the dependence of our results on the redshift upper limit and we have verified that our main conclusions are not affected.

The archetypical relation demonstrating this idea is plotted in Figure~\ref{fig:gammarsp}, where we have also included the relations found in \cite{More_2015}, \cite{Diemer_2017b}, and \cite{diemer2020part}, to provide additional context. We find good agreement, even though a perfect match is not necessarily expected. The Hydrangea clusters represent a biased sample, selected to be mostly isolated at low redshift \citep{Bahe2017}. \edit{While the effect of this selection on the accretion rate distribution is not fully known, we} show below that a connection between cluster environment and this quantity exists, and the presence of mergers might therefore influence it. \edit{This is not surprising since a connection between accretion and large-scale bias is already known \citep[e.g., ][]{2010MNRAS.401.2245F}.}

We show this relationship explicitly in Figure~\ref{fig:slope} by using one of the parameters of the profile model. As visible in the figure, the power-law index of material outside of splashback correlates with the accretion rate. \edit{We find that this is true for both subhaloes and galaxies and that the difference between the two is consistent with sample variance}.

To try and explain this behaviour, we use a fully consistent model of spherical collapse introduced by \citet{1985ApJS...58...39B}, which was also used in \citet{Contigiani_2019}. The set-up of this toy model is the same as what is shown in Equation~\eqref{eq:motion}, but with a mass profile that also needs to be solved for. Starting from an initial guess for $M(r, t) = \mathcal{M}(r/R(t))$, orbits are integrated and their mass distribution is calculated. Iterating this process multiple times returns a self-similar density profile and orbits consistent with each other.

The result of this calculation is also shown in Figure~\ref{fig:slope}. Because the mass-profile prediction is not a power law, we plot a filled line displaying the range of logarithmic slopes allowed between $r_\text{sp}^\text{c}$ and $2r_\text{sp}^\text{c}$. The fact that this prediction is not a function of accretion rate implies that the correlation between the slope and the accretion rate seen in the simulations is not purely dynamical, and suggests a connection between the cluster environment and accretion rate.

We stress here that previous splashback works have mostly focused on stacked halo profiles, for which the expectation of the spherically symmetric calculation shown above is roughly verified, even in the presence of dark energy \citep{Shi2016}. We also recover this result for our sample (see the star symbol in Figure~\ref{fig:slope}), but we point out that this is a simple conclusion. Because Newtonian gravity is additive, stacking enough clusters should always recover the spherically symmetric result. Despite this, we also note that results from the literature do not always agree with this prediction. However, we do not linger on these discrepancies since 1) this was never the focus of previous articles, and 2) different methods to extract the power-law have been employed.

\subsection{Anisotropy}

This departure from the spherical case implies that anisotropies play a role in shaping the accretion rate $\Gamma$. To study the impact of accretion flows on the cluster boundary, we study $72$ sky projections of the Hydrangea clusters ($3$ each, perpendicular to the $x, y,$ and $z$ axes of the simulation boxes) and rotate them to align the preferred accretion axes in these planes. For each projection, we define this direction $\theta\in (-\pi/2, \pi/2)$ in two ways: 1) to capture the filamentary structure around the cluster between $r_\text{200m}$ and $5r_\text{200m}$, we divide the subhalo distribution in $20$ azimuthal bins and mark the direction of the most populated one, and 2) to capture the major axis of the BCG, we use unweighted quadrupole moments of the central galaxy's stellar profile within $10$ kpc from its centre. The mean projected distributions according to these two methods are presented in the left and right top panels of Figure~\ref{fig:filaments}, respectively.

\begin{figure}
    \centering
    \includegraphics[width=0.47\textwidth]{5power.pdf}
    \caption{The distribution of subhaloes and galaxies outside the cluster edge as a function of accretion rate. Faster growing haloes display a more concentrated distribution of satellites outside of their boundary. This behaviour seen in individual clusters is not explained by simple models of spherical collapse (blue shaded area), but the average profile (marked by a star) matches the expectation. This suggests that non-isotropic processes shape this relation.}
    \label{fig:slope}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=0.47\textwidth]{6filaments.pdf}
    \caption{The impact of filaments and accretion flows on the cluster's edge. We rotate the $2$D subhalo distributions of different clusters to align their accretion axes. The top panels show the resulting mean distributions in a square region of size $5r_\text{200m}$ obtained with two definitions of this direction: one based on the presence of filaments outside $r_\text{200m}$ (left), and one based on the central galaxy's major axis (right). The first one better identifies the filamentary structures around the clusters, but the second one is closer to what can be observed. In the bottom panels, we show how the inferred $3$-dimensional logarithmic slope inside the quadrants aligned with the accretion direction (darker shade) differs from the profile outside (lighter shade). The results from the bottom panel imply that the central galaxy's major axis traces the direction of infalling material.}
    \label{fig:filaments}
\end{figure}

Looking at the top-left panel of the figure, it is not surprising that filamentary structures of the cosmic web are visible around the central cluster -- this is by construction. \edit{Because of the higher contrast between outside and inside regions, the subhalo distribution exhibits a sharper feature in the directions pointing towards voids (see central panel of Figure~\ref{fig:filaments}).} More surprisingly, however, these same traits are also noticeable in the mean distributions aligned according to the central galaxy's axis (see lower panel).

This result implies that the distribution of stellar mass within the central $10$ kpc of the cluster contains information about the distribution of matter at radii which are a factor $10^2$ larger. In fact, the connection between the shape of the dark matter halo and the ellipticity of the brightest cluster galaxy (BCG, which is also the central galaxy for massive galaxy clusters) is known \citep{Okumura_2009, Herbonnet_2019, 2020MNRAS.495.2436R}. And, similarly to other results \citep{Conroy_2007, 2007MNRAS.375....2D}, the Hydrangea simulations predict that the stellar-mass buildup of the BCG is driven by the stripping of a few massive satellites after their first few pericentre passages (Bah\'{e} et al., in prep.). Because these galaxies quickly sink to the centre, the material they leave behind is, therefore, a tracer of their infalling direction.

\section{The mass--size relation}
\label{sec:msrelation}

In our sample, we find that the splashback feature seen in the galaxy, subhalo, and total matter profiles are all at the same location. The mean fractional difference between any two of $r_\text{sp}^g, r_\text{sp}^s$, or $r_\text{sp}$ is consistent with zero, with a mean standard deviation of $3$ per cent. We also verified that this statement is unaffected by cuts in subhalo mass or galaxy stellar mass. \edit{Due to the limited size of our sample, the effects of dynamical friction on the distribution of high-mass subhaloes are not visible \citep{2016JCAP...07..022A}.}

We emphasize, however that this does not mean that galaxy selection effects have no impact on these quantities. For example, it is an established result, both in the Hydrangea simulations \citep{oman2020homogeneous} and in observations \citep{2020arXiv200811663A}, that the location of a galaxy in projected phase-space correlates with its colour and star-formation rate. This is because a red colour preferentially selects quenched galaxies that have been orbiting the halo for a longer time.

Until their first apocentre after turnaround, galaxies act as test particles orbiting the overdensity as the halo grows in mass. In the standard cold dark matter paradigm, based on a non-interacting particle, it is not surprising then that the edge formed in their distribution is identical to the one seen in the dark matter profile. It should be noted, however, that this is not necessarily true in extended models in which dark matter does not act as a collisionless fluid.  Due to their infalling trajectories, the distribution of galaxies will always display a splashback feature, even if the dark matter profile does not exhibit one.

In the cold dark matter scenario, our result implies that galaxies can be used to trace the edge of clusters. We note, in particular, that this measurement has already been performed several times using photometric surveys \citep{Baxter_2017, Nishizawa_2017, Chang_2018, Z_rcher_2019, Shin_2019}. Furthermore, due to the large number of objects detected, galaxy distributions obtained through this method offer the most precise measurements of splashback. The accuracy of the results, however, depends heavily on the details of the cluster finding algorithm \citep{Busch_2017, Shin_2019}.

With this in mind, we build an observational mass-size relation between the location of this feature in the galaxy distribution ($r_\text{sp}^g$) and the mass enclosed within it ($M_\text{sp}^g$). In Figure~\ref{fig:msscaling} we present the correlation between the two for the Hydrangea clusters. Because the splashback radius is roughly located at $r_\text{200m}$ (see Figure~\ref{fig:gammarsp}), this relationship can be understood as a generalization of the virial mass-radius relation, where we have introduced a dependence on accretion rate. Surprisingly, we find that the dependence on $\Gamma_{0.3}$ is well captured by a simple form:
\begin{equation}
    \label{eq:scaling}
    \frac{M_\text{sp}}{r_\text{sp}^3 } \propto (1+\Gamma_\text{0.3})^{\beta}.
\end{equation}
While we do not constrain $\beta$ precisely, we find that $\beta = 1.5$ provides an adequate fit by reducing the total scatter from $0.25$ dex to about half of this value. This choice of exponent and functional form is supported by the model of self-similar collapse used for Figure~\ref{fig:slope}, where we find that a power-law $\beta=1.45$ fits this relationship with the same precision as the exponential functions calibrated to numerical simulations shown in Figure~\ref{fig:gammarsp}. For a more extensive comparison with these predictions, we refer the reader to Section~\ref{sec:conclusions}.

The virial relation is a trivial connection between the mass and size of haloes based on an overdensity factor, but its observational power is limited by the fact that these masses are usually extracted from parametric fits to weak-lensing profiles that do not extend to the respective overdensity radii. Because of this, the overdensity masses have a strong dependence on the assumed mass-concentration relation \citep[see, e.g.][]{Umetsu_2020}. The splashback feature, on the other hand, naturally predicts a mass-size relation for galaxy clusters and does so without the need for external calibrations.

In Figure~\ref{fig:msscaling} we also plot the expected change in this relation to due modifications of gravity. We use the symmetron gravity model of \cite{Contigiani_2019} with parameters $f=1$ and $z_\mathrm{ssb}=1.5$, and assume that the change affects only the splashback radius and not the mass contained within it. The exact result depends on the theory parameters, but the expected change in this relation is around $0.15$ dex.

Experimentally, we argue that this relation can be probed using a combination of galaxy density profiles (to extract $r_\text{sp}^g$) and weak lensing measurements. Aperture masses \citep{2000ApJ...539..540C}, in particular, can be used to extract in a model-independent fashion the average projected mass within a large enough radius. If necessary, the aperture mass can also be deprojected to obtain a low-bias estimate \citep{Herbonnet_2020}.

\begin{figure}
    \centering
    \includegraphics[width=0.47\textwidth]{7scaling.pdf}
    \caption{ The mass size relation of galaxy clusters. In the top panel, we show how the size of the cluster boundary seen in the galaxy distribution, $r_\text{sp}^{g}$, scales with its enclosed mass, $M_\text{sp}^g$. In the same panel we also show the median relation in Equation~\eqref{eq:scaling} obtained for $\beta = 0$ and how modifications of gravity are expected to affect this relation (blue dashed line, see text for more details). In this relation, a secondary dependence on the accretion rate $\Gamma_{0.3}$ is a source of scatter that can be captured if $\beta\neq 0$. As visible in the residuals in the bottom panel, a simple power-law form well reproduces this dependence. In the considered sample, we find that half of the total scatter ($0.25$ dex) is due to the mass accretion rate distribution.}
    \label{fig:msscaling}
\end{figure}

\iffalse
\begin{figure}
    \centering
    \includegraphics[width=0.47\textwidth]{scaling.pdf}
    \caption{The observed mass-size relation of galaxy clusters. For each cluster in our simulation, we take three projections and  }
    \label{fig:projection}
\end{figure}
\fi

\section{Redshift evolution}
\label{sec:redshift}
So far, we have only considered the simulation predictions at $z=0$. In this section, we extend our analysis to higher redshifts by exploring two other snapshots of the Hydrangea simulations at $z=0.474$ and $z=1.017$.

At these higher redshifts, we find that the scatter in the splashback relation for individual haloes is large. This is visible in Figure~\ref{fig:gammarspz}, where we plot the equivalent of Figure~\ref{fig:gammarsp} for these two snapshots. We recover the general result of \cite{diemer2020part} that the average values of $r_\mathrm{sp}/r_\mathrm{200m}$ and $\Gamma$ should be higher at early times, but the correlation between the two is completely washed out by $z=1$. We connect this to three causes: 1) The fixed time interval between the snapshots does not allow us to reliably estimate $\Gamma$ at higher redshift when the crossing times are smaller. 2) The lower number of resolved galaxies and subhaloes means that the residuals of the individual profile fits are larger around the virial radius. And finally, 3) the higher frequency of mergers at high redshift means that the numbers of haloes with profiles not displaying a clear splashback feature increases.

We find that Equation~\ref{eq:scaling} is still valid, even if our ability to constrain the scatter at high redshift is impeded by the sample variance. Furthermore, we report that the splashback overdensity $M_\text{sp}/r_\text{sp}^3$ has a redshift dependence. Or, in other words, that the logarithmic zero-point that was not specified in Equation~\ref{eq:scaling} is a function of redshift. Not accounting for the $\Gamma$ dependence, our best fit values for the logarithm of the average overdensity $\log_{10}(M_\text{sp}/\text{M}_\odot)-3\log_{10}( r_\text{sp}/\text{Mpc})$ are $[13.3, 13.8, 14.1]\pm 0.3$ at redshifts $[0, 0.5, 1]$.

Regarding the anisotropy in the splashback feature due to filamentary structures, we report that this phenomenon exists also at high redshift. In Figure~\ref{fig:filamentsz} we compare the sky-projected subhalo profiles $\Sigma_s(R)$ towards different directions, similarly to what we have done for Figure~\ref{fig:filaments}. In this case, however, we explicitly discuss the connection with observations by plotting directly the ratio of the density profiles inside quadrants oriented towards and perpendicular to the two accretion directions, instead of focusing on the result of the profile fits. The mean and variance of these ratios are calculated assuming that the different projections are independent. We find that the orientation of the major axis of the brightest cluster galaxy does not correlate with a splashback anisotropy at $z=1$.  This is because, in most cases, the identification of a central, brightest galaxy is not straightforward at this redshift. At early times, the future central galaxy is still in the process of being created from the mergers of multiple bright satellites located close to the cluster's centre of potential.

To conclude this section, we point out that in the region around $r_\mathrm{200m}$, the difference between the profiles perpendicular and parallel to the central galaxy's major axis is about $10$ per cent at redshift $z\lesssim0.5$. This departure is well within the precision of galaxy profiles extracted from large surveys \cite[e.g.][]{2020arXiv200811663A}. Therefore, this measurement might already be possible using such catalogues.

\begin{figure}
    \centering
    \includegraphics[width=0.46\textwidth]{8gammacorrz.pdf}
    \caption{The splashback radius and its correlation with the accretion rate as a function of redshift. This plot is an extension of Figure~\ref{fig:gammarsp} for redshifts $z=0.474$ (orange crosses) and $z=1.017$ (light blue plus symbols). The ratio between the splashback radius and the $200$m overdensity radius should correlate with the accretion rate $\Gamma$, but for the Hydrangea snapshot at $z=1.017$ the large sample variance washes out this correlation. Despite this, we still recover the expectation of previous results (plotted lines), a larger average $r_\mathrm{sp}/r_\mathrm{200m}$ at higher redshift. }
    \label{fig:gammarspz}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=0.46\textwidth]{9filamentsz.pdf}
    \caption{The impact of filaments and accretion flows on the outer density profile of massive haloes as a function of redshift. We plot the mean value and variance of the ratio between the $2$D subhalo distributions in quadrants perpendicular ($\Sigma_\perp$) and parallel ($\Sigma_\parallel$) to the accretion direction defined through two tracers. This ratio is closely related to what can be measured in observations. While the difference in profile towards and away from filamentary structures is visible at all redshifts, the orientation of the central galaxy is not a good tracer of the splashback anisotropy at high redshift. This is because the central BCG is still forming and its orientation is not yet finalized.}
    \label{fig:filamentsz}
\end{figure}

\section{Discussion and conclusions}
\label{sec:conclusions}

On its largest scales, the cosmic web of the Universe is not formed by isolated objects, but by continuously flowing matter distributed in sheets, filaments and nodes. For accreting (and hence non-virialized) structures such as galaxy clusters, the splashback radius $r_
\text{sp}$ represents a physical boundary motivated by their phase-space distributions. \edit{To exploit the information content of this feature, in this paper we have introduced and studied two observable quantities related to it. }

\edit{First, we have shown that the full galaxy profile can be used to define a cluster mass, i.e. the mass within $r_\mathrm{sp}$. This is an extension of the traditional approach of using richness as mass proxy \citep[see, e.g., ][]{Simet_2016}.} Because of the dynamical nature of the equivalent feature in the dark matter profile, we conclude that, observationally, the splashback feature in the galaxy profile \emph{defines} the physical halo mass. Moreover, we have shown here that the natural relation between the mass and size of haloes according to this definition (see Fig.~\ref{fig:msscaling}) can be used to constrain new physics at cluster scales. Because this boundary is delimited by recently accreted material, we found that a majority of the scatter in this mass-size relation can be explained through a secondary dependence on accretion rate $\Gamma$.

\edit{Second, we have explored how this connection to the accretion rate might be interpreted as a connection between the geometry of the cosmic web and how clusters are embedded in it.} The relation between the two is made explicit in Figure~\ref{fig:slope}, Figure~\ref{fig:filaments}, and Figure~\ref{fig:filamentsz}. In these figures, we have investigated how the cluster environment affects both the halo growth and the stellar distribution of the central galaxy. This information, combined with the scatter of the mass-size relation, can therefore be used as a consistency check for any property that claims to select for accretion rate.

\subsection{The role of simulations}
\label{sec:sims}

In the last few years, the study of the splashback feature has evolved into a mature field both observationally and theoretically. We use this section to discuss explicitly the connection between the two, in light of this work and its connection to previous endeavours.

\edit{In the context of splashback, simulations have guided the formulation of theoretical principles and hypotheses. However, as more measurements become viable, it becomes necessary to provide clear and powerful observables.} Following this spirit, we used high-resolution hydrodynamical simulations to explore directly the connection between measurements based on sky-projected galaxy distributions and theoretical predictions.

Our conclusions regarding the mass-size relation and its redshift evolution are similar to the results of \cite{Diemer_2017b} and \cite{diemer2020part}, which are based on more extensive N-body simulations. For the sake of completeness, it important to note that, in the same papers, it was also found that the splashback overdensity is not universal, but has both a mild dependency on $M_\mathrm{200m}$ and a strong dependency on cosmology, especially at low redshift ($z\approx0.2$). Due to our limited sample, we are clearly unable to model these effects in this work. Nonetheless, we point out that our goal here is to construct a pure splashback scaling-relation based on galaxy profiles and weak lensing mass measurements. Every other dependency, if present, should be captured either as additional scatter or through different parameter values.

We also point out that that these previous works are based on the apocentre definition of splashback (see Section~\ref{sec:definition}). In contrast, we defined the splashback radius as the point of steepest slope according to a model fit to the density profiles of galaxy clusters. While we do not necessarily expect the two definitions to differ, our choice is based on its connection to observations, and the desire to highlight the fact that the splashback radius is not only some abstract halo property but can be defined as a characteristic of individual profiles, such as, e.g., the concentration parameter \citep{Navarro_1997}.

\edit{An alternative method employed by other studies \citep{Diemer_2014, Mansfield_2017, xhakaj2019accurately} to obtain a measure of $r_\mathrm{sp}$ makes use of the minimum of the logarithmic slope in smoothed profiles. While this approach is much faster than profile fitting when only $r_\mathrm{sp}$ is of interest, it does not describe the full shape visible in Figure~\ref{fig:comparison}. In particular, a model that captures the width of this feature is necessary to define the slope of the outer region without an arbitrary choice of which radial scales to consider. Because the model used here contains an asymptotic outer slope, this definition is unique.}

Our decision has, of course, its drawbacks. The versatility of the fitted model is necessary to capture the variance of the individual profiles, but the resulting intrinsic scatter is large and not the best suited to study tight splashback correlations (such as Figure~\ref{fig:gammarsp}). At the same time, the large parameter space might also be seen by some as a chance to study a multitude of correlations between different model parameters. However, we resist this temptation, as inferences based on such correlations might say more about the particular model employed than provide any physical information.

A subtler difference between our method to characterize splashback and other ones present in the literature is related to the definition of spherical density profiles. \cite{Mansfield_2017} and \cite{2020MNRAS.tmp.3386D} found that the most successful method to achieve a clear splashback feature for individual haloes is to measure the median profile along multiple angular directions. In light of the results of Figure~\ref{fig:filaments}, we argue that the distribution of splashback as a function of direction is skewed by the presence of a few dense filaments and hence the difference between a median and mean splashback can be substantial. Therefore, we stress that future works should exercise caution when employing such methods. The use of median profiles smooths substructure by focusing on the halo boundary in the proximity of voids, but because this process is itself correlated with the halo growth rates (see Figure~\ref{fig:slope}), the connection with observations is not as simple as one might expect.

\subsection{Next steps}

Because in this work we have focused only on high-mass objects ($M\sim10^{14.5}~M_\odot$), a natural future step is to investigate if the results apply also in other regimes. For example, a larger sample over a wide range of masses and redshifts is required to confirm the simple form of Equation~\eqref{eq:scaling} and verify if it applies to lower mass groups ($M\sim10^{13.5}~M_\odot$).

Exploring a wider range in mass, both in observations and simulations, can also be used to confirm a key prediction: because the median accretion rate is expected to be a function of mass and redshift \citep{More_2015}, we expect the mass-size relation for an observed halo sample to not necessarily follow a simple form.

Finally, we point out that our results encourage a concentrated effort towards understanding the relationship between cluster environment and splashback. What is discussed in Figure~\ref{fig:slope}, \ref{fig:filaments} and \ref{fig:filamentsz} suggests that the connection between accretion-flows, filaments, and cluster boundary is not a simple one. To better understand this process, it will become necessary to complement the usual \emph{inside-out} theoretical approaches to splashback, that look at haloes to define their boundaries, with \emph{outside-in} approaches, that connect the cosmic web to its nodes. In this context, the amount of splashback data gathered by projects such as the Kilo Degree Survey \citep{2013ExA....35...25D}, Dark Energy Survey \citep{2005astro.ph.10346T}, and, in the future, LSST \citep{2009arXiv0912.0201L} and \emph{Euclid} \citep{laureijs2011euclid}  will provide a powerful probe for the study of structure formation.

\section*{Acknowledgements}
We thank Benedikt Diemer for providing valuable feedback on the manuscript.
OC is supported by a de Sitter Fellowship of the Netherlands Organization for Scientific Research (NWO). HH acknowledges support from the VICI grant
639.043.512 from NWO. YMB acknowledges funding from the EU Horizon 2020 research and innovation programme under Marie Sk{\l}odowska-Curie grant agreement 747645 (ClusterGal) and the NWO through VENI grant 639.041.751. The Hydrangea simulations were in part performed on the German federal maximum performance computer ``HazelHen'' at the maximum performance computing centre Stuttgart (HLRS), under project GCS-HYDA / ID 44067 financed through the large-scale project ``Hydrangea'' of the Gauss Center for Supercomputing. Further simulations were performed at the Max Planck Computing and Data Facility in Garching, Germany.
\iffalse
\section*{Data Availability}

The inclusion of a Data Availability Statement is a requirement for articles published in MNRAS. Data Availability Statements provide a standardised format for readers to understand the availability of data underlying the research results described in the article. The statement may refer to original data generated in the course of the study or to third-party data analysed in the article. The statement should describe and provide means of access, where possible, by linking to the data or providing the required accession numbers for the relevant databases or DOIs.

\fi

\bibliographystyle{mnras}


\subsection{I. Introduction}
Complex networks, especially networks with hubs, have been studied
extensively in the last two decades. The degree distribution in the networks
is important, because it is related with the features that define
  scale-free networks. Scale-free networks have hubs with a
  large number of links,
  and they are crucial in many applications.
  Complex networks can be applied in multidisciplinary areas such as sociology,
  social psychology, ethnology, and economics \cite{nw,MN}.
  In statistical physics, complex systems can be described in terms of such
  networks. Studies on such topics have led to the development of
new fields of research such as sociophysics \cite{galam} and
econophysics \cite{Cont,Egu,Stau,Curty,nuno}, where financial markets
and opinion dynamics on networks are studied \cite{net, SF, DN}.
As a model for scale-free network, Barab\'{a}si and Albert (BA)
  model is well known \cite{BA}. It is an evolving network model and new nodes are
  linked with the network based on the preferential attachment process, where
  the probability that a node is selected by a new node is proportional
  to its degree. In addition, by introducing the fitness of each node
  to modify the probability, the relation with Bose statistics has been found \cite{BA1,BA2}.
  The relation between networks and Fermi statistics has also beed studied \cite{Bi,Ci,Hal}.

Herein, we discuss several networks from a unified perspective.
We propose a model that connects the random graph, Barab\'{a}si-Albert (BA) model,
 and lattice model with a parameter $\omega$, which represents the weights of
  the incoming and outgoing links. Our model is an evolving network model,
   and newly added nodes are more likely to be linked to other nodes
   with high ``popularity.’’. The popularity is defined as
   the sum of the number of incoming and the product of the number of
   outgoing links with $\omega$.
 When $\omega$ is positive, the evolving mechanism is the preferential attachment
process in the BA model, and the model exhibits positive feedback.
When $\omega=1$, the weights for the incoming and outgoing
links are equal, and the model takes the form of the BA model \cite{BA}.
When $\omega$ is negative, there is a limit on the number of the
  links from newly added nodes in the evolution of the network.
In this case, the model incorporates negative feedback. In addition, $\omega$ is
related to the length of the memory, and indicates the age of the node to which a new node connects.
Long and short memories indicate that new nodes tend to connect to older and newer nodes, respectively.
When $\omega$ is positive, the network has long memory, and this corresponds to scale-free networks.
Networks with negative $\omega$ have intermediate memories, and the limit $\omega=-1$ corresponds
to extended lattices.

To understand $\omega$ more deeply, we introduce the fitness model, which is an extension of
the BA model\cite{BA1,BA2}.
Fitness is introduced in each node to adjust popularity, and
it can be transformed to the ``energy’’ using the Boltzmann weight.
The fitness decreases as the energy increases, and the degree distribution can thus be
  interpreted in terms of quantum statistics. When $\omega$ is positive, the number of
  outgoing links are not limited, as in the case of Bose statistics. When $\omega$ is negative,
the number of outgoing links are limited, as in the case of Fermi statistics,
  where an energy level can accommodate one particle at most.
We obtain a distribution that connects the Bose and Fermi statistics.
When $\omega$ is negative, Fermi degeneracy can be observed in the network at
low temperatures. Nodes with high fitness have the maximum number of links, and nodes
with low fitness have no links except for the incoming links. When $\omega$ is
positive, Bose-Einstein condensation (BEC) was realized in the network, and
a small number of hubs with low energy (high popularity) is connected to most of the nodes in the network.
In other words, the ``winner takes all’’ scenario occurs.
Furthermore, $\omega$ is related to the threshold of BEC, and the BEC phase becomes
wider as $\omega$ increases. When $\omega=0$, which corresponds to the random network,
we obtain the Boltzmann distribution.

The remainder  of this paper is organized  as follows.
In section 2, we introduce  the networks which have the parameter $\omega$.
In the section 3,   we discuss the relation between the quantum statistics and networks.
Finally, the conclusions are presented in section 4.

\subsection{II. Networks}

We begin this section by discussing the generation of the networks.
The proposed model is an evolving network model, and we
consider the case where a newly added node selects at most
$r$ different nodes based on the popularity. The process is sequential
like in the BA model \cite{BA, BA1.5}.
When the $i$-th node enter the network
at time $t=i$, it selects $r$
nodes for the connections, as shown in Fig\ref{BA2}.
For node $i$, the links from the selected $r$ nodes are
incoming links,
and these are shown by the incoming arrows in Fig.\ref{BA2}.
For the $r$ nodes selected by node $i$, the links to the node $i$
are outgoing links.
In Fig.\ref{BA2}, we draw the links with the outgoing arrows.
We denote the number of incoming and outgoing links of
node $i$ as $k_{i}^{IN}$ and $k_{i}^{OUT}$, respectively.
The degree of node $i$ is given by $k_{i}=k_{i}^{IN}+k_{i}^{OUT}$.
The popularity $l_i$ of node $i$ is defined as
\begin{equation}
l_{i}=\bar{\omega}\cdot k_{i}^{IN}+\omega \cdot k_{i}^{OUT}. \label{eq:pop}
\end{equation}
Here, $\bar{\omega}$ and $\omega$ are the weights for $k_{i}^{IN}$
and $k_{i}^{OUT}$, respectively. When $\omega=\bar{\omega}=1$,
the model attains the form of the BA model. Hereafter,
we set $\bar{\omega}=1$ for normalization without loss of
generality.

\begin{figure}[h]
\includegraphics[width=110mm]{BA.eps}
\caption{Sample graph with $r=3$. Node $i$ selects three nodes and
  is selected by three nodes. The arrows represent the directions
  from the selected nodes to the selecting nodes.
  $k_{i}^{IN}=k_{i}^{OUT}=3$, and $l_i=3\omega+3$.
Node $j$ selects three nodes and is selected by one node.
$k_{j}^{IN}=3,k_{j}^{OUT}=1$, and $l_j=3+\omega $. $l_j$
changes from $3+\omega$ to $3+2\omega$ when a new node
  joins to the network and selects node $j$.
  $k_{j}^{OUT}$ increases from $1$ to $2$.}
\label{BA2}
\end{figure}

The initial state of the model is the complete network with $r+1$ nodes,
and we label them as node $i,i=1,\cdots,r+1$.
We set $k_{i}^{IN}=i-1$ and $k_{i}^{OUT}=r+1-i$.
We add node $t\ge r+2$ and link it with $r$ nodes sequentially.
The probability that node $i<t$ is selected by node $t$ is
\begin{equation}
P(\mbox{node} \,\, i \,\, \mbox{is selected by node}\,\,t)
=p_{i}(t)=\frac{\mbox{Max}(l_i,0)}{\sum_{s=1}^{t-1}\mbox{Max}(l_s,0)} \label{eq:prob}
\end{equation}
We note that when $l_i$ becomes negative or zero,
node $i$ cannot be selected. This occurs for negative values of
$\omega$ and $k_{i}^{OUT}\ge \frac{k_{i}^{IN}}{|\omega|}$.
The maximum value of $k_{i}^{OUT}$ is $\lceil k_{i}^{IN}/|\omega| \rceil$,
where $\lceil x\rceil$ is the ceiling function, and $k_{i}^{OUT}$ cannot
exceed $\lceil \frac{k_{i}^{IN}}{|\omega|} \rceil$. We denote the number
of the nodes with positive $p_{i}(t),i=1,\cdots,t-1$ by $r'$. If $r'$
is less than or equal to $r$, all nodes are linked by node $t$, and we set
$k_{t}^{IN}=r'$. If $r'>r$, $r$ different nodes are selected
from $r'$ nodes with probability $p_{i}(t)$, and we set $k_{t}^{IN}=r$.
$k_{i}^{OUT}$ for the selected nodes increases by unity.

$\omega=1$ and $\omega=0$ correspond to the BA model \cite{BA} and
the random network, respectively. This is because when
$\omega=0$, all nodes are selected with equal probability.
We consider that the range of negative $\omega$ is $-1\leq\omega<0$.
When $\omega>0$, the popularity and the probability of
selection increase during node selection.
This is called positive feedback or preferential attachment.
When $\omega<0$, the popularity and the probability selection
decrease during node selection.

Samples of the network are presented in Fig.\ref{sample}.
In Fig.\ref{sample} (a), we show the case of $r=1$,
which corresponds to tree networks, and in (b), we show the case of $r=3$.
We upload the R script that perform the numerical studies of
this paper \cite{git}.

\begin{figure}[h]
\begin{center}
\begin{tabular}{c}
\begin{minipage}{0.5\hsize}
\begin{center}
\includegraphics[clip, width=8cm]{N100r1.eps}
\hspace{1.6cm} (a)
\end{center}
\end{minipage}
\begin{minipage}{0.5\hsize}
\begin{center}
\includegraphics[clip, width=8cm]{N100r3.eps}
 \hspace{1.6cm} (b)
\end{center}
\end{minipage}
 \end{tabular}
\caption{Sample networks with $\omega \in \{-1,0,1,100\}$ in 100 steps.
  (a) $r=1$. In this case, the network becomes a tree network. (b) $r=3$.
  $\omega=1$ corresponds to the BA model. When $\omega>0$, the network is a scale-free network.
  When $\omega=0$, the network is a random network.
  When $\omega=-1$, the network corresponds to an extended lattice.}
\label{sample}
\end{center}
\end{figure}

\subsection{A. Degree distribution}
We consider the degree distribution of the networks. The number of links in the
node $i$ is $k_i$. When $\omega=1$ for the BA model, $l_i=k_i$. The sum of the
popularities is $\sum_j l_j=(\omega+1)r[t-1/2(r+1)]\sim (\omega+1)rt$ for large $t$.
Note that when $\omega$ is negative, the popularity decreases as the number of links to
  the node increases. However, the total popularity does not decrease when $\omega>-1$,
  and new nodes can join. We can obtain the differential equation for the expected
value of $k_i$ as
\begin{equation}
\frac{d k_i}{dt}=\frac{(k_i-r)\omega+r}{(1+\omega)t}.
\label{df}
\end{equation}

Note that when $\omega$ is negative, if $k_{i}^{OUT}>k_{i}^{IN}/|\omega|$,
$l_i=k_{i}^{IN}-|\omega| k_{i}^{OUT}<0$ and $p_{i}=0$.
We discuss this effect in more detail later.

\subsubsection{i. Case when $\omega>0$ }
Eq.(\ref{df}) can be solved with the initial condition $k_i(t_i)=r$ and we
obtain,
\[
k_i=\frac{r}{\omega}\left(\frac{t}{t_i}\right)^{\omega/(1+\omega)}-r\frac{1-\omega}{\omega}.
\]
We estimate the average probability that node $i$ is selected as
\[
\frac{l_i}{(1+\omega)t}=\frac{r}{1+\omega}\left(\frac{1}{t_i}\right)^{\frac{\omega}{1+\omega}}t^{-1/(1+\omega)}.
\]
The average frequency of selection of node $i$,  called ``memory,’’
decreases as a power law with $t$, and the index is $1/(1+\omega)$.
We summarize the power index of the memory $1/(1+\omega)$ in Fig.\ref{memory}(a).
As $\omega$ increases, the power index becomes smaller, and the memory becomes longer.
In the limit $\omega\to \infty$, the memory covers all history equally,
and the power index of the memory becomes 0.
We term the cases where the power index is larger than 1
and smaller than 1 as the cases of intermediate memory and
long memory, respectively.
This classification depends on whether the integral of the memory over $t$
is finite\cite{Long,Hisakado6,Hisakado7}.

We can obtain the cumulative degree distribution of the network as

\begin{eqnarray}
  P[k_i(t)<k]&=&P\left[(\frac{r}{\omega})^{(1+\omega)/\omega}
    \frac{t}{(k+r\frac{1+\omega}{\omega})^{(1+\omega)/\omega}}<t_i\right]
\nonumber \\
&=&
1-(\frac{r}{\omega})^{(1+\omega)/\omega}
\frac{1}{(k+r\frac{1-\omega}{\omega})^{(1+\omega)/\omega}}\nonumber.
\end{eqnarray}

Note that $k_i\geq r$.
Hence, the degree distribution is
\begin{equation}
p(k)=\frac{\partial P[k_i(t)<k]}{\partial k}
\propto k^{-\frac{1+2\omega}{\omega}}.
\end{equation}
In this region, the degree distribution exhibits power-law decay.
When $\omega=1$, the network corresponds to the BA model with distribution $k^{-3}$.
When $\omega\rightarrow\infty$ and $\omega=1/2$, the networks
have distributions $k^{-2}$ and $k^{-4}$, respectively.
We show the power index in Fig.\ref{memory} (b).
We also show the numerical estimates of $\gamma$ with 100 sample networks
of $10^4$ nodes with $r=3$. We estimate $\gamma$ and its standard error
using the maximum likelihood estimator\cite{newman}. This is similar
to the extended BA model \cite{beta,Krap,Krap2, Krap3,Doro,Doro2,Basu},
which has a discount factor for each node. If the discount factor is introduced,
  the transition from the lattice to the scale-free network can be observed.
  However, these networks do not include trees, whereas our networks include the tree network.
The minimum power index is 2 when the network is graphical. This means that a degree sequence
can be converted into a simple graph\cite{nw}. Fig.\ref{memory} (c) shows the empirical
  degree distribution for the data in Fig.\ref{memory} (b) with $\omega=0.5,1,3,10$.
  When $\omega=0.5$ or $10$, the discrepancy from the power-law behavior becomes apparent.
  This may be attributed to the finite size effect. It becomes difficult to observe the
  power-law behavior for a small $\omega<<1$ and a large $\omega>>1$.

\begin{figure}[htbp]
\begin{center}
  \begin{tabular}{cc}

      \begin{minipage}{0.5\hsize}
        \begin{center}
          \includegraphics[clip, width=6cm]{t_vs_memory.eps}
          \hspace{-3cm}(a)
        \end{center}
      \end{minipage}
      &

      \begin{minipage}{0.5\hsize}
        \begin{center}
          \includegraphics[clip, width=6cm]{omega_vs_gamma.eps}
          \hspace{-3cm}(b)
        \end{center}
      \end{minipage}
      \\

      \begin{minipage}{0.5\hsize}
        \begin{center}
          \includegraphics[clip, width=6cm]{omega_vs_dist_k.eps}
          \hspace{-3cm}(c)
        \end{center}
      \end{minipage}
      &
      \begin{minipage}{0.5\hsize}
        \begin{center}
          \includegraphics[clip, width=6cm]{omega_vs_dist_k_2.eps}
          \hspace{-3cm}(d)
        \end{center}
      \end{minipage}
\end{tabular}
  \caption{(a)Plot of the average probability of the choice of node $i$ (memory) vs. $t$.
    The power index of the decay of the memory is $1/(1+\omega)$.
(b) Plot of $\gamma$ vs. $\omega>0$, which is the power index of the degree distribution in Eq.(4).
    We show the estimates of $\gamma$ and its standard errors using symbols ($\circ$) and error-bars. We sample 100 networks of $10^4$ nodes for each $\omega\in \{0.1,0.5,1,2,3,4,5,10\}$,
      and estimate $\gamma$ using the empirical degree distributions.
      (c) Plot of the empirical degree distributions for $\omega=0.5,1,3,10$ for the data in (b). The straight lines are the power-law equations with the theoretical power index $\gamma=2+1/\omega$.
(d) Plot of the empirical distribution for $\omega=0,-1$.
}

\label{memory}
\end{center}
\end{figure}

\subsubsection{ii. Case when $\omega=0$}
We solve Eq.(\ref{df}) with the initial condition $k_i(t_i)=r$ when $\omega=0$ to obtain

\[
k_i(t)=r \log(\frac{t}{t_i})+r.
\]

We estimate the probability that node $i$ is selected as
\[
\frac{l_i}{t}=\frac{r}{t}.
\]
The selected frequency decreases as a power law, and the index is $1$.
Therefore, the memory is long, but shorter than the case when $\omega>0$.
We can obtain the cumulative degree distribution of the network as

\[
P[k_i(t)<k]=P[t_i>t e^{-\frac{k}{r}}]=
1- e^{-\frac{k}{r}+1}.
\]

The degree distribution is
\begin{equation}
p(k)\propto    e^{-\frac{k}{r}}.
\end{equation}
The degree distribution decays exponentially (Fig.\ref{memory}(d)), and
the model corresponds to a random network\cite{nw}.

\subsubsection{iii. Case when $-1<\omega<0$}

In this case, there is a probability that the popularity $l_i$ becomes negative and $p_i=0$.
Eq.(\ref{df}) becomes an approximated equation when $1/|\omega|$  is not an integer.
When $1/|\omega|$ is an integer, $p_i\geq 0$ for all $i$.

We solve Eq.(\ref{df}) with the initial condition $k_i(t_i)=r$ for $\omega<0$ as
\[
k_i(t)= r+r\frac{1}{|\omega|}\left(1-
\left(\frac{t_i}{t}\right)^{\frac{|\omega|}{1-|\omega|}}\right).
\]
Note that in this case, $r\leq k_i  \leq r+\lceil r/|\omega| \rceil$.
We estimate the probability that node $i$ is selected as
\[
\frac{l_i}{(1+\omega)t}=\frac{r}{1+\omega}t_i^{-\frac{\omega}{1+\omega}}\cdot t^{-\frac{1}{1+\omega}}.
\]
The selected frequency of node $i$ decreases as per the power law with $t$, and the index is
$1/(1+\omega)$.
When $-1<\omega<0$, the power index is larger than 1, and the memory
becomes intermediate. We can obtain the cumulative degree distribution of the network as

\begin{eqnarray}
P[k_i(t)<k]&=&
P\left[(\frac{|\omega|}{r})^{\frac{1-|\omega|}{|\omega|}}
  (\frac{1+|\omega|}{|\omega|}r-k)^{\frac{1-|\omega|}{|\omega|}}t>t_i\right]
\nonumber \\
&=&
\left(1-\frac{|\omega|}{r}(k-r)\right)^{\frac{1-|\omega|}{|\omega|}}.
\end{eqnarray}

In this region, the distribution is the power decay in the range
$r$ to $r+\lceil r/|\omega| \rceil$.
The power index is $2-1/|\omega|$.
When $\omega=-1/2$, the power index is zero, and the distribution
becomes a uniform distribution between $r$ and $3r$.
We show the empirical degree distribution in Fig.\ref{degree} for $r=3$,
and confirm the uniform distribution for $\omega=-1/2$.
When $\omega=-0.9$ and $\omega=-0.999$, $k=6=3+\lceil 3/|\omega|\rceil$.
After the node $i$ is selected by four new nodes,
$l_i$ becomes negative, and $p_{i}=0$.
We can observe this effect to a larger extent for $\omega=-0.9$
than for $\omega=-0.999$.

\begin{figure}[h]
\includegraphics[width=120mm]{neg_omega_vs_deg_dis.eps}
\caption{Empirical degree distribution for $r=3$,
  $\omega=-1/3, -1/2, -0.9$, and $-0.999$. We sample $10^3$ networks of $10^4$ nodes.}
\label{degree}
\end{figure}

\subsubsection{iv. Limit $\omega\to -1$}

Before we discuss the limit $\omega\to -1$, we make a comment about the initial
configuration of the network.
If the initial condition of the network before its evolution is a complete graph,
the limit $\omega\to -1$ is trivial. The number of linkable nodes $r'$ is $r$ or $r-1$
and the network becomes an extended lattice. Here, we consider a more general
initial configuration of the network.
In the limit $\omega\to -1$, Eq.(\ref{df}) is undefined.
The differential equation for the expected value of $k_i$ becomes
\begin{equation}
\frac{d k_i}{dt}=\frac{-(k_i-k_{i}^{IN})+k_{i}^{IN}}{C}=\frac{-(k_i-2k_{i}^{IN})}{C},
\label{df3}
\end{equation}
where $C$ is a constant that depends on the initial condition.
For example, $C$ is the number of nodes that have no links at time $0$.
We note that node $t$ cannot necessarily be linked with $r$ nodes when it
enters the network, as the number $r'$ of nodes with $p_i>0$ might be less than $r$.
The solution of Eq.(\ref{df3}) with $k(t_i)=k_{i}^{IN}$ is
\[
k_i=2k_{i}^{IN}-k_{i}^{IN} e^{-\frac{1}{C}(t-t_i)}.
\]
$k_i$ increases exponentially from $k_{i}^{IN}$ to $2k_{i}^{IN}$.
We estimate the probability that node $i$ is selected as
\[
\frac{l_i}{C}\sim \frac{r}{C}e^{-\frac{1}{C}(t-t_i)}.
\]
Here, we assume $k_{i}^{IN}=r$ for most nodes, which is true for $\omega>-1$.
The selected frequency decreases exponentially, and the memory becomes short.
We can confirm this in Fig.\ref{memory} (a).
In the case $C\rightarrow 0$, the network becomes an extended lattice,
which can be confirmed in Fig.\ref{memory} (d).
The degree distribution
is the superposition of delta functions at $k=2r$ and $2(r-1)$.
Most of the nodes are selected by the previous $r$ or $r-1$ nodes and
the numbers of the incoming and outgoing arrows are equal.

\subsection{B. Network properties}

We study the average distance of the network $\L(\omega)$
and the cluster coefficient $C(\omega)$ in Fig.\ref{cluster}(a).
We sample $10^3$ networks with $10^3$ nodes and $r=3$, and
  estimate the average distance and the cluster coefficients.
The average distance $\L(\omega)$ decreases rapidly as $\omega$ increases
from $\omega=-1$, which corresponds to an extended lattice.
The cluster coefficient $C(\omega)$ decreases more slowly than
$L(\omega)$, and it becomes minimum at $\omega\sim 0$.
The regions with large cluster coefficients and small average distances are
observed for $\omega\sim -1$ and $\omega>0$.
These properties are observed in several real-world networks \cite{nw}.
We compared the average distance and the cluster coefficient with those
in the Watts--Strogatz (WS) model, which uses a parameter $p$ instead of $\omega$ \cite{WS}.
$p$ is known as the probability of rewiring the extended lattice.
  When $p=0$ and $p=1$, the network becomes the extended lattice and the random network, respectively.
There is a difference between our model and the WS model in the regions $p\sim 0$ and
$\omega\sim -1$ near the extended lattice where the cluster coefficient is large, and
the average distance is small.
The region with this property is wider in the WS model than in our model.
  We also study the system size dependence of the region in Fig.\ref{cluster}(b).
  As the number of nodes $N$ increases from $10^3$ t0 $10^4$,
  the region moves leftward. The result suggests that the region disappears in the limit $N\to\infty$.
  About the WS mode, the curve of the cluster coefficient does not move with $N$ and only the
  curve for the average distance moves leftward. The region becomes wider the the increase of $N$.
Near the random network $p\sim 1$ and $\omega\sim 0$, we can confirm a small cluster
coefficient and a small average distance.

\begin{figure}[h]
\begin{center}
\begin{tabular}{c}
\begin{minipage}{0.5\hsize}
\begin{center}
\includegraphics[clip, width=8cm]{omega_vs_CL.eps}
\hspace{1.6cm} (a)
\end{center}
\end{minipage}
\begin{minipage}{0.5\hsize}
\begin{center}
\includegraphics[clip, width=8cm]{omega_vs_CL2.eps}
 \hspace{1.6cm} (b)
\end{center}
\end{minipage}
 \end{tabular}
\caption{
  Average distance $L(\omega)$ and cluster coefficient $C(\omega)$ of the networks.
  (a) $N=10^3,r=3$ and we collect $10^3$ sample networks. WS is the Watts--Strogatz model.
  (b)$N=10^2,10^3,10^4$ and we collect $10^3,10^3$ and $10^2$ networks, respectively.}
\label{cluster}
\end{center}
\end{figure}

\subsection{III. Quantum statistics and networks}
In this section, we introduce the fitness model for the nodes \cite{BA1,BA2}
to establish the relation between the degree distribution and the quantum statistics.
The fitness is the weight for popularity $l_i$, and
the popularity of node $i$ is transformed as $\hat{l}_i=e^{-\beta\epsilon_i} l_i$.
Here, $e^{-\beta\epsilon_i}$ is the fitness, $\beta$ is the inverse temperature,
  and $\epsilon_i$ is the associated energy level of node $i$ in the terminology of
  quantum statistical mechanics.
The equation for the expected number of $k_i$ in Eq.(\ref{df}) is modified as

\begin{equation}
\frac{d l_i}{dt}=
\omega r \frac{e^{-\beta\epsilon_i} l_i}{\sum_i e^{-\beta\epsilon_i} l_i}.
\label{def}
\end{equation}

Here, we assume the solution of this equation as
\begin{equation}
l_i=r\left(\frac{t}{t_i}\right)^{f(\epsilon_i)},
\label{sol}
\end{equation}
which satisfies the initial condition, $l_i=r$ at $t=t_i$.
We substitute Eq.(\ref{sol}) into Eq.(\ref{def}) and
obtain

\begin{equation}
f(\epsilon_i)=\omega \frac{rt e^{-\beta\epsilon_i}}{\sum_i e^{-\beta\epsilon_i} l_i},
\label{f}
\end{equation}

\begin{equation}
 \sum_i e^{-\beta\epsilon_i}l_i\sim
 \int e^{-\beta\epsilon} g(\epsilon)d\epsilon
 \int_1^t l_idt_i
 = \int e^{-\beta\epsilon} g(\epsilon)r\frac{t-t^f}{1-f(\epsilon)}d\epsilon,
 \label{norm}
\end{equation}
where $g(\epsilon)$ is the distribution of the energy level $\epsilon_i$.

In the large-$t$ limit, we calculate Eq.(\ref{norm}) using $f(\epsilon)<1$,
\begin{equation}
  \sum_i e^{-\beta\epsilon_i} l_i\sim
rt\int
  g(\epsilon)
 \frac{e^{-\beta\epsilon} }{1-f(\epsilon)}d\epsilon
 =rte^{-\beta\mu},
 \label{Z2}
\end{equation}
where
\begin{equation}
 e^{-\beta\mu}=\int
  g(\epsilon)
 \frac{ e^{-\beta\epsilon}}{1-f(\epsilon)}d\epsilon.
 \label{17}
\end{equation}
$\mu$ becomes the chemical potential in quantum statistics.
Thus, we can rewrite Eq.(\ref{f}) using Eq.(\ref{Z2}) as
\begin{equation}
    f(\epsilon)=\omega e^{-\beta(\epsilon-\mu)}.
\end{equation}
Then, we use Eq.(\ref{17}) to obtain
\begin{equation}
1=  \int
g(\epsilon)  \frac{1}{e^{\beta(\epsilon-\mu)}-\omega}d\epsilon=I(\omega, \beta, \mu).
\label{fb}
\end{equation}
Eq.(\ref{fb}) is the normalization.
$1/(e^{\beta(\epsilon-\mu)}-\omega)$ corresponds to the distribution function in quantum statistics.
It is an extension of the Bose and Fermi distributions and is similar to fractional exclusion statistics \cite{Hal,Bi,Ci}.
When $\omega=1(-1)$, it corresponds to the Bose (Fermi) distribution that
connects the Bose distribution with the Fermi distribution.

\subsubsection{A. Boson-type Network: $\omega>0$}

We rewrite the distribution function in Eq.(\ref{fb}) as follows:
\begin{equation}
\frac{1}{e^{\beta(\epsilon-\mu)}-\omega}
=
\frac{1}{\omega(e^{\beta(\epsilon-\mu+\delta_{\mu})}-1)}=\frac{1}{\omega(e^{\beta(\epsilon-\mu ')}-1)},
\label{bo}
\end{equation}
where
$\omega=e^{-\beta\delta_{\mu}}$ and $\mu'=\mu-\delta_{\mu}$.
When $\omega=1$, which corresponds to the BA model,
$\delta_{\mu}=0$, and we can obtain the Bose distribution.

We consider the BEC of the network \cite{BA1,BA2}
under the condition of large $\beta$.
We can rewrite Eq.(\ref{fb}) using Eq.(\ref{bo}) as
\begin{equation}
1=\frac{1}{\omega}\int
g(\epsilon)  \frac{1}{e^{\beta(\epsilon-\mu')}-1}d\epsilon=\frac{I'(\omega,\beta,\mu')}{\omega}=   I(\omega, \beta, \mu),
\end{equation}
where
$I'(\omega,\beta,\mu')=\int
g(\epsilon) /(e^{\beta(\epsilon-\mu')}-1)d\epsilon$.
We can obtain the condition of the BEC in Eq.(\ref{fb}) as
\begin{equation}
I'(\omega,\beta,\mu')=\omega.
\label{bec}
\end{equation}

$I'(\omega, \beta, \mu')$ increases as $\mu'$.
If $I'(\omega, \beta, 0)>\omega$, we obtain the solution $\mu'(<0)$ for Eq.(\ref{bec}).
If $I'(\omega, \beta, 0)\leq \omega$, no solution exists.
$\omega-I'(\omega, \beta, \mu')$ of nodes are connected to the maximum fitness node corresponds to the hub in the network.
This corresponds to the ``the winner takes all’’ or BEC of the network.

To confirm this condition, we consider the case where the
distribution of the energy level $\epsilon_i$ is
\begin{equation}
g(\epsilon)=A \epsilon^{\theta},
\end{equation}
where $\theta$ is a parameter for the power distribution, and $A=(\theta+1)/\epsilon_{max}^{\theta+1}$.
$\epsilon_{max}$ is the maximum energy for the fitness model.

The condition of the BEC in Eq.(\ref{bec}) is
\begin{equation}
I'(\omega, \beta,0)=\frac{\theta+1}{(\beta \epsilon_{max})^{\theta+1}}
\int_0^{\beta_{max}}dx\frac{x^{\theta}}{e^x-1}<\omega
\end{equation}
In the large-time limit, we can obtain the lower bound of the transition temperature as
\begin{equation}
T_{BEC}>\epsilon_{max}\left(\frac{\omega}{\zeta(\theta+1)\Gamma(\theta+2)}\right)^{1/(\theta+1)},
\label{tbec}
\end{equation}
where $\Gamma(x)$ is the gamma function and $\zeta(x)$ is the
zeta function.
We can observe the BEC when $\theta>0$.

As $\omega$ increases, the lower bound of the transition temperature increases.
In the limit $\omega\rightarrow \infty$, the phase is the BEC phase.
In the limit $\omega\rightarrow 0$, corresponding to the random network limit, the BEC phase disappears.
In summary, $\omega$ is related to the threshold of the BEC phase by Eq.(\ref{tbec}).
As $\omega$ is increased, lower energy levels are occupied.
In Fig.\ref{bf} (a), we show the depiction of the Bose statistics.

\subsubsection{B. Classical network: $\omega=0$ }
When $\omega=0$, which corresponds to the random network case,
the density function becomes $e^{-\beta(\epsilon-\mu)}$. This corresponds to the
Boltzmann distribution.

\subsubsection{C. Fermion-type network: $-1\leq\omega<0$}
We revise the distribution function in Eq.(\ref{fb}) as
\begin{equation}
\frac{1}{e^{\beta(\epsilon-\mu)}+|\omega|}
=
\frac{1}{|\omega|(e^{\beta(\epsilon-\mu+\delta_{\mu})}+1)}=\frac{1}{|\omega|(e^{\beta(\epsilon-\mu ')}+1)},
\label{fe}
\end{equation}
where
$|\omega|=e^{-\beta\delta_{\mu}}$ and
$\mu'=\mu-\delta_{\mu}$.
When $0<|\omega|\leq 1$, $\delta_{\mu}\geq 0$.

We consider the case when $\beta>>1$, i.e., the low-temperature limit.
Eq.(\ref{fe}) becomes the step function at $\epsilon=\mu'$.
When $\epsilon<\mu' (\epsilon>\mu')$, Eq.(\ref{fe}) becomes $1/|\omega| (0)$.
The maximum occupied number, which corresponds to the number of outgoing arrows,
is $-/|\omega|$ per incoming arrow.
$\omega$ corresponds to the change in the chemical potential and
the maximum occupied number in Eq.(\ref{fe}).
On one hand, the lower $\epsilon$ nodes are connected to other $r+\lceil r/|\omega| \rceil$ nodes.
When a node joins the network, the node is connected to $r$ nodes,
and then the node is connected to $\lceil r/|\omega| \rceil$ nodes with the maximum number of links.
On the other hand, the high $\epsilon$ nodes are connected to other $r$ nodes
that are connected only when the node joins.
This is similar to the Fermi degeneracy of the network.
As $\omega<0$ is increased, the maximum occupied number increases,
and the Fermi energy $\mu'$ decreases.
In the limit $\omega\rightarrow 0$, corresponding to the limit of the random
network, there is no maximum occupied number.
In the limit $\omega=-1$, $\delta_{\mu}=0$, we can obtain the Fermi distribution.
The maximum occupied number obtained after the node joins is one per incoming arrow.
In Fig.\ref{bf} (b), we depict the Fermi statistics.

\begin{figure}[h]
\begin{center}
\begin{tabular}{c}
\begin{minipage}{0.5\hsize}
\begin{center}
\includegraphics[clip, width=8cm]{boson.eps}
\hspace{1.6cm} (a)
\end{center}
\end{minipage}
\begin{minipage}{0.5\hsize}
\begin{center}
\includegraphics[clip, width=8cm]{fermion.eps}
 \hspace{1.6cm} (b)
\end{center}
\end{minipage}
 \end{tabular}
\caption{(a) Bose statistics. (b) Fermi statistics. The horizontal axis corresponds to $\epsilon$.
  The dotted (real) circles correspond to the incoming (outgoing) arrows.
  In the Fermion case, the occupied number is limited.}
\label{bf}
\end{center}
\end{figure}

\subsection{VI. Concluding Remarks }

In this article, we discuss several types of networks:
the random graph, Barab\'{a}si-Albert(BA) model, and lattice networks, using a parameter $\omega$.
$\omega=1,0, and -1$ correspond to the networks above, respectively. The parameter
$\omega$ is related to the preferential attachment of the nodes in the networks.
In our model, we set different weights for the incoming and outgoing arrows.
The parameter $\omega$ is also related to the memory length, which indicates how old a
node is with respect to the new node it connects to.
When $\omega>0$ newly added nodes tend to be connected with older nodes,
and the degree distribution exhibits power-law behavior. $\omega=0$ corresponds to a
random network and long memory with larger power index than the positive $\omega$ case.
Negative $\omega$ corresponds to intermediate memory. The case of $\omega=-1$
corresponds to an extended lattice and short memory.

In addition, we discuss the correspondence between quantum statistics
and the degree distribution of the network model. Positive (negative)
$\omega$ corresponds to Bose (Fermi)-like statistics, and we obtain a
distribution that connects the Bose and Fermi statistics.
When $\omega$ is positive, it corresponds to the threshold of BEC.
As $\omega$ decreases, the area of the BEC phase becomes narrower.
When $\omega$ is negative, $r+\lceil r/|\omega| \rceil$ is the maximum occupied number.
At low temperatures, we can observe the Fermi degeneracy of the network.
The lower energy nodes have the maximum number of links, and the higher energy
nodes have no links except for the initial links.
As $\omega$ decreases, the maximum occupied number decreases.

We can observe directed arrow graphs (DAGs) in several fields.
If we consider an arrow from the selected node to the selecting node, our networks become DAGs.
DAGs have recently been used in some cryptocurrencies for scalability.
The scalability is the property of a system that handles the growing confirmation of transactions.
For example, the Bitcoin network is a 1D lattice \cite{bit} without scalability, and
a Tangle (network) of IOTA, which is the network of authentications, is a DAG \cite{iota}.
Tangle has two incoming arrows and two outgoing arrows.
This is an example of a fermionic network with $r=2$ and $\omega=-1$.
For the case of IOTA, the node corresponds to the transaction, and the network
corresponds to the directed links.
Short memory is a necessary property of cryptocurrency networks, because
the confirmations of recent transactions have an incentive.