\section{Introduction}
Networks are a convenient way to represent real-world complex systems, such as social interactions \cite{newman_2003},\cite{onnela_2007}, metabolic interactions \cite{boccaletti_2006}, the Internet/world wide web \cite{vespignani_2001}, transportation systems \cite{barrat_2006,barthelemy_2011}, etc. For several systems it is interesting to find and describe areas of the network which are more densely connected, i.e. the communities of the network. In 20 years of complex networks history extensive work was conducted on community detection in static - non evolving - networks, see \cite{newman_2006,blondel_2008,fortunato_2007} and the review \cite{fortunato_2016} for an overview on community detection in static graphs.

However, many networks have a temporal dimension and need a dynamic mesoscopic description at risk of non-negligible information losses if studied as static networks. Therefore the description of large temporal graphs has been a hot topic of research for the last decade, see the excellent reviews \cite{HolmeSaramaki_2012} and \cite{Holme_2015} for a complete description of temporal networks. Most recently the detection of dynamic communities, that is communities on temporal networks, has become one of the main interests in network science, as temporal networks require to adapt the methods of static community detection. So far no consensual method was found and around 60 methods have been proposed to try to detect dynamic communities evolving with temporal networks. A total of four published reviews try to classify and summarize them \cite{Aynaud_2013}, \cite{Hartmann_2016}, \cite{Masuda_2016} and \cite{Rossetti_2018}.

In \cite{Rossetti_2018}, these methods are classified into 3 main categories: (a) \emph{instant optimal}, (b) \emph{temporal trade-off} and (c) \emph{cross-time}. These methods aim to detect clusters at different times t, i.e. for many snapshots of the temporal network. As these clusters are only dependent on the state of the network at time t, it is then necessary to match the communities at different t with some similarity measures, e.g. Jaccard based \cite{morini_2017, lorenz_2017, greene_2010}, core-node \cite{wang_2008}. Methods in category (b) define clusters at t depending on current and past states of the network. Clusters are incrementally temporally smooth. However such methods are subject to drift as clusters are added up to each other locally. There is no compromise between temporal smoothness and 'optimal' partition at time $t$, see for example \cite{rossetti_2017,guo_2014,gorke_2010,gorke_2013}. Finally, in category (c) clusters at $t$ depend on both past and future states of the network, see \cite{duan_2009,mucha_2010,matias_2016,ghasemian_2016}. Clusters are completely temporally smooth and not subject to drift, but they do not respect causality as communities at $t$ are determined using network's information at $t+n$, i.e. communities at time $t$ can change depending on what comes next, which makes these methods inappropriate for use on-the-fly.

In this article, we present a new tool to achieve a mesoscopic description of dynamic networks, which tries to find a good compromise between 'global' and 'local' methods. We apply our method on two data sets of scientific articles, to show how it can describe the emergence and evolution of scientific disciplines. The main difficulties for meta-community detection methods are twofold: Finding the right temporal smoothing and quantifying the `stability' of communities. It is difficult to distinguish if changes between snapshots are due to structural evolution of the community or algorithm instability, as static community detection methods used at each time $t$ can find different communities for a same topology (see \cite{Rossetti_2018} for a complete description of pros and cons of each clustering category). Here, we propose an algorithm which aims to find a good balance between temporal inertia (smoothness) and `optimal' partition at any given time $t$. We compare this method to the most basic approach, which optimizes the modularity of the aggregated network using the Louvain algorithm \cite{blondel_2008}. The latter can be assimilated to a category (c) method in \cite{Rossetti_2018}. We then describe the methods to analyze differences between partitions: mutual information (MI) measures and bipartite network (BN) representations. We show that MI based measures are interesting but give a limited amount of information on how different two partitions are, whereas bipartite network representation allows to see how streams split between partitions. We used the methods on two bibliographic data sets: (1) the scientific publications of a scientific institution, ENS Lyon and (2) publications related to the emergence of a new mathematical tool, the 'wavelets'. We show that the global approach represents a good approximation when the dynamics is simple, i.e. when there are mainly parallel streams without much interaction, as in the ENS Lyon case. However, when the dynamics is more complex (and interesting), when many communities are born, die, split or merge, one needs a more sophisticated approach, and we show that our algorithm performs well compared to an expert-based dynamics.

\section{Methods}
We start by presenting the two building blocks used in the algorithms we want to compare: how we define and partition a Bibliographic Coupling (BC) network and how we match clusters from successive time periods to create `streams' (temporal meta-clusters). We then introduce the \emph{BiblioMaps} platform used for visualizing dynamic communities and finally we present two standard methods to compare partitions derived from the different methods : Normalized Mutual Information and Bipartite Networks.

\subsection{Bibliographic Coupling partitioning}
Given a set of publications on a given period, a Bibliographic Coupling (BC) network can be defined based on the relative overlap between the references of each pair of publications. More specifically, we compute Kessler's similarities $\omega_{ij}=R_{ij}/\sqrt{R_{i}R_{j}}$, where $R_{ij}$ is the number of shared references between publications $i$ and $j$ and $R_{i}$ is the number of references of publication $i$. In the BC network, each publication corresponds to a node and two publications $i$ and $j$ share a link of weight $\omega_{ij}$. If they don't share any reference, they are not linked ($\omega_{ij}=0$); if they have an identical set of references, their connexion has a maximal weight ($\omega_{ij}=1$). Here, we consider that the link between two publications is only meaningful if they share at least two references and we impose $\omega_{ij}=0$ if they share only one reference.\\

We use weighted links to reinforce the dense (in terms of links per publication) regions of the BC networks. This reinforcement facilitates the partition of the network into meaningful groups of cohesive publications, or communities. We measure the quality of the partition with the \emph{modularity} Q (eq. \ref{eq1}), a quantity that roughly compares the weight of the edges inside the communities to the expected weight of these edges if the network were randomly produced:

\begin{equation}
Q = \frac{1}{2\Omega} \sum_{i,j} \left[ \omega_{i,j} - \frac{\omega_i \omega_j}{2\Omega} \right]\delta\left(c_i, c_j \right),
\label{eq1}
\end{equation}

\noindent where $\omega_i = \sum_j \omega_{ij}$ is the sum of the weights of the edges linked to node i, $c_i$ and $c_j$ are the communities containing respectively nodes $i$ and $j$, $\delta$ is the Kronecker function ($\delta (u,v)$ is 1 if $u = v$ and 0 otherwise) and $\Omega = \frac{1}{2}\sum_{ij} \omega_{ij}$ is the total weight of edges. We compute the graph partition using the efficient heuristic algorithm presented in \cite{blondel_2008}.

\subsection{Matching communities from successive time periods}
Given the sets of communities $\{C_1^t, ...C_{k_t}^t \}$ in each time window $t$, the problem at hand is to identify a set of relevant historical communities, or \emph{streams}, that correspond to a chain of communities from successive time periods (at most one per period). In order to decide which community of a given period should be added to a chain of communities from previous periods, we need to use some measure of the similarity between communities from different time periods. A standard measure is the Jaccard index, which computes the proportion of shared nodes between clusters of successive and overlapping periods (see e.g. \cite{gingras_2016,morini_2017}). One drawback of this measure is the need to use overlapping periods which implies that there is no bijection between the publications and the streams (a given publication can belong to several streams).\\

Here, we take advantage of the fact that links can be computed between nodes from different time periods (publications from different periods can have common references). We could for example define a similarity measure between two clusters $C_a$ and $C_b$ from different periods either by the total sum of the links between pairs of publications from these clusters $\Omega_{a,b} = \sum_{i \in C_a, j \in C_b} \omega_{i,j}$ or by a normalized version of this sum  $\omega_{a,b}=\Omega_{a,b} /|C_a||C_b|$, which is comprised between 0 and 1. While these two measures may appear quite intuitive, each of them has some drawbacks as well: using $\Omega_{a,b}$ may bias the construction of the streams by linking two `large' (in terms of publications) but rather dissimilar (in terms of shared references) clusters. On the opposite, using $\omega_{a,b}$ may create some biases by linking two very similar clusters of very different sizes rather than the two clusters that have the second-best similarity and have similar sizes. To be coherent with our construction of clusters (maximizing the modularities within each time period), we propose here to use a modularity-based concept to match clusters from successive time periods. The similarity measure we use is thus $\delta Q = \Omega_{a,b} - \Omega_{a}\Omega_{b}/2\Omega_{A,B}$ which corresponds to an increase in the modularity of the BC network built from the two periods $A$ and $B$ when, starting from partitions defined on each period, clusters $a$ and $b$ are merged.\\

{\bf Matching Algorithm}\\
\begin{algorithm}[H]
  Only compare pairs of communities $(a,b)$ with a minimum similarity $\omega_{a,b}> \Theta=10^{-6}$.\\
  Define the best match of each cluster by the one maximizing $\delta Q$.\\
 \For{each temporal window}{
  Define the {\bf predecessor} of each cluster as its best match from the previous time period.\\
  Define the {\bf successor} of each cluster as its best match from the next time period.
  }
 Two clusters are said to be \emph{paired} if they are predecessors / successors of each other.\\
 If a cluster is not the successor of its predecessor, we have a \emph{split}.\\
 If a cluster is not the predecessor of its successor, we have a \emph{merge}.\\
 Streams are defined as chains of paired clusters.\\
\end{algorithm}

An illustration of this algorithm is given in Figure \ref{fig:matching}.

\begin{figure}
\centering
\resizebox*{12cm}{!}{\includegraphics{Figures/fig_matching.png}}
\caption[Matching clusters from successive time periods.]{Matching clusters from successive time periods. We start with the BC network built from clusters of publications detected independently in periods A and B. For each cluster $a$ of period A and cluster $b$ of period B, we compute the modularity change $\delta Q_{ab}$ obtained by grouping these two clusters in the 2-periods BC network. The `successor' of cluster $a$ is defined as the cluster $b$ from period B maximizing this quantity and the `predecessor' of cluster $b$ is defined as the cluster $a$ from period B maximizing this quantity. In the stream visualisation, `paired' clusters (successor / predecessor of each other) are represented on the same y-position, and we only show the BC links between successors or predecessors, which highlight dynamical events such as merge and splits. In these figures, the size of the nodes are proportional to the number of publications in the corresponding clusters and the thickness of the links represent the average weights of the BC links between publications from two clusters.}
\label{fig:matching}
\end{figure}

\subsection{Different algorithms used to define historical streams}
We compare the results of two types of algorithms which build historical communities, or \emph{streams}, starting from publications data sets. The first method is `global', as it considers the whole data set to compute the communities. The second is `local', as it starts from successive windows of $\Delta T$ years and starts by building a mesoscopic description adapted to that specific window. Hereafter, we present the results from the two main variants and refer the reader to the Annex \ref{annex1} for more detailed results.\\

\subsubsection{Global Algorithm (GA)}
The Global Algorithm builds a $global$ BC network by taking into account all the publications in the data set. Streams are defined as time evolution of these (static) communities maximizing the global modularity. Since we are working in a single (large) time period, this approach does not yield any dynamical events such as splitting / merging of communities, but it provides a simple reference.

\subsubsection{Best-Combination Local Communities (BCLC)}
This $local$ algorithm starts by running, for each time period, $N$ independent runs (we used $N=100$) of the Louvain algorithm. Because of the noise inherent to the Louvain algorithm, the best modularity partitions in each time period are not necessarily the ones that best match each other across successive time periods. We thus optimize the inter-period combination by the following algorithm: \\

{\bf BCLC Algorithm}\\
\begin{algorithm}[H]
Compute the Bibliographic Coupling Graph \;
Split the data set into temporal windows $\Delta t$ \;
\For{each of the $N=100$ partitions of the first period and each of the $N=100$ partitions of the second period}{
 Run the matching algorithm to define the 2-periods streams\;
}
Among the $N*N$ defined streams, select the ones maximizing the modularity of the BC network on the first 2 periods. \;
Define the `best combination' partitions of the first 2 periods as those corresponding to those streams\;

\For{each pair of successive temporal windows $A$ and $B$, starting from the second one}{
   \For{each of the $N=100$ partitions of the period $B$}{
      Run the matching algorithm between these partitions and the `best combination' partition of period $A$ (known from a previous step) to define 2-periods streams\;}
 Among the $N$ defined streams, select the ones maximizing the modularity of the BC network on periods $A$ and $B$\;
 Define the `best combination' partition of period $B$ as the one corresponding to those streams
}
  \
\end{algorithm}

Note that maximizing a global indicator over the $T$ periods with $N$ runs would take too long as there would be $N^T$ possibilities to explore. For this reason, we choose the best combinations between the first two periods ($N^2$ checks) and then one period at a time ($N(T-2)$ checks).\\

This algorithm returns temporal streams we call \emph{BCLC-streams}. These streams still maximize the modularity at each time $t$ while using some cross-time information to improve the global modularity. Figure \ref{fig:bestcombinaison} is an illustration of different runs of this algorithm. Choosing the value of the period $T$ is a trade-off. It needs to be long enough so that communities within each period have enough articles to be meaningful and limit clustering variability. But it also needs to be small enough to follow scientific dynamics. For example, the mean cited half-life of scientific articles is close to 7 years \cite{half_life}. After trying different periods, we chose, for our databases, a period {$T = 5$} years.\\

\begin{figure}
\centering
\resizebox*{12cm}{!}{\includegraphics{Figures/fig_bestcombination.png}}
\caption[Choosing the best combination.]{Choosing the best combination. Given a set of clusters from period A, we perform the matching algorithm between these clusters and N=100 sets of clusters from period B obtained by independent runs of the Louvain algorithm. In each case, we compute the modularity $Q$ obtained by grouping paired clusters in a single clusters in the 2-periods BC network. The `best combination' is defined as the set of clusters from period B maximizing this quantity.}
\label{fig:bestcombinaison}
\end{figure}

\subsection{BiblioTools / BiblioMaps}\label{sec:bibliotools}
All the data sets were extracted from the ISI Web of Knowledge Core Collection database\footnote{\textit{http://apps.isiknowledge.com/}}. The bibliographic records were parsed and analyzed using Bibliotools, a Python-based open-source software and the historical streams figures were generated using the web-based visualisation platform BiblioMaps \cite{Grauwin2011,lund2017carte,grauwin2018bibliomaps}. Bibliotools and its extension BiblioMaps were developed by one of us and are available online, as all the data analysis presented in this paper\footnote{\textit{http://www.sebastian-grauwin.com/bibliomaps/}}.

\subsection{Normalized Mutual Information}
The mutual information (MI) is a widely used measure for comparing community detection algorithms. It is defined as a measure of the statistical independence between two random variables (see eq. \ref{eq2}). In other words, if $H(P_{X})$ is the entropy associated with partition $X$ and $H(P_{Y})$ is the entropy associated with partition $Y$ (the entropy is a measure of how partitioned is our network, the more communities - here temporal streams - the higher the entropy), then $MI(P_{X},P_{Y})$ represents the overlap of the two partitions. In layman's terms, it tells us how much we know about the partition $P_{X}$ when the partition $P_{Y}$ is given. You may refer to  \cite{Wagner_2007,Kvalseth_2017} for a deeper description on mutual information. In particular, note that the mutual information is a symmetrical measure

\begin{equation}
MI(P_{X},P_{Y}) = H(P_{Y}) - H(P_{Y}\vert P_{X}) = MI(P_{Y},P_{X})
\label{eq2}
\end{equation}
\

The MI is defined on $[0,+\infty]$, therefore it is difficult to make sense of it without an upper-bound. There exists different ways to normalize the mutual information. The idea is to take into account the entropies of the partitions we consider to gauge the proportion of mutual information between the partitions. Normalizing by the entropy of one of the partition, e.g. $H(P_{X})$ (see eq. \ref{eq3}) measures how much of the partition $P_{X}$ is included in the partition $P_{Y}$. We call this normalized mutual information $NMI_{X}$. If it reaches its maximum value 1, it means that it is possible to retrieve all the information (the partition) of $P_{X}$ from the partition $P_{Y}$. However this measure does not take into account the size of the other partition, $P_{Y}$. A partition $P_{Y}$ where each node would be its own community would make $NMI_{X}$ equals to 1 even though both partitions are very different. This measure then needs to be combined with at least another NMI which takes into account the relative size of both partitions (see eq. \ref{eq4}). Here the mutual information is normalized by $\sqrt{H(P_{X})*H(P_{Y})}$, which shows how much of the two entropies overlap on a scale between 0 and 1. This expresses how similar the partitions are. It is equal to 1 when the partitions are the same. Moreover, this last NMI is symmetrical, so it takes into account both retrieval of $P_{X}$ from $P_{Y}$ and retrieval of $P_{Y}$ from $P_{X}$.

\begin{equation}
NMI_{X}(P_{X},P_{Y}) = \frac{MI(P_{X},P_{Y})}{H(P_{X})}
\label{eq3}
\end{equation}

\begin{equation}
NMI(P_{X},P_{Y}) = \frac{MI(P_{X},P_{Y})}{\sqrt{H(P_{X})*H(P_{Y})}} = NMI(P_{Y},P_{X})
\label{eq4}
\end{equation}
\
While Mutual information based measures give a value of similarity between two partitions, it is not straightforward to analyze. For example, it does not allow to track where the (dis)similarities come from. To allow in depth comparison, we represent pairs of partitions as bipartite networks.

\subsection{Bipartite Network of streams}

To track and quantify differences between partitions $X$ and $Y$, we compute a bipartite network where the $n_{X}^i \in N_{X}$ are the first kind of nodes. They represent the streams $s_{X}^i \in P_{X}$ (hence $\vert N_{X}\vert = \vert P_{X}\vert$). It follows that the second kind of nodes $n_{Y}^j \in N_{Y}$ represent the streams $s_{Y}^j \in P_{Y}$. A weighted directed edge is drawn between $n_{X}^i$ and $n_{Y}^j$ only if their corresponding streams $s_{GA}^i$ and $s_{BCLC}^j$ share articles. For a given pair of nodes ($n_{X}^i$,$n_{Y}^j$) the weights of the two edges between them (one in each direction) are defined in eq.\ref{eq5}. We quantify differences between streams of two partitions from this graph, quantities are given in table \ref{tab:flows}.

\begin{equation}
  \left\{
      \begin{aligned}
			w_{n_{X}^i\rightarrow n_{Y}^j} = \frac{\vert s_{X}^i \cap s_{Y}^j\vert}{\vert s_{X}^i\vert}\\
			w_{n_{Y}^j\rightarrow n_{X}^i} = \frac{\vert s_{Y}^j \cap s_{X}^i\vert}{\vert s_{Y}^j\vert}
      \end{aligned}
    \right.
\label{eq5}
\end{equation}

\section{Data sets}\label{sec:dyncom_data set}
In this section, we present the specificity of each data set and the motivations for studying them. Key statistical data are summarized in table \ref{tab:data sets}.

\begin{table}
\centering
\resizebox{\textwidth}{!}{
\begin{tabular}{rrrrrrrrrr}
\hline
data set & Type & Period & N & $N_{BC}$ & $\rho_{links}$ & $< d >$ & $< w >$ & $Q_{GA}$\\
\hline
Wavelets & Thematic & 1963-2012 & 6,582 & 5,568 & 0.0065 & 35.98 & 0.000719 & 0.677\\
ENS-Lyon & Institution & 1988-2017 & 16,679 & 14,389 & 0.0019 & 27.04 & 0.000175 & 0.919\\
\hline
\end{tabular}}
\caption[Statistics on data sets investigated]{Statistics on data sets investigated in the paper. \emph{Type} is the type of organization data come from. \emph{Period} is the period over which spans the data set. \emph{N} is the number of publications in the data set. \emph{$N_{BC}$} is the number of articles in the BC table. $\rho_{links}$ is the density of links in the BC network. $< d > = (N_{BC}-1)*\rho_{links}$, it indicates the average number of publications a given publication shares references with. $< w >$ indicates the average link weight. \emph{Q} is the modularity of the network using global partitioning (GA).}
\label{tab:data sets}
\end{table}

\subsection{ENS-Lyon Publications data set}
The ENS-Lyon Publications data set contains all publications produced by researchers affiliated to the \emph{\'Ecole Normale Sup\'erieure de Lyon} in natural science fields. It spans the 1988-2017 period and contains 16,679 publications. As for many scientific institutions, its publication records is highly structured by disciplinary academic departments. Here, we compare our temporal clustering methods to a partition that clusters articles according to their authors' laboratories (reference partition, $P_{REF}$).

\begin{sidewaysfigure}
\begin{subfigure}[b]{\textwidth}
\centering
\resizebox*{20cm}{!}{\includegraphics{Figures/avecliens_ENSL_GA_mostfreqAU.png}}
\caption[Historical streams computed from the ENS Lyon natural sciences publications]{Historical streams computed from the ENS Lyon natural sciences publications. Streams were determined using the global algorithm (GA).}
\label{fig:streams_ensl_a}
\end{subfigure}
\end{sidewaysfigure}

\begin{sidewaysfigure}\ContinuedFloat
\begin{subfigure}[b]{\textwidth}
\centering
\resizebox*{20cm}{!}{\includegraphics{Figures/avecliens_ENSL_BCLC_mostfreqAU.png}}
\caption[Historical streams computed from the ENS Lyon natural sciences publications]{Historical streams computed from the ENS Lyon natural sciences publications. Streams were computed using our local method (BCLC).}
\label{fig:streams_ensl_b}
\end{subfigure}
\caption[Historical streams computed from the ENS Lyon natural sciences publications]{Labels on each stream correspond to the most frequent author name in that stream during a given period. Streams with the same color have close research topics (here the proximity of streams to each other is computed from the weight of BC network links between clusters of a same period). Bar height is proportional to the number of publications in a given year. Links between streams show the streams that are preceding/following each other.}
\label{fig:streams_ensl}
\end{sidewaysfigure}

\subsection{Wavelets Publications data set}
The Wavelets Publications data set contains all publications related to wavelets and spans from 1910 to 2012 (however the period before 1960 contains only a few publications). This data set contains 6,582 publications, corresponding to all the publications of a list of 83 key actors in the field of wavelets selected by expert advice and bibliographic searches (for more details, see \cite{morini_2017}). The study of this data set represents a difficult task because it emerged from the collaboration of several research fields, constituted by many entangled subfields. Based on the knowledge of one of the authors (PF), a field's expert, we built manually a temporal partition drawing the history of wavelets. We refer to this partition as $P_{REF}$ and compare our automatically generated partitions to this partition of reference. We acknowledge that this partition is not an absolute ground truth as it relies on the subjectivity of an expert. However, we assume that this reference gives a reasonable picture of the field’s evolution.

\begin{sidewaysfigure}
\begin{subfigure}[b]{\textwidth}
\centering
\resizebox*{23cm}{!}{\includegraphics{Figures/avecliens_wav_GA_mostfreqAU.png}}
\caption[Historical streams computed from the wavelets field of research publications]{Historical streams computed from the wavelets field of research publications. Streams determined using the global algorithm (GA).}
\label{fig:streams_wavelets_a}
\end{subfigure}
\end{sidewaysfigure}

\begin{sidewaysfigure}\ContinuedFloat
\begin{subfigure}[b]{\textwidth}
\centering
\resizebox*{23cm}{!}{\includegraphics{Figures/avecliens_wav_BCLC_mostfreqAU.png}}
\caption[Historical streams computed from the wavelets field of research publications]{Historical streams computed from the wavelets field of research publications. Streams were computed using our local method (BCLC).}
\label{fig:streams_wavelets_b}
\end{subfigure}
\end{sidewaysfigure}

\begin{sidewaysfigure}\ContinuedFloat
\begin{subfigure}[b]{\textwidth}
\centering
\resizebox*{23cm}{!}{\includegraphics{Figures/avecliens_wav_REF_mostfreqAU_thr05pc.png}}
\caption[Historical streams computed from the wavelets field of research publications]{Historical streams computed from the wavelets field of research publications. Streams are the reference streams determined manually.}
\label{fig:streams_wavelets_c}
\end{subfigure}
\caption[Historical streams computed from the wavelets field of research publications]{Labels on each stream correspond to the most frequent author name in that stream during a given period. Streams with the same color have close research topics (here the proximity of streams to each other is computed from the weight of BC network links between clusters of a same period). Bar height is proportional to the number of publications in a given year. Links between streams show the streams that are preceding/following each other.}
\label{fig:streams_wavelets}
\end{sidewaysfigure}

\section{Results}

It is difficult to represent the richness of the information conveyed by streams in paper figures. To be able to attribute scientific meaning to each of the streams, and characterize them through their main authors, references, keywords... an interactive stream visualization is available at {\textit{ http://www.sebastian-grauwin.com/streams/BCstreams.html}}.

\subsection{General features}

As illustrated in Figures \ref{fig:streams_ensl} and \ref{fig:streams_wavelets}, the global method cannot lead to a rich dynamics description. By construction, GA streams are well separated from each other (Figure \ref{fig:streams_ensl_a}) and show only a few links in Figure \ref{fig:streams_wavelets_a}, which could be interpreted as splits or merges of subfields. On the opposite, BCLC streams lead to a more dynamical history for both data sets. There are only a few links Figure \ref{fig:streams_ensl_b}, because different streams correspond to different scientific (sub)disciplines, which are known to be only marginally connected. However, our method rightly spots teams that split to focus on different research topics (as streams `Blichert-Toft' and `Lecuyer', 6$^{th}$ and 7$^{th}$ from the top). Similarly, many splits and merges occur in Figure \ref{fig:streams_wavelets_b}. To analyze these differences, we compare the partitions for each data set using the two measures defined above : Mutual information (table \ref{tab:NMIs}) and community similarity from the bipartite network representation.

\begin{table}
\centering
\begin{tabular}{rrr}
\hline
Measures & ENS-Lyon & Wavelets \\
\hline
$\vert P_{GA}\vert$ & 57 & 27 \\
$\vert P_{BCLC}\vert$ & 97 & 36 \\
$\vert P_{REF}\vert$ & 17 & 36 \\
\hline
$H(P_{GA})$ & 3.63 & 2.87 \\
$H(P_{BCLC})$ & 4.05 & 3.04 \\
$H(P_{REF})$ & 2.37 & 3.18 \\
\hline
$MI(GA,REF)$ & 1.93 & 2.03 \\
$MI(BCLC,REF)$ & 1.93 & 2.49 \\
$MI(GA,BCLC)$ & 3.10 & 1.90 \\
\hline
$NMI_{REF}(GA,REF)$ & 0.82 & 0.64 \\
$NMI_{REF}(BCLC,REF)$ & 0.81 & 0.80 \\
$NMI(GA,BCLC)$ & 0.81 & 0.64 \\
\hline
\end{tabular}
\caption[Entropies and Mutual Information measures]{$\vert P_{X}\vert$ is the number of streams in partition $X$. $H(P_{X})$ is the entropy of partition $X$. $MI(P_{X},P_{Y})$ is the mutual information between the partitions $X$ and $Y$. $NMI_{REF}$ is the mutual information MI normalized by $H(P_{REF})$. $NMI(GA,BCLC)$ is the symmetrical normalized mutual information (normalized by $\sqrt{H(GA)*H(BCLC)}$).}
\label{tab:NMIs}
\end{table}

\subsection{Results on ENS Lyon data set}
Table \ref{tab:NMIs} shows the highly different number of streams of each partition : 57 streams for the global method, 97 for the local one and only 17 for the reference partition (the 17 laboratories of the ENS Lyon). The high values of $NMI_{REF}(GA,REF)$ (0.82) and $NMI_{REF}(BCLC,REF)$ (0.81) suggest that the extra streams in both $P_{GA}$ and $P_{BCLC}$ are mostly hierarchical subdivisions of the laboratory streams from $P_{REF}$. A partition being a subdivision of another does not result in a decrease of MI between them. The MI decreases only if communities of a partition need to be mixed to become communities of another. These results suggest that $P_{GA}$ and $P_{BCLC}$ are merely a smaller-scale division of $P_{REF}$. Similarly, the high value of $NMI(GA,BCLC)$ (0.81) suggests that $P_{BCLC}$ and $P_{GA}$ convey the same information.

The measures from Table \ref{tab:flows} confirm this analysis. $\overline{1^{st}E}(GA,REF)$ shows that streams from $P_{GA}$ share on average $86 \pm 17\%$ of their articles with a stream from $P_{REF}$ and an average of $3.37\pm 1.76$ streams from $P_{GA}$ are needed to retrieve 80\% of streams from $P_{REF}$. Similar observations can be made for $P_{BCLC}$. Moreover, $\overline{Sum_{80}}(GA,BCLC)$ shows that it takes on average two streams from $P_{BCLC}$ to reach 80\% of streams from $P_{GA}$.

Figure \ref{fig:ensl_bipartite_net} shows a part of the bipartite network between $P_{GA}$ (left) and $P_{BCLC}$ (right) on the ENS Lyon publications data set. The part of the network is centered on nine streams from $P_{GA}$ equivalent to 17 streams from $P_{BCLC}$. It suggests that streams from $P_{REF}$ are not a mix of different streams from $P_{GA}$ or $P_{BCLC}$. They are rather unions of (almost) entire streams, which means that GA and BCLC yield almost the same partitions, but at different scales.

\begin{figure}
\centering
\resizebox*{12cm}{!}{\includegraphics{Figures/ensl_bipartite_ga_BCLC.png}}
\caption[Bipartite network representation of ENS Lyon data set between $P_{GA}$ and $P_{BCLC}$]{Part of the bipartite network representation of ENS Lyon data set. This network shows the links between temporal communities from $P_{GA}$ (left in red) and $P_{BCLC}$ (right in blue). On each node is given the stream ID and the most frequent author name of the temporal community. Size of nodes accounts for the size of the streams, each stream contains at least 20 articles.}
\label{fig:ensl_bipartite_net}
\end{figure}

\subsection{Results on Wavelets data set}

\subsubsection{Overall comparison}
Describing the history of the wavelets research field is a complicated task as it was born from the collaboration of multiple fields and sub-fields. The values from table \ref{tab:NMIs} show that, even though partitions have a similar number of streams (27 for $P_{GA}$ and 36 for $P_{BCLC}$), there are significant differences between the local and global method. In this case, $NMI(BCLC,REF)$ is significantly higher than $NMI(GA,REF)$ (0.82 vs. 0.68). Moreover $NMI(GA,BCLC)$ is rather low (0.64) which suggests that differences do not only arise from differences of scale. We visualize some of these differences in section \ref{example:wavelets}. From Table \ref{tab:flows} we see that most similar streams between $P_{GA}$ and $P_{REF}$ share $75\% \pm20\%$ of articles on average, whereas the corresponding figure for $P_{BCLC}$ and $P_{REF}$ is $87\% \pm17\%$.

\begin{table}
\centering
\begin{tabular}{rrr}
\hline
Measures & ENS-Lyon & Wavelets \\
\hline
\multirow{2}{*}{$\overline{1^{st}E}(GA,REF)$} & $0.86\pm0.17$ & $0.75\pm0.20$ \\ & $0.49\pm0.20$ & $0.81\pm0.17$ \\\cline{2-3}
\multirow{2}{*}{$\overline{Sum_{80}}(GA,REF)$} & $1.26\pm0.54$ & $1.88\pm0.93$ \\ & $3.37\pm1.76$ & $1.5\pm0.73$ \\\cline{2-3}
\hline
\multirow{2}{*}{$\overline{1^{st}E}(BCLC,REF)$} & $0.89\pm0.14$ & $0.87\pm0.17$ \\ & $0.49\pm0.26$ & $0.87\pm0.15$ \\\cline{2-3}
\multirow{2}{*}{$\overline{Sum_{80}}(BCLC,REF)$} & $1.23\pm0.44$ & $1.26\pm0.50$ \\ & $4.87\pm3.35$ & $1.31\pm0.57$ \\\cline{2-3}
\hline
\multirow{2}{*}{$\overline{1^{st}E}(GA,BCLC)$} & $0.74\pm0.23$ & $0.72\pm0.23$ \\ & $0.85\pm0.16$ & $0.83\pm0.19$ \\\cline{2-3}
\multirow{2}{*}{$\overline{Sum_{80}}(GA,BCLC)$} & $1.96\pm1.14$ & $1.88\pm0.96$ \\ & $1.34\pm0.51$ & $1.61\pm0.83$ \\\cline{2-3}
\hline
\end{tabular}
\caption[Bipartite graph measures]{In this table each cell contains two lines. Each measure $M(X,Y)$ is made on edges. The first line correspond to $M$ measured on edges from $n_X$ to $n_Y$ and the second line corresponds to $M$ being measured on edges from $n_Y$ to $n_X$. So, the first row in $\overline{1^{st}E}(X,Y)$ is the average proportion of articles $n_X$ shares with $n_Y$ $\pm$ its standard deviation. The second row is the average proportion of articles $n_Y$ shares with $n_X$ $\pm$ its standard deviation. For instance, for the ENS-Lyon, this means that streams of $P_{GA}$ share on average 86\% of their articles with their most similar stream in $P_{REF}$, whereas streams from $P_{REF}$ only share on average 49\% of their articles with their most similar stream in $P_{GA}$. $\overline{Sum_{80}}(X,Y)$ is the average number of streams from $P_{Y}$ it takes to retrieve 80\% of the streams' articles from $P_{X}$. For example in the case of the Wavelet data set, on average $1.88\pm 0.96$ streams from $P_{BCLC}$ are needed to retrieve 80\% of a stream from $P_{GA}$.}
\label{tab:flows}
\end{table}

\subsubsection{Examples of major differences}\label{example:wavelets}
We now show some major differences between $P_{GA}$ and $P_{BCLC}$ for the wavelets data set. From Figure \ref{fig:wavelet_bipartite_ga_BCLC}, we can see two types of differences between partitions: scale differences (e.g. g\_9 with s\_51 and s\_166) as in the ENS Lyon case; and more significant differences, when fractions of $P_{GA}$ streams have to be combined to retrieve $P_{BCLC}$ streams (for example, the group of streams around g\_5 and g\_10). Interestingly, g\_7 combines scale and mixing differences. Looking at the BN representation of these $P_{BCLC}$ streams with corresponding $P_{REF}$ streams (Figure \ref{fig:wavelet_bipartite_BCLC_ref}), we see that our $P_{BCLC}$ description is quite similar to the reference description. There are more `stream-to-stream' equivalences, represented by the double arrow on each side of the edge linking streams. Note that, though $P_{BCLC}$ is closer to $P_{REF}$, there are still scale differences (e.g. s\_21, s\_111) and mixing differences (e.g. s\_85, s\_52).

To understand the origin of the better match of $P_{BCLC}$ to the reference, it is instructive to inspect some of the differences between the local ($P_{BCLC}$) and the global partition $P_{GA}$. Let's look first at the difference between the global method stream g\_7, which corresponds to a merger of four local streams, among which s\_21 and s\_53 (see Fig. \ref{fig:wavelet_bipartite_ga_BCLC}). Figure \ref{fig:streams_wavelets_b} shows that these streams do not belong to the same time period. Stream s\_53 corresponds to the bottom stream (labelled 'Frisch'), and represents early works on wavelets, from 1963 to 1994, focusing on multi fractal analysis and turbulence. The second stream, s\_21 (1987 - 2014, 6th from the top, labelled 'Arneodo'), addressed similar issues in a first period and then, since the early 90's, enlarged the subject matter to include mathematical formalization, together with new applications beyond turbulence, such as genome characterization. Our method appropriately distinguishes these two streams, which correspond to different subfields. The second difference relates to the evolution of one of the authors’ (PF) activities. In the global approach, most PF articles belong to a single cluster that gathers papers in signal representations, and especially time-frequency representations that have been at the heart of his works over the years (third stream starting from the bottom in Fig\ref{fig:streams_wavelets_a}). This is a good approximation, but a finer description of the subjects addressed by PF during his career include three topics: (a) time-frequency methods per se, (b) relations of these methods with wavelets and (c) wavelet methods related to self-similarity, in domains such as turbulence. Using the interactive stream visualization, it is possible to look for $BCLC$ streams containing PF's publications. One finds three streams, addressing the three topics described above, and corresponding to (the stream label refers to its position in Fig\ref{fig:streams_wavelets_b}, starting from the top) respectively streams 3, 8 and 6. These two examples suggest that $BCLC$ is able to capture the complexity of a field dynamics’, including relevant subfields, while the global approach tends to merge streams that represent different fields of inquiry.

\begin{figure}[h]
\centering
\centering
\resizebox*{10.5cm}{!}{\includegraphics{Figures/wavelets_bipartite_ga_BCLC.png}}
\caption[Bipartite network representation of the Wavelets data set ($P_{GA}$ and $P_{BCLC}$)]{Part of the bipartite network representation of Wavelets data set. This network shows the links between temporal communities from $P_{GA}$ (left in red) and $P_{BCLC}$ (right in blue). Each node is labelled by the stream ID and the most frequent author name of the temporal community. Node size accounts for stream size.}
\label{fig:wavelet_bipartite_ga_BCLC}
\end{figure}

\begin{figure}\ContinuedFloat
\centering
\centering
\resizebox*{11cm}{!}{\includegraphics{Figures/wavelets_bipartite_ref_BCLC.png}}
\caption[Bipartite network representation of Wavelets data set between $P_{BCLC}$ and $P_{REF}$]{Part of the bipartite network representation of Wavelets data set. This network shows the links between temporal communities from $P_{BCLC}$ (left in blue) and $P_{REF}$ (right in green). Each node is labelled by the stream ID and the most frequent author name of the temporal community. Node size accounts for stream size.}
\label{fig:wavelet_bipartite_BCLC_ref}
\label{fig:wavelet_bipartite}
\end{figure}

\section{Discussion}

We have presented a coherent approach to create a dynamic mesoscopic description of a temporal network. As the standard method to used to create static communities, our method only uses modularity to build the dynamic communities.  We have compared our method to the static (global) approach. We first showed that both methods give the same result for networks with well-separated streams (high modularity), as in the case of ENS-Lyon publications. However, when analyzing data sets with more complex dynamics, as for the birth of wavelets (section \ref{example:wavelets}), our method can generate a more satisfactory dynamics, as compared to an expert reference partition.

Clearly, much more work is needed to develop a standard approach for describing dynamical networks at a mesoscopic scale. The stochastic character of many partitioning algorithms (as Louvain's \cite{blondel_2008}), and the different scales generated by each method make comparisons difficult. Moreover, the dynamical character of the communities renders the definition of an acceptable reference partition even trickier than for static networks.

\newpage
\bibliographystyle{agsm}


\section{Introduction}
Einstein's general relativity is the most elegant and successful gravitational theory subjected to current experimental tests. However, this theory is still imperfect. One of the most important problems of Einstein's gravity is that singularity is inevitable in almost all of its solutions. Since the existence of the singularity in spacetime, causality will break down and predicability of the theory will be lost at this point thus to find the solution for the singularity problem becomes theoretically essential. There are many proposals devoted to solving this problem, one attempt is to construct modified gravitational theories which have the properties of conformal invariance~\cite{Englert:1976ep,narlikar:1977nf,Fradkin:1978yw,
Kaku:1982xu,tHooft:2009wdx,Mannheim:2011ds,Bars:2013yba,Prester:2013fia}.
A gravity theory is conformal invariant if its action is invariant under following transformation of metric
\begin{equation} \label{metrictransformation}
{g_{\mu \nu }} \to {g'_{\mu \nu }} = S \left( x \right){g_{\mu \nu }},
\end{equation}
where $S \left( x \right)$ is a regular conformal factor depending on spacetime coordinate. Obviously Einstein's gravity is not conformal invariant. There are several ways to construct a conformal invariant gravity theory, for example, by introducing an auxiliary field~\cite{Dirac:1973gk} or the conformal invariant Weyl curvature tensor in the action~\cite{Mannheim:2011ds,Mannheim:2016lnx}.

Einstein's gravity is invariant under coordinate transformation and the singularity at the event horizon can be removed by coordinate transformation therefore this singularity is not intrinsic as to Einstein's gravity. Similarly in the context of conformal gravity theory, since the singularity at $r=0$ can in principle be removed by suitable choice of conformal transformation, hence conformal gravity is free from the singularity.
However, the universe we are living in is not conformal invariant and if the conformal gravity theory can describe the real universe then there must be a phase transition at which the conformal symmetry is spontaneously broken. In the broken phase, there are singular and regular metrics with different physical properties due to that the conformal symmetry no longer exists and it was assumed that our nature has a kind of unknown mechanism to select among these different conformal invariant spacetimes~\cite{Toshmatov:2017kmw}. We believe conformal gravity theory is viable and actually, a class of non-singular black hole spacetime solutions from conformal gravity are compatible with the data from the X-ray observation \cite{Bambi:2017yoz,Zhou:2018bxk}.

In this paper, we will study the nonrotating singularity-free black hole which was found in~\cite{Modesto:2016max,Bambi:2016wdn}. The formation and evaporation process of this black hole has been analyzed in~\cite{Bambi:2016yne}, the quasinormal modes of various kinds of perturbations such as scalar, electromagnetic, axial, and gravitational perturbations in this black hole have been studied in~\cite{Toshmatov:2017bpx,Chen:2019iuo,Liu:2020ddo}. Besides, the violation of energy condition for this black hole has been investigated in~\cite{Toshmatov:2017kmw}.
However, the greybody factor of this black hole is barely explored. So far only the greybody factor of the massless scalar field from this black hole has been numerically studied in \cite{Lin:2020rvv} recently. In this paper, we give a rigorous lower bound for the greybody factor of the massless scalar field radiated from the black hole and study the sparsity of the Hawking radiation. Additionally, we investigate the effects of the variation of conformal factor parameters $N$ and $L$ (see eq.(\ref{S(r)})) on the greybody factor and the sparsity of the Hawking radiation and obtain values for conformal parameters which minimize the sparsity of the Hawking radiation numerically. Our study shows that the Hawking radiation of the massless scalar field is extremely sparse from the black hole and become even sparser for the large values of conformal parameters $N$ and $L$.

The paper is organized as follows. In section.\ref{Review} we briefly review the nonrotating singularity-free black hole solution obtained in~\cite{Modesto:2016max,Bambi:2016wdn} and the Regge-Wheeler equation for massless scalar field from this black hole derived in~\cite{Toshmatov:2017bpx}. In section.\ref{HG} we study the power spectrum of Hawking radiation and calculate the rigorous lower bound of the greybody factor for the massless scalar field. Moreover, we give the values for conformal parameters which maximize the lower bound on the greybody factor numerically. Various numerical results for the upper bound of dimensionless parameter $\eta$ used to measure the sparsity of Hawking radiation are given in section.\ref{Sparsity}. Finally, in section.\ref{Summary} we summarize the results of this paper and give a discussion. The geometric units, $\hbar {\rm{ = }}c = {k_{\rm{B}}} = G = 1$, are adopted in this paper.

\section{Non-singular black holes in conformal gravity}\label{Review}
In this section, we briefly review the nonrotating singularity-free black hole solution in conformal gravity~\cite{Modesto:2016max,Bambi:2016wdn}.

The construction of the black hole solution is straightforward and is obtained by simply rescaling the Schwarzschild metric with a conformal factor as follows
\begin{equation} \label{metric}
d{s^2} = S(r)ds_{Schw}^2,
\end{equation}
where $ds_{Schw}^2$ is the line element of Schwarzschild black hole
\begin{equation} \label{Schwarzschild black hole}
ds_{Schw}^2 =  - \left(1-\frac{r_{\text{H}}}{r}\right)d{t^2} + \frac{dr^2}{\left(1-\frac{r_{\text{H}}}{r}\right)} + {r^2}d{\Omega_2 ^2},
\end{equation}
in which $r_{\text{H}}$ is the horizon radius and $S(r)$ is the conformal factor
\begin{equation} \label{S(r)}
S(r) = {\left( {1 + \frac{{{L^2}}}{{{r^2}}}} \right)^{2N}}
\end{equation}
and $N$ is an arbitrary positive integer and $L$ is a new length scale\footnote{Evidently when $N=0$ or $L=0$, the metric (\ref{metric}) reduces to Schwarzschild metric.} which can be either of the order of Planck scale or of the order of the mass $M$ of the black hole. Due to the metric is conformal to the Schwarzschild black hole, the Hawking temperature of the black hole is still $T={1}/{(4\pi r_\text{H})}$. This metric can be the solution of a family of conformally invariant theories since usually solutions of Einstein's gravity theory are a subset of the wider class solutions of conformally invariant theories \cite{Grumiller:2013mxa,Sultana:2012qp,Said:2014lua} and indeed, metric in eq.(\ref{metric}) is obtained by the conformal transformation of Schwarzschild metric. It can be verified that the square of the Riemann tensor ${R_{\alpha \beta \gamma \delta }}{R^{\alpha \beta \gamma \delta }}$ of this spacetime is regular everywhere including $r=0$. Besides the spacetime is geodesically complete in the sense that it will take infinite proper time for a particle to reach $r=0$ (for more details see Ref.\cite{Bambi:2016wdn}).

The Klein-Gordon equation for a massless scalar field is
\begin{equation} \label{KG eq}
\frac{1}{{\sqrt { - g} }}{\partial _\mu }({g^{\mu \nu }}\sqrt { - g} {\partial _\nu }\Phi ) = 0.
\end{equation}
Decomposing $\Phi {\rm{ = }}{e^{ - i\omega t}}\frac{{\psi \left( r \right)}}{{r\sqrt {S(r)} }}Y_m^\ell (\theta ,\phi )$ and using the tortoise coordinate $d{r_*} = \frac{{dr}}{{f(r)}}$, the resulting Regge-Wheeler equation is~\cite{Toshmatov:2017bpx}
\begin{equation} \label{Master equation}
\left( {\frac{{{d^2}}}{{d{r_*}^2}} + {\omega ^2} - {V_{{\rm{eff}}}}} \right)\psi \left( r \right) = 0,
\end{equation}
where the effective potential ${{V_{{\rm{eff}}}}}$ is defined as
\begin{equation} \label{effective potential}
{V_{{\rm{eff}}}} \equiv f(r)\left[ {\frac{{\ell \left( {\ell  + 1} \right)}}{{{r^2}}} + \frac{{\frac{d}{{dr}}\left( {f(r)\frac{d}{{dr}}\left( {r\sqrt {S(r)} } \right)} \right)}}{{r\sqrt {S(r)} }}} \right].
\end{equation}

\section{Hawking radiation and Greybody factor}\label{HG}
As discussed in \cite{Gray:2015pma,Miao:2017jtr,Chowdhury:2020bdi}, the energy emitted per unit time for a 4-dimensional black hole is
\begin{equation} \label{energy emit}
\frac{{dE\left( \omega  \right)}}{{dt}} = \sum\limits_\ell  {{\Gamma _\ell }\left( \omega  \right)} \frac{\omega }{{{e^{{\omega  \mathord{\left/
 {\vphantom {\omega  T}} \right.
 \kern-\nulldelimiterspace} T}}} - 1}}\hat{\vec k} \cdot \hat{\vec n}\frac{{{d^3}kdA}}{{{{\left( {2\pi } \right)}^3}}},
\end{equation}
where $dA$ is the infinitesimal surface area, $\hat{\vec k}$ is the unit tangent vector of the momentum of the scalar field, $\hat{\vec{n}}$ is the unit normal vector of $dA$ and ${{\Gamma _\ell }\left( \omega  \right)}$ is the greybody factor. Since the scalar field we consider is massless thus $| {\vec k} | = \omega$, then the energy emitted per unite time on finite surface $A$ is

\begin{equation} \label{energy emit2}
P = \sum\limits_\ell  {\int {{P_\ell }\left( \omega  \right)d\omega } }
\end{equation}
with
\begin{equation} \label{energy emit3}
{P_\ell }\left( \omega  \right) = \frac{A}{{8{\pi ^2}}}\frac{{{\Gamma _\ell }\left( \omega  \right){\omega ^3}}}{{{e^{{\omega  \mathord{\left/
 {\vphantom {\omega  T}} \right.
 \kern-\nulldelimiterspace} T}}} - 1}}.
\end{equation}
It was shown that $A$ can be taken as 27/4 times the horizon area for Schwarzschild black hole to match high frequency results \cite{Page:1976df,Gray:2015pma} and here we can also take this value since the cross section in geometrical optical limit for massless particle in this model is the same as that of Schwarzschild black hole~\cite{Lin:2020rvv}. Besides, the Hawking temperature $T$ is independent of conformal parameters $N$ and $L$ (see section.\ref{Review}). Therefore, the dependence of power spectrum ${P_\ell }\left( \omega  \right)$ for Hawking radiation on conformal parameters $N$ and $L$ is as same as that of greybody factor ${\Gamma _\ell }\left( \omega  \right)$.

Calculating the analytical expression for greybody factor is generally impractical and usually, it is obtained from numerical computations~\cite{Harris:2004mf,Harris:2005jx,Kanti:2005ja,Duffy:2005ns}. The analytical calculations are available only in some special cases, such as the low frequency limit~\cite{Unruh:1976fm,Rocha:2009xy,Harmark:2007jy,Kanti:2014dxa,
Anderson:2015bza,Sporea:2015hla}, the near extremal limit~\cite{Cvetic:2009jn,Bredberg:2009pv,Chen:2010bsa,Chen:2010as,Chen:2010ywa,Chen:2010yu}, or the large spacetime dimension limit~\cite{Emparan:2013xia,Emparan:2013moa,Guo:2015swu,Emparan:2020inr}. In this paper we follow \cite{Visser:1998ke,Boonserm:2008zg,Boonserm:2014rma} to give a rigorous lower bound for the greybody factor. The bound for greybody factor can be written as follow
\begin{equation} \label{Bound 0}
{{\Gamma _\ell }\left( \omega  \right)} \geqslant {\rm{sec}}{{\rm{h}}^2}\left\{ {\int_{ - \infty }^\infty  {\vartheta d{r_*}} } \right\}
\end{equation}
with
\begin{equation} \label{vartheta1 }
\vartheta  = \frac{{\sqrt {{{\left( {h'\left( r \right)} \right)}^2} + {{\left( {{\omega ^2} - {V_{{\rm{eff}}}} - h{{\left( r \right)}^2}} \right)}^2}} }}{{2h\left( r \right)}},
\end{equation}
where $h(r)$ is an arbitrary positive definite function and satisfy the limits $h( - \infty ) = h( + \infty ) = \omega $. A particularly simple
choice of $h(r)$ is $h(r) = \omega$,
which leads to
\begin{equation} \label{Integrate of vartheta 1 }
\int_{ - \infty }^\infty  {\vartheta d{r_*}}  = \int_{2M}^\infty  {\frac{{{V_{{\rm{eff}}}}}}{{2\omega f(r)}}dr}.
\end{equation}
Combining eq.(\ref{effective potential}) and eq.(\ref{Integrate of vartheta 1 }) we have
\begin{equation} \label{Integrate of vartheta 2}
\int_{ - \infty }^\infty  {\vartheta d{r_*}}  = \frac{{2NM(2N - 1)}}{{\omega {L^2}}}\ln \left( {1 + \frac{{{L^2}}}{{4{M^2}}}} \right) + \frac{{N\left( {3N - 2} \right)}}{{{\rm{2}}\omega L}}\left[ {2\arctan \left( {\frac{{2M}}{L}} \right) - \pi } \right] + \frac{{{{\left( {2N - 1} \right)}^2}}}{{8\omega M}} + \frac{{\ell \left( {\ell  + 1} \right)}}{{4\omega M}},
\end{equation}
where $M$ is the mass of the black hole. Substituting the above equation into eq.(\ref{Bound 0}), we obtain the bound as
\begin{equation} \label{Bound 1}
{{\Gamma _\ell }\left( \omega  \right)} \geqslant{\rm{sec}}{{\rm{h}}^2}\left\{ {\frac{{2NM(2N - 1)}}{{\omega {L^2}}}\ln \left( {1 + \frac{{{L^2}}}{{4{M^2}}}} \right) + \frac{{N\left( {3N - 2} \right)}}{{{\rm{2}}\omega L}}\left[ {2\arctan \left( {\frac{{2M}}{L}} \right) - \pi } \right] + \frac{{{{\left( {2N - 1} \right)}^2}}}{{8\omega M}} + \frac{{\ell \left( {\ell  + 1} \right)}}{{4\omega M}}} \right\}.
\end{equation}
Now let us analyze the dependence of this bound on the parameters $N$ and $L$ in various situations. At first,
for $N=0$ we have
\begin{equation} \label{Schwarzschild bound}
{{\Gamma _\ell }\left( \omega  \right)}\geqslant{\rm{sec}}{{\rm{h}}^2}\left\{ {\frac{{\ell \left( {\ell  + 1} \right)}}{{4\omega M}} + \frac{1}{{8\omega M}}} \right\},
\end{equation}
which is exactly the result for Schwarzschild black hole as it should be. Note that this result has been obtained in \cite{Boonserm:2008zg}.
In addition, when $L$ is large, the bound in eq.(\ref{Bound 1}) can be simplified as follows
\begin{equation} \label{Schwarzschild bound1}
{{\Gamma _\ell }\left( \omega  \right)}\geqslant{\text{sec}}{{\text{h}}^2}\left\{ {\frac{{\ell \left( {\ell  + 1} \right)}}{{4\omega M}} + \frac{{{{\left( {2N - 1} \right)}^2}}}{{8\omega M}}} \right\},
\end{equation}
while in small $L$ limit eq.(\ref{Bound 1}) reduces to eq.(\ref{Schwarzschild bound}).

Furthermore, it is not hard to find out specific values for the conformal parameters $N$ and $L$ that can maximize the lower bound obtained in eq.(\ref{Bound 1}). Specifically it can be done as follows. We first rewrite eq.(\ref{Bound 1}) as
\begin{equation} \label{Bound 2}
{\Gamma _\ell }\left( \omega  \right)\geqslant{\text{sec}}{{\text{h}}^2}\left\{ {{{J\left( {L,N,\ell } \right)} \mathord{\left/
 {\vphantom {{h\left( {L,N,\ell } \right)} {\left( {\omega {r_{\text{H}}}} \right)}}} \right.
 \kern-\nulldelimiterspace} {\left( {\omega {r_{\text{H}}}} \right)}}} \right\},
\end{equation}
where we have defined
\begin{equation} \label{h function}
J\left( {L,N,\ell } \right) \equiv \frac{{N(2N - 1)}}{{{{\left( {L/{r_{\text{H}}}} \right)}^2}}}\ln \left( {1 + \frac{{{L^2}}}{{r_{\text{H}}^{\text{2}}}}} \right) + \frac{{N\left( {3N - 2} \right)}}{{{\text{2}}\left( {L/{r_{\text{H}}}} \right)}}\left[ {2\arctan \left( {\frac{{{r_{\text{H}}}}}{L}} \right) - \pi } \right] + \frac{{{{\left( {2N - 1} \right)}^2}}}{4} + \frac{{\ell \left( {\ell  + 1} \right)}}{2}.
\end{equation}
Note the function $y(x) = {\operatorname{sech} ^2}(x)$ is even, monotonously decreasing when $x>0$, and reaches its maximum at $x=0$, thus at fixed $\omega$, the values of $N$ and $L$ that minimize  $\left| {J\left( {L,N,\ell } \right)} \right|$ will maximize ${\text{sec}}{{\text{h}}^2}\left\{ {{{J\left( {L,N,\ell } \right)} \mathord{\left/
 {\vphantom {{h\left( {L,N,\ell } \right)} {\left( {\omega {r_{\text{H}}}} \right)}}} \right.
 \kern-\nulldelimiterspace} {\left( {\omega {r_{\text{H}}}} \right)}}} \right\}$. Therefore to obtain the values for $N$ and $L$ that can maximize the greybody factor we only need to find out the values that can minimize $\left| {J\left( {L,N,\ell } \right)} \right|$.

In TABLE.\ref{tab:table1} and TABLE.\ref{tab:table2}, for given values of ${L \mathord{\left/
 {\vphantom {L {{r_{\text{H}}}}}} \right.
 \kern-\nulldelimiterspace} {{r_{\text{H}}}}}$($N$), we find out numerically the corresponding $N$(${L \mathord{\left/
 {\vphantom {L {{r_{\text{H}}}}}} \right.
 \kern-\nulldelimiterspace} {{r_{\text{H}}}}}$) which can maximize the lower bound obtained in eq.(\ref{Bound 1}).
\begin{table}[hbp]
  \begin{center}
 \caption{For different values of $L/r_{\text{H}}$ the corresponding values of $N$ which maximize the lower bound on greybody factor with $\ell=1$.}
    \label{tab:table1}
    \begin{tabular}{c|c|c|c|c|c|c|c|c|c}
     \hline \hline
     ${L \mathord{\left/
 {\vphantom {L {{r_{\text{H}}}}}} \right.
 \kern-\nulldelimiterspace} {{r_{\text{H}}}}}$&$0.1$&$0.3$&$0.5$&$0.7$&$0.9$&$1.5$&$2$&$2.5$&$3$\\
      \hline
      $N$&8&8&6&3&2&1&1&1&1\\
      \hline\hline
    \end{tabular}
  \end{center}
\end{table}

\begin{table}[h!]
  \begin{center}
 \caption{For different values of $N$ the corresponding values of $L/{r_{\text{H}}}$ which maximize the lower bound on greybody factor with $\ell=1$.}
    \label{tab:table2}
    \begin{tabular}{c|c|c|c|c|c|c|c|c}
     \hline \hline
     $N$&$1$&$3$&$5$&$7$&$9$&$12$&$15$&$20$\\
      \hline
      ${L \mathord{\left/
 {\vphantom {L {{r_{\text{H}}}}}} \right.
 \kern-\nulldelimiterspace} {{r_{\text{H}}}}}$&2.893&0.785&0.558&0.337&0.456&0.395&0.299&0.257\\
 \hline\hline
    \end{tabular}
  \end{center}
\end{table}

To see the dependence of greybody factor and power spectrum of Hawking radiation on conformal parameters $N$, $L$ and the frequency $\omega$, intuitively, we plot the lower bound of greybody factor and power spectrum of Hawking radiation as a function of frequency for different values of $N$ and $L$. FIG.\ref{figure.1} shows that, when the value of ${L \mathord{\left/{\vphantom {L {{r_{\rm{H}}}}}} \right.\kern-\nulldelimiterspace} {{r_{\rm{H}}}}}$ is fixed, the increase of $N$ leads greybody factor to increase at first (top left) and then to decrease (top right). Similarly the increase of  ${L \mathord{\left/{\vphantom {L {{r_{\rm{H}}}}}} \right.\kern-\nulldelimiterspace} {{r_{\rm{H}}}}}$ will enhance the greybody factor at the beginning (bottom left) and suppress the greybody factor later (bottom right) for a given value of $N$. Obviously FIG.\ref{figure.1} is consistent with TABLE.\ref{tab:table1} and TABLE.\ref{tab:table2}.

Furthermore, as shown in FIG.\ref{figure.2}, for a fixed value of $L/{r_{\rm H}}$, the peak value of power spectrum of Hawking radiation increases as the increase of $N$ in the first place and decreases if $N$ continue to grow. While from the bottom of FIG.\ref{figure.2} one can find that for a relatively small value of $L/{r_{\rm H}}$ the growth of  $L/{r_{\rm H}}$ will cause the peak value of Hawking radiation power spectrum to increase. Once $L/{r_{\rm H}}$ becomes relatively larger the increase of $L/{r_{\rm H}}$ will give rise to the decrease of the peak value of power spectrum when $N$ is fixed. Also, we find that the values of $L/{r_{\rm H}}$ and $N$ which maximize the power spectrum of Hawking radiation are exactly those values for $L/{r_{\rm H}}$ and $N$ that maximize the greybody factor, which is in agreement with our analysis about the dependence of Hawking radiation's power spectrum on $N$ and $L$ above.
\begin{figure}
\centering
\subfigure{
\includegraphics[width=5cm]{ImBDN.eps}
}
\subfigure{
\includegraphics[width=5cm]{ImBDN2.eps}
}
\\
\subfigure{
\includegraphics[width=5cm]{ImBDL.eps}
}
\subfigure{
\includegraphics[width=5cm]{ImBDL2.eps}
}
\DeclareGraphicsExtensions.
\caption{(Top) Plot of the bound of Greybody factor ${{\Gamma _\ell }\left( \omega  \right)}$ as a function of $\omega r_\text{H}$ for different values of $N$ with $\ell=1$ and ${L}/{r_\text{H}}=0.7$; (Bottom) Plot of the bound of Greybody factor ${{\Gamma _\ell }\left( \omega  \right)}$ as a function of $\omega r_\text{H}$ for different values of ${L}/{r_\text{H}}$ with $\ell=1$ and $N=20$.}
\label{figure.1}
\end{figure}
\begin{figure}
\centering
\subfigure{
\includegraphics[width=5.5cm]{ImPSN.eps}
}\subfigure{
\includegraphics[width=5.5cm]{ImPSN2.eps}
}
\subfigure{
\includegraphics[width=5.5cm]{ImPSN3.eps}
}
\\
\subfigure{
\includegraphics[width=5.5cm]{ImPSL.eps}
}
\subfigure{
\includegraphics[width=5.5cm]{ImPSL2.eps}
}
\subfigure{
\includegraphics[width=5.5cm]{ImPSL3.eps}
}
\DeclareGraphicsExtensions.
\caption{(Top) Plots of power spectrum  ${r_{\rm{H}}}{P_\ell }\left( \omega  \right)$ as a function of $\omega r_\text{H}$ for different values of $N$ with $\ell=1$ and ${L}/{r_\text{H}}=0.7$; (Bottom) Plots of power spectrum ${r_{\rm{H}}}{P_\ell }\left( \omega  \right)$ as a function of $\omega r_\text{H}$ for different values of ${L}/{r_\text{H}}$ with $\ell=1$ and $N=20$.}
\label{figure.2}
\end{figure}

\section{The sparsity of Hawking radiation}\label{Sparsity}
The sparsity of Hawking radiation can be measured by the following quantity~\cite{Gray:2015pma}
\begin{equation} \label{Sparsity parameter}
\eta  = \frac{{{\tau _{{\rm{gap}}}}}}{{{\tau _{{\rm{emission}}}}}},
\end{equation}
where ${{\tau _{{\rm{gap}}}}}$ is the average time interval between the emission of two Hawking quanta, which is proportional to inverse of the total Hawking radiation power $P$ in eq.(\ref{energy emit2}). More specifically, it can be computed by~\cite{Gray:2015pma,Paul:2016xvb}
\begin{equation} \label{Average time}
{\tau _{{\rm{gap}}}} = \frac{{{\omega _{{\rm{peak}}}}}}{P}.
\end{equation}
Here ${{\omega _{{\rm{peak}}}}}$ is the value at which ${{{\omega ^3}} \mathord{\left/
 {\vphantom {{{\omega ^3}} {\left( {{e^{\omega T}} - 1} \right)}}} \right.
 \kern-\nulldelimiterspace} {\left( {{e^{\omega T}} - 1} \right)}}$ is maximized, i.e.
\begin{equation} \label{Peak value}
{\omega _{{\rm{peak}}}} = T\left[ {3 + W\left( { - 3{e^{ - 3}}} \right)} \right],
\end{equation}
where $W\left( x \right)$ is Lambert W-function which is the inverse function of $F\left( W \right) = W{e^W}$. ${{\tau _{{\rm{emission}}}}}$ is the characteristic time for an individual quantum emitted from black hole, which can be taken equal to the period of the quantum emitted, since the quanta emitted at peak value can only be temporarily localized within a few oscillation periods. Therefore $\tau _{{\text{emission}}}$ obeys the bound
\begin{equation} \label{characteristic time}
{\tau _{{\text{emission}}}} \geqslant \frac{{2\pi }}{{{\omega _{{\text{peak}}}}}},
\end{equation}
which leads to an upper bound for sparsity
\begin{equation} \label{etamax}
\eta  \leqslant \frac{{\omega _{{\text{peak}}}^2}}{{2\pi P}} \equiv {\eta _{\max }}.
\end{equation}
In summary, the sparsity of radiation in eq.(\ref{Sparsity parameter}) tells us that if  $\eta  \gg 1$, i.e. the average interval time $\tau _{{\rm{gap}}}$ between the emission of two quanta is much larger than the time $\tau _{{\text{emission}}}$  taken by an individual quantum emitted from the black hole, then naturally the Hawking radiation would be very sparse. On the contrary, if $\eta  \ll 1$  we will have almost continuous Hawking radiation.

In what follows we will study the dependence of sparsity on the conformal parameters $N$ and $L$ numerically.
Firstly, let us investigate the behavior of $\eta_{max}$ with different $N$ while keeping $L$ fixed. In TABLE.\ref{tab:table3}  the values of ${\eta _{\max }}$ for different values of $N$ with $L$ fixed are listed.
The TABLE.\ref{tab:table3} shows that given ${L \mathord{\left/
 {\vphantom {L {{r_{\text{H}}}}}} \right.
 \kern-\nulldelimiterspace} {{r_{\text{H}}}}}=0.7$ and $\ell=1$, as $N$ increasing, the values of ${\eta _{\max }}$ will decreases at first, and begin to increase at $N=3$.
In the next step, we go on to study the case with fixed $N=20$ (also $\ell=1$) and letting $L/r_H$ vary. The results are presented in TABLE.\ref{tab:table4}, from which it is clear that the situation in TABLE.\ref{tab:table4} is similar to that of TABLE.\ref{tab:table3}.  TABLE.\ref{tab:table4} shows that ${\eta _{\max }}$ takes its minimum value at ${L \mathord{\left/
 {\vphantom {L {{r_{\text{H}}}}}} \right.
 \kern-\nulldelimiterspace} {{r_{\text{H}}}}}=0.257$, and as ${L \mathord{\left/
 {\vphantom {L {{r_{\text{H}}}}}} \right.
 \kern-\nulldelimiterspace} {{r_{\text{H}}}}}$ getting larger than this values, ${\eta _{\max }}$  will increase as the increase of ${L \mathord{\left/
 {\vphantom {L {{r_{\text{H}}}}}} \right.
 \kern-\nulldelimiterspace} {{r_{\text{H}}}}}$.

In TABLE.\ref{tab:table3} and TABLE.\ref{tab:table4} the sparsity of $N=0$ or ${L \mathord{\left/
 {\vphantom {L {{r_{\text{H}}}}}} \right.
 \kern-\nulldelimiterspace} {{r_{\text{H}}}}}=0$ corresponds to Schwarzschild black hole, we show that when the $N$ and $L$ is relatively small (comparing with according values which minimize the sparsity) one can have nonrotating singularity-free black holes whose Hawking radiation are less sparse than that of Schwarzschild black hole with the same mass.

Also note that, as shown in TABLE.\ref{tab:table3} and TABLE.\ref{tab:table4}, the values for conformal parameters that minimize the sparsity of Hawking radiation are those values which maximize the greybody factor listed in TABLE.\ref{tab:table1} and TABLE.\ref{tab:table2}.
 This can be understood as follows,
as mentioned in section.\ref{HG} the dependence of the power spectrum (see eq.(\ref{energy emit3})) on $N$ and $L$ are totally coming from the greybody factor, hence, the values of $N$ and $L$ that maximize the greybody factor as presented in TABLE.\ref{tab:table1} and TABLE.\ref{tab:table2} will also maximize the power spectrum for Hawking radiation, which in turns minimize the upper bound ${\eta _{\max }}$ in eq.(\ref{etamax}).

The large values of ${\eta _{\max }}$ from TABLE.\ref{tab:table3} and TABLE.\ref{tab:table4} implies that the Hawking radiation from this black hole is extremely sparse. Moreover, the Hawking radiation becomes even more sparse for the increase of relatively large values of conformal parameters $N$ and $L$.

\begin{table}[h!]
  \begin{center}
 \caption{Numerical values of dimensionless parameter ${\eta _{\max }}$ for different values of $N$  with $\ell=1$ and ${L}/{{{r_{\rm{H}}}}}=0.7$.}
    \label{tab:table3}
    \begin{tabular}{c|c|c|c|c|c|c|c|c|c|c|c|c|c}
     \hline \hline
     $N$&$0$&$1$&$2$&$3$&$4$&$5$&$6$&$7$&$8$&$9$&$10$&$15$&$20$\\
      \hline
      ${\eta _{\max }}$&2506.84&2068.86&1823.34&1719.36&1736.27&1877.39&2170.99&2679.09&3517.53&4895.76&7193.68&95611.3&273617\\
      \hline\hline
    \end{tabular}
  \end{center}
\end{table}
\begin{table}[h!]
  \begin{center}
 \caption{Numerical values of dimensionless parameter ${\eta _{\max }}$ for different values of ${L}/{{{r_{\rm{H}}}}}$ with $\ell=1$ and $N=20$.}
    \label{tab:table4}
    \begin{tabular}{c|c|c|c|c|c|c|c|c|c|c|c|c|c}
     \hline \hline
     ${L \mathord{\left/
 {\vphantom {L {{r_{\rm{H}}}}}} \right.
 \kern-\nulldelimiterspace} {{r_{\rm{H}}}}}$&0&0.05&0.1&0.15&0.2&0.257&0.3&0.35&0.4&0.5&0.7&1&1.5\\
      \hline
      ${\eta _{\max }}$&2506.84&2441.15&2268.87&2051.76&1862.66&1768.19&1846.84&2225.15&3233.5&13293.1&$2.7{\text{4}} \times {\text{1}}{{\text{0}}^{\text{6}}}$&$2.4{\text{8}} \times {10^{11}}$&$1.2{\text{8}} \times {10^{20}}$\\
      \hline\hline
    \end{tabular}
  \end{center}
\end{table}

\section{Summary}\label{Summary}
We studied the sparsity of Hawking radiation from a nonrotating singularity-free black hole in conformal gravity proposed in \cite{Modesto:2016max,Bambi:2016wdn}. We mainly focused on analyzing the dependence of sparsity of Hawking radiation on the conformal parameters $N$ and $L$ (see eq.(\ref{S(r)})).

We first obtained a rigorous lower bound on the greybody factor (see eq.(\ref{Bound 1})) for the massless scalar field radiated from the black hole and analyzed the dependence of the greybody factor on the conformal parameters $N$ and $L$ via the bound. We found that there exist values for conformal parameters $N$ and $L$ that can maximize the lower bound on the greybody factor. For fixed $L$ (or $N$), the values of $N$ (or $L$)  maximizing this bound are given numerically (see TABLE.\ref{tab:table1} and TABLE.\ref{tab:table2}). Furthermore, we investigated the power spectrum of Hawking radiation (see eq.(\ref{energy emit3})) and showed that the dependence of power spectrum for Hawking radiation on conformal parameters $N$ and $L$ is the same as that of the greybody factor.

Finally, with the rigorous lower bound on greybody factor and power spectrum for Hawking radiation in hand, we studied the sparsity of Hawking radiation by the upper bound ${\eta _{\max }}$ (see eq.(\ref{etamax})) of dimensionless parameter $\eta$ (see eq.(\ref{Sparsity parameter})). Since the dependence of  ${\eta _{\max }}$ on conformal parameters $N$ and $L$ is completely coming from the greybody factor thus those values for conformal parameters $N$ and $L$ that maximize the lower bound on greybody factor in TABLE.\ref{tab:table1} and TABLE.\ref{tab:table2} can also minimize  ${\eta _{\max }}$ (sparsity). We discovered that when the conformal parameters $N$ and $L$ are relatively small comparing with the values that minimize the sparsity one can have nonrotating singularity-free black holes with Hawking radiation less sparse than that of Schwarzschild black hole with the same mass (see TABLE.\ref{tab:table3} and TABLE.\ref{tab:table4}). Also, we found that generally ${\eta _{\max }} \gg {\text{1}}$ as listed in TABLE.\ref{tab:table3} and TABLE.\ref{tab:table4}, which means the Hawking radiation is extremely sparse.  When the values of conformal parameters surpassing the values that minimize the sparsity of Hawking radiation, the further increase of conformal parameters will enhance ${\eta _{\max }}$ rapidly, which implies the Hawking radiation becomes even more sparse.

In future work, we plan to extend this study to the case of rotating singularity-free black hole which has also been proposed in \cite{Modesto:2016max,Bambi:2016wdn} since rotating black holes are more likely to exist in our universe.

\section*{Acknowledgement}
We would like to thank Jia-Rui Sun, Bin Wu, and Zhen-Ming Xu for useful discussions. This work was supported by the National Natural Science Foundation of China (No.~11675272).