# Latency Aware Semi-synchronous Client Selection and Model Aggregation for
Wireless Federated Learning

Liangkun Yu, Xiang Sun, Rana Albelaihi, and Chen Yi L. Yu, X. Sun, and R.
Albelaihi are with the SECNet Lab., University of New Mexico, Albuquerque, NM
87131, USA. E-mail: $\\{$liangkun, sunxiang<EMAIL_ADDRESS>Yi is
with the Chongqing Key Laboratory of Signal and Information Processing, School
of Communication and Information Engineering, Chongqing University of Posts
and Telecommunications, Chongqing 400065, China. E-mail:
<EMAIL_ADDRESS>work was supported by the National Science Foundation
under Award CNS-2148178.

###### Abstract

Federated learning (FL) is a collaborative machine learning framework that
requires different clients (e.g., Internet of Things devices) to participate
in the machine learning model training process by training and uploading their
local models to an FL server in each global iteration. Upon receiving the
local models from all the clients, the FL server generates a global model by
aggregating the received local models. This traditional FL process may suffer
from the straggler problem in heterogeneous client settings, where the FL
server has to wait for slow clients to upload their local models in each
global iteration, thus increasing the overall training time. One of the
solutions is to set up a deadline and only the clients that can upload their
local models before the deadline would be selected in the FL process. This
solution may lead to a slow convergence rate and global model overfitting
issues due to the limited client selection. In this paper, we propose the
Latency awarE Semi-synchronous client Selection and mOdel aggregation for
federated learNing (LESSON) method that allows all the clients to participate
in the whole FL process but with different frequencies. That is, faster
clients would be scheduled to upload their models more frequently than slow
clients, thus resolving the straggler problem and accelerating the convergence
speed, while avoiding model overfitting. Also, LESSON is capable of adjusting
the tradeoff between the model accuracy and convergence rate by varying the
deadline. Extensive simulations have been conducted to compare the performance
of LESSON with the other two baseline methods, i.e., FedAvg and FedCS. The
simulation results demonstrate that LESSON achieves faster convergence speed
than FedAvg and FedCS, and higher model accuracy than FedCS.

###### Index Terms:

Federated learning, client selection, model aggregation, semi-synchronous

##  I Introduction

With the development of Internet of Things (IoT), numerous smart devices, such
as smartphones, smartwatches, and virtual reality headsets, are widely used to
digitize people’s daily lives. Traditionally, a huge volume of data generated
by these IoT devices is uploaded to and analyzed by a centralized data center
that generates high-level knowledge and provides corresponding services to
users, thus facilitating their lives [1]. A typical example is smart home,
where various IoT devices, such as smart meters, thermostats, motion
detectors, and humidity sensors, are deployed to monitor the status of the
smart homes. The data generated by the IoT devices would be uploaded to a
centralized data center, which applies deep reinforcement learning model to
intelligently and autonomously control, for example, the smart bulbs and air
conditioners in the smart homes, to improve the quality of experience and
reduce the energy usage of smart homes [2].

On the other hand, sharing data with third-party data centers may raise
privacy concerns as data generated by the IoT devices may contain personal
information, such as users’ locations and personal preferences [3, 4]. As a
result, various policies have been made, such as General Data Protection
Regulation (GDPR) made by the European Union [5], to regulate and hinder data
sharing. In order to fully utilize these personal data while preserving
privacy (i.e., without sharing the data), federated learning (FL) is proposed
to distributively train machine learning models by enabling different IoT
devices to analyze their data locally without uploading them to a central
facility [6]. A typical example to use FL is to train a next word prediction
model, which is used to predict what word comes next based on the existing
text information [7]. Basically, as shown in Fig. 1, an FL server would first
initialize the parameters of the global model and then broadcast the global
model to all the clients via wireless networks. Each client would train the
received global model based on its local data sets (i.e., their text messages)
and upload the updated model to the FL server via wireless networks. The FL
server then aggregates the received models from the clients to generate a new
global model, and then starts a new iteration by broadcasting the new global
model to the clients. The iteration continues until the model is converged.

Figure 1: Wireless federated learning.

The traditional FL method can resolve the privacy issue to allow the clients
to train the model locally, and it applies the synchronous strategy, where the
FL server has to wait until it receives the models from all the clients in
each global iteration. This may lead to the straggler problem when the
configurations of clients are heterogeneous, meaning that they have different
computing and communications capabilities. Hence, some stragglers take much
longer time to train and upload their models in a global iteration because of
their lower computing and communications capabilities, and thus significantly
prolong the model training process. To resolve the straggler problem, many
client selection methods have been proposed [8, 9, 10, 11, 12], which would
select the qualified clients that can finish their model training and
uploading before a predefined deadline. Normally, client selection and
resource allocation are jointly optimized to maximize the number of the
selected qualified clients. Selecting the qualified clients can resolve the
straggler problem, but may raise other issues. First, the proposed client
selection may significantly reduce the number of participated clients, which
may slow down the convergence speed [13], thus leading to longer training
latency (which equals to the sum of the latency for all the global
iterations). Second, the proposed client selection may result in the model
overfitting issue caused by the reduction of data diversity. That is, if the
FL server only selects the qualified clients to participate in the model
training, then the generated model can only fit the data samples in these
qualified clients, but not the non-qualified clients. The model overfitting
issue would be compounded if the data samples of the qualified clients are not
sufficient [14]. The other solution to solve the straggler problem is to apply
asynchronous FL, where the FL server does not need to wait until the deadline
expires for each global iteration, but would update the global model once it
receives a local model from a client [15]. However, the asynchronous FL may
suffer from 1) the high communications cost since both the FL server and
clients will more frequently exchange their models, and 2) the stale issue,
where some slow clients are training based on an outdated global model, which
may lead to slow convergence rate or even global model divergence [16, 17].

In order to solve the slow convergence and model overfitting issues in the
synchronous FL while avoiding model divergence in the asynchronous FL, we
propose a semi-synchronous FL method, i.e., Latency awarE Semi-synchronous
client Selection and mOdel aggregation for federated learNing (LESSON). The
basic idea of LESSON is to allow all the clients to participate in the whole
learning process with different frequencies. Specifically, the clients are
clustered into different tiers based on their model training and uploading
latency. The clients in a lower tier (i.e., lower model training and uploading
latency) would participate in the learning process more frequently than those
in a higher tier (with higher model training and uploading latency). As a
result, the straggler problem can be resolved (since the FL server does not
need to wait for stragglers in each global iteration) and the model
overfitting problem can be fixed (since all the clients join the learning
process to provide high data diversity). The main contributions of the paper
are summarized as follows.

  1. 1.

In LESSON, we propose a latency-aware client clustering method to group
different clients into different tiers based on their model training and
uploading latency.

  2. 2.

In LESSON, we design a model aggregation based on the proposed client
clustering. The designed model aggregation method determines the weight and
time for the clients in each tier to aggregate their local models. Also,
LESSON can flexibly adjust its performance to change the tradeoff between
model accuracy and model convergence rate.

  3. 3.

The performance of LESSON is evaluated via extensive experiments. The results
demonstrate that LESSON achieves a faster convergence speed than FedAvg and
FedCS, and higher model accuracy than FedCS.

The rest of this paper is organized as follows. The related work is summarized
in Section II. System models are described in Section III. Section IV
elaborates the proposed LESSON algorithm, which comprises client clustering
and model aggregation. In Section V, the performance of LESSON is compared
with the other two baseline algorithms via extensive simulations, and
simulation results are analyzed. Finally, Section VI concludes this paper.

##  II Related Work

How to solve the straggler issue is one of the main challenges in synchronous
FL. The existing solutions mainly focus on jointly optimizing client selection
and resource allocation. Nishio and Yonetani [18] aimed to maximize the number
of selected clients that can finish their model training and uploading before
a predefined deadline in each global iteration. By assuming that the selected
clients have to iteratively upload their models to the FL server, they
designed FedCS that jointly optimizes the uploading schedule and client
selection to achieve the objective. Abdulrahman and Tout [19] designed a
similar client selection method FedMCCS. The goal of FedMCCS is to maximize
the number of selected clients, who can not only finish the model training and
uploading before a predefined deadline but also guarantee the resource
utilization less than the threshold to avoid device dropout. Albelaihi _et
al._ [12] proposed a client selection method that tries to achieve the same
objective as FedCS, but they argued that the latency of a client in waiting
for the wireless channel to be available for model uploading should be
considered; otherwise, the selected clients may not upload their local models
before the deadline. Yu _et al._ [20] proposed to dynamically adjust and
optimize the tradeoff between maximizing the number of selected clients and
minimizing the total energy consumption of the selected clients by picking
suitable clients and allocating appropriate resources (in terms of CPU
frequency and transmission power) in each global iteration. Shi _et al._ [21]
also jointly optimized the client selection and resources allocation for FL.
However, the objective is to minimize the overall learning latency (which
equals the product of the average latency of one global iteration and the
number of global iterations), while achieving a certain model accuracy. They
formulated a system model to estimate the number of global iterations given
the global model accuracy requirement. All the mentioned client selection
method for synchronous FL can alleviate the straggler problem, but may lead to
the slow convergence rate and model overfitting issues as we illustrated
previously.

Instead of selecting qualified clients to avoid stragglers, Li _et al._ [22]
proposed to let fast clients train their local models by running more gradient
descent iterations in each global iteration. In this way, the fast and slow
clients could upload their models at the similar time, but the fast clients
may provide better models to fit their local data samples, thus potentially
speedup the model convergence. This method, however, may lead to local model
overfitting issue when some fast clients run too many gradient descent
iterations over the limited data samples. The overfitted local models would
significantly slow down the global model convergence. Wu _et al._[23] proposed
to perform spilt learning, where a global neural network is divided into two
parts. The parameters in the former and latter layers are trained in the
clients and the FL server, respectively. As a result, the computational
complexity of the clients is reduced, thus potentially reducing the training
time and energy consumption of the clients.

Other works aim to design asynchronous FL, where the FL server does not need
to wait for the selected clients to upload their local models; instead, once
the FL server receives a local model from a client, it would aggregate the
received local model to update the global model, and then send the updated
global model to the client [24, 25, 26, 27]. However, as mentioned before, the
fast clients and the FL server have to more frequently exchange their models,
thus leading to higher communications cost for both the fast clients and FL
server [28, 29]. Meanwhile, in asynchronous FL, the slower clients may train
their local models based on an outdated global model, which results in slow
convergence or even leads to model divergence [30, 31].

To overcome the drawbacks in synchronous and asynchronous FL, we propose the
semi-synchronous FL to allow all the clients to participate in the whole
learning process with different frequencies. Although the term ”semi-
synchronous FL” has been used by the existing works, the definitions are
different from what we defined in this paper. For example, Stripelis and
Ambite [32] defined semi-synchronous FL as the clients to train their local
models over different sizes of local data sets, depending on their computing
capabilities. The semi-synchronous FL proposed in [33] periodically re-selects
a number of clients, and follows the same method as asynchronous FL to
aggregate and update the global model.

##  III System Models

###  III-A Federated learning preliminary

The idea of FL is to enable distributed clients to cooperatively train a
global model such that the global loss function, denoted as
${\bm{\mathcal{F}}\left(\bm{\omega}\right)}$, can be minimized. That is,

$\mathop{\arg\min}\limits_{\bm{\omega}}{\bm{\mathcal{F}}\left(\bm{\omega}\right)}=\mathop{\arg\min}\limits_{\bm{\omega}}\sum\limits_{i\in\bm{\mathcal{I}}}\frac{{\left|{\bm{\mathcal{D}}_{i}}\right|}}{{\left|{\bm{\mathcal{D}}}\right|}}{{f_{i}}\left(\bm{\omega}\right)},$
(1)

where $\bm{\omega}$ is the set of the parameters for the global model,
$\bm{\mathcal{I}}$ is the set of the selected clients,
${\left|{\bm{\mathcal{D}}}\right|}$ is the number of the training data samples
of all the clients, ${\left|{\bm{\mathcal{D}}_{i}}\right|}$ is the number of
the training data samples at client $i$ (where
$\bm{\mathcal{D}}=\bigcup\limits_{i\in\bm{\mathcal{I}}}\bm{\mathcal{D}}_{i}$),
and ${f_{i}}\left(\bm{\omega}\right)$ is the local loss function of client
$i$, i.e.,

${f_{i}}\left(\bm{\omega}\right)=\frac{1}{{\left|{{{\bm{\mathcal{D}}_{i}}}}\right|}}\sum\limits_{n\in{{\bm{\mathcal{D}}_{i}}}}{f\left({\bm{\omega},{\bm{a}_{i,n}},{b_{i,n}}}\right)}.$
(2)

Here, $\left({{\bm{a}_{i,n}},{b_{i,n}}}\right)$ is the input-output pair for
the $n^{th}$ data sample in user $i$’s data set, and
$f\left({\bm{\omega},{\bm{a}_{i,n}},{b_{i,n}}}\right)$ captures the error of
the local model (with parameter $\bm{\omega}$) over
$\left({{\bm{a}_{i,n}},{b_{i,n}}}\right)$.

In each global iteration, FL comprises four steps.

  1. 1.

In the $k$-th global iteration, the FL server broadcasts the global model
generated in the previous global iteration, denoted as
${{\bm{\omega}}^{\left({k-1}\right)}}$, to all the selected clients.

  2. 2.

Each client $i$ trains its local model over its local data set
${\bm{\mathcal{D}}_{i}}$, i.e.,
$\bm{\omega}_{i}^{\left({k}\right)}=\bm{\omega}_{i}^{\left(k-1\right)}-\delta\nabla{f_{i}}\left({\bm{\omega}_{i}^{\left(k-1\right)}}\right)$,
where $\delta$ is the learning rate.

  3. 3.

After deriving the local model ${\bm{\omega}_{i}^{\left(k\right)}}$, client
$i$ uploads its local model to the FL server.

  4. 4.

The FL server aggregates the local models from the clients and updates the
global model based on, for example, FedAvg [6], i.e.,
${{\bm{\omega}}^{\left({k}\right)}}=\sum\limits_{i\in\bm{\mathcal{I}}}{\frac{{\left|{\bm{\mathcal{D}}_{i}}\right|}}{{\left|{\bm{\mathcal{D}}}\right|}}{\bm{\omega}}_{i}^{\left({k}\right)}}$.

The global iteration keeps executed to update the global model
${{\bm{\omega}}^{\left({k}\right)}}$ until the global model converges.

Figure 2: Illustration of client scheduling in LESSON.

###  III-B Latency models of a client

There are four steps in each global iteration for FL, and so the latency of a
global iteration equals the sum of the latency among these four steps. The
local model training latency in Step 2) and local model uploading latency in
Step 3) are different among the clients, depending on their computing and
communications capacities. In addition, the broadcast latency in Step 1) and
model aggression latency in Step 4) are the same for all the clients, and they
are negligible as compared to local model training and uploading latency.
Thus, we define the latency of client $i$ in a global iteration as follows.

$t_{i}=t_{i}^{comp}+t_{i}^{upload},$ (3)

where $t_{i}^{comp}$ is the computing latency of client $i$ in training its
local model over its local data samples, and $t_{i}^{upload}$ is the uploading
latency of client $i$ in uploading its local model to the FL server.

####  III-B1 Computing latency

The computing latency of client $i$ in a global iteration can be estimated by
[34]

$t^{comp}_{i}=\theta\log_{2}\left(\frac{1}{\epsilon}\right)\frac{C_{i}{\left|{\bm{\mathcal{D}}}_{i}\right|}}{f_{i}},$
(4)

where $\theta$ is a constant determined by the structure of the desired model;
$\theta\log_{2}\left(\frac{1}{\epsilon}\right)$ is the number of iterations
required for a desired accuracy $\epsilon$; $C_{i}$ is the number of CPU
cycles required for training one data sample of the model ${\bm{\omega}}$,
varied by the complexity of a sample and CPU design;
${\left|{\bm{\mathcal{D}}}_{i}\right|}$ is number of samples on each client;
and $f_{i}$ is CPU frequency for training the model.

####  III-B2 Uploading latency

The achievable data rate of client $i$ can be estimated by

${r_{i}}=b\log_{2}\left({1+\frac{{{p_{i}}g_{i}}}{{{N_{0}}}}}\right),$ (5)

where $b$ is the amount of bandwidth allocated to each participate client,
${p_{i}}$ is the transmission power of client $i$, ${g_{i}}$ is the channel
gain from client $i$ to the BS, and ${N_{0}}$ is the average background noise
and inter-cell interference power density. Assume that the size of the local
model is $s$, and so the latency of client $i$ in uploading its local model to
the BS is

$t^{upload}_{i}=\frac{s}{r_{i}}=\frac{s}{b\log_{2}\left(1+\frac{pg_{i}}{N_{0}}\right)}.$
(6)

##  IV Latency awarE Semi-synchronous client Selection and mOdel aggregation
for federated learNing (LESSON)

Different from synchronous FL, the proposed LESSON method aims to allow all
the clients to participate in the whole learning process while avoiding the
straggler problem. The basic idea of LESSON is to cluster the clients into
different tiers based on their latency and the deadline. The clients in
different tiers would train and upload their local models at different
frequencies.

###  IV-A Latency-aware client clustering

Denote $\tau$ as the deadline of a global iteration. The FL server would
accept all the models uploaded from the clients before the deadline $\tau$ and
reject the rest in each global iteration. Hence, we cluster the clients into a
number of tiers, and $x_{ij}$ is used to indicate client $i$ is in tier $j$
(i.e., $x_{ij}=1$) or not (i.e., $x_{ij}=0$). Basically, if client $i$ can
finish its local model training and uploading before deadline $\tau$, i.e.,
$t_{i}\leq\tau$, then client $i$ is in tier-1, i.e., $x_{i1}=1$. Similarly, if
client $i$ can finish its local model training and uploading between $\tau$
and $2\times\tau$, i.e., $\tau<t_{i}\leq 2\times\tau$, then client $i$ is in
tier-2, i.e., $x_{i2}=1$. The following equation provides a general
mathematical expression to cluster client $i$ into a specific tier.

$\displaystyle x_{ij}=\left\\{{\begin{array}[]{*{20}{l}}{1,}&{{\rm{if}}\
\tau\times(j-1)<t_{i}\leq\tau\times j,}\\\
{0,}&{{\rm{otherwise,}}}\end{array}}\right.$ (9)

where $j$ is the index of tiers.

###  IV-B Semi-synchronized model aggregation

In each global iteration, clients from different tiers would upload their
local models, and the FL server is able to estimate when a client may upload
its local model according to its associated tier. Fig. 2 provides one example
to illustrate the scheduling of the clients from four tiers in LESSON. For
example, the clients in tier-1 are expected to upload their local models in
each global iteration, and the clients in tier-2 are expected to upload their
local models in every two global iterations. Denote $k$ as the index of the
global iterations, and let $y_{jk}$ be the binary variable to indicate whether
the clients in tier-$j$ are expected to upload their local models by the end
of $k^{th}$ global iteration, where

$\displaystyle y_{jk}=\left\\{{\begin{array}[]{*{20}{l}}{1,}&{{\rm{if}}\
k\%j=0,}\\\ {0,}&{{\rm{otherwise,}}}\end{array}}\right.$ (12)

where $\%$ is the modulo operation, and so $k\%j=0$ indicates $k$ is divisible
by $j$. Meanwhile, let $z_{ik}$ be the binary variable to indicate whether
client $i$ is expected to upload its local model by the end of $k^{th}$ global
iteration ($z_{ik}=1$) or not ($z_{ik}=0$) , where

$z_{ik}=x_{ij}y_{jk}.$ (13)

Based on the value of $z_{ik}$, the FL server would expect which clients will
upload their local models in global iteration $k$, and then aggregate all the
received local models based on

${{\bm{\omega}}^{\left({k}\right)}}=\sum\limits_{i\in\bm{\mathcal{I}}}{\frac{{\left|{\bm{\mathcal{D}}_{i}}\right|}}{\sum\limits_{i\in\bm{\mathcal{I}}}{\left|{\bm{\mathcal{D}}_{i}}\right|}z_{ik}}{\bm{\omega}}_{i}^{\left({k}\right)}z_{ik}},$
(14)

where
$\sum\limits_{i\in\bm{\mathcal{I}}}{\left|{\bm{\mathcal{D}}_{i}}\right|}z_{ik}$
indicates the total number of the data samples among all the clients. who
would upload their local models in global iteration $k$.

Note that, in synchronous FL (e.g., FedAvg), each selected client would update
its local model $\bm{\omega}_{i}^{\left({k}\right)}$ based on Eq. (15) in each
global iteration.

$\bm{\omega}_{i}^{\left({k}\right)}=\bm{\omega}_{i}^{\left(k-1\right)}-\delta\nabla{f_{i}}\left({\bm{\omega}_{i}^{\left(k-1\right)}}\right),$
(15)

where $\delta$ is the step size, which is the same for all the selected
clients. In LESSON, different clients update their local models in different
frequencies, depending on their associated tiers, and so it is reasonable to
adopt different step sizes for the clients in different tiers [24]. Here, we
adjust the step size of a client proportional to its tier (i.e., $\delta\times
j$). Thus, when $k\%j=0$, client $i$ in tier-$j$ would update its local model
based on

$\bm{\omega}_{i}^{\left({k}\right)}=\bm{\omega}_{i}^{\left(k-j\right)}-j\delta\nabla{f_{i}}\left({\bm{\omega}_{i}^{\left(k-j\right)}}\right)$
(16)

###  IV-C Summary of LESSON

Algorithm 1 summarizes the LESSON algorithm. Initially, the FL server
estimates the latency of all the clients and cluster the clients into
different tiers based on Eq. (9), i.e., Steps 1-2 in Algorithm 1. Then, the FL
server broadcasts the initial global model ${\bm{\omega}}^{\left(0\right)}$ to
start the collaborative model training process, which comprises many global
iterations. In any global iteration $k$, each client trains and updates its
local model ${\bm{\omega}}_{i}^{\left(k\right)}$ based on Eq. (16). If
$z_{ik}=1$, client $i$ should upload its local model to the FL server by the
end of global iteration $k$. Then, client $i$ would wait until it receives the
updated global model ${\bm{\omega}}^{\left(k\right)}$ from the FL server to
start the next round of local model training. On the other hand, the FL server
keeps receiving the local models from the clients in global iteration $k$.
Once the deadline expires, the FL server updates the global model
${\bm{\omega}}^{\left(k\right)}$ based on Eq. (14), and then broadcasts the
new global model ${\bm{\omega}}^{\left(k\right)}$ to the clients, who just
uploaded their local models in global iteration $k$.

Note that the deadline of a global iteration, i.e., $\tau$, is a very
important parameter to adjust the performance of LESSON. Specifically, if
$\tau\rightarrow+\infty$, there is only one tier and all the clients in the
network would be in this tier, and so LESSON acts as FedAvg, where the FL
server has to wait until it receives the local models from all the clients in
each global iteration. On the other hand, if $\tau\rightarrow 0$, LESSON acts
as asynchronous FL, where clients with different latency (i.e., $t_{i}$) will
be in different tiers. The FL server will immediately aggregate the local
models from the clients with the same latency and then update and broadcast
the global model. Thus, one of the advantages of LESSON is that it can change
$\tau$ to adjust its performance, and we will discuss how $\tau$ affects the
performance of LESSON in Section V.

1 Estimate the latency of all the clients based on Eq. (3).

2Cluster the clients into different tiers based on Eq. (9).

3The FL server initializes the global model ${\bm{\omega}}^{\left(0\right)}$
and broadcasts to all the clients.

4$k=1$.

5for  _each global iteration $k$_ do

6

7 Client side:{

8

9 Update its local model based on Eq. (16);

10 if  _$z_{ik}=1$_ then

11 Upload its local model ${\bm{\omega}}_{i}^{\left(k\right)}$ to the FL
server;

12 Wait until it receives the global model ${\bm{\omega}}^{\left(k\right)}$
broadcast by the FL server;

13 end if

14

15 $k:=k+1$;

16

17

18 }

19 FL server side:{

20 Receive all the local models uploaded from the clients during global
iteration $k$;

21 Update the global model ${\bm{\omega}}^{\left(k\right)}$ based on Eq. (14);

22 Broadcast the updated global model ${\bm{\omega}}^{\left(k\right)}$ to the
corresponding clients;

23 $k:=k+1$;

24

25

26 }

27 end for

Algorithm 1 LESSON algorithm

##  V Simulation

In this section, we conduct extensive simulations to evaluate the performance
of LESSON.

Figure 3: Clients’ latency distribution. Figure 4: Probability distribution of
10 categories samples for 5 clients with different $\beta$.

###  V-A Simulation setup

####  V-A1 Configuration of clients

We assume that there are $50$ clients that are uniformly distributed in a 2 km
$\\!\\!\times\\!\\!$ 2 km area, which is covered by a BS located at the center
of the area. All the clients upload their local models to the FL server via
the BS. The pathloss between the BS and client $i$ is calculated based on
$128.1+37.6d_{i}$, where $d_{i}$ is the distance in kilometer between the BS
and client $i$. Meanwhile, the transmission power $p_{i}$ is set to be $1$
Watt for all the clients, and the amount of available bandwidth for each
client in uploading its local model is $30$ kHz. In addition, each client has
around ${\left|{\bm{\mathcal{D}}}_{i}\right|}=1,000$ data samples (i.e.,
$50,000$ combined for all clients) to train its local model. The number of CPU
cycles required for training one data sample (i.e., $C_{i}$) among clients are
randomly selected from a uniform distribution, i.e.,
$C_{i}\sim\mathcal{U}\left(3,5\right)\\!\times\\!10^{5}$ CPU cycles/sample.
The CPU frequency of a client $f_{i}$ is also randomly selected from a uniform
distribution, i.e., $f_{i}\sim\mathcal{U}\left(0.8,3\right)$ GHz. Other
simulation parameters are listed in Table I.

TABLE I: Simulation parameters Parameter | Value
---|---
Noise and inter-cell interference ($N_{0}$) |  $-94$ dBm
Size of the local model ($s$) |  $100$ kbit
Number of local iterations $(\theta\log_{2}(\frac{1}{\epsilon}))$ | $1\\!\times\\!\log_{2}(\frac{1}{0.05})$
Number of local samples ${\left|{\bm{\mathcal{D}}}_{i}\right|}$ | $1000$
Number of local epochs | $1$
Number of local batch size | $20$
Non-IID Dirichlet distribution parameter $\beta$ | $1$

Fig. 3 shows the probability density function of the latency (i.e., $t_{i}$)
among all the clients in a specific time instance. The median value in the
clients’ latency is around 10 seconds, and so we initially set up the deadline
$\tau$ of a global iteration to be 10 seconds, and we will change the value of
$\tau$, later on, to see how it affects the performance of LESSON.

####  V-A2 Machine learning model and training datasets

We will use two benchmark datasets to train the corresponding machine learning
model.

  1. 1.

CIFAR-10[35] is an image classification dataset containing 10 labels/classes
of images, each of which has 6,000 images. Among the 60,000 images, 50,000 are
used for model training and 10,000 for model testing.

  2. 2.

MNIST [36] is a handwritten digits dataset that includes many $28\times 28$
pixel grayscale images of handwritten single digits between 0 and 9. The whole
dataset has a training set of 60,000 examples and a test set of 10,000
examples.

We apply the convolutional neural network (CNN) to classify the CIFAR-10
images. The CNN model has four 3x3 convolution layers (where the first layer
has 32 channels, and each of the following three layers has 64 channels. Also,
only the first two layers are followed with $2\times 2$ max pooling), followed
with a dropout layer with rate of 75$\%$, a fully connected 256 units ReLu
layer, and a 10 units softmax output layer. There are totally 1,144,650
parameters in this CNN model.

With respect to MNIST, which is a much simpler image dataset than CIFAR-10, a
smaller CNN model has been used. Specifically, the CNN model has two 5x5
convolution layers (where the two layers have 6 and 16 channels, respectively,
each of which is followed with a $2\times 2$ max pooling), followed with two
fully connected ReLu layers with 120 and 84 units, respectively, and a 10
units softmax output layer. There are totally 61,706 parameters in this CNN
model.

In addition, we partition the MNIST/CIFAR-10 dataset among the 50 clients
based on non-independent and identical distribution (Non-IID), and the
probability of having $\eta_{m}$ images in label class $m$ at client $i$ is
assumed to follow a Dirichlet distribution [37], i.e.,

$f(\eta_{1},\eta_{2},\dots,\eta_{M};\beta)=\frac{\Gamma(\beta
M)}{\Gamma(\beta)^{M}}\prod_{m=1}^{M}\eta_{m}^{\beta-1},$ (17)

where $M$ is the total number of label classes (i.e., $M=10$ for both CIFAR-10
and MNIST), $\Gamma()$ is the gamma
function111$\Gamma(z)=\int_{0}^{\infty}x^{z-1}e^{-x}dx$, and $\beta$ is the
concentration parameter that determines the level of label imbalance. A larger
$\beta$ results in more balanced data partition among different labels within
a client (i.e., much closer to IID), and vice versa. Fig. 4 shows how
different labels of images are distributed by varying $\beta$. In addition, we
assume each client would train its CNN model over
${\left|{\bm{\mathcal{D}}}_{i}\right|}=20$ data samples locally based on
stochastic gradient descent (SGD) with base learning rate $\delta=0.1$ [38]
and epoch equal to 1.

####  V-A3 Baseline comparison methods

####  V-A4 Baseline algorithms

The performance of LESSON will be compared with the other two baseline client
selection algorithms, i.e., FedCS [18] and FedAvg [13]. FedCS only selects the
clients that can finish their model uploading before the deadline $\tau$ in
each global iteration, i.e., only the clients in tier-1 will be selected to
participate in the training process. In FedAvg, all the clients in the network
will be selected to participate in the training process for each global
iteration. That is, the FL server would wait until it receives the local
models from all the clients, and then update the global model for the next
global iteration.

###  V-B Simulation results

Assume that $\beta=1$ and $\tau=20$ seconds. Fig. 5 shows the test accuracy of
the three algorithms over the global iterations and simulation time for
CIFAR-10 and MNIST. From Figs 5(a) and 5(b), we can find that LESSON and
FedAvg have similar test accuracy, i.e., $\sim 70\%$ for CIFAR-10 and $\sim
95\%$ for MNIST. However, the test accuracy achieved by FedCS is lower than
LESSON and FedAvg, i.e., $\sim 60\%$ for CIFAR-10 and $\sim 90\%$ for MNIST.
This is because FedCS only selects fast clients to participate in the training
process, and so the derived global model can only well fit the data samples
for fast clients, not slow clients, thus reducing the model accuracy.
Meanwhile, the convergence rate with respect to the number of global
iterations for LESSON and FedAvg is also very similar, which is slightly
faster than FedCS. However, by evaluating the convergence rate with respect to
the time, as shown in Figs 5(c) and 5(d), we find out that LESSON is faster
than FedAvg. For example, the global model in LESSON has already converged at
20,000 seconds for MNIST, but the global model in FedAvg is still under-
trained. This is because the FL server in FedAvg has to wait until the local
models from all the clients have been received in each global iteration, and
thus the latency of a global iteration incurred by FedAvg is much higher than
that incurred by LESSON. Table II shows the average delay of a global
iteration incurred by different algorithms, where the average latency of a
global iteration incurred by LESSON is 48 seconds faster than FedAvg. As a
result, FedAvg only runs around $588$ global iterations at 20,000 seconds in
Figs and 5(d), respectively, but LESSON runs $2,000$ global iterations.

TABLE II: Average latency per global iteration for different algorithms. Algorithms | Average latency of a global iteration
---|---
FedAvg | 68 seconds
FedCS | 20 seconds
LESSON | 20 seconds
Figure 5: Test accuracy of different algorithms for CIFAR-10 and MNIST, where
(a) test accuracy vs number of global iterations for CIFAR-10, (b) test
accuracy vs number of global iterations for MNIST, (c) test accuracy vs time
in CIFAR-10, and (d) test accuracy vs time in MNIST. Figure 6: Test accuracy
over number of global iterations for CIFAR-10, where (a) $\beta=0.1$, (b)
$\beta=1.0$, and (c) $\beta=10.0$. Figure 7: Test accuracy over the time for
CIFAR-10, where (a) $\beta=0.1$, (b) $\beta=1.0$, and (c) $\beta=10.0$.

We further investigate how the data sample distribution affects the
performance of the algorithms based on CIFAR-10. As mentioned before, $\beta$
is used to change the data sample distribution, i.e., a larger $\beta$ implies
a more balanced data partition among the labels in a client or data sample
distribution much closer to IID, and vice versa. Assume $\tau=20$ seconds, and
Fig. 6 shows the test accuracy of different algorithms over the number of
global iterations by selecting different values of $\beta$. From the figures,
we can see that if the data sample distribution exhibits non-IID, i.e.,
$\beta=0.1$, FedAvg has higher test accuracy than LESSON and FedCS. As $\beta$
increases, i.e., the data sample distribution is getting closer to IID, the
test accuracy gap between FedAvg and LESSON is getting smaller, while the test
accuracy of FedCS remains unchanged. On the other hand, Fig. 7 shows the test
accuracy of different algorithms over the time by selecting different values
of $\beta$. From the figure, we can see that LESSON achieves much faster
convergence rate than FedAvg and FedCS under different values of $\beta$.
Therefore, we conclude that LESSON achieves a faster convergence rate at the
cost of reducing the model accuracy, especially when the data distribution
exhibits non-IID.

Figure 8: Test accuracy over the number of global iterations for CIFAR-10,
where (a) $\beta=0.1$ and (b) $\beta=1.0$. Figure 9: Test accuracy over the
time for CIFAR-10, where (a) $\beta=0.1$ and (b) $\beta=1.0$.

As mentioned in Section III.C, the deadline of a global iteration, i.e.,
$\tau$, is a very important parameter to adjust the performance of LESSON.
Basically, if $\tau\rightarrow+\infty$, LESSON acts as FedAvg, and if
$\tau\rightarrow 0$, LESSON acts as asynchronous FL. Fig. 8 shows the test
accuracy of LESSON by having different values of $\tau$ and $\beta$ for
CIFAR-10 over global iterations, where we can find that $\tau=60$ seconds
incurs the highest test accuracy for both $\beta=0.1$ and $\beta=1$. This is
because a larger $\tau$ 1) reduces the number of tiers in the system, thus
alleviating the performance degradation caused by stale issue, and 2)
increases the average number of clients in uploading their local models in a
global iteration, which can mitigate the impact caused by non-IID. Also, we
can see that the blue curve, i.e., $\tau=10$ seconds, is more sensitive to the
change of $\beta$ than the other two curves. This is because as $\tau$
reduces, LESSON acts more like asynchronous FL, which has the convergence
issue under non-IID. Figs. 9 show the test accuracy of LESSON by having
different values of $\tau$ and $\beta$ for CIFAR-10 over the time. From Figs.
9(a), we can find that $\tau=60$ exhibits the slowest convergence rate with
respect to the time because a larger $\tau$ implies a longer latency of a
global latency, i.e., $\tau=60$ runs the fewest global iterations than
$\tau=10$ and $\tau=20$ within a time period. Therefore, we conclude that
changing $\tau$ can adjust the tradeoff between the model accuracy and model
convergence rate with respect to the time. A large $\tau$ can increase the
model accuracy but reduce the model convergence rate, and vice versa.
Meanwhile, as the data distribution is closer to IID (i.e., as $\beta$
increases), the difference of the convergence rate among the three algorithms
increases as shown in Fig. 9(b).

##  VI Conclusion

In order to solve the data diversity and straggler issues in the synchronous
FL while avoiding model divergence in the asynchronous FL, we propose a semi-
synchronous client selection and model aggregation algorithm, i.e., LESSON,
which allows all the clients to participate in the FL process but with
different frequencies. Simulation results demonstrate that, in the IID
setting, the model accuracy incurred by LESSON and FedAvg is similar but high
than FedCS, and the convergence rate incurred by LESSON is the highest. In the
non-IID setting, FedAvg achieves the highest model accuracy at the cost of
lowering the convergence rate, and LESSON still incurs the highest convergence
rate. Also, the other advantage of LESSON is that we can adjust the deadline
$\tau$ in LESSON to change the tradeoff between the model accuracy and model
convergence rate with respect to the time. In the future, we will investigate
on how to dynamically adjust $\tau$ such that overall performance of LESSON
can be maximized.

## References

  * [1] X. Sun and N. Ansari, “EdgeIoT: Mobile edge computing for the Internet of Things,” _IEEE Commun. Mag._ , vol. 54, no. 12, pp. 22–29, 2016.
  * [2] Y. Liu, C. Yang, L. Jiang, S. Xie, and Y. Zhang, “Intelligent edge computing for iot-based energy management in smart cities,” _IEEE Network_ , vol. 33, no. 2, pp. 111–117, 2019.
  * [3] M. M. Dhanvijay and S. C. Patil, “Internet of things: A survey of enabling technologies in healthcare and its applications,” _Computer Networks_ , vol. 153, pp. 113–131, 2019. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S1389128619302695
  * [4] J. Granjal, E. Monteiro, and J. Sá Silva, “Security for the Internet of Things: A survey of existing protocols and open research issues,” _IEEE Commun. Surveys Tuts._ , vol. 17, no. 3, pp. 1294–1312, 2015.
  * [5] M. Goddard, “The eu general data protection regulation (gdpr): European regulation that has a global impact,” _International Journal of Market Research_ , vol. 59, no. 6, pp. 703–705, 2017.
  * [6] B. McMahan, E. Moore, D. Ramage, S. Hampson, and B. A. y Arcas, “Communication-Efficient Learning of Deep Networks from Decentralized Data,” in _Proc. 20th Intl. Conf. Artif. Intell. Stat._ PMLR, 20–22 Apr 2017, pp. 1273–1282.
  * [7] A. Hard, K. Rao, R. Mathews, S. Ramaswamy, F. Beaufays, S. Augenstein, H. Eichner, C. Kiddon, and D. Ramage, “Federated learning for mobile keyboard prediction,” _arXiv preprint arXiv:1811.03604_ , 2018.
  * [8] A. Imteaj and M. H. Amini, “Fedar: Activity and resource-aware federated learning model for distributed mobile robots,” in _2020 19th IEEE International Conference on Machine Learning and Applications (ICMLA)_. IEEE, 2020, pp. 1153–1160.
  * [9] W. Wu, L. He, W. Lin, R. Mao, C. Maple, and S. Jarvis, “Safa: A semi-asynchronous protocol for fast federated learning with low overhead,” _IEEE Transactions on Computers_ , vol. 70, no. 5, pp. 655–668, 2020.
  * [10] A. Reisizadeh, I. Tziotis, H. Hassani, A. Mokhtari, and R. Pedarsani, “Straggler-resilient federated learning: Leveraging the interplay between statistical accuracy and system heterogeneity,” _arXiv preprint arXiv:2012.14453_ , 2020.
  * [11] Z. Xu, Z. Yang, J. Xiong, J. Yang, and X. Chen, “Elfish: Resource-aware federated learning on heterogeneous edge devices,” _Ratio_ , vol. 2, no. r1, p. r2, 2019.
  * [12] R. Albelaihi, X. Sun, W. D. Craft, L. Yu, and C. Wang, “Adaptive participant selection in heterogeneous federated learning,” in _2021 IEEE Global Communications Conference (GLOBECOM)_ , 2021, pp. 1–6.
  * [13] J. Wang and G. Joshi, “Cooperative sgd: A unified framework for the design and analysis of communication-efficient sgd algorithms,” _arXiv preprint arXiv:1808.07576_ , 2018.
  * [14] C. Shorten and T. M. Khoshgoftaar, “A survey on image data augmentation for deep learning,” _Journal of Big Data_ , vol. 6, no. 1, pp. 1–48, 2019.
  * [15] L. Li, Y. Fan, and K.-Y. Lin, “A survey on federated learning,” in _2020 IEEE 16th International Conference on Control & Automation (ICCA)_. IEEE, 2020, pp. 791–796.
  * [16] C. Xu, Y. Qu, Y. Xiang, and L. Gao, “Asynchronous federated learning on heterogeneous devices: A survey,” _arXiv preprint arXiv:2109.04269_ , 2021\.
  * [17] G. Damaskinos, R. Guerraoui, A.-M. Kermarrec, V. Nitu, R. Patra, and F. Taiani, “Fleet: Online federated learning via staleness awareness and performance prediction,” in _Proceedings of the 21st International Middleware Conference_ , ser. Middleware ’20. New York, NY, USA: Association for Computing Machinery, 2020, p. 163–177. [Online]. Available: https://doi.org/10.1145/3423211.3425685
  * [18] T. Nishio and R. Yonetani, “Client selection for federated learning with heterogeneous resources in mobile edge,” in _2019 IEEE International Conference on Communications (ICC)_ , 2019, pp. 1–7.
  * [19] S. Abdulrahman, H. Tout, A. Mourad, and C. Talhi, “Fedmccs: Multicriteria client selection model for optimal iot federated learning,” _IEEE Internet of Things Journal_ , vol. 8, no. 6, pp. 4723–4735, 2021.
  * [20] L. Yu, R. Albelaihi, X. Sun, N. Ansari, and M. Devetsikiotis, “Jointly optimizing client selection and resource management in wireless federated learning for internet of things,” _IEEE Internet of Things Journal_ , vol. 9, no. 6, pp. 4385–4395, 2022.
  * [21] W. Shi, S. Zhou, and Z. Niu, “Device scheduling with fast convergence for wireless federated learning,” in _ICC 2020 - 2020 IEEE International Conference on Communications (ICC)_ , 2020, pp. 1–6.
  * [22] T. Li, A. K. Sahu, M. Zaheer, M. Sanjabi, A. Talwalkar, and V. Smith, “Federated optimization in heterogeneous networks,” in _Proceedings of Machine Learning and Systems_ , I. Dhillon, D. Papailiopoulos, and V. Sze, Eds., vol. 2, 2020, pp. 429–450. [Online]. Available: https://proceedings.mlsys.org/paper/2020/file/38af86134b65d0f10fe33d30dd76442e-Paper.pdf
  * [23] D. Wu, R. Ullah, P. Harvey, P. Kilpatrick, I. Spence, and B. Varghese, “Fedadapt: Adaptive offloading for iot devices in federated learning,” _arXiv preprint arXiv:2107.04271_ , 2021.
  * [24] Y. Chen, Y. Ning, M. Slawski, and H. Rangwala, “Asynchronous online federated learning for edge devices with non-iid data,” in _2020 IEEE International Conference on Big Data (Big Data)_. IEEE, 2020, pp. 15–24.
  * [25] Y. Lu, X. Huang, Y. Dai, S. Maharjan, and Y. Zhang, “Differentially private asynchronous federated learning for mobile edge computing in urban informatics,” _IEEE Transactions on Industrial Informatics_ , vol. 16, no. 3, pp. 2134–2143, 2019.
  * [26] B. Gu, A. Xu, Z. Huo, C. Deng, and H. Huang, “Privacy-preserving asynchronous federated learning algorithms for multi-party vertically collaborative learning,” _arXiv preprint arXiv:2008.06233_ , 2020.
  * [27] X. Lian, W. Zhang, C. Zhang, and J. Liu, “Asynchronous decentralized parallel stochastic gradient descent,” in _International Conference on Machine Learning_. PMLR, 2018, pp. 3043–3052.
  * [28] Z. Chai, Y. Chen, L. Zhao, Y. Cheng, and H. Rangwala, “Fedat: A communication-efficient federated learning method with asynchronous tiers under non-iid data,” _ArXivorg_ , 2020.
  * [29] H. R. Feyzmahdavian, A. Aytekin, and M. Johansson, “A delayed proximal gradient method with linear convergence rate,” in _2014 IEEE International Workshop on Machine Learning for Signal Processing (MLSP)_. IEEE, 2014, pp. 1–6.
  * [30] J. Jiang, B. Cui, C. Zhang, and L. Yu, “Heterogeneity-aware distributed parameter servers,” in _Proceedings of the 2017 ACM International Conference on Management of Data_ , ser. SIGMOD ’17. New York, NY, USA: Association for Computing Machinery, 2017, p. 463–478. [Online]. Available: https://doi.org/10.1145/3035918.3035933
  * [31] W. Zhang, S. Gupta, X. Lian, and J. Liu, “Staleness-aware async-sgd for distributed deep learning,” in _Proceedings of the Twenty-Fifth International Joint Conference on Artificial Intelligence_ , ser. IJCAI’16. AAAI Press, 2016, p. 2350–2356.
  * [32] D. Stripelis and J. L. Ambite, “Semi-synchronous federated learning,” _arXiv preprint arXiv:2102.02849_ , 2021.
  * [33] J. Hao, Y. Zhao, and J. Zhang, “Time efficient federated learning with semi-asynchronous communication,” in _2020 IEEE 26th International Conference on Parallel and Distributed Systems (ICPADS)_. IEEE, 2020, pp. 156–163.
  * [34] Z. Yang, M. Chen, W. Saad, C. S. Hong, M. Shikh-Bahaei, H. V. Poor, and S. Cui, “Delay Minimization for Federated Learning Over Wireless Communication Networks,” _arXiv e-prints_ , p. arXiv:2007.03462, Jul. 2020\.
  * [35] A. Krizhevsky, “Learning multiple layers of features from tiny images,” 2009.
  * [36] Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner, “Gradient-based learning applied to document recognition,” _Proceedings of the IEEE_ , vol. 86, no. 11, pp. 2278–2324, 1998.
  * [37] Q. Li, Y. Diao, Q. Chen, and B. He, “Federated learning on non-iid data silos: An experimental study,” _arXiv preprint arXiv:2102.02079_ , 2021.
  * [38] P. Kairouz, H. B. McMahan, B. Avent, A. Bellet, M. Bennis, A. N. Bhagoji, K. Bonawitz, Z. Charles, G. Cormode, R. Cummings, R. G. L. D’Oliveira, H. Eichner, S. E. Rouayheb, D. Evans, J. Gardner, Z. Garrett, A. Gascón, B. Ghazi, P. B. Gibbons, M. Gruteser, Z. Harchaoui, C. He, L. He, Z. Huo, B. Hutchinson, J. Hsu, M. Jaggi, T. Javidi, G. Joshi, M. Khodak, J. Konecný, A. Korolova, F. Koushanfar, S. Koyejo, T. Lepoint, Y. Liu, P. Mittal, M. Mohri, R. Nock, A. Özgür, R. Pagh, H. Qi, D. Ramage, R. Raskar, M. Raykova, D. Song, W. Song, S. U. Stich, Z. Sun, A. T. Suresh, F. Tramèr, P. Vepakomma, J. Wang, L. Xiong, Z. Xu, Q. Yang, F. X. Yu, H. Yu, and S. Zhao, “Advances and open problems in federated learning,” _Foundations and Trends® in Machine Learning_ , vol. 14, no. 1–2, pp. 1–210, 2021. [Online]. Available: http://dx.doi.org/10.1561/2200000083


# The Splashback Mass Function in the Presence of Massive Neutrinos

Suho Ryu and Jounghun Lee<EMAIL_ADDRESS><EMAIL_ADDRESS>Astronomy program, Department of Physics and Astronomy, Seoul National
University, Seoul 08826, Republic of Korea

###### Abstract

We present a complementary methodology to constrain the total neutrino mass,
$\sum m_{\nu}$, based on the diffusion coefficient of the splashback mass
function of dark matter halos. Analyzing the snapshot data from the Massive
Neutrino Simulations, we numerically obtain the number densities of distinct
halos identified via the SPARTA code as a function of their splashback masses
at various redshifts for two different cases of $\sum m_{\nu}=0.0\,{\rm eV}$
and $0.1\,{\rm eV}$. Then, we fit the numerical results to the recently
developed analytic formula characterized by the diffusion coefficient that
quantifies the degree of ambiguity in the identification of the splashback
boundaries. Our analysis confirms that the analytic formula works excellently
even in the presence of neutrinos and that the decrement of its diffusion
coefficient with redshift is well described by a linear fit, $B(z-z_{c})$, in
the redshift range of $0.2\leq z\leq 2$. It turns out that the massive
neutrino case yields significantly lower value of $B$ and substantially higher
value of $z_{c}$ than the massless neutrino case, which indicates that the
higher masses the neutrinos have, the more severely the splashback boundaries
become disturbed by the surroundings. Given our result, we conclude that the
total neutrino mass can in principle be constrained by measuring how rapidly
the diffusion coefficient of the splashback mass function diminishes with
redshifts at $z\geq 0.2$. We also discuss the anomalous behavior of the
diffusion coefficient found at lower redshifts for both of the $\sum m_{\nu}$
cases, and ascribe it to the fundamental limitation of the SPARTA code at
$z\leq 0.13$.

Unified Astronomy Thesaurus concepts: Large-scale structure of the universe
(902); Cosmological models (337)

##  1 Introduction

The dark matter halos form through gravitational collapse of the initially
overdense sites. As the gravitational collapse is a highly nonlinear process
occurring after the shell crossing, the linear or even higher-order
perturbation theory is not eligible for the analytical description of the halo
formation. Only for the extreme case that the collapse of an isolated
overdense site occurs completely in a spherically symmetry way, the halo
formation can be analytically described by the top-hat spherical dynamics,
which basically predicts that the virialized halos form from the sites whose
linearly extrapolated density contrast exceeds a certain value, called the
virial density threshold (Gunn & Gott, 1972). A unique aspect of the top-hat
spherical dynamics is that the virial density threshold is independent of halo
mass and quite insensitive to the initial conditions of the universe (e.g.,
Press & Schechter, 1974; Lahav et al., 1991; Eke et al., 1996; Pace et al.,
2010).

In realistic cases where the overdense sites are usually not isolated and the
gravitational collapse does not proceed in a spherically symmetric way, the
simple analytical description of the halo formation in terms of a constant
value of the virial density threshold is no longer attainable (Bond & Myers,
1996; Sheth & Tormen, 1999; Chiueh & Lee, 2001; Sheth et al., 2001; Tinker et
al., 2008). Nevertheless, several theoretical works based on $N$-body
simulations indicated that the condition for the realistic formation of dark
matter halos can still be expressed in terms of the virial density threshold
by treating it as a mass-dependent stochastic variable rather than a
deterministic constant value (e.g., Robertson et al., 2009; Maggiore & Riotto,
2010a, b; Corasaniti & Achitouv, 2011a; Corasaniti, & Achitouv, 2011b). The
mass-dependence of the mean value of the virial density threshold stems from
the departure of the gravitational collapse from the top-hat spherical
dynamics, while its stochastic nature reflects the ambiguity in the
identification of non-isolated halos.

In the high-mass limit corresponding to the massive cluster halos, the virial
density threshold becomes less non-spherical but more stochastic (e.g.,
Robertson et al., 2009; Maggiore & Riotto, 2010a, b). Since the gravitational
collapse of highest density peaks occur more or less in a spherically
symmetric way (Bernardeau, 1994), the mean value of the virial density
threshold converges to the spherical value, $\delta_{sc}=1.686$, (Gunn & Gott,
1972). On the other hand, the massive cluster halos often reside at the nodes
of cosmic filaments and thus are more vulnerable to the disturbance from the
surroundings. In consequence, it becomes more ambiguous to identify their
clear-cut boundaries, which results in a higher degree of stochasticity of the
virial density threshold in this limit (Maggiore & Riotto, 2010b).

Very recently, it has been suggested that the ambiguity in the cluster
identification can be considerably reduced by using the splashback radius as a
physical boundary instead of the virial counterpart (Ryu & Lee, 2021). The
splashback radius of a DM halos was originally found in $N$-body simulations
as its outer limit where the density profile exhibits an abrupt drop (Wang et
al., 2009; Wetzel et al., 2014; Diemer & Kravtsov, 2014; Adhikari et al.,
2014; Diemer, 2021). It corresponds to the radial distance to the apocenter of
the elliptical orbit of the latest infalling particle (Adhikari et al., 2014;
More et al., 2015), distinguishing between infalling and orbiting particles,
and naturally taking into account the fly-by objects (More et al., 2015). If
the halo boundary is defined by its splashback radius, which is in fact a
solution to the self-similar spherical infall model (Fillmore & Goldreich,
1984; Bertschinger, 1985), the unrealistic assumption of isolated overdense
region is no longer required to describe the halo formation, since the
splashback boundary generically incorporates the pseudo-evolution of halo mass
(Diemer et al., 2013).

Noting that the self-similar spherical infall model gives a lower mean value
of the splashback density threshold (Shapiro et al., 1999), Ryu & Lee (2021)
substituted it for the virial counterpart in the the generalized excursion set
formalism (Maggiore & Riotto, 2010a, b; Corasaniti & Achitouv, 2011a;
Corasaniti, & Achitouv, 2011b) and proposed an analytic formula for the
splashback mass function characterized by a single parameter, diffusion
coefficient, which measures the degree of the stochasticity of the splashback
density threshold. Testing the analytically derived splashback mass functions
against the $N$-body results for the Planck and WMAP7 cosmologies (Planck
Collaboration et al., 2014; Komatsu et al., 2011), Ryu & Lee (2021) verified
its validity in a wide mass range up to the redshift $z\sim 3$. Showing that
the redshift evolution of the diffusion coefficient sensitively depends on the
matter density parameter $\Omega_{m}$, Ryu & Lee (2021) suggested that the
diffusion coefficient of the splashback mass function of dark matter halos be
a sensitive probe of the initial conditions of the universe.

In this paper, we explore if the diffusion coefficient of the splashback mass
function can effectively constrain the total neutrino mass, $\sum m_{\nu}$,
speculating that the degree of the stochasticity of the splashback density
threshold may be sensitive not only to the amount of dark matter but also to
its nature. The upcoming Sections contain the following contents: a brief
review of the analytic formula for the splashback mass function (Section 2); a
comparison between the analytic formula and the numerically obtained
splashback mass function in the presence of neutrinos and a difference in the
redshift evolution of the diffusion coefficient between the massless and
massive neutrino cases (Section 3); a summary and conclusion (Section 4).
Throughout this paper, we will assume a cosmology where the neutrino ($\nu$)
is present, and the cosmological constant ($\Lambda$) and cold dark matter
(CDM) are dominant, i.e., $\nu\Lambda$CDM cosmology.

##  2 A brief review of the analytic model

In the original excursion set formalism (Press & Schechter, 1974; Bond et al.,
1991), the virial mass function, $dN(M,z)/d\ln M$, defined as the number
densities of DM halos as a function of its logarithmic virial mass, $\ln M$,
at redshift $z$, can be analytically written as

$\frac{dN(M,z)}{d\ln
M}=\frac{\bar{\rho}}{M}\Bigg{|}\frac{d\ln\sigma^{-1}}{d\ln
M}\Bigg{|}f(\sigma)\,,$ (1)

where $\bar{\rho}$ is the mean mass density of the universe, $\sigma(M,z)$ is
the rms fluctuation of the linear density field on the mass scale of $M$ at
redshift $z$, and $f(\sigma)$ is the multiplicity function that counts the
number of the initial density peaks whose linearly extrapolated density
contrast exceeds some virial density threshold when the rms density
fluctuation has the value of $\sigma$. The excursion set theory relates
$dN/d\ln M$ to the background cosmology through $\sigma(M,z)$ that is often
computed as
$\sigma^{2}(M,z)\equiv(2\pi^{2})^{-1}\int\,dk\,k^{2}P(k,z)W^{2}(k,M)$ where
$P(k,z)$ is the linear density power spectrum and $W(k,M)$ is the top-hat
window function on the virial mass scale $M$ (Press & Schechter, 1974; Bond et
al., 1991).

This original excursion set theory was generalized by Maggiore & Riotto
(2010a, 2010b) (hereafter, Maggiore-Riotto) who treated the virial density
threshold as a stochastic random variable whose mean value is
$\delta_{sc}=1.686$ according to the top-hat spherical model (Gunn & Gott,
1972). The resulting Maggiore-Riotto model turned out to describe very well
the virial mass function of cluster halos, despite that it has only one
fitting parameter, called the diffusion coefficient, which characterizes the
degree of the stochasticity of the virial density threshold. Noting the
limited validity of the Maggiore-Riotto formula to the cluster mass scale,
Corasaniti & Achitouv (2011a, 2011b) modified it by taking into account
possible variation of the mean value of the virial density threshold with mass
scale. The modified formula, which has one additional parameter called the
drifting average coefficient, was shown to improve the agreements with the
numerical results on the galactic mass scale.

It was Ryu & Lee (2021) who for the first time proposed that the generalized
excursion set theory can also be used for the evaluation of the splashback
mass function. From here on, we let $M$ denote not the virial but the
splashback mass of a DM halo. They put forth an idea that the splashback
density threshold must be much less stochastic than the virial counterpart and
that its mean value would not change with the mass scale. Then, they claimed
that the splashback multiplicity function could be described by the following
single-parameter formula, which is the same as the Maggiore-Riotto model but
with $\delta_{sc}=1.686$ replaced by $\delta_{sp}=1.52$, the mean splashback
density threshold predicted in the self-similar spherical infall model
(Shapiro et al., 1999):

$\displaystyle f(\sigma;D_{B})$ $\displaystyle=$
$\displaystyle\frac{\delta_{sp}}{\sigma}\sqrt{\frac{2}{\pi(1+D_{B})}}\Biggl{\\{}\left[1-\frac{\kappa}{(1+D_{B})}\right]\exp\left[-\frac{\delta_{sp}^{2}}{2\sigma^{2}(1+D_{B})}\right]$
(2)
$\displaystyle+\frac{\kappa}{2(1+D_{B})}\Gamma\left[0,\frac{\delta_{sp}^{2}}{2\sigma^{2}(1+D_{B})}\right]\Biggl{\\}}\,.$

Here, $\Gamma$ is the incomplete gamma function, $D_{\rm B}$ is the diffusion
coefficient defined as the ratio of the variance of the splashback density
threshold to the mass variance, and $\kappa$ is related to the cross
correlation of the linear density field on two different splashback mass
scales $M$ and $M^{\prime}$ as111In the work of Ryu & Lee (2021), they adopted
the same approximation of $\kappa\approx 0.475$ as used in Corasaniti, &
Achitouv (2011b), assuming a $\Lambda$CDM cosmology. However, in the current
work, we rigorously calculate $\kappa$ by Equation (3).:

$\kappa\approx\frac{\sigma^{2}(M)}{\sigma^{2}(M^{\prime})\left[\sigma^{2}(M)-\sigma^{2}(M^{\prime})\right]}\left[\langle\delta(M)\delta(M^{\prime})\rangle-\sigma^{2}(M^{\prime})\right]\,.$
(3)

Comparing Equations (1)-(3) with the numerical results from the Erebos
$N$-body experiment (Diemer & Kravtsov, 2014, 2015) and determining $D_{\rm
B}$ as a function of redshift, Ryu & Lee (2021) found the following. First,
despite that it has only one parameter, the analytic formula excellently works
in the splashback mass range of $5\times 10^{12}\leq M/(h^{-1}\,M_{\odot})\leq
10^{15}$ up to $z\sim 3$ for two different $\Lambda$CDM cosmologies. Second,
the splashback mass function is characterized by a much lower value of $D_{\rm
B}$ than the virial mass function, which justifies the assumption that the
splashback density threshold is less stochastic than the virial counterpart.

Third, the diffusion coefficient almost monotonically decreases with $z$, well
approximated by a linear fit, and vanishing at some critical redshift,
$z_{c}$. This result implied that the splashback mass function at $z\geq
z_{c}$ can be described by the following purely analytic formula with no free
parameter:

$\frac{dN(M,z)}{d\ln
M}=\frac{\bar{\rho}}{M}\Bigg{|}\frac{d\ln\sigma^{-1}}{d\ln
M}\Bigg{|}\frac{\delta_{sp}}{\sigma}\sqrt{\frac{2}{\pi}}\Bigg{\\{}\left(1-\kappa\right)\exp\left(-\frac{\delta_{sp}^{2}}{2\sigma^{2}}\right)+\frac{\kappa}{2}\Gamma\left[0,\frac{\delta_{sp}^{2}}{2\sigma^{2}}\right]\Bigg{\\}}\,.$
(4)

It was also found by Ryu & Lee (2021) that the critical redshift, $z_{c}$,
appears to sensitively depend on the initial conditions, exhibiting a tendency
to have a lower value for the case of the less amount of dark matter.

##  3 The effect of massive neutrinos on the splashback mass function

We are going to numerically tackle the following two issues. First, is the
analytic model reviewed in Section 2 still valid for the $\nu\Lambda{\rm CDM}$
cosmology? Second, what effect does the presence of neutrinos has on the
evolution of the diffusion coefficient, $D_{\rm B}(z)$? To address these
issues, we utilize the snapshot data from the Cosmological Massive Neutrino
Simulations (MassiveNuS) (Liu et al., 2018) for which the effect of neutrinos
was incorporated into the background with the help of the analytic linear
response approximation (Ali-Haïmoud & Bird, 2013). In the periodic box of
linear size $512\,h^{-1}$Mpc, the MassiveNuS contained a total of $N_{\rm
dm}=1024^{3}$ DM particles with individual mass $10^{10}\,h^{-1}\,M_{\odot}$
(Liu et al., 2018).

The snapshot data from the MassiveNuS are available only for two
$\nu\Lambda{\rm CDM}$ models which share the same matter density parameter
($\Omega_{m}$) and same large-scale amplitude of the linear density power
spectrum ($A_{s}$) but different total neutrino mass, $\sum m_{\nu}/\,{\rm
eV}=0.0$ and $0.1$, respectively. Table 1 lists the values of $\Omega_{m}$,
$A_{s}$ and $\sigma_{8}$, for the two models. Note that $\sigma_{8}$ has a
slightly lower value in the presence of massive neutrinos due to the
suppression of the small-scale powers caused by the neutrino free streaming
(Lesgourgues & Pastor, 2014). The catalogs of distinct halos and their
subhalos identified by the Rockstar algorithm (Behroozi et al., 2013) in the
snapshot data are also publicly available at the MassiveNuS website
222http://astronomy.nmsu.edu/aklypin/SUsimulations/MassiveNuS/.

We apply the publicly available SPARTA (Subhalo and PARticle Trajectory
Analysis) algorithm (Diemer, 2017) to a total of $45$ snapshots and rockstar
catalogs in the redshift range of $0\leq z\leq 5.285$ from the MassiveNuS for
the two $\nu\Lambda{\rm CDM}$ cosmologies. The redshift interval, $\Delta z$,
between the adjacent snapshots in the redshift range of $z\leq 1$ is
approximately $\delta z\approx 0.04$. The SPARTA basically inspects the orbits
of all particles in each rockstar halo over redshifts to compute the distances
to the apocenters of their first orbits, $r_{sp}$, and takes the average over
them, $R^{mn}_{sp}$, to define the splashback radius of each halo. The SPARTA
algorithm also provides several different definitions of the splashback
radius, $R^{50}_{sp}$, $R^{70}_{sp}$, and $R^{90}_{sp}$, which correspond to
the $50\%$, $70\%$, and $90\%$ percentiles of $r_{sp}$, respectively. The
analytic formula reviewed in Section 2 was proven to describe well the
splashback mass function for the case of the $\Lambda$CDM cosmology, no matter
which definition of $R_{sp}$ was used. In the current work, we exclusively use
$R^{mn}_{sp}$ for the halo splashback radius, since the linear fit
approximation to the evolution of the diffusion coefficient turned out to work
best for the case of $R^{mn}_{sp}$ (Ryu & Lee, 2021).

To determine the splashback mass function, we consider only the distinct halos
that are not embedded within the splashback radii of any larger halos. Note
that the criterion for a distinct halo depends on the halo finder algorithm.
For example, some rockstar object classified as a distinct halo by the
spherical overdensity (SO) code could be classified as a subhalo by the MORIA
code, a subrouine contained in the SPARTA algorithm, as its splashback radius
can be overlapped with that of a neighbor larger halo. Selecting the distinct
halos with splashback masses $M\geq 5\times 10^{12}\,\,h^{-1}\\!M_{\odot}$, we
split the range of $\ln M$ into short bins with equal length of $d\ln M$ and
count the number of the selected halos, $dN$, whose logarithmic masses are in
the range of $\left[\ln M,\ln M+d\ln M\right]$. Then, the splashback mass
function is numerically determined as the ratio, $dN/d\ln M$, at each
redshift, to which the analytic model, Equations (1)-(3), is adjusted by
varying $D_{\rm B}$. Dividing the halo sample into eight subsamples of equal
size, we also repeat the same calculation separately for each of the eight
subsamples to obtain eight different splashback mass functions,
$\\{dN_{i}/d\ln M\\}_{i=1}^{8}$. The one standard deviation scatter among
$\\{dN_{i}/d\ln M\\}_{i=1}^{8}$ is then determined as the Jackknife errors in
the original splashback mass function. As done in Ryu & Lee (2021), the
standard $\chi^{2}$-statistics is employed to find the best-fit value of
$D_{\rm B}$, while the error in $D_{\rm B}$ is computed as the associated
Fisher information.

Figure 1 plots the numerically obtained splashback mass function (black closed
circles) with the Jackknife errors at $z=0$, and compares them with the
analytic formula (red solid line) with the best-fit value of $D_{\rm B}$, for
the two $\nu\Lambda{\rm CDM}$ models in the first and third from the top
panels. The ratios of the numerical result to the analytical formula for each
model (red solid line) with one standard deviation scatter (gray area) are
also shown in the third and first from the bottom panels. As can be seen, for
both of the $\nu\Lambda{\rm CDM}$ models, despite that it has only one
parameter, the analytic formula excellently matches the numerical results in
the wide mass range of $5\times 10^{12}\leq M/(h^{-1}\,M_{\odot})\leq 10^{15}$
at $z=0$.

Repeating the same calculation but at higher redshifts, we trace $D_{\rm
B}(z)$. Figures 2-3 plot the same as Figure 1 but at $z=0.47$ and $1.05$,
respectively, demonstrating how well the analytic formula, Equations (1)-(3)
with fixed $\delta_{sp}=1.52$, describes the splashback mass function even in
the presence of neutrinos over a broad range of $z$. Figure 4 shows how
$D_{\rm B}(z)$ (black filled circles) changes with redshifts for the two cases
of $\sum m_{\nu}$. As can be seen, for both of the cases, $D_{\rm B}(z)$
decreases almost linearly with redshifts at $z\geq 0.2$, showing a trend of
converging to zero at some critical redshift, $z_{c}$. If the diffusion
coefficient vanishes, the analytic formula of the splashback mass function
becomes parameter free. Figure 5 compares the numerical results at three
edshifts beyond $z_{c}$ with the parameter-free model, Equation (4). Although
the numerical results suffer from large Jackknife errors at these high
redshifts due to the low number densities of massive cluster halos, the
splashback mass functions agree quite well with the purely analytic model for
both of the $\nu\Lambda{\rm CDM}$ models at $z>z_{c}$. These result are
consistent with what Ryu & Lee (2021) found for the case of the $\Lambda$CDM
cosmology, confirming the validity of Equations (1)-(4) even in the presence
of neutrinos.

As done in Ryu & Lee (2021), approximating $D_{\rm B}(z)$ in the range of
$0.2\leq z\leq z_{c}$ to a linear fit, $B(z_{c}-z)$, we determine the best
values of $B$ and $z_{c}$ via the standard $\chi^{2}$-statistics, which are
listed in Table 1 . As can be read, the two models differ in the best-fit
values of $B$ and $z_{c}$, the massive neutrino case yielding a substantially
larger value of $z_{c}$ and significantly lower value of $B$ than the massless
neutrino case. This result indicates that the presence of massive neutrinos
has an effect of making $D_{\rm B}(z)$ more diffusive in this redshift range,
slowing down the decrease of $D_{\rm B}(z)$ with $z$. In other words, in the
presence of hot DM particles like massive neutrinos that have large free
streaming lengths, the DM halos experience larger amount of disturbance from
the surrounding, which make it more ambiguous to precisely identify their
physical boundaries.

Unlike the result of Ryu & Lee (2021) that $D_{\rm B}(z)$ linearly decreases
with redshifts in the whole range from $z=0$ to $z_{c}$, however, we find an
anomalous behavior of $D_{\rm B}(z)$ at $z<0.2$ in the presence of neutrinos:
it increases with $z$, deviating from the linear fit at $z<0.2$, as can be
seen in Figure 4. We also note that the massive neutrino case exhibits a much
more conspicuous deviation of $D_{\rm B}(z)$ from the linear fit than the
massless case at these low redshifts. We suspect that it should be caused or
at least linked closely with the inherent limitation of the SPARTA code in
locating the splashback boundaries of DM halos at $z\leq 0.13$ where the
unknown future orbits of infalling particles make it impossible to determine
particle apocenters (private communication with B. Diemer 2022).

##  4 Summary and Conclusion

We have numerically investigated the effect of neutrinos on the splashback
mass functions of DM halos in the mass range of $0.5\leq
M/(10^{13}\,h^{-1}\,M_{\odot})\leq 10^{2}$ by applying the SPARTA code
(Diemer, 2017) to the snapshots of MassiveNuS (Liu et al., 2018) for two
different $\nu\Lambda{\rm CDM}$ cosmologies with total neutrino mass of $\sum
m_{\nu}=0.0\,{\rm eV}$ and $0.1\,{\rm eV}$. Comparing the numerically obtained
splashback mass functions to the analytic single parameter formula proposed by
Ryu & Lee (2021) at various redshifts, we have trailed the evolution of the
single parameter, dubbed the diffusion coefficient $D_{\rm B}(z)$, in the
range of $0\leq z\leq 3$. It has turned out that the analytic formula with
fixed splashback density threshold $\delta_{sp}=1.52$ validly describes the
splashback mass function for both of the $\nu\Lambda{\rm CDM}$ models at
redshifts up to $z\approx 3$ and that the diffusion coefficient decreases
almost linearly with redshifts at $z\geq 0.2$, evanescing at some critical
redshift, which are all in consistent with the claim of Ryu & Lee (2021).

Fitting $D_{\rm B}(z)$ at $z\geq 0.2$ to a linear scaling relation of
$B(z-z_{c})$ for both of the $\nu\Lambda{\rm CDM}$ models, we have determined
the best-fit values of the slope, $B$, and critical redshift, $z_{c}$, and
find that they sensitively depend on the total neutrino mass. The
$\nu\Lambda{\rm CDM}$ model with $\sum m_{\nu}=0.1\,{\rm eV}$ yields a
significantly lower value of $B$ and a substantially higher value of $z_{c}$
than the other model with $\sum m_{\nu}=0.1\,{\rm eV}$, which indicates that
the diffusion coefficient decreases more slowly with $z$ in the presence of
massive neutrinos. Since the diffusion coefficient of the analytic formula is
a measure of the stochasticity of the splashback density threshold, this
result implies that the presence of hot DM particles like massive neutrinos
has an effect of disturbing more severely the splashback boundaries and making
the halo identification more ambiguous. The dependence of $B$ and $z_{c}$ on
$\sum m_{\nu}$ also implies that the diffusion coefficient of the splashback
mass function can in principle be used as a probe of the total neutrino mass.

A failure of the linear fit approximation to $D_{\rm B}(z)$ has been witnessed
for both of the models at $z<0.2$ when the diffusion coefficient exhibit
increment rather than decrement with $z$, an anomalous tendency never detected
in the absence of neutrinos for which case the linear fit always matches the
diffusion coefficient in the entire range of $0\leq z<3$. Given the
fundamental limitation of the SPARTA code at $z\leq 0.13$ that it is incapable
of locating the apocenters of infalling DM particles whose future orbits are
known (private communication with B. Diemer 2022), we suspect that the
anomalous behavior of $D_{\rm B}(z)$ witnessed at $z<0.2$ should not be a real
one caused by the presence of massive neutrinos but a spurious one caused by
the inherent flaw of the SPARTA code at low redshifts.

It is, however, worth discussing a caveat that our conclusion is subject to.
The MassiveNuS employed the analytic linear response approximation (Ali-
Haïmoud & Bird, 2013) to mimic the presence of massive neutrinos rather than
directly including the massive neutrinos as component particles. Although this
approximation has been known to work quite well provided that $\sum
m_{\nu}\leq 0.3\,{\rm eV}$ (Bird et al., 2018), it is not guaranteed to
accommodate all possible effects that the massive neutrinos would have on the
nonlinear evolution of cosmic web. It will require a more rigorous and
accurate treatment of the neutrino effects to investigate the true behavior
$D_{\rm B}(z)$ in the nonlinear regime and its dependence on the total
neutrino mass. It will be also quite desirable to have an analytic model for
$D_{\rm B}(z)$ derived from the first principles, with which one can
quantitatively probe the nature of dark matter including massive neutrinos.
Our future work is in this direction.

We are very grateful to B. Diemer for his many helps with the application of
the SPARTA code to the MassiveNuSdata from the Columbia Lensing group as well
as his enlightening comments, which helped us avoid a misinterpretation of our
results and significantly improve the original manuscript. We thank the
Columbia Lensing group for making their suite of simulated maps available at
the website (http://columbialensing.org), and NSF for supporting the creation
of those maps through grant AST-1210877 and XSEDE allocation AST-140041. We
thank the New Mexico State University (USA) and Instituto de Astrofisica de
Andalucia CSIC (Spain) for hosting the Skies & Universes site for cosmological
simulation products. We thank J.Liu for providing us with the MassiveNuS
snapshot data. We acknowledge the support by Basic Science Research Program
through the National Research Foundation (NRF) of Korea funded by the Ministry
of Education (No.2019R1A2C1083855).

## References

  * Adhikari et al. (2014) Adhikari, S., Dalal, N., & Chamberlain, R. T. 2014, JCAP, 2014, 019
  * Ali-Haïmoud & Bird (2013) Ali-Haïmoud, Y. & Bird, S. 2013, MNRAS, 428, 3375
  * Bond & Myers (1996) Bond, J. R., & Myers, S. T. 1996, ApJS, 103, 1
  * Behroozi et al. (2013) Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013, ApJ, 762, 109
  * Bernardeau (1994) Bernardeau, F. 1994, ApJ, 427, 51
  * Bertschinger (1985) Bertschinger, E. 1985, ApJS, 58, 39
  * Bird et al. (2018) Bird, S., Ali-Haïmoud, Y., Feng, Y., et al. 2018, MNRAS, 481, 1486
  * Bond et al. (1991) Bond, J. R., Cole, S., Efstathiou, G., et al. 1991, ApJ, 379, 440
  * Bond et al. (1996) Bond, J. R., Kofman, L., & Pogosyan, D. 1996, Nature, 380, 603
  * Borzyszkowski et al. (2017) Borzyszkowski, M., Porciani, C., Romano-Díaz, E., et al. 2017, MNRAS, 469, 594
  * Chiueh & Lee (2001) Chiueh, T. & Lee, J. 2001, ApJ, 555, 83
  * Corasaniti & Achitouv (2011a) Corasaniti, P. S. & Achitouv, I. 2011a, Phys. Rev. Lett., 106, 241302
  * Corasaniti, & Achitouv (2011b) Corasaniti, P. S. & Achitouv, I. 2011b, Phys. Rev. D, 84, 023009
  * Diemer et al. (2013) Diemer, B., More, S., & Kravtsov, A. V. 2013, ApJ, 766, 25
  * Diemer & Kravtsov (2014) Diemer, B. & Kravtsov, A. V. 2014, ApJ, 789, 1
  * Diemer & Kravtsov (2015) Diemer, B. & Kravtsov, A. V. 2015, ApJ, 799, 108
  * Diemer (2017) Diemer, B. 2017, ApJS, 231, 5
  * Diemer et al. (2017) Diemer, B., Mansfield, P., Kravtsov, A. V., et al. 2017, ApJ, 843, 140
  * Diemer (2020a) Diemer, B. 2020a, ApJ, 903, 87
  * Diemer (2020b) Diemer, B. 2020b, ApJS, 251, 17
  * Diemer (2021) Diemer, B. 2021, ApJ, 909, 112
  * Eke et al. (1996) Eke, V. R., Cole, S., & Frenk, C. S. 1996, MNRAS, 282, 263
  * Fillmore & Goldreich (1984) Fillmore, J. A. & Goldreich, P. 1984, ApJ, 281, 1
  * García et al. (2021) García, R., Rozo, E., Becker, M. R., et al. 2021, MNRAS, 505, 1195
  * Gao & Theuns (2007) Gao, L. & Theuns, T. 2007, Science, 317, 1527. doi:10.1126/science.1146676
  * Gunn & Gott (1972) Gunn, J. E. & Gott, J. R. 1972, ApJ, 176, 1
  * Jenkins et al. (2001) Jenkins, A., Frenk, C. S., White, S. D. M., et al. 2001, MNRAS, 321, 372
  * Komatsu et al. (2011) Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18
  * Lahav et al. (1991) Lahav, O., Lilje, P. B., Primack, J. R., et al. 1991, MNRAS, 251, 128
  * Lee (2019) Lee, J. 2019, ApJ, 872, L6
  * Lesgourgues & Pastor (2012) Lesgourgues, J., & Pastor, S. 2012, Invited Review article for the Special Issue on Neutrino Physics, 608515; doi:10.1155/2012/608515
  * Lesgourgues & Pastor (2014) Lesgourgues, J. & Pastor, S. 2014, New Journal of Physics, 16, 065002
  * Lewis et al. (2000) Lewis, A., Challinor, A., & Lasenby, A. 2000, ApJ, 538, 473
  * Liu et al. (2018) Liu, J., Bird, S., Zorrilla Matilla, J. M., et al. 2018, JCAP, 2018, 049
  * Maggiore & Riotto (2010a) Maggiore, M., & Riotto, A. 2010a, ApJ, 711, 907
  * Maggiore & Riotto (2010b) Maggiore, M., & Riotto, A. 2010b, ApJ, 717, 515
  * More et al. (2015) More, S., Diemer, B., & Kravtsov, A. V. 2015, ApJ, 810, 36
  * Pace et al. (2010) Pace, F., Waizmann, J. C., & Bartelmann, M. 2010, MNRAS, 1865, 1874
  * Planck Collaboration et al. (2014) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2014, A&A, 571, A16
  * Press & Schechter (1974) Press, W. H., & Schechter, P. 1974, ApJ, 187, 425
  * Robertson et al. (2009) Robertson, B. E., Kravtsov, A. V., Tinker, J., et al. 2009, ApJ, 696, 636
  * Ryu & Lee (2021) Ryu, S. & Lee, J. 2021, ApJ, 917, 98
  * Shapiro et al. (1999) Shapiro, P. R., Iliev, I. T., & Raga, A. C. 1999, MNRAS, 307, 203
  * Sheth & Tormen (1999) Sheth, R. K. & Tormen, G. 1999, MNRAS, 308, 119
  * Sheth et al. (2001) Sheth, R. K., Mo, H. J., & Tormen, G. 2001, MNRAS, 323, 1
  * Sheth & Tormen (2002) Sheth, R. K., & Tormen, G. 2002, MNRAS, 329, 61
  * Tinker et al. (2008) Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709
  * Wang et al. (2009) Wang, H., Mo, H. J., & Jing, Y. P. 2009, MNRAS, 396, 2249
  * Warren et al. (2006) Warren, M. S., Abazajian, K., Holz, D. E., et al. 2006, ApJ, 646, 881
  * Wetzel et al. (2014) Wetzel, A. R., Tinker, J. L., Conroy, C., et al. 2014, MNRAS, 439, 2687

Figure 1: Analytical formula of the splashback mass function (red solid line)
with the diffusion coefficient $D_{\rm B}$ determined by adjusting the formula
to the numerical results (black closed circles with errors) for the two models
at $z=0$.

Figure 2: Same as Figure 1 but at $z=0.47$.

Figure 3: Same as Figure 1 but at $z=1.05$.

Figure 4: Best-fit values of the diffusion coefficient of the splashback mass
function versus redshifts.

Figure 5: Purely analytical parameter free formula of the the splashback mass
function (red solid line) compared with the numerical results (black closed
circles with errors) for the two models at three different redshifts higher
than the critical redshift, $z>z_{c}$.

Table 1: Best-fit values of the critical redshift for the massless and massive neutrino cases. $\sum m_{\nu}$ | $\Omega_{m}$ | $A_{s}$ | $\sigma_{8}$ | $B$ | $z_{c}$
---|---|---|---|---|---
$[\,{\rm eV}]$ |  | $(10^{-9})$ |  |  |
$0.0$ | $0.3$ | $2.1$ | $0.85$ | $-0.032\pm 0.001$ | $2.097\pm 0.047$
$0.1$ | $0.3$ | $2.1$ | $0.83$ | $-0.029\pm 0.001$ | $2.198\pm 0.044$