# Fast unfolding of communities in large networks: 15 years later

Vincent Blondel1, Jean-Loup Guillaume2 and Renaud Lambiotte3

( 1UCLouvain, Belgium<EMAIL_ADDRESS>
2L3i, La Rochelle Université, France<EMAIL_ADDRESS>
3Turing Institute, University of Oxford, London, UK,
<EMAIL_ADDRESS>
###### Abstract

The Louvain method was proposed 15 years ago as a heuristic method for the
fast detection of communities in large networks. During this period, it has
emerged as one of the most popular methods for community detection, the task
of partitioning vertices of a network into dense groups, usually called
communities or clusters. Here, after a short introduction to the method, we
give an overview of the different generalizations and modifications that have
been proposed in the literature, and also survey the quality functions, beyond
modularity, for which it has been implemented.

Keywords: Networks, Community Detection, Optimization, Louvain method

##  1 Introduction

Networks provide a powerful language to model and analyse interacting systems,
in a broad range of disciplines [1]. Facebook and food-webs, the brain and the
job market are all systems where connectivity is critical, and where elements
– here users, species, neurons and workers – are represented by vertices,
while pairwise interactions – here friendship, predation, synchronization and
job mobility – are represented by edges. Over the last two decades, several
tools have been developed to extract information from the myriads of
interactions composing large networks. Among those, community detection has
played plays a central role. Community detection aims at uncovering the
community structure of a network, that is its organization into groups where
vertices in the same group are strongly connected with each other, and weakly
connected with vertices outside the group [2, 3]. Community detection is
usually considered as an unsupervised clustering problem where the purpose is
to find the best – according to a given definition – division into groups
together with the most appropriate number of groups. Community detection is a
rich and challenging problem, that has been considered from a variety of
perspectives and principles [4] and finds many applications, from network
visualization [5] to vertex classification [6]. In its most popular form, the
communities are assumed to form a partition of the network, so that each
vertex belongs to one single community. Partitioning methods have the
advantage that their outcome is easier to interpret than that of overlapping
methods, but also to involve less parameters to estimate and to have sound
mathematical foundations. Even within the family of partitioning methods,
several complementary approaches have been developed, leading to different
functions to measure the quality of a partition defined in terms of e.g. the
density of links inside communities [7], the number of links between
communities [8], the confinement of random walkers inside communities [9, 10]
or its likelihood of having generated the observed network [11].

A common challenge across methods is the need for efficient algorithms to
optimize their respective quality function, that is finding the partition of
the network with an optimal value. This problem is not new and finds its root
in graph partitioning, where different algorithms have been proposed since the
1970s. It is more recently, arguably since the seminal works of Newman [7]
that the research community started focusing on the related problem of
community detection. It rapidly gained popularity due to the increasing access
to large network data whose analysis required efficient clustering methods.
The most popular quality function is the so-called Newman-Girvan modularity,
and community detection is often considered - in practice - synonymous to
modularity optimization, despite its limitations and the availability of
alternative approaches. Until 2007, the most efficient methods for community
detection were either spectral, i.e. defining modules based on dominant
eigenvectors of matrices associated to the graph, or greedy, with the early
works of Clauset et al. [12] and Wakita and Tsurumi [13]. However, the former
are by construction limited to networks of the order of $10^{4}$ vertices,
orders of magnitude smaller than large-scale social and information networks,
whereas the proposed greedy methods did not perform sufficiently well in terms
of accuracy, due to the ruggedness of the optimization landscape.

On March 4 2008, we posted on the arXiv an article entitled “Fast unfolding of
communities in large networks ” proposing a greedy method for modularity
optimization, and made available an efficient C++ code. A workshop was
organized a few days later, on March 8, entitled “Detection and visualization
of communities in large complex networks”, where we invited the leading
researchers in the field to our research group in Louvain-la-Neuve. After a
desk rejection from Nature and a swift rejection from PNAS, the paper was
submitted to JSTAT, accepted on September 3 and published on October 9 [14].
15 years later, the paper has been cited just over 20,000 times on Google
Scholar, this is one new citation every 6 hours for more than 15 years. The
method had no name in the original article, and was for a time named BGLL
after the authors’ initials. The term ”Louvain method” appears for the first
time in an article published in September 2009 [15]. The Louvain method,
implemented in standard libraries for graph mining (see Table 1), is at the
core of the network science ecosystem, with applications in several fields of
science (see Figure 1). The purpose of this perspective paper is take a step
back and to consider the legacy of the paper, by giving an overview of the
research surrounding it, from improvements of the original method to its
applications to other quality functions.

Figure 1:  (Top left) Collaboration of researchers having cited the Louvain
method. All articles citing [14] have been downloaded from Scopus on the 6th
of April 2023. A bipartite network has been constructed relating author IDs,
as identified by Scopus, and papers. The collaboration network was then
constructed by first projecting the network according to [16], removing edges
with a weight smaller than $0.2$ and keeping the largest connected component.
The resulting network is composed of 2490 vertices. Running the Louvain method
identifies 39 communities, with a modularity of 0.926. The size of the
vertices is proportional to their weighted degree. (Top right) Originally
published in a Physics journal, the Louvain method has found applications in a
variety of scientific domains, as shown by the list of the top 20 journals
citing [14]. Further down in the list, one can find journals as diverse as
Poetics (9 articles), Journal of Transport Geography (8 articles) and Nature
Immunology (5 articles). (Bottom) Histogram of the total number of citations
per year according to Google Scholar (20,054 citations in total on July 6
2023).

Libraries
---
Language | Link | Objective function
Python | https://networkx.org/ | mutiresolution $Q$, directed
Python | https://networkit.github.io/ | mutiresolution $Q$, parallel
Python | https://github.com/vtraag/leidenalg | several functions, Leiden
Julia | https://github.com/afternone/CommunityDetection.jl | several functions
Javascript | https://graphology.github.io/ | mutiresolution $Q$
Matlab | https://github.com/GenLouvain/GenLouvain | several functions, multilayer
Java | https://github.com/jcnguyen/research-2015 | several functions
Multiple | https://igraph.org/ | mutiresolution $Q$

Softwares
---
Name | Objective | Link
Cytoscape (plugins) | Analysis and visualisation | https://cytoscape.org/
Gephi | Analysis and visualisation | https://gephi.org/
Tulip | Analysis and visualisation | https://tulip.labri.fr/site/
Pajek | Analysis and visualisation | http://mrvar.fdv.uni-lj.si/pajek/
SAS Viya | Analytics | https://www.sas.com/en_us/home.html
Neo4j | Database | https://neo4j.com/fr/
Memgraph | Database | https://memgraph.com/

Table 1: List of implementations of the Louvain method.

##  2 Basics on networks and modularity

A network is represented mathematically by a graph $\mathcal{G}(V,E)$, made of
a set of vertices $V$ of cardinality $n\equiv|V|$, and a set of edges
$E=\\{\\{i,j\\}\;|\;i,j\in V\\}$ of cardinality $m\equiv|E|$. Without loss of
generality we will identify the vertex set $V$ with the set
$\\{0,\ldots,n-1\\}$. In this introductory section, we will assume for the
sake of simplicity that the network is undirected, so that edges do not have a
direction, and unweighted, so that the graph only encodes the presence or
absence of an edge between two vertices. Note that the concepts described in
this section, and in particular expressions for modularity, have been
generalized to more general settings, as we discuss in detail in Section 4.3.
A useful way to encode and investigate the structure of a graph is its
adjacency matrix. The adjacency matrix $\bm{A}$ is an $n\times n$ matrix
encoding the presence or absence of an edge between any pair of vertices in
the graph. Its entries $A_{ij}$ are equal to 1 if $i$ is adjacent to $j$, and
0 otherwise. For undirected graphs we have $A_{ij}=A_{ji}$, i.e., the
adjacency matrix is symmetric ($\bm{A}=\bm{A}^{\top}$). Several properties of
a graph can be readily calculated from the adjacency matrix, including the
degree of each vertex $i$ as $k_{i}=\sum_{j}A_{ij}$, but above all, this
representation makes it possible to exploit concepts and tools from linear
algebra.

Since the seminal works of Watts and Strogatz [17], it has been recognised
that real-world networks exhibit a mix of random and non-random patterns.
Important types of structures, found in a wide range of empirical data,
include clustering, here understood in terms of the density of triangles, fat-
tailed degree distribution, power-law or not [18], and the presence of
communities in a network [19], that is the presence of dense clusters that are
weakly connected with one another. Community detection is an algorithmic
problem aiming at finding the community structure of a network. As such, it
can be seen as an unsupervised clustering problem, and some methods actually
consist in first embedding the vertices in a space, before applying classical
clustering methods, such as $k$-means [20]. In this section, we describe one
of the most popular methods for community detection, despite some limitations
described below, based on the notion of Newman-Girvan modularity [7],
modularity in short.

Modularity is a measure of the quality of the partition of network into
groups. Let us first consider a group of vertices defined by a set $\alpha$.
The intuitive idea behind modularity is to compare the number of links inside
$\alpha$ with an expectation of this number under a random null model. The
choice of null model is in principle open and should depend on the mechanisms
that constrain the formation of edges [21]. The most popular choice is, by
far, the soft configuration or Chung-Lu model [22], a random graph model
fixing, on expectation, the degrees of each vertex. For this choice, the
expected number of edges between two vertices $i$ and $j$ is given by
$\frac{k_{i}k_{j}}{2m}$, and the contribution of community $\alpha$ to
modularity is given by

$\sum_{i,j\in\alpha}\left(A_{ij}-\frac{k_{i}k_{j}}{2m}\right).$ (1)

Note that each term in the sum may either be positive or negative depending on
the presence or absence of a link between two vertices. Intuitively, a good,
dense community will be one with a high value for this expression. Note also
that the null model penalises more pairs of vertices with a higher degree,
which tends to favor configurations where high degree vertices are placed in
different communities [10].

The Newman-Girvan modularity does not focus on one single community, but on
the partition of the vertices into $C$ communities
$\mathcal{P}=\\{\alpha_{1},\alpha_{2},\ldots,\alpha_{C}\\}$. The quality of a
partition is obtained from the sum of Eq.(1) over the communities, leading to
the expression

$Q=\frac{1}{2m}\sum_{\alpha\in\mathcal{P}}\sum_{i,j\in\alpha}\left(A_{ij}-\frac{k_{i}k_{j}}{2m}\right),$
(2)

where the additional prefactor $1/(2m)$ ensures that modularity is in the
interval $[-1,1]$. Note that modularity can also be interpreted in terms of
the assortativity coefficient, a correlation measure of vertex attributes,
here the community labels, at the endpoints of edges in the network [23], in
terms of a covariance defined on the graph [24] and in terms of the clustered
covariance matrix from the random walk exploration of the graph [25].
Intuitively, modularity gives a high value to a partition if its communities
concentrate an unexpectedly large number of edges inside its communities. In
contrast with classical concepts such as the cut size, modularity allows to
compare partitions made of different numbers of communities, and the partition
maximising modularity is usually non-trivial, i.e. neither made of one large
community or of $n$ singletons. Modularity has become an essential element of
several community detection methods. Modularity optimization methods aim at
optimising modularity, that is finding the “best” partition of a network
having the highest value of modularity. As modularity optimization is NP-hard
[26], several heuristics have been proposed for modularity optimization.
Popular techniques include spectral methods [27] but large networks require
faster, greedy methods proceeding by agglomerating groups of vertices into
larger ones [28].

Before going further, note that the contribution of vertex $i$, assigned to a
certain modularity $\alpha$, to the modularity of a partition is given by

$Q_{i}=\frac{1}{2m}\left(\sum_{j\in\alpha}A_{ij}-\frac{k_{i}k_{\alpha}}{2m}\right),$
(3)

where $k_{\alpha}$ is the degree of community $\alpha$, defined as the sum of
the degrees of its vertices. This expression shows that the contribution of
vertex $i$ requires only local information, its number of neighbors inside the
community, its degree and the degrees inside the community. Locality is
essential to design efficient greedy algorithm for modularity optimization.

Modularity, as a quality function, suffers from some limitations, which are
inherited by any modularity optimization algorithm. Among those, let us note
its tendency to over-fit, finding communities even in the case of random
graphs [29], and to have, implicitely, a characteristic scale making it search
for communities with a size that might not be compatible with the underlying
structure [30]. Generalizations of modularity have been proposed to alleviate
these limitations, for instance by incorporating a resolution parameter
allowing to tune the characteristic size of the communities [31, 25, 32].

##  3 The Louvain method

The Louvain method is a greedy optimization heuristic originally proposed to
optimize modularity, but its flexibility allows to optimize other quality
functions as well, as described below. The approach consists in two phases
that are iteratively repeated, until a maximum of modularity is obtained. The
first phase consists in moving vertices from their community to that of their
neighbors until a local maximum is reached. The second phase then consists in
forming a new graph whose vertices are obtained by aggregating the vertices
attached to the same community in the previous phase. A combination of the two
phases is called a ”pass”. Starting from a network composed of $n$ vertices,
the number of vertices decreases at each pass, until no improvement of
modularity can be obtained.

First phase: vertex mover (VM). The first phase begins with an undirected
weighted graph having $C$ vertices to which an index between $0$ and $C-1$ has
been randomly assigned. $C$ is the number of vertices in the first pass, and
the number of communities found in previous pass otherwise. One starts by
placing each vertex into its own community, that is we start with a partition
made of $C$ communities. We then consider the first vertex, i.e. with index
$0$, and evaluate the change of modularity by removing $0$ from its community
and placing it in the community of each of its neighbors. Importantly for the
speed of the method, this step only requires local information – See Eq.(3) –,
and its number of operations is proportional to the degree of the vertex.
Vertex $0$ is then assigned to the community where the increase is maximum, if
this maximum increase is positive. Otherwise, vertex $0$ is left in its
original community – See Algorithm 1 (all algorithms are adapted from [33]).
This step is repeated sequentially to all the vertices based on their index.
When reaching vertex $C-1$, the process restarts at vertex $0$, and the
iteration continues until no vertex is moved during a complete iteration over
the $C$ vertices – See Algorithm 2. The phase ends in a local maximum, where
greedy move of the vertices between community do not allow for an increase of
the modularity. At this moment, a new partition of the vertices into
$C^{{}^{\prime}}$ communities has been produced. If $C^{{}^{\prime}}\neq C$,
one proceeds to the next phase, and update the value of the number of
communities. Else, the algorithm is finished and returns the partition as its
estimate for optimum of modularity.

Figure 2:  Each pass of the Louvain method is made of two phases. The first
phase (Vertex Mover) consists in sequentially moving vertices to neighboring
communities leading to a maximum increase of modularity. The second phase
aggregates vertices and constructs a meta-graph whose vertices are the
communities found after the first phase. The two phases are repeated until no
improvement of modularity is observed. In this illustration, the algorithm
converges after two passes, uncovering an optimal partition made of 2
communities, and a value of $Q=0.45$. In this example, it is easy to show that
aggregating further the two communities into one big community would result in
a decline of modularity, as $Q=0$ in that case. Figure adapted from [34].

Second phase: aggregation. The second phase consists in building a new
weighted graph whose vertices are the $C$ communities found during the first
phase. The weight of the link between two communities is given by the sum of
the weights of the links between the vertices of these two communities.
Moreover, links between the vertices of a same community lead to self-loops in
the aggregated network. This operation is performed to ensure that the
modularity of a partition in the aggregated graph has the same value as that
of the same partition on the original graph. Once the second phase is
finished, the first phase is reapplied on the aggregated network.

The output of the algorithm is the partition uncovered in the last pass, as it
is the best partition found in terms of modularity – See Figure 2 and
Algorithm 3. The method can also provide intermediary steps, that is the set
of partitions of each pass, which provides a hierarchical representation of
the network, each community at a pass being composed of a combination of
combinations at the previous pass. Since its inception, the method has
consistantly been shown to provide a good balance between accuracy and speed
of execution, and to have a complexity close to linear in benchmarks and
empirical networks. The speed of the method is ensured by the locality of each
operation, as vertices only move to neighboring communities. Its accuracy
arises from its flexibility, as vertices may be removed from their community
and re-assigned to others in the process, and its multi-scale nature, as the
community exploration is first done locally, and then over longer distance as
the vertices are aggregated. Combined with its simplicity, these strengths
have made the Louvain method a popular choice to optimize not only modularity,
but also other quality functions, including the Map Equation [9] and Markov
stability [32], as we discuss more in detail below. When applied on
modularity, as is usually the case, it is important to remember that the
Louvain method inherits the limitations of modularity, including its
resolution limit.

Require: $G=(V,E,w)$ a weighted graph

Require: $\mathcal{P}$ a partition of $V$ ($\mathcal{P}[i]=$ community of
vertex $i$)

Ensure: $\Delta Q$, the gain of modularity

begin

$c_{old}\leftarrow\mathcal{P}[i]$

REMOVE$(i,c_{old})$

$\mathcal{C}\leftarrow\\{\mathcal{P}[j]\mid(i,j)\in E\\}\cup\\{c_{old}\\}$

$c_{new}\leftarrow\arg\max_{c\in\mathcal{C}}\\{$GAIN$(i,c)\\}$

if  _$GAIN(i,c_{new})-LOSS(i,c_{old}) >0$_ then

INSERT$(i,c_{new})$

return $GAIN(i,c_{new})-LOSS(i,c_{old})$

else

INSERT$(i,c_{old})$

return 0

Algorithm 1 SingleVertexMover

Require: $G=(V,E,w)$ a weighted graph

Require: $\epsilon$ a stopping criterion

Ensure: a partition $\mathcal{P}$ of $V$

begin

INIT$(\mathcal{P})$

do

$increase\leftarrow 0$

forall  _vertices $i$_ do

$increase\leftarrow increase+$SingleVertexMover(i)

while  _$increase >\epsilon$_

return $\mathcal{P}$

Algorithm 2 VertexMoverIteration

Require: $G=(V,E,w)$ a weighted graph

Ensure: a partition $\mathcal{P}$ of $V$

begin

repeat

$\mathcal{P}\leftarrow$VertexMoverIteration$(G)$

$G\leftarrow$PartitionToGraph$(\mathcal{P},G)$

until  _no improvement is possible_

Algorithm 3 Louvain algorithm

##  4 Louvain++, a selection of enhancements

Because of its effectiveness and simplicity, the Louvain method has been
extensively studied, and numerous modifications have been proposed, whether to
improve results, further increase calculation speed or generalize its scope of
use. We present these modifications in four different directions, each
detailed in its respective subsection:

  * •

the modification of the algorithm’s core, including its initialization,
stopping criterion, vertex mover, etc;

  * •

the use of quality functions other than modularity, functions which may be
similar to modularity but with other parameters, or functions which may be
completely different;

  * •

generalizations of Louvain to deal with more than simple graphs. These may
include weighted, oriented, spatially constrained, temporal networks, etc;

  * •

and finally the use of parallelism to speed up computation time.

###  4.1 Modification of the algorithm

The original article and an extended version of it [34] (published in French
in 2011 and in English in 2013) already described several improvements that
have been proposed independently and studied in greater details by other
authors:

  * •

a partial optimization that consists in stopping the optimization when the
modularity gain is below a given threshold (both for iterations or
aggregation). This very simple optimization was already implemented in the
initial version of Louvain;

  * •

removing leaf vertices, which consists of deleting the leaves since they can
only be placed in the community of their single neighbor;

  * •

a VM that consists in not considering all vertices at each iteration;

  * •

returning to lower levels to move again vertices;

  * •

a simulated annealing-inspired VM allowing for non-optimal moves.

In the following, we describe modifications proposed in the different steps of
the algorithm, that is its initialization, the vertex mover step and the
aggregation step.

Initialization: Several types of optimization have targeted the initialization
phase. The first one is to reduce the size of the graph before computation by
removing a number of vertices that have little or no impact on the results.
Based on the observation that a leaf (any vertex with a single neighbor) will
inevitably be placed in its neighbor’s community, several authors propose
deleting these leaves [34, 35] or, more generally, all tree structures, i.e.
the initial leaves and all those that will appear when deleting leaves [36].
The resulting graph corresponds to the 2-core of the original graph and the
leaves are added back once the communities are found. A generalization has
been proposed in [37] where a k-core decomposition is performed and
communities are only computed on the k-core. All the removed vertices are then
added back and assigned to a community. The authors propose to use label
propagation or any modularity maximization algorithm for this task.

Other approaches have modified the initial partition by not starting with
communities containing a single vertex. The authors of [38] identify cliques
containing at least 3 vertices that serve as initial communities. However, the
way cliques are identified and the management of overlaps between these
cliques are not made explicit.

Vertex mover: In the Louvain method, the best choice is always preferred
during the vertex mover operation. Based on the principle of simulated
annealing, it was proposed to choose a non optimal move, even with a negative
gain, with a low probability rather than the optimal one [39].

Several articles propose to avoid considering all vertices at each iteration.
The simplest version is defined in [40]. If a vertex does not move for several
iterations (with a fixed parameter) then it will not be considered anymore. In
[36], the authors show that if two communities $A$ and $B$ have not changed in
the last iteration, then no vertex can move from $A$ to $B$ or from $B$ to $A$
in the current iteration. A similar idea is used in [41] where the author
identify vertices that are more likely to move in the next iteration and only
consider these vertices. They propose four situations, the simplest of which
considers that if a vertex $i$ has just moved from the community $A$ to the
community $B$, then all $i$’s neighbors who are not in $B$ are likely to join
$i$ at the next iteration. The other three situations are a little more
complex and are not used in practice. Limiting oneself to considering only a
few vertices at the next iteration could lead to iterations being stopped
prematurely and thus to a loss of quality, but in practice the impact is very
slight.

In a different fashion, authors of [35] use a concept of seed vertices which
are vertices more central than the average (for a given centrality measure
defined in the paper). Non-seed vertices are then moved in priority to the
neighboring community that contains a seed-vertex and maximizes the (positive)
gain of modularity. If there is no such neighboring community or if no gain is
positive, then the classical vertex mover is used.

In [42], instead of considering all neighbor communities and pick the best,
one neighbor is picked at random and only its community is considered – other
random selection procedure have also been considered, e.g. a community is
selected proportionally to its size. Even though random selection will not
always make the best choice, if a vertex has many neighbors in a given
community then this community is more likely to be picked. The author notes
that this operation accelerates convergence. This principle has been reused in
the Leiden method [43] with other elements.

In most implementations of Louvain, the order in which vertices are considered
is random. Although this can be seen as a problem, it does allow us to explore
more possible solutions. If the aim is to obtain an identical partition, any
fixed order can be used. In [44] the authors treat the vertices by decreasing
centrality. This choice appears to improve results slightly, but above all, it
speeds up the algorithm.

Other improvements have been proposed to simplify VM in a distributed
computing context. We describe them in Section 4.4

Aggregation and refinement: A key element of the Louvain method is that the
vertex mover operation is not greedy since any given vertex can be moved
several times until convergence is reached. As this step is not irreversible,
it is less likely to be penalized by poor initial choices. However, the
aggregation phases are irreversible, and the original method does not allow to
go back once the graph has been aggregated. For this reason, a number of
papers have proposed to challenge this principle and allow to move lower-level
vertices even after aggregation. This procedure is called partition
refinement.

In the simpler versions [45], a single refinement of the partition is
performed on the first level, i.e., the original graph where all vertices are
subject to the vertex mover procedure until convergence. In [46, 47], a multi-
level refinement is used where, starting from the top level of the partition,
the graph is unaggregated once and vertices are allowed to move. This process
is repeated moving down the hierarchy until they are back on the initial
graph. A different refinement is proposed in [48, 49] where the subgraph
induced by each community is partitionned using one level or the complete
Louvain method. All these parts of each communities are used to build the
aggregated network. A more complex scheme is used in the Leiden method [43]
where the vertex mover during refinement is limited and a move is accepted if
it increases the modularity even though it is not the best move. Finally, in
[50], if small (less than $4$ vertices) or weak (that score low using their
quality function) communities are found, these communities are broken into
individual vertices and a VM step is applied once.

The Louvain method and its evolutions do not impose a priori constraints on
the number or size of communities. If the number of communities is known to be
$k$, then [51] propose to stop the aggregation in Louvain at the highest level
that contains more than $k$ communities (i.e. the next aggregation would
result in fewer than $k$ communities). They then use a greedy method to merge
the communities that maximize the modularity gain to obtain exactly $k$
communities. Finally, VM iterations are done until convergence.

###  4.2 Modification of the quality function

Despite the shortcomings of modularity, several quality functions have been
inspired by it and used as a replacement in the Louvain method. For those, the
algorithm remains unchanged after adapting the expression for the modularity
change during the VM step. These quality functions mostly consider a positive
contribution for intra-community edges and a negative contribution for inter-
community edges or in the absence of edges.

  * •

Several variations of modularity are using a resolution parameter, and can be
shown to be equivalent, up to a trivial map, to the so-called Potts modularity
[31]

$Q_{\gamma}=\frac{1}{2m}\sum_{\alpha\in\mathcal{P}}\sum_{i,j\in\alpha}\left(A_{ij}-\gamma\frac{k_{i}k_{j}}{2m}\right),$
(4)

where the resolution parameter $\gamma$ allows to give more or less importance
to the negative term, and hence to modulate the size of the communities. This
expression has the same optimum as the so-called Markov stability, where the
quality of a partition is estimated by its tendency to capture random walkers
for long times inside communities [25, 32] and where time plays the role of
the resolution parameter. Most implementations of Louvain have been
generalized for one (or both) of these equivalent variations.

  * •

The balanced modularity [52, 33] uses a symmetric version that also considers
the absence of links

$\frac{1}{2m}\sum_{i,j\in
V}\left[\left(a_{ij}-\frac{k_{i}k_{j}}{2m}\right)x_{ij}-\left(\bar{a}_{ij}-\frac{(n-k_{i}).(n-k_{j})}{n^{2}-2m}\right)\bar{x}_{ij}\right].$

  * •

The authors of [39] use a function that considers both edges within
communities and edges between communities. The linear version of their
function is

$\frac{1}{2m}\sum_{i,j\in
V}\left[\left(a_{ij}-\frac{k_{i}k_{j}}{2m}\right)x_{ij}-\beta\left(a_{ij}-\frac{k_{i}k_{j}}{2m}\right)^{\alpha}\bar{x}_{ij}\right],$

where $\alpha$ and $\beta$ are resolution parameters that can be tuned to
adjust the size of communities.

  * •

In [53], the authors propose Jaccard Cosine Share Measure (JCSM). For two
vertices $i$ and $j$, the main term combines Jaccard (size of the intersection
of neighborhoods of $i$ and $j$ divided by the size of their union), Cosine
(normalized dot product of rows $i$ and $j$ of the adjacency matrix) and
classical modularity ($a_{ij}-k_{i}k_{j}/2m$).

  * •

Other functions mostly based on the number of connections within or between
communities have been proposed. For instance [50] propose a new function
defined as the sum for each community of

$F2(C_{i})=\frac{[d_{in}(C_{i})]^{2}}{[d_{in}(C_{i})+d_{out}(C_{i})]^{2}}.$

  * •

Based on other principles than modularity, the map equation introduced in [9]
uses the length of description of random walks in a network. The equation uses
a two levels description to encode both the modules and the vertices. vertices
from different modules can use the same code (the authors use the analogy of
street names, which may be similar in different cities). The optimization of
the map equation follows closely the steps of the Louvain method. Note that
the map equation has also been generalized to uncover a hierarchy of
communities [49].

More generally, works have explored the conditions that need to be satisfied
by a quality function so that it can be optimized by the Louvain method. In
that direction, the authors of [33] have used relational coding to describe
several quality function using the same formalism. A quality function $F$ is
said to be linear if $F(X,A)=\sum_{i,j\in V}\phi(a_{ij})x_{ij}+K$, where
$a_{ij}=1$ if $i$ and $j$ are connected ($0$ otherwise), $x_{ij}=1$ if $i$ and
$j$ are in the same community ($0$ otherwise), $\phi$ and $K$ are respectively
a function and a constant that only depends on the original data. For example,
modularity can be expressed as $Q=\sum_{i,j\in
V}\left(a_{ij}-k_{i}k_{j}/2m\right)x_{ij}$. If a quality function is linear
then the gain obtained using VM can be computed locally, henceforth without
loss of efficiency. The authors also study separable quality function, i.e.
that can be written as $F(X,A)=\sum_{i,j\in V}\phi(a_{ij})\psi(x_{ij})+K$,
where $\psi$ only depends on $X$. If the function $\psi$ only depends on the
community of interest then the corresponding quality function can also be used
efficiently in Louvain. The authors study five non-linear but separable
function and clarify why only some can be optimized efficiently. A different
formulation is given in [54]. The Louvain method can be used to optimize any
quality function of the form $trace\ H^{\top}[F-ab^{\top}]H$, where $H$ is the
partition indicator matrix, $F$ is a general matrix derived from the network,
and $a$ and $b$ are two $n$ dimensional vectors.

###  4.3 Louvain for non-simple graphs

Modularity has originally been defined for undirected, unweighted networks.
Since then, several generalizations have been proposed to allow for more
complex graph structures. As compared to the expression Eq.(3), the first term
often remains unchanged, to count the total number, or total weight, of edges
inside communities, but the second term is modified to take into account the
additional constrains induced by the graph structure. For instance, in the
case of directed networks, the null model is now chosen according to the
directed configuration model, to properly estimate the expected number of
edges between two vertices, given their respective in-degree and out-degree
[55]. The Louvain method can then be directly generalized with an adapted
modularity gain during VM [56]. Other examples include the optimization of
modularity for spatial networks where, again, the null model of the quality
function can be modified in order to account for spatial constraints [21].
Alternative approaches keep standard expressions for modularity but modify
instead the Louvain method to impose a spatial contiguity of the clusters in
the uncovered partition [57].

Community detection is also an active field of study for temporal networks
[58, 59]. Compared to its static counterpart, communities are now dynamical
objects that may encounter different types of change, such as growth or
splitting, and algorithms have the additional task to characterize these
events. Two different types of approach build on community optimization.
First, communities can be uncovered in different snapshots. In that case, some
constraints on the communities in adjacent snapshots may be incorporated in
order to ensure their continuity in time, such as in modularity maximization
under estrangement constraint [60]. After the communities have been found in
each snapshot, additional methods determine, with a certain statistical
accuracy, when and how a reorganization takes place [61]. A second family of
methods aims at finding a decomposition of the network in an extended network
in which snapshots are concatenated, e.g. by connecting the same vertex in
consecutive snapshots to favor continuity of communities across time. Here,
instead of searching communities in each time window, the community detection
is performed in one single operation, which can be formulated as a modularity
maximization, on a larger, multilayer network [62] and again optimized by the
Louvain method [63].

Modularity has also been generalized to quantify the quality of the partition
of signed networks. Here, the intuition is place together vertices such that
positive edges are concentrated within the modules and negative edges between
the modules. In its simplest setting, the signed modularity consists in
dividing the signed graph into two graphs, one defined by the positive
interactions and another one defined by the negative interactions, and to
search to maximize the difference of the modularities of the positive and
negative graph [64]. The Louvain method can be directly generalized for this
quality function, and the only step that needs to be considered with care is
VM. Consider a vertex $i$. One expects to find other vertices of the community
of $i$ in the neighborhood of $i$ only for the positive edges but its negative
neighbors should instead be placed in other communities. Ways to address this
complication are either to consider all the communities in the network in the
VM step - and not only the neighboring ones - , but this involves a massive
slow down of the algorithm, or to consider second neighbors, that is neighbors
of neighbors, in the case of negative edges [65].

Finally, some works have focused on attributed networks where nodes carry
labels. The aim is to create well-connected communities that are also
homogeneous with regard to labels, in order to meet both structural and
homophilic criteria. Three families of approaches are proposed in [66]: early
fusion where a new graph is built to encode both topology and attributes,
simultaneous fusion where the clustering is performed simultaneously on both
domains and late fusion where both clusterings are computed independently
before being joined. Among those, SAC2 [67] is an early fusion method where
each node retains only its k-nearest neighbors based on a weighted combination
of topology and similarity between attributes. Louvain is then used on this
construction. Simultaneous fusion generally involves a composite function of
modularity for the structural part and a similarity function for the
attributes. This is the case in [68] using purity for attributes. Relatedly,
SAC1 [67] uses different similarities for attributes, and [69] uses inertia.
In all cases Louvain is used to optimize these functions. A late fusion is
proposed in [70] where structural communities are found with Louvain and
attributed ones with k-means. The structural communities are only used if they
are sufficiently good, otherwise the attribute ones are used. More generally
any consensus-based solution can be used to merge several partitions.

###  4.4 Distributed/parallel/GPU solutions

Many efforts have been made to parallelize the Louvain method, the main
challenge being to perform the VM step on several CPUs simultaneously. There
are two main paradigms: shared-memory implementations, where processes share
the same memory and must therefore synchronize for writing, and distributed-
memory implementations, where processors communicate via messages to maintain
a consistent view. The former is limited by the number of physical processing
units that can be used, and the latter by the quantity of messages sent. In
most distributed implementations, if not all, the input network is partitioned
so that vertices and edges are distributed evenly among the processing units.
Evaluating the correctness of a parallel implementation can be complex.
Because of the impact of the order in which vertices are processed, one does
not expect to obtain the same results from one run to the next, or when
compared with a sequential run.

Shared-memory implementations: In shared-memory implementations several
vertices evaluate simultaneously and independently the best VM. They all base
their choice on a structure that may no longer be valid at the time of the
decision. The assignment of vertices to communities is protected and can only
happen in a sequential way so that no inconsistency can occur but a vertex $i$
may have been assigned between the time a vertex $j$ made its choice and the
time it validated it. Hence convergence might be slower or even not
attainable. Several approaches have been considered.

  * •

To keep the VM sequential, [71] propose an implementation by identifying the
loops in the code that can be executed in parallel. This includes
initialization of the algorithm, identifying all the neighboring communities
of a vertex and finally identifying the neighboring community with the best
gain. Community merging can also be partially parallelized.

  * •

To avoid conflicts that may results in a slower convergence one solution would
be to find vertices that do not conflict with each other and parallelize them.
While the extra cost is generally considered too high [72], some authors
propose to identify isolate sets of vertices, i.e. vertices that are not
neighbors of each other, and to process them in parallel [73]. The authors of
[72] propose to neglect this risk of non-convergence but limit the number of
iterations. They also compare asynchronous moves with synchronous moves, where
all parallely processed vertices are moved simultaneously.

  * •

The authors of [74] propose to adapt the number of processes dynamically
depending on the number of neighbors of a vertex and available resources.

  * •

In [75], the authors use techniques already presented to avoid considering all
vertices at every iteration. In the first iterations they use a pull method,
which is the classic method where neighboring communities are recalculated at
each iteration for every vertex. In the later iterations there are fewer
vertex moves, they switch to a push method where neighboring communities are
simply updated (pushed) with each vertex move. This optimization is combined
with classic code parallelization, using locks for updates of neighboring
communities.

Distributed memory implementations: In distributed memory implementations,
vertices are distributed among the various processes. Each process will
therefore store a given number of vertices. A number of problems arise that
require algorithmic adaptation of the Louvain method. First, vertices or edges
must be distributed evenly between processes to improve overall calculation
time. In complex networks, the presence of very large degree vertices (or
hubs) makes vertices/edges distribution more complex. Second, if the ends of
an edge are not on the same process, messages will have to be sent to
synchronize them and the total number of messages must be kept as low as
possible. Finally, since calculations are made independently, it is necessary
to be able to question local choices. Several works have looked at, and
addressed, these issues.

  * •

In [76], the authors use three levels of parallelism: (1) the network is
divided in subgraphs and edges between these subgraphs are ignored. Louvain is
applied on each subgraph. Then the individual results are joined,
reintroducing the missing edges, and Louvain is reapplied. (2) vertices of
each subnetwork are considered in batches and vertices of a given batch are
processed in parallel. Potential conflicts if two vertices of a batch are
neighbors are not considered but the update of vertex’s community are still
performed sequentially to avoid inconsistency. (3) Modularity gain for
neighboring communities is computed in parallel.

  * •

Processes can store knowledge of the neighbors of the vertices it owns even
though such vertices are stored on another device. Such vertices are called
”ghost vertices”. Similarly, a process can store information about ”ghost”
communities adjacent to the ones it owns. Ghost vertices and communities are
used in [77]. Initially, contiguous block of vertices are distributed trying
to balance the number of edges received by each process. At each iteration,
some messages are exchanged between processes owning ghost vertices to update
the information about community assignment. Graph coarsening is also performed
in a distributed way so as to determine new global ids of communities and to
exchange partial edge lists of neighbors communities so that each process can
have a complete edge list. A similar procedure adapted for GPUs is detailled
in [78]

  * •

The authors of [79] use two hash tables to store the graph and also propose to
identify the vertices that will provide the largest improvement of modularity
and to only consider such vertices to accelerate the VM step. They show that
the amount of such vertices decays exponentially fast so that only the first
few iterations are costly.

  * •

The initial distribution of vertices can have an impact on the workload but
also on the number of ghost vertices and therefore on the quantity of messages
to be transmitted. This is particularly the case for hubs so it is worth
choosing the initial partition carefully. In [76], the authors distribute
vertices evenly based on their id. This is studied in more details in [80]
where the authors show that the degree distribution of such a partition is
close to the degree distribution of the whole network, which ensures a fair
workload. Then they propose an enhancement that also considers the degree of
ghost vertices to improve workload fairness.

  * •

To account for hubs it is also possible to duplicate them. Each duplicate
(called a delegate) retains the neighbors belonging to its partition by
default, but this can be modified to ensure a better balance of edges between
partitions [81]. Then, delegates communicate with each other to ensure
synchronization. This general principal is used in [82] as a preprocessing
step for a distributed Louvain implementation. The authors also take delegates
into account when calculating the modularity gain.

  * •

The authors of [83] use PMETIS, a parallel implementation of METIS [84] to
obtain a partition of the network and distribute theses parts. Each one is
then processed independently, then aggregation is performed and the rest of
the algorithm is executed sequentially.

  * •

Independent calculations on each process can result in incorrect communities
if vertices are incorrectly assigned to a process. This problem is studied in
[85] with the identification of doubtful vertices. The relative commitment for
a vertex is calculated by computing the ratio between its internal degree and
the maximum internal degree in its community, and a vertex’s affinity for a
community is calculated by taking into account the relative commitment of the
vertex and its neighbors. Vertices with low commitment are considered doubtful
and can be moved to another device if they have many connections with vertices
on that device.

##  5 Perspectives

In this article, we have overviewed the different ways in which the Louvain
method has been used, generalised and improved over the last 15 years. As we
have shown, a wide variety of works have built on its original design,
exploring how to use it for other quality functions than modularity, and
adapting its different steps for faster, more accurate and/or distributed
solutions. A great deal is known about the strengths and weaknesses of the
methods as well as its theoretical foundations, yet it remains an active
source of inspiration as well as a standard method to cluster large-scale
networks in practice. Never did we imagine the impact of the method when it
was originally released and it is equally difficult for us to foresee now what
will become of it in the years to come. The method has accompanied the
different trends that have energed in network science, from multiplex networks
to embeddings and temporal networks. For this reason, it is likely that the
Louvain method will find applications and generalisations in the increasingly
popular field of higher-order networks [86, 87, 88, 89], but also in
continuous generalisations of networks like graphons [90] and in extensions of
graphs to non-real weights [91, 92]. But we are even more sure that we will be
surprised, and we look forward to sitting down in 15 years, discussing
together, discovering in awe the many ways the research community has modified
our original idea and, who knows, sharing it with you [93].

## Acknowledgments

Our special thanks go to Etienne Lefebvre, without whom the Louvain method
would never have existed. Only one publication as a master student, but what a
success! We also thank John Pougué Biyong for his help to analyse the Scopus
data and the many colleagues with whom we have collaborated on this topic over
the year, in particular Thomas Aynaud, Mauricio Barahona, Romain Campigotto,
Patricia Conde-Céspedes, Jean-Charles Delvenne, Mason Porter, Martin Rosvall
and Michael Schaub. Finally, many thanks to Sébastien Amoury and Shazia Babul
for carefully reading this manuscript.

The work of RL was supported by the EPSRC grants EP/V013068/1 and
EP/V03474X/1. The work of JLG was supported by the ANR MITIK project, French
National Research Agency (ANR), PRC AAPG2019.

## References

  * [1] Newman M 2018 Networks (Oxford university press)
  * [2] Fortunato S 2010 Physics reports 486 75–174
  * [3] Fortunato S and Newman M E 2022 Nature Physics 18 848–850
  * [4] Schaub M T, Delvenne J C, Rosvall M and Lambiotte R 2017 Applied network science 2 1–13
  * [5] Garza S E and Schaeffer S E 2019 Physica A: Statistical Mechanics and its Applications 534 122058
  * [6] Devooght R, Mantrach A, Kivimäki I, Bersini H, Jaimes A and Saerens M 2014 Random walks based modularity: application to semi-supervised learning Proceedings of the 23rd international conference on World wide web pp 213–224
  * [7] Newman M E and Girvan M 2004 Physical review E 69 026113
  * [8] Von Luxburg U 2007 Statistics and computing 17 395–416
  * [9] Rosvall M and Bergstrom C 2008 Proc. Natl. Acad. Sci. USA 105 1118–1123
  * [10] Lambiotte R, Delvenne J C and Barahona M 2014 IEEE Transactions on Network Science and Engineering 1 76–90
  * [11] Peixoto T P 2013 Physical review letters 110 148701
  * [12] Clauset A, Newman M E and Moore C 2004 Physical review E 70 066111
  * [13] Wakita K and Tsurumi T 2007 Finding community structure in mega-scale social networks Proceedings of the 16th international conference on World Wide Web pp 1275–1276
  * [14] Blondel V D, Guillaume J L, Lambiotte R and Lefebvre E 2008 Journal of statistical mechanics: theory and experiment 2008 P10008
  * [15] Lambiotte R and Panzarasa P 2009 Journal of Informetrics 3 180–190 ISSN 1751-1577 science of Science: Conceptualizations and Models of Science URL https://www.sciencedirect.com/science/article/pii/S1751157709000248
  * [16] Newman M E 2001 Physical review E 64 016132
  * [17] Watts D J and Strogatz S H 1998 nature 393 440–442
  * [18] Broido A D and Clauset A 2019 Nature communications 10 1017
  * [19] Simon H A 1977 Models of discovery: And other topics in the methods of science 245–261
  * [20] MacQueen J B 1967 Some methods for classification and analysis of multivariate observations Proc. of the fifth Berkeley Symposium on Mathematical Statistics and Probability vol 1 ed Cam L M L and Neyman J (University of California Press) pp 281–297
  * [21] Expert P, Evans T S, Blondel V D and Lambiotte R 2011 Proceedings of the National Academy of Sciences 108 7663–7668
  * [22] Chung F and Lu L 2002 Annals of combinatorics 6 125–145
  * [23] Peel L, Delvenne J C and Lambiotte R 2018 Proceedings of the National Academy of Sciences 115 4057–4062
  * [24] Devriendt K, Martin-Gutierrez S and Lambiotte R 2022 SIAM Review 64 343–359
  * [25] Delvenne J C, Yaliraki S N and Barahona M 2010 Proceedings of the national academy of sciences 107 12755–12760
  * [26] Brandes U, Delling D, Gaertler M, Gorke R, Hoefer M, Nikoloski Z and Wagner D 2007 IEEE transactions on knowledge and data engineering 20 172–188
  * [27] Newman M E J 2013 Physical Review E 88 042822
  * [28] Fortunato S 2010 Physics Reports 486 75–174 ISSN 0370-1573
  * [29] Guimera R, Sales-Pardo M and Amaral L A N 2004 Physical Review E 70 025101
  * [30] Ghasemian A, Hosseinmardi H and Clauset A 2019 IEEE Transactions on Knowledge and Data Engineering 32 1722–1735
  * [31] Reichardt J and Bornholdt S 2006 Physical review E 74 016110
  * [32] Lambiotte R, Delvenne J C and Barahona M 2008 Laplacian dynamics and multiscale modular structure in networks arXiv 0812.1770
  * [33] Campigotto R, Conde Céspedes P and Guillaume J L 2014 arXiv preprint arXiv:1406.2518
  * [34] Aynaud T, Blondel V D, Guillaume J L and Lambiotte R 2013 Multilevel Local Optimization of Modularity (John Wiley & Sons, Ltd) chap 13, pp 315–345 ISBN 9781118601181 (Preprint https://onlinelibrary.wiley.com/doi/pdf/10.1002/9781118601181.ch13) URL https://onlinelibrary.wiley.com/doi/abs/10.1002/9781118601181.ch13
  * [35] Lui D, Huang K, Zhang C, Wu D and Wu S 2021 Scientific Programming
  * [36] Zhang J, Fei J, Song X and Feng J 2021 Mathematical Problems in Engineering
  * [37] Peng C, Kolda T G and Pinar A 2014 Accelerating community detection by using k-core subgraphs (Preprint 1403.2226)
  * [38] Abbas E A and Nawaf H N 2020 Improving louvain algorithm by leveraging cliques for community detection 2020 International Conference on Computer Science and Software Engineering (CSASE) pp 244–248
  * [39] Zhou Z, Wang W and Wang L 2012 Community detection based on an improved modularity Pattern Recognition ed Liu C L, Zhang C and Wang L (Berlin, Heidelberg: Springer Berlin Heidelberg) pp 638–645 ISBN 978-3-642-33506-8
  * [40] Ryu S and Kim D 2016 Quick community detection of big graph data using modified louvain algorithm 2016 IEEE 18th International Conference on High Performance Computing and Communications; IEEE 14th International Conference on Smart City; IEEE 2nd International Conference on Data Science and Systems (HPCC/SmartCity/DSS) pp 1442–1445
  * [41] Ozaki N, Tezuka H and Inaba M 2016 International Journal of Computer and Electrical Engineering 8 207–218
  * [42] Traag V A 2015 Phys. Rev. E 92(3) 032801 URL https://link.aps.org/doi/10.1103/PhysRevE.92.032801
  * [43] Traag V A, Waltman L and van Eck N J 2019 Scientific Reports 9 5233 ISSN 2045-2322 URL https://doi.org/10.1038/s41598-019-41695-z
  * [44] Aldabobo A, Sharieh A and Riad J 2022 Journal of Theoretical and Applied Information Technology 100
  * [45] Du W and He X 2016 A common strategy to improve community detection performance based on the nodes’ property Bio-inspired Computing – Theories and Applications ed Gong M, Pan L, Song T and Zhang G (Singapore: Springer Singapore) pp 355–361 ISBN 978-981-10-3614-9
  * [46] Gach O and Hao J K 2014 Improving the louvain algorithm for community detection with modularity maximization Artificial Evolution: 11th International Conference, Evolution Artificielle, EA 2013, Bordeaux, France, October 21-23, 2013. Revised Selected Papers 11 (Springer) pp 145–156
  * [47] Rotta R and Noack A 2011 ACM J. Exp. Algorithmics 16 ISSN 1084-6654 URL https://doi.org/10.1145/1963190.1970376
  * [48] Waltman L and van Eck N J 2013 The European Physical Journal B 86
  * [49] Rosvall M and Bergstrom C T 2011 PLoS ONE 6 e18209 (Preprint 1010.0431)
  * [50] Yao B, Zhu J, Ma P, Gao K and Ren X 2023 Applied Sciences 13 ISSN 2076-3417 URL https://www.mdpi.com/2076-3417/13/6/4045
  * [51] Darst R K, Nussinov Z and Fortunato S 2014 Phys. Rev. E 89(3) 032809 URL https://link.aps.org/doi/10.1103/PhysRevE.89.032809
  * [52] Céspedes P C and Marcotorchino J F 2013 Comparing different modularization criteria using relational metric Geometric Science of Information ed Nielsen F and Barbaresco F (Berlin, Heidelberg: Springer Berlin Heidelberg) pp 180–187 ISBN 978-3-642-40020-9
  * [53] Chaudhary L and Singh B 2019 Community detection using an enhanced louvain method in complex networks Distributed Computing and Internet Technology ed Fahrnberger G, Gopinathan S and Parida L (Cham: Springer International Publishing) pp 243–250 ISBN 978-3-030-05366-6
  * [54] Schaub M, Delvenne J, Lambiotte R and Barahona M 2019 Physical Review E 99 062308–1–062308–18 URL http://dx.doi.org/10.1103/PhysRevE.99.062308
  * [55] Leicht E A and Newman M E 2008 Physical review letters 100 118703
  * [56] Li L, He X and Yan G 2018 Improved louvain method for directed networks Intelligent Information Processing IX: 10th IFIP TC 12 International Conference, IIP 2018, Nanning, China, October 19-22, 2018, Proceedings 10 (Springer) pp 192–203
  * [57] Wang C, Wang F and Onega T 2021 Transactions in GIS 25 1065–1081
  * [58] Cazabet R and Rossetti G 2019 Temporal network theory 181–197
  * [59] Masuda N and Lambiotte R 2016 A guide to temporal networks (World Scientific)
  * [60] Kawadia V and Sreenivasan S 2012 Scientific reports 2 794
  * [61] Greene D, Doyle D and Cunningham P 2010 Tracking the evolution of communities in dynamic social networks 2010 international conference on advances in social networks analysis and mining (IEEE) pp 176–183
  * [62] Mucha P J, Richardson T, Macon K, Porter M A and Onnela J P 2010 science 328 876–878
  * [63] Jutla I S, Jeub L G, Mucha P J et al. 2011 URL http://netwiki. amath. unc. edu/GenLouvain
  * [64] Traag V A and Bruggeman J 2009 Physical Review E 80 036115
  * [65] Pougué-Biyong J and Lambiotte R In preparation
  * [66] Chunaev P 2020 Computer Science Review 37 100286 ISSN 1574-0137 URL https://www.sciencedirect.com/science/article/pii/S1574013720303865
  * [67] Dang T and Viennet E 2012 Community detection based on structural and attribute similarities International Conference on digital society pp 7–12
  * [68] Citraro S and Rossetti G 2020 Applied Network Science 5 1–20 URL https://appliednetsci.springeropen.com/articles/10.1007/s41109-020-00302-1
  * [69] Combe D, Largeron C, Géry M and Egyed-Zsigmond E 2015 I-Louvain: An Attributed Graph Clustering Method Intelligent Data Analysis (Saint-Etienne, France: LaHC, University of Saint-Etienne, France) URL https://ujm.hal.science/ujm-01219447
  * [70] Elhadi H and Agam G 2013 Structure and attributes community detection: Comparative analysis of composite, ensemble and selection methods Proceedings of the 7th Workshop on Social Network Mining and Analysis SNAKDD ’13 (New York, NY, USA: Association for Computing Machinery) ISBN 9781450323307 URL https://doi.org/10.1145/2501025.2501034
  * [71] Bhowmick S and Srinivasan S 2013 A template for parallelizing the louvain method for modularity maximization (Springer Basel) pp 111–124 Modeling and Simulation in Science, Engineering and Technology funding Information: This work was completed with the support of the College of Information Science and Technology, University of Nebraska at Omaha (UNO), and the FIRE grant from the UNO Office of Research and Creative Activity. Publisher Copyright: © 2013, Springer Science+Business Media New York.
  * [72] Shi J, Dhulipala L, Eisenstat D, Łăcki J and Mirrokni V 2021 Proc. VLDB Endow. 14 2305–2313 ISSN 2150-8097 URL https://doi.org/10.14778/3476249.3476282
  * [73] Qie H, Li S, Dou Y, Xu J, Xiong Y and Gao Z 2022 Scientific Reports 12
  * [74] Fazlali M, Moradi E and Tabatabaee Malazi H 2017 Microprocessors and Microsystems 54 26–34 ISSN 0141-9331 URL https://www.sciencedirect.com/science/article/pii/S014193311630240X
  * [75] Tithi J J, Stasiak A, Aananthakrishnan S and Petrini F 2020 Prune the unnecessary: Parallel pull-push louvain algorithms with automatic edge pruning Proceedings of the 49th International Conference on Parallel Processing ICPP ’20 (New York, NY, USA: Association for Computing Machinery) ISBN 9781450388160 URL https://doi.org/10.1145/3404397.3404455
  * [76] Cheong C Y, Huynh H P, Lo D and Goh R S M 2013 Hierarchical parallel algorithm for modularity-based community detection using gpus Proceedings of the 19th International Conference on Parallel Processing Euro-Par’13 (Berlin, Heidelberg: Springer-Verlag) p 775–787 ISBN 9783642400469
  * [77] Ghosh S, Halappanavar M, Tumeo A, Kalyanaraman A, Lu H, Chavarrià-Miranda D, Khan A and Gebremedhin A 2018 Distributed louvain algorithm for graph community detection 2018 IEEE International Parallel and Distributed Processing Symposium (IPDPS) pp 885–895
  * [78] Gawande N, Ghosh S, Halappanavar M, Tumeo A and Kalyanaraman A 2022 Parallel Computing 111 102898 ISSN 0167-8191 URL https://www.sciencedirect.com/science/article/pii/S0167819122000060
  * [79] Que X, Checconi F, Petrini F and Gunnels J A 2015 Scalable community detection with the louvain algorithm 2015 IEEE International Parallel and Distributed Processing Symposium pp 28–37
  * [80] Zeng J and Yu H 2015 Parallel modularity-based community detection on large-scale graphs Proceedings of the 2015 IEEE International Conference on Cluster Computing CLUSTER ’15 (USA: IEEE Computer Society) p 1–10 ISBN 9781467365987 URL https://doi.org/10.1109/CLUSTER.2015.11
  * [81] Pearce R, Gokhale M and Amato N M 2014 Faster parallel traversal of scale free graphs at extreme scale with vertex delegates SC ’14: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis pp 549–559
  * [82] Zeng J and Yu H 2018 A scalable distributed louvain algorithm for large-scale graph community detection 2018 IEEE International Conference on Cluster Computing (CLUSTER) pp 268–278
  * [83] Wickramaarachchi C, Frîncu M, Small P E and Prasanna V K 2014 2014 IEEE High Performance Extreme Computing Conference (HPEC) 1–6 URL https://api.semanticscholar.org/CorpusID:11970878
  * [84] Karypis G and Kumar V 1998 SIAM Journal on Scientific Computing 20 359–392
  * [85] Bhowmick A, Vadhiyar S and PV V 2022 Concurrency and Computation: Practice and Experience 34 e6987 (Preprint https://onlinelibrary.wiley.com/doi/pdf/10.1002/cpe.6987) URL https://onlinelibrary.wiley.com/doi/abs/10.1002/cpe.6987
  * [86] Lambiotte R, Rosvall M and Scholtes I 2019 Nature physics 15 313–320
  * [87] Battiston F, Cencetti G, Iacopini I, Latora V, Lucas M, Patania A, Young J G and Petri G 2020 Physics Reports 874 1–92
  * [88] Bianconi G 2021 Higher-order networks (Cambridge University Press)
  * [89] Bick C, Gross E, Harrington H A and Schaub M T 2023 SIAM Review 65 686–731
  * [90] Caron F and Fox E B 2017 Journal of the Royal Statistical Society Series B: Statistical Methodology 79 1295–1366
  * [91] Böttcher L and Porter M A 2022 arXiv preprint arXiv:2212.06257
  * [92] Tian Y and Lambiotte R 2023 arXiv preprint arXiv:2307.01813
  * [93] Blondel V, Guillaume J L and Lambiotte R 2038 In progress


# Two-photon Hong-Ou-Mandel interference and quantum entanglement between the
frequency-converted idler photon and the signal photon

Jiaxuan Wang  1,2,* Alexei V. Sokolov1,2  and Girish S. Agarwal  1,2,3
1Department of Physics and Astronomy, Texas A&M University, College Station,
Texas 77843, USA
2Institute for Quantum Science and Engineering, Texas A&M University, College
Station, Texas 77843, USA
3Department of Biological and Agricultural Engineering, Texas A&M University,
College Station, Texas 77843, USA<EMAIL_ADDRESS>
††journal: opticajournal

Quantum frequency up-conversion is a cutting-edge technique that leverages the
interaction between photons and quantum systems to shift the frequency of
single photons from a lower frequency to a higher frequency. If the photon
before up-conversion was one of the entangled pair, then it is important to
understand how much entanglement is preserved after up-conversion. In this
study, we present a theoretical analysis of the transformation of the time-
dependent second-order quantum correlations in photon pairs and find the
preservation of such correlations under fairly general conditions. We also
analyze the two-photon Hong-Ou-Mandel interference between the frequency-
converted idler photon and the signal photon. The visibility of the two-photon
interference is sensitive to the magnitude of the frequency conversion, and it
improves when the frequency separation between two photons goes down.

##  1 Introduction

Optical frequency up-conversion has been the subject of much interest and
investigation since its discovery in 1967 by Midwinter and Warner [1]. The
driving factor behind the importance of up-conversion [2] has been the
possibility of detecting infrared radiation [3, 4, 5, 6, 7] by converting it
to the optical domain as the detection technology in the optical domain is
much better. Thus, techniques for improving the efficiency of the up-
conversion process have been developed [8, 9]. The improved efficiency is
important for many applications, for example, in imaging [10, 11]. Of
particular importance in the context of quantum technologies is the up-
conversion of a single photon. The quantum theory of the up-conversion process
has been developed [12, 13]. With the growing interest in quantum
entanglement, Kumar and coworkers [14, 15, 8] demonstrated preservation of the
quantum correlations between up-converted idler photons and signal photons
produced by an optical parametric amplifier. In particular, they showed the
survival of the non-classical intensity correlation between two photons at
1064 nm when one of these was converted to 532 nm. Since these early
experiments, the interest shifted to single photons, and several experiments
studied the quality of the up-converted single photons [16, 17] and verified
the preservation of quantum correlations during the process [18, 19, 20]. As
known, a versatile source of single photons is the spontaneous parametric
down-conversion process where one has entangled pairs, i.e., signal and idler
photons whose frequencies can be quite different, still exhibiting strong
second-order quantum correlations represented by $g^{(2)}(\tau)$. Direct
measurement of $g^{(2)}(\tau)$ is difficult as the correlation time could be
100 fs or less, and one uses the second-order Hong-Ou-Mandel interference [21,
22, 23, 24, 25, 26, 27, 28] for getting information at these time scales. In a
recent work [29, 30] the correlation time was about 80 ps [31], and the
frequency of the idler photon was up-converted by about 120 THz. Tyumenev et
al [29] demonstrated the preservation of such correlation time between the up-
converted idler photon and the signal photon. Given the considerable interest
in the quantum properties of the up-converted photons, we present a first
principle theoretical calculation of the non-classical time correlations of
the up-converted and down-converted photon and the signal photon. We find the
preservation of such correlations under fairly general conditions. We also
analyze the two-photon Hong-Ou-Mandel interference between the frequency-
converted idler photon and the signal photon. The visibility of the two-photon
interference is sensitive to the magnitude of the frequency conversion, and it
improves when the frequency separation between two photons goes down. Our
theoretical results on HOM interference are quite relevant to several
experiments, for example, [18, 19, 20] on the quality of frequency-converted
single photons. Our results also apply to cases when the frequency is down-
converted though we focus on up-conversion processes.

##  2 Second-order Correlation of the Up-shifted Idler and Original Signal
Biphotons

As shown in Fig. 1, we start with the entangled biphoton state [32] comprised
of a idler photon of frequency $\omega_{s}$ and an signal photon of frequency
$\omega_{i}$

$\displaystyle\ket{\Phi}=\int\int
d\omega_{s}d\omega_{i}f(\omega_{s},\omega_{i})a_{s}^{\dagger}(\omega_{s})a_{i}^{\dagger}(\omega_{i})\ket{0,0},$
(1)

where $a_{s(i)}^{\dagger}(\omega_{s(i)})$ is the creation operator for the
idler (signal) mode that satisfies the
$[a_{s(i)}(\omega),a^{\dagger}_{s(i)}(\omega^{\prime})]=\delta(\omega-\omega^{\prime})$,
and
$[a_{s(i)}(\omega),a_{i(s)}(\omega^{\prime})]=[a_{s(i)}(\omega),a^{\dagger}_{i(s)}(\omega^{\prime})]=0$.
Function $f(\omega_{s},\omega_{i})$ is the two-mode frequency distribution. We
consider the input photon statistics to be an entangled Gaussian distribution

$\displaystyle
f(\omega_{s},\omega_{i})=\frac{1}{\sqrt{2\pi\sigma_{P}\sigma_{-}}}e^{-(\omega_{s}+\omega_{i}-\omega_{p})^{2}/(16\sigma_{P}^{2})}e^{-(\omega_{s}-\omega_{i}-\triangle)^{2}/(4\sigma_{-}^{2})},$
(2)

where $\omega_{p}$ and $\sigma_{P}$ are the frequency and bandwidth of the
pump; $\triangle$ is the frequency difference between the central frequency
$\omega_{i0}$ and $\omega_{s0}$; $\sigma_{-}$ is the bandwidth of the photon
pairs. The distribution is normalized, i.e.,
$\int\int|f(\omega_{s},\omega_{i})|^{2}d\omega_{s}d\omega_{i}=1$.

Figure 1: Experimental setup for the detection of second-order correlation
$g^{(2)}(\tau)$ of the up-shifted idler and original signal biphotons.

As for the frequency up-conversion process, The idler photon travels through a
medium that up-shifts its frequency

$\displaystyle
b_{i}(\omega)=T(\omega,\Omega)a_{i}(\omega-\Omega)+\eta(\omega,\Omega)a_{v}(\omega),$
(3)

where $T(\omega,\Omega)$ is the up-conversion rate; $a_{v}(\omega)$ is the
vacuum noise creation operator and
$\eta(\omega,\Omega)=\sqrt{1-|T(\omega,\Omega)|^{2}}$. The operators
$b_{i}(\omega)$ and $b^{\dagger}_{i}(\omega^{\prime})$ satisfy the commutation
relations
$[b_{i}(\omega),b^{\dagger}_{i}(\omega^{\prime})]=\delta(\omega-\omega^{\prime})$.
This formula has the same format as a transmission output-input relationship
in an ordinary sample, while now the only difference is that the output
frequency is up-converted, and the conversion rate should be dependent on both
the original frequency and the frequency shift. Many experiments on frequency
conversion use modulation at frequency $\Omega$. Thus the input photons at
$\omega$ would be converted to $\omega+\Omega$ and $\omega-\Omega$. The
modulation was in the GHz domain for some of the experiments, as noted in the
introduction. The signal photon frequency remains the same. A phase shift is
introduced due to the different paths taken by the photon pair

$\displaystyle b_{s}(\omega)=a_{s}(\omega)e^{-i\omega\tau_{0}}.$ (4)

The second order correlation $g^{(2)}(t,t+\tau)$ is defined as

$\displaystyle
g^{(2)}(t,t+\tau)=\expectationvalue{b^{\dagger}_{i}(t)b^{\dagger}_{s}(t+\tau)b_{s}(t+\tau)b_{i}(t)}.$
(5)

Considering Eq. (1), (3), (4), we obtain

$\displaystyle b_{s}(t+\tau)b_{i}(t)\ket{\Phi}$
$\displaystyle=\frac{1}{2\pi}\int\int
d\omega_{1}d\omega_{2}T(\omega_{1},\Omega)a_{s}(\omega_{2})a_{i}(\omega_{1}-\Omega)e^{i\omega_{1}t}e^{i\omega_{2}(t+\tau-\tau_{0})}$
$\displaystyle\times\int\int
d\omega_{s}d\omega_{i}f(\omega_{s},\omega_{i})a_{s}^{\dagger}(\omega_{s})a_{i}^{\dagger}(\omega_{i})\ket{0,0},$
(6)

which gives the result for Eq. (5) as

$\displaystyle g^{(2)}(t,t+\tau)$ $\displaystyle=\frac{1}{4\pi^{2}}|\int\int
d\omega_{1}d\omega_{2}T(\omega_{1},\Omega)\times
e^{i\omega_{1}t}e^{i\omega_{2}(t+\tau-\tau_{0})}f(\omega_{1}-\Omega,\omega_{2})|^{2}.$
(7)

Note that the vacuum noise term does not contribute to the normal-ordered
correlation. Assuming that the phase-matching is satisfied over the entire
frequency width of the idler photon, then the up-conversion rate is a constant
$T(\omega,\Omega)=T$ in the range of interest. Applying the parameter
transformation $\Omega_{1}=\omega_{1}-\Omega$, Eq. (7) can be simplified to

$\displaystyle g^{(2)}(t,t+\tau)$
$\displaystyle=\frac{T^{2}}{4\pi^{2}}|\int\int
d\Omega_{1}d\omega_{2}e^{i\Omega_{1}t}e^{i\omega_{2}(t+\tau-\tau_{0})}f(\Omega_{1},\omega_{2})|^{2}.$
(8)

Considering a photon pair following Gaussian distribution as in Eq. (2), we
obtain

$\displaystyle g^{(2)}(t,t+\tau)$
$\displaystyle=\frac{T^{2}}{(2\pi)^{3}\sigma_{P}\sigma_{-}}|\int\int
d\Omega_{1}d\omega_{2}e^{i(\omega_{2}+\Omega_{1})\frac{(t+\tau^{\prime}-\tau_{0})+t}{2}}$
$\displaystyle\times
e^{i(\omega_{2}-\Omega_{1})\frac{(t+\tau-\tau_{0})-t}{2}}e^{-\frac{(\omega_{2}+\Omega_{1}-\omega_{p})^{2}}{16\sigma_{P}^{2}}}e^{-\frac{(\omega_{2}-\Omega_{1}-\triangle)^{2}}{4\sigma_{-}^{2}}}|^{2}.$
(9)

By orthogonally transforming parameters to
$u=\frac{\omega_{2}+\Omega_{1}}{\sqrt{2}}$ and
$v=\frac{\omega_{2}-\Omega_{1}}{\sqrt{2}}$, the double integral reduces to a
product of two single-parameter integrals.

$\displaystyle g^{(2)}(t,t+\tau)$
$\displaystyle=\frac{T^{2}}{(2\pi)^{3}\sigma_{P}\sigma_{-}}|\int
due^{iu\frac{2t+\tau-\tau_{0}}{\sqrt{2}}}e^{-(u-\frac{\omega_{p}}{\sqrt{2}})^{2}/(8\sigma_{P}^{2})}\int
dve^{iv\frac{\tau-\tau_{0}}{\sqrt{2}}}e^{-(v-\frac{\triangle}{\sqrt{2}})^{2}/(2\sigma_{-}^{2})}|^{2}.$

Notably, the two integrals present in the expression are the Fourier
transforms of Gaussian functions. The shifts $\frac{\omega_{p}}{\sqrt{2}}$ and
$\frac{\triangle}{\sqrt{2}}$ thus don’t contribute to the modulus square. We
then obtain

$\displaystyle
g^{(2)}(t,t+\tau)=\frac{2\sigma_{P}\sigma_{-}T^{2}}{\pi}e^{-2\sigma_{P}^{2}(2t+\tau-\tau_{0})^{2}}e^{-\frac{1}{2}\sigma_{-}^{2}(\tau-\tau_{0})^{2}}.$
(11)

Note that $\omega_{p}$ and $\triangle$ do not explicitly contribute to the
formula. This observation is due to the fact that the photon pair in question
follows a standard bivariate normal distribution $f_{XY}(x,y)$. The
correlation between the two variables, denoted by $\rho$, is described by the
given formula [33, 34]

$\displaystyle
f_{XY}(x,y)=\frac{1}{2\pi\sqrt{1-\rho^{2}}}e^{-\frac{x^{2}-2\rho
xy+y^{2}}{2(1-\rho^{2})}}.$ (12)

In our case the correlation
$\frac{4\sigma_{p}^{2}-\sigma_{-}^{2}}{4\sigma_{p}^{2}+\sigma_{-}^{2}}$ for
the photon statistics $|f(\omega_{s},\omega_{i})|^{2}$ is a function of
$\sigma_{-}$ and $\sigma_{p}$ only.

Figure 2: Normalized second-order correlation $g^{(2)}(\tau)$ in Eq. (14) of
the up-shifted idler and original signal photon pair (yellow), compared to
$g^{(2)}(\tau)$ of the original photon pair (blue), for $\sigma_{-}=2\pi\times
1\>THz$, $\sigma_{p}=\sigma_{-}/10$, $\tau_{0}=0.2\>ps$ and $T_{R}=100\>ps$.

The two-time correlation in Eq. (11) is to be averaged over the resolution
time $T_{R}$ of the detector.

$\displaystyle g^{(2)}(\tau)=\int_{-T_{R}/2}^{+T_{R}/2}dtg_{2}(t,t+\tau),$
(13)

which gives

$\displaystyle
g^{(2)}(\tau)=\frac{1}{2\sqrt{2\pi}}\sigma_{-}T^{2}e^{-\frac{1}{2}\sigma_{-}^{2}(\tau-\tau_{0})^{2}}[erf(\sqrt{2}\sigma_{P}(\tau-\tau_{0}+T_{R}))-erf(\sqrt{2}\sigma_{P}(\tau-\tau_{0}-T_{R}))],$

where $erf(x)$ stands for the error function. We consider two scenarios. In
the first scenario, when the large-detection time condition is satisfied,
i.e., $\sigma_{P}T_{R}\rightarrow\infty$ and $T_{R}\gg\tau-\tau_{0}$, say with
$T_{R}\sim 10\>ps$ or $100\>ps$, we obtain

$\displaystyle
g^{(2)}(\tau)=\frac{\sigma_{-}T^{2}}{\sqrt{2\pi}}e^{-\frac{1}{2}\sigma_{-}^{2}(\tau-\tau_{0})^{2}}.$
(15)

The FWHM of $g^{(2)}(\tau)$ between the up-shifted idler and original signal
photon pair is the same as that between the original photon pairs. This is
illustrated in Fig. 2, where the normalized $g^{(2)}(\tau)$ (by dropping the
factor $T^{2}$) is plotted for the up-shifted idler and original signal photon
pair, and demonstrates the presence of only a peak shift arising from the
difference in path length. Note that Fig. 2 is obtained by using the exact
result in Eq. (14). In the second scenario, $\sigma_{P}T_{R}\rightarrow 0$,
and $\sigma_{P}(\tau-\tau_{0})\rightarrow 0$ we obtain

$\displaystyle
g^{(2)}(\tau)=\frac{T_{R}\sigma_{P}\sigma_{-}T^{2}}{\sqrt{\pi}}e^{-\frac{1}{2}\sigma_{-}^{2}(\tau-\tau_{0})^{2}}.$
(16)

In both limits, the width of the $g^{(2)}(\tau)$ function is unchanged except
for the phase difference imparted by the up-conversion medium from the
correlation when there is no up-conversion (and no $\tau_{0}$).

$\displaystyle
g^{(2)}_{0}(\tau)\propto\frac{\sigma_{-}}{\sqrt{2\pi}}e^{-\frac{1}{2}\sigma_{-}^{2}(\tau)^{2}}.$
(17)

Next, we examine the effect of phase-matching on the two-photon correlation.
For this purpose, we now take the factor $T(\omega,\Omega)$ as

$\displaystyle
T(\omega,\Omega)=Te^{-\frac{(\omega-\omega_{i0})^{2}}{2\beta^{2}}},$ (18)

where the conversion rate has a maximum at $\omega_{i0}$ to achieve maximum
conversion. It decreases significantly outside the range indicated by $\beta$.
We update Eq. (11) and (14) to

$\displaystyle
g^{(2)}(t,t+\tau)=\frac{4T^{2}\beta^{2}\sigma_{P}\sigma_{-}}{\pi(2\beta^{2}+4\sigma_{P}^{2}+\sigma_{-}^{2})}e^{-\frac{\Omega^{2}+16\sigma_{P}^{2}\sigma_{-}^{2}(t+\tau-\tau_{0})^{2}+8\beta^{2}\sigma_{P}^{2}(2t+\tau-\tau_{0})^{2}+2\beta^{2}\sigma_{-}^{2}(\tau-\tau_{0})^{2}}{2(2\beta^{2}+4\sigma_{P}^{2}+\sigma_{-}^{2})}},$
(19)

and

$\displaystyle
g^{(2)}(\tau)=\frac{T^{2}\beta^{2}\sigma_{-}}{\sqrt{2\pi(2\beta^{2}+4\sigma_{P}^{2}+\sigma_{-}^{2})(\sigma_{-}^{2}+2\beta^{2})}}e^{-\frac{\Omega^{2}}{2(2\beta^{2}+4\sigma_{P}^{2}+\sigma_{-}^{2})}}e^{-\frac{4(\tau-\tau_{0})^{2}\beta^{2}\sigma_{-}^{2}\sigma_{P}^{2}}{(2\beta^{2}+4\sigma_{P}^{2}+\sigma_{-}^{2})(\sigma_{-}^{2}+2\beta^{2})}}e^{-\frac{\beta^{2}\sigma_{-}^{2}(\tau-\tau_{0})^{2}}{(2\beta^{2}+4\sigma_{P}^{2}+\sigma_{-}^{2})}}$
$\displaystyle\times\\{erf[\sqrt{\frac{2(\sigma_{-}^{2}+2\beta^{2})}{2\beta^{2}+4\sigma_{P}^{2}+\sigma_{-}^{2}}}(2\tau-2\tau_{0}+T_{R})]-erf[\sqrt{\frac{2(\sigma_{-}^{2}+2\beta^{2})}{2\beta^{2}+4\sigma_{P}^{2}+\sigma_{-}^{2}}}(2\tau-2\tau_{0}-T_{R})]\\}.$
(20)

The expressions (19) and (20) reduce to equations (11) and (14) for
$\beta^{2}\rightarrow\infty$. If the condition $T_{R}\gg(\tau-\tau_{0})$ is
satisfied, then

$\displaystyle g^{(2)}(\tau)\propto
e^{-\frac{\Omega^{2}}{2(2\beta^{2}+4\sigma_{P}^{2}+\sigma_{-}^{2})}}e^{-\frac{\beta^{2}\sigma_{-}^{2}}{(\sigma_{-}^{2}+2\beta^{2})}(\tau-\tau_{0})^{2}},$
(21)

showing explicit dependence on the parameter $\beta$ and the frequency change
$\Omega$.

##  3 Hong–Ou–Mandel Interference of the Upshifted Signal and Original Idler
Biphotons

Figure 3: Experimental setup for the Hong–Ou–Mandel (HOM) measurement of the
upshifted idler and original signal biphotons, where $\tau_{T}$ is a tunable
phase.

The Hong–Ou–Mandel (HOM) interferometer [21, 35, 36], a widely used tool in
measuring biphoton joint frequency distribution, especially for time scales of
10’s of fs [37, 38, 39, 40]. We, therefore, study the two-photon interference
between the up-converted idler and the signal photon. By introducing a tunable
phase $\tau_{T}$ in the path of the signal photon, as depicted in Fig. 3, it
is possible to measure both the frequency shift $\Omega$ and phase shift
$\tau_{0}$ in a single setup. For the signal photon,

Figure 4: (a) Normalized $R_{c}(\tau_{T})$ function of the up-shifted idler
and original signal photons (yellow), compared to the $R_{c}(\tau_{T})$ of the
original photon pair (blue), for $\sigma_{-}=2\pi\times 1\>THz$,
$\tau_{0}=0.2\>ps$, $\triangle=2\pi\times 2\>THz$ and
$\triangle-\Omega=2\pi\times 0.05\>THz$. (b) HOM dip visibility as a function
of the shifted frequency $\Omega$ when $\sigma_{-}=2\pi\times 1\>THz$ and
$\triangle=2\pi\times 2\>THz$.

$\displaystyle b_{s}(\omega)=a_{s}(\omega)e^{-i\omega(\tau_{T}+\tau_{0})},$
(22)

while the idler photon travels through the same up-conversion medium. After
frequency conversion, the two photons are recombined at a 50:50 beam-splitter.
The output of the beam-splitter is expressed as

$\displaystyle c(\omega)=\frac{b_{i}(\omega)+ib_{s}(\omega)}{\sqrt{2}},$
$\displaystyle d(\omega)=\frac{ib_{i}(\omega)+b_{s}(\omega)}{\sqrt{2}}.$ (23)

As a result, the averaged coincidence rate of the output fields

$\displaystyle R_{c}(\tau_{T})=\frac{1}{4\pi^{2}}\int dt\int
d\tau\expectationvalue{d^{\dagger}(t)c^{\dagger}(t+\tau)c(t+\tau)d(t)}$ (24)

is measured, where we assume the detection time is much larger than the pump
correlation time $1/\sigma_{p}$ and the entanglement time $1/\sigma_{-}$.
Changing the modes to the frequency domain $d(t)=\frac{1}{\sqrt{2\pi}}\int
d\omega_{1}d(\omega_{1})e^{-i\omega_{1}t}$ and
$c(t+\tau)=\frac{1}{\sqrt{2\pi}}\int
d\omega_{2}c(\omega_{2})e^{-i\omega_{2}(t+\tau)}$, we obtain

$\displaystyle R_{c}(\tau_{T})=\frac{1}{8\pi^{2}}$ $\displaystyle[\int\int
d\omega_{1}d\omega_{2}f^{*}(\omega_{1}-\Omega,\omega_{2})f(\omega_{1}-\Omega,\omega_{2})T^{*}(\omega_{1},\Omega)T(\omega_{1},\Omega)$
$\displaystyle-\int\int$ $\displaystyle
d\omega_{1}d\omega_{2}f^{*}(\omega_{1}-\Omega,\omega_{2})f(\omega_{2}-\Omega,\omega_{1})T^{*}(\omega_{1},\Omega)T(\omega_{2},\Omega)e^{-i(\omega_{1}-\omega_{2})(\tau_{T}+\tau_{0})}].$

For a system with the same constant conversion rate assumption and photon
distribution as in Eq. (2), we obtain

$\displaystyle
R_{c}(\tau_{T})=\frac{T^{2}}{8\pi^{2}}(1-e^{-\frac{(\Omega-\triangle)^{2}}{2\sigma_{-}^{2}}}e^{-\frac{1}{2}\sigma_{-}^{2}(\tau_{T}+\tau_{0})^{2}}).$
(26)

The FWHM of $R_{c}(\tau_{T})$ is the same as the simple $g^{(2)}$ function as
in Eq. (15), (16). The position of the peak in the HOM dip provides
information on the phase difference, while the visibility
$e^{-\frac{(\Omega-\triangle)^{2}}{2\sigma_{-}^{2}}}$ of the dip reveals the
frequency shift $\Omega$. As shown in Fig. 4, the visibility of the HOM dip
increases if the frequency difference between the signal and idler photons
decreases due to the frequency conversion. Note that Eq. (26) shows that the
visibility of the HOM dip remains sufficient as long as the parameter
$f=\frac{\Omega-\triangle}{\sigma_{-}}\lesssim 1$, which is the case for
reported experiments [18, 19, 20].

Figure 5: (a) HOM dip visibility, and (b) FWHM for both $R_{c}(\tau_{T})$ and
$g^{{2}}(\tau)$, as a function of the phase-matching factor $\beta$ when
$\sigma_{-}=2\pi\times 1\>THz$, $\sigma_{p}=\sigma_{-}/10$,
$\triangle=2\pi\times 2\>THz$ and $\triangle-\Omega=2\pi\times 0.05\>THz$. The
dashed lines indicate the limits when $\beta\gg\sigma_{-}$.

Next, we examine the effects of imperfect phase-matching in up-conversion
using Eq. (18) and Eq. (25); we obtain

$\displaystyle R_{c}(\tau_{T})$
$\displaystyle=\frac{\sqrt{2}T^{2}\beta}{8\pi^{2}}[\frac{1}{\sqrt{4\sigma_{P}^{2}+\sigma_{-}^{2}+2\beta^{2}}}e^{-\frac{2\Omega^{2}}{4\sigma_{P}^{2}+\sigma_{-}^{2}+2\beta^{2}}}$
$\displaystyle-\frac{\beta}{\sqrt{(\beta^{2}+2\sigma_{P}^{2})(2\beta^{2}+\sigma_{-}^{2})}}e^{-\frac{(\Omega+\triangle)^{2}}{4(\beta^{2}+2\sigma_{P}^{2})}}e^{-\frac{(\Omega-\triangle)^{2}}{2\sigma_{-}^{2}}}e^{-\frac{\sigma_{-}^{2}\beta^{2}(\tau_{T}+\tau_{0})^{2}}{2\beta^{2}+\sigma_{-}^{2}}}],$
(27)

which reduces to Eq. (26) for $\beta^{2}\rightarrow\infty$. Note that Eq. (27)
shares the same FWHM as Eq. (21). From Eq. (27), we obtain the visibility of
the HOM dip

$\displaystyle
V=\frac{\beta\sqrt{4\sigma_{P}^{2}+\sigma_{-}^{2}+2\beta^{2}}}{\sqrt{(\beta^{2}+2\sigma_{P}^{2})(2\beta^{2}+\sigma_{-}^{2})}}e^{-\frac{(\Omega+\triangle)^{2}}{4(\beta^{2}+2\sigma_{P}^{2})}}e^{-\frac{(\Omega-\triangle)^{2}}{2\sigma_{-}^{2}}}e^{\frac{2\Omega^{2}}{4\sigma_{P}^{2}+\sigma_{-}^{2}+2\beta^{2}}}.$
(28)

The visibility is maximized when $\Omega=\triangle$ as

$\displaystyle
V=\frac{\beta\sqrt{4\sigma_{P}^{2}+\sigma_{-}^{2}+2\beta^{2}}}{\sqrt{(\beta^{2}+2\sigma_{P}^{2})(2\beta^{2}+\sigma_{-}^{2})}}e^{-\frac{(8\sigma_{P}^{2}+\sigma_{-}^{2}+4\beta^{2})\Omega^{2}}{(4\sigma_{P}^{2}+\sigma_{-}^{2}+2\beta^{2})(\beta^{2}+2\sigma_{P}^{2})}}.$
(29)

When $\beta\gg\sigma_{-}$ and $\sigma_{P}$, the visibility $V\rightarrow 1$.
However, a very satisfactory visibility of 0.96 can be obtained when
$\beta\sim 2\sigma_{-}$ for a system as shown in Fig. 5(a). The FWHM depends
both on $\sigma_{-}$ and $\beta$. As shown in Fig. 5(b), when $\beta\sim
2\sigma_{-}$ the peak is only expanded by 50% compared to the value when
$\beta\gg\sigma_{-}$.

##  4 Conclusions

The results of our theoretical study provide a comprehensive insight into the
effect of frequency up-conversion on the time-dependent quantum correlation of
a photon pair. By deriving the second-order correlation function, we
demonstrate that the full width at half maximum (FWHM) remains unchanged,
while the peak height and position shift after the up-conversion procedure.
This finding agrees with the recent experiment [29], where the correlation
time was within reach of the detector resolution time. This is true if the
phase matching is satisfied over the spectral width of the idler photon. To
gain a deeper understanding of the impact of the frequency up-conversion on
the quantum correlation, we examined the two-photon Hong-Ou-Mandel
interferometry. We showed how the visibility of the two-photon interference
was sensitive to the magnitude of the frequency change in the up-conversion
process, the bandwidth of the signal photon, and the phase matching factor
$\beta$. These results are fairly general and are expected to apply to other
types of frequency conversion.

## Funding

JW and GSA thank the support of Air Force Office of Scientific Research (Award
No FA-9550-20-1-0366) and the Robert A Welch Foundation (A-1943-20210327). AVS
thanks Welch Foundation (Grant No A-1547) for support.

## Disclosures

The authors declare no conflicts of interest.

## Data Availability

No data were generated or analyzed in the presented research.

## References

  * [1] J. Midwinter and J. Warner, “Up-conversion of near infrared to visible radiation in lithium-meta-niobate,” Journal of Applied Physics 38, 519–523 (1967).
  * [2] L. Ma, O. Slattery, and X. Tang, “Single photon frequency up-conversion and its applications,” Physics Reports 521, 69–94 (2012).
  * [3] R. T. Thew, S. Tanzilli, L. Krainer, S. C. Zeller, A. Rochas, I. Rech, S. Cova, H. Zbinden, and N. Gisin, “Low jitter up-conversion detectors for telecom wavelength ghz qkd,” New Journal of Physics 8, 32 (2006).
  * [4] A. P. Vandevender and P. G. Kwiat, “High efficiency single photon detection via frequency up-conversion,” Journal of Modern Optics 51, 1433–1445 (2004).
  * [5] M. A. Albota and F. N. Wong, “Efficient single-photon counting at 1.55 $\mu$m by means of frequency upconversion,” Optics Letters 29, 1449–1451 (2004).
  * [6] C. Langrock, E. Diamanti, R. V. Roussev, Y. Yamamoto, M. M. Fejer, and H. Takesue, “Highly efficient single-photon detection at communication wavelengths by use of upconversion in reverse-proton-exchanged periodically poled linbo 3 waveguides,” Optics Letters 30, 1725–1727 (2005).
  * [7] L. Ma, O. Slattery, and X. Tang, “Experimental study of high sensitivity infrared spectrometer with waveguide-based up-conversion detector,” Optics Express 17, 14395–14404 (2009).
  * [8] Y.-P. Huang, V. Velev, and P. Kumar, “Quantum frequency conversion in nonlinear microcavities,” Optics Letters 38, 2119–2121 (2013).
  * [9] S. Lefrancois, A. S. Clark, and B. J. Eggleton, “Optimizing optical bragg scattering for single-photon frequency conversion,” Phys. Rev. A 91, 013837 (2015).
  * [10] A. Barh, P. J. Rodrigo, L. Meng, C. Pedersen, and P. Tidemand-Lichtenberg, “Parametric upconversion imaging and its applications,” Advances in Optics and Photonics 11, 952–1019 (2019).
  * [11] A. Ashik, C. F. O’Donnell, S. C. Kumar, M. Ebrahim-Zadeh, P. Tidemand-Lichtenberg, and C. Pedersen, “Mid-infrared upconversion imaging using femtosecond pulses,” Photonics Research 7, 783–791 (2019).
  * [12] W. H. Louisell, A. Yariv, and A. E. Siegman, “Quantum fluctuations and noise in parametric processes. i.” Phys. Rev. 124, 1646–1654 (1961).
  * [13] J. Tucker and D. F. Walls, “Quantum theory of parametric frequency conversion,” Annals of Physics 52, 1–15 (1969).
  * [14] P. Kumar, “Quantum frequency conversion,” Optics Letters 15, 1476–1478 (1990).
  * [15] J. Huang and P. Kumar, “Observation of quantum frequency conversion,” Physical Review Letters 68, 2153 (1992).
  * [16] H. J. McGuinness, M. G. Raymer, C. J. McKinstrie, and S. Radic, “Quantum frequency translation of single-photon states in a photonic crystal fiber,” Phys. Rev. Lett. 105, 093604 (2010).
  * [17] S. Ates, I. Agha, A. Gulinatti, I. Rech, M. T. Rakher, A. Badolato, and K. Srinivasan, “Two-photon interference using background-free quantum frequency conversion of single photons emitted by an inas quantum dot,” Physical Review Letters 109, 147405 (2012).
  * [18] L. Fan, C.-L. Zou, M. Poot, R. Cheng, X. Guo, X. Han, and H. X. Tang, “Integrated optomechanical single-photon frequency shifter,” Nature Photonics 10, 766–770 (2016).
  * [19] C. Chen, J. E. Heyes, J. H. Shapiro, and F. N. Wong, “Single-photon frequency shifting with a quadrature phase-shift keying modulator,” Scientific Reports 11, 1–7 (2021).
  * [20] D. Zhu, C. Chen, M. Yu, L. Shao, Y. Hu, C. Xin, M. Yeh, S. Ghosh, L. He, C. Reimer _et al._ , “Spectral control of nonclassical light pulses using an integrated thin-film lithium niobate modulator,” Light: Science & Applications 11, 1–9 (2022).
  * [21] C.-K. Hong, Z.-Y. Ou, and L. Mandel, “Measurement of subpicosecond time intervals between two photons by interference,” Physical Review Letters 59, 2044 (1987).
  * [22] M. Barbieri, E. Roccia, L. Mancino, M. Sbroscia, I. Gianani, and F. Sciarrino, “What hong-ou-mandel interference says on two-photon frequency entanglement,” Scientific Reports 7, 1–6 (2017).
  * [23] S. Walborn, A. De Oliveira, S. Pádua, and C. Monken, “Multimode hong-ou-mandel interference,” Physical Review Letters 90, 143601 (2003).
  * [24] R.-B. Jin, R. Shiina, and R. Shimizu, “Quantum manipulation of biphoton spectral distributions in a 2d frequency space toward arbitrary shaping of a biphoton wave packet,” Optics Express 26, 21153–21158 (2018).
  * [25] N. Fabre and S. Felicetti, “Parameter estimation of time and frequency shifts with generalized hong-ou-mandel interferometry,” Physical Review A 104, 022208 (2021).
  * [26] T. Douce, A. Eckstein, S. P. Walborn, A. Z. Khoury, S. Ducci, A. Keller, T. Coudreau, and P. Milman, “Direct measurement of the biphoton wigner function through two-photon interference,” Scientific Reports 3, 1–6 (2013).
  * [27] G. Boucher, T. Douce, D. Bresteau, S. P. Walborn, A. Keller, T. Coudreau, S. Ducci, and P. Milman, “Toolbox for continuous-variable entanglement production and measurement using spontaneous parametric down-conversion,” Physical Review A 92, 023804 (2015).
  * [28] X. Liu, T. Li, J. Wang, M. R. Kamble, A. M. Zheltikov, and G. S. Agarwal, “Probing ultra-fast dephasing via entangled photon pairs,” Optics Express 30, 47463–47474 (2022).
  * [29] R. Tyumenev, J. Hammer, N. Joly, P. S. J. Russell, and D. Novoa, “Tunable and state-preserving frequency conversion of single photons in hydrogen,” Science 376, 621–624 (2022).
  * [30] A. V. Sokolov, “Giving entangled photons new colors,” Science 376, 575–576 (2022).
  * [31] J. Hammer, M. V. Chekhova, D. R. Häupl, R. Pennetta, and N. Y. Joly, “Broadly tunable photon-pair generation in a suspended-core fiber,” Physical Review Research 2, 012079 (2020).
  * [32] Y. Shih, _An introduction to quantum optics: photon and biphoton physics_ (CRC press, 2020).
  * [33] J. F. Kenney and E. S. Keeping, _Mayhematics of Statistics_ (D. van Nostrand, 1939).
  * [34] E. T. Whittaker and G. Robinson, _The Calculus of Observations: a Treatise on Numerical Mathematics_ (Blackie and Son limited, 1924).
  * [35] F. Bouchard, A. Sit, Y. Zhang, R. Fickler, F. M. Miatto, Y. Yao, F. Sciarrino, and E. Karimi, “Two-photon interference: the hong–ou–mandel effect,” Reports on Progress in Physics 84, 012402 (2020).
  * [36] A. Hochrainer, M. Lahiri, M. Erhard, M. Krenn, and A. Zeilinger, “Quantum indistinguishability by path identity and with undetected photons,” Reviews of Modern Physics 94, 025007 (2022).
  * [37] A. Lyons, G. C. Knee, E. Bolduc, T. Roger, J. Leach, E. M. Gauger, and D. Faccio, “Attosecond-resolution hong-ou-mandel interferometry,” Science advances 4, eaap9416 (2018).
  * [38] G. Di Martino, Y. Sonnefraud, M. S. Tame, S. Kéna-Cohen, F. Dieleman, i. m. c. K. Özdemir, M. S. Kim, and S. A. Maier, “Observation of quantum interference in the plasmonic hong-ou-mandel effect,” Phys. Rev. Appl. 1, 034004 (2014).
  * [39] B. Ndagano, H. Defienne, D. Branford, Y. D. Shah, A. Lyons, N. Westerberg, E. M. Gauger, and D. Faccio, “Quantum microscopy based on hong–ou–mandel interference,” Nature Photonics 16, 384–389 (2022).
  * [40] A. Eshun, B. Gu, O. Varnavski, S. Asban, K. E. Dorfman, S. Mukamel, and T. Goodson III, “Investigations of molecular optical properties using quantum light and hong–ou–mandel interferometry,” Journal of the American Chemical Society 143, 9070–9081 (2021).