{'fraction_non_alphanumeric': 0.03797550854883445, 'fraction_numerical': 0.009761820275329678, 'mean_word_length': 4.65178242540826, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A framework for the construction of generative models for mesoscale structure in multilayer networks

Multilayer networks allow one to represent diverse and coupled connectivity patterns—such as time-dependence, multiple subsystems, or both—that arise in many applications and which are difﬁcult or awkward to incorporate into standard network representations. In the study of multilayer networks, it is important to investigate mesoscale (i.e., intermediate-scale) structures, such as dense sets of nodes known as communities, to discover network features that are not apparent at the microscale or the macroscale. The ill-deﬁned nature of mesoscale structure and its ubiquity in empirical networks make it crucial to develop generative models that can produce the features that one encounters in empirical networks. Key purposes of such models include generating synthetic networks with empirical properties of interest, benchmarking mesoscale-detection methods and algorithms, and inferring structure in empirical multilayer networks. In this paper, we introduce a framework for the construction of generative models for mesoscale structures in multilayer networks. Our framework provides a standardized set of generative models, together with an associated set of principles from which they are derived, for studies of mesoscale structures in multilayer networks. It uniﬁes and generalizes many existing models for mesoscale structures in fully ordered (e.g., temporal) and unordered (e.g., multiplex) multilayer networks. One can also use it to construct generative models for mesoscale structures in partially ordered multilayer networks (e.g., networks that are both temporal and multiplex). Our framework has the ability to produce many features of empirical multilayer networks, and it explicitly incorporates a user-speciﬁed dependency structure between layers. We discuss the parameters and properties of our framework, and we illustrate examples of its use with benchmark models for community-detection methods and algorithms in multilayer networks.


I. INTRODUCTION
One can model many physical, technological, biological, financial, and social systems as networks. The simplest type of network is a graph [1], which consists of a set of nodes (which represent entities) and a set of edges that encode interactions between those nodes. One can consider either unweighted graphs or weighted graphs, in which each edge has a weight that quantifies the strength of its associated interaction. Edges can also incorporate directions to represent asymmetric interactions or signs to differentiate between positive and negative interactions.
patterns in somebody's Facebook friendship network today may depend both on the connectivity patterns in that person's Facebook friendship network last year (temporal) and on the connectivity patterns in that person's Twitter followership network today (multiplex). Data sets that have multilayer structures are increasingly available (e.g., see Table 2 of [4]). A natural type of multilayer network consists of a sequence of dependent single-layer networks, where layers may correspond to different temporal snapshots, different types of related interactions that occur during a given time interval, and so on. Following existing terminology [8,9], we refer to an instance of a node in one layer as a "state node" (see Sec. II A).
There is a diverse set of models for multilayer networks. (We overview them in Sec. II C.) Many of these models take a specific type of dependency (e.g., a temporal one) as their starting point. In this paper, we introduce a framework for the construction of generative models for multilayer networks that incorporate a wide variety of structures and dependencies. It is broad enough to unify many existing, more restrictive interlayer specifications, but it is also easy to customize to yield multilayer network models for many specific cases of interest. Key purposes of such generative models include (1) generating synthetic networks with empirical features of interest, (2) benchmarking methods and algorithms for detecting mesoscale structures, and (3) inferring structure in empirical multilayer networks.

A. A unifying framework
A key feature of multilayer networks is their flexibility, which allows one to incorporate many different types of data as part of a single structure. In this spirit, our goal is to provide a general, unifying framework that enables users to construct generative models of multilayer networks with a large variety of features of interest in empirical multilayer networks by appropriately constraining the parameter space of our framework. We accomplish this in two consecutive steps. First, we partition the set of state nodes of a multilayer network. Second, we allocate edges, given a multilayer partition. We focus on modeling dependency at the level of partitions (as was done in Ref. [10]), rather than with respect to edges. Additionally, we treat the process of generating a multilayer partition separately from that of generating edges for a given multilayer partition. This modular approach yields random structures that can capture a wide variety of interlayer-dependency structures, including temporal and/or multiplex networks, appearance and/or disappearance of entities, uniform or nonuniform dependencies between state nodes from different layers, and others. For a specified interlayer-dependency structure, one can then use any (single-layer or multilayer) network model with a planted partition (i.e., a planted-partition network model) to generate a wide variety of network features, including weighted edges, directed edges, and spatially embedded layers.
The flexibility of our framework to generate multilayer networks with a specified dependency structure between different layers makes it possible to (1) gain insight into whether, when, and how to build interlayer dependencies into methods for studying many different types of multilayer networks and (2) generate tunable "benchmark models" that allow a principled comparison for community-detection (and, more generally, "meso-set-detection") tools for multilayer networks, including for complicated situations that arise in many applications (such as networks that are both temporal and multiplex) but thus far have seldom or never been studied. In many benchmark models, one plants a partition of a network into well-separated "meso-sets" (e.g., communities), and one thereby imposes a so-called "ground truth" (should one wish to use such a notion) [11] that a properly deployed meso-set-detection method ought to be able to recover. Benchmark networks with known structural properties can be important for (1) analyzing and comparing the performance of different meso-set-detection tools, (2) achieving a better understanding of the inner workings of meso-set-detection tools, and (3) determining which tool(s) may be most appropriate in a given situation.
One can also use our framework to generate synthetic networks with desired empirical properties, to generate null networks, and to explore "detectability limits" of mesoscale structures. (See, for example, the study of detectability thresholds in Ref. [10] for a model that is a special case of our framework.) With some further work, it is also possible to use our framework to develop models for statistical inference. In our concluding discussion (see Sec. VII), we suggest directions for how to apply our framework to the task of statistical inference. Our intention in designing such a flexible framework is to ensure that the generative models that it provides remain useful as researchers consider progressively more general multilayer networks in the coming years.

B. Paper outline
Our paper proceeds as follows. In Sec. II, we give an introduction to multilayer networks, overview mesoscale structure in networks, and review related work on generative models for mesoscale structure. In Sec. III, we formally introduce the notation that we use throughout the paper. In Sec. IV, we explain how we generate a multilayer partition with a specified dependency structure between layers. We also give examples of how to constrain the parameter space of our framework to generate several different types of multilayer networks and discuss what we expect to be common use cases (including temporal structure [10,12,13] and multiplex structure [9,13]) in detail. In Sec. V, we describe how we generate edges that are consistent with the planted partition. In Sec. VI, we illustrate the use of our framework as a way to construct benchmark models for multilayer community detection. In Sec. VII, we summarize our main results and briefly discuss possible extensions of our work to enable statistical inference on empirical multilayer networks with our framework.
Along with this paper, we include code [14] that users can modify to readily incorporate different types of "null distributions" (see Sec. IV A), "interlayer-dependency structures" (see Sec. IV B), and planted-partition network models (see Sec. V). The model instantiations that one needs for generating the figures in Sec. VI are available online. 1

A. Multilayer networks
The simplest type of network is a graph G = (V, E ), where V = {1, . . . , n} is a set of nodes (which correspond to entities) and E ⊆ V × V is a set of edges. Using a graph, one can encode the presence or absence of connections (the edges) between entities (the nodes). However, in many situations, it is desirable to include more detailed information about connections between entities. A common extension is to allow each edge to have a weight, which one can use to represent the strength of a connection. We assign a weight to each edge using a weight function w : E → R.
One can also generalize the notion of a graph to encode different aspects of connections between entities, such as connections at different points in time or multiple types of relationships. We adopt the framework of multilayer networks [4,[6][7][8]15] to encode such connections in a network. In a multilayer network, a node can be present in a variety of different states, where each state is also characterized by a variety of different aspects. In this setting, edges connect state nodes, each of which is the instantiation of a given node in a particular state, to each other. The set of all state nodes of a given entity corresponds to a single physical node (which represents the entity). The set of all state nodes in a given state and the edges between those state nodes constitute one layer of a multilayer network. In the remainder of this paper, we use the terms "physical node" and "node" interchangeably and the terms "layer" and "state" interchangeably. One can think of the aspects as features that one needs to specify to identify the state of a node. In other words, a state is a collection of exactly one element from each aspect. For convenience, we introduce a mapping that assigns an integer label to each element of an aspect. That is, we map the l a elements of aspect a to the elements of a set {1, . . . , l a } of integers. Aspects can be unordered (e.g., social-media platform) or ordered (e.g., time). Most empirical investigations of multilayer networks focus on a single aspect (e.g., temporal [16][17][18] or multiplex [19][20][21]). However, many real-world situations include more than one aspect (e.g., a multiplex network that changes over time). For an ordered aspect, we require the mapping to respect the ordering of the aspect. For example, t i → i for time, where t 1 · · · t l a is a set of discrete time points. A multilayer network can include an arbitrary number of ordered aspects and an arbitrary number of unordered aspects, and one can generalize these ideas further (e.g., by introducing a time horizon) [4].
To illustrate the above ideas, consider a hypothetical social network with connections on multiple social-media platforms (Facebook, Twitter, and LinkedIn) between the same set of people at different points in time. In this example network, there are two aspects: social-media platform and time. [One can consider "type of connection" (e.g., "friendship" and "following") as a third aspect, but we restrict our example to two aspects for simplicity.] The first aspect is ordered, and its number of elements is equal to the number of time points or time intervals. (For simplicity, we often refer simply to "time points" in our discussions.) The second aspect is unordered and consists of three elements: Facebook (which we label with the integer "1"), Twitter (which we label with "2"), and LinkedIn (which we label with "3"). If we assume that the FIG. 1. Toy example of a multilayer network with three physical nodes and two aspects. We represent undirected intralayer edges using solid black lines, directed interlayer edges using dotted red lines with an arrowhead, and undirected interlayer edges using dashed blue arcs. The first aspect is ordered and corresponds to time. It has two time points (which we label as "1" and "2"). That is, the set of labels of the 1st aspect is L 1 = {1, 2}. The second aspect is unordered and represents a social-media platform. It has three elements: Facebook (labeled "1"), Twitter (labeled "2"), and LinkedIn (labeled "3"). That is, the set of labels of the 2nd aspect is L 2 = {1, 2, 3}. In this example, a state node takes the form (i, (α 1 , α 2 )), with physical node i ∈ {1, 2, 3}, aspect α 1 ∈ {1, 2}, and aspect α 2 ∈ {1, 2, 3}. The total number of layers is l = |L 1 | × |L 2 | = 6. time resolution is daily and spans the year 2010, an example of a state is the tuple (01-Jan-2010, Twitter), which is (1,2) in our shorthand notation. We show a multilayer representation of this example social network in Fig. 1.
Edges in a multilayer network can occur either between nodes in the same state (i.e., intralayer edges) or between nodes in different states (i.e., interlayer edges). An example of an intralayer edge in Fig. 1 is (1, (1, 1)) ↔ (2, (1, 1)), indicating that entity 1 is friends with entity 2 on Facebook at time 1. All interlayer edges in Fig. 1 are diagonal (because they connect state nodes that correspond to the same physical node). Interlayer edges between layers at different times for a given social-media platform are ordinal (because they connect state nodes with successive time labels 2 ) and directed. Interlayer edges between concurrent layers that correspond to different social-media platforms are categorical and undirected. (In this example, such edges occur between state nodes from all pairs of concurrent layers, although one can envision situations in which interlayer edges occur only between state nodes from a subset of such layers.) An example of an ordinal intralayer edge in Fig. 1 is (1, (1, 1)) → (1, (2, 1)), indicating a connection from entity 1 in Facebook at time 1 to entity 1 in Facebook at time 2. An example of a categorical intralayer edge in Fig. 1 is (1, (1, 1)) ↔ (1, (1, 2)), indicating a connection between node 1 in Facebook at time 1 and node 1 in Twitter at time 1. One can also imagine other types of interlayer edges in an example like the one in Fig. 1, as there may be edges between state nodes at 2 It is also possible to have ordinal interlayer edges in aspects other than time.

B. Mesoscale structures in networks
Given a (single-layer or multilayer) network representation of a system, it is often useful to apply a coarse-graining technique to investigate features that lie between those at the microscale (e.g., nodes, pairwise interactions between nodes, or local properties of nodes) and those at the macroscale (e.g., total edge weight, degree distribution, or mean clustering coefficient) [1]. One thereby studies mesoscale features, such as community structure [22][23][24][25], core-periphery structure [26,27], role structure [28], and others.
Our framework can produce multilayer networks with any of the above mesoscale structures. Notwithstanding this flexibility, an important situation is multilayer networks with "community structure," which is the most commonly studied type of mesoscale structure [1]. Community detection is part of the standard toolkit for studying single-layer networks [22,23,29], and efforts at community detection in multilayer networks have led to insights in applications such as brain and behavioral networks in neuroscience [30], financial correlation networks [16], committee and voting networks in political science [31,32], networks of interactions between bacterial species [21], disease-spreading networks [12], social networks in Indian villages [33], and much more.
Loosely speaking, a community in a network is a set of nodes that are "more densely" connected to each other than they are to nodes in the rest of the network [22,23,29,34,35]. Typically, a "good community" should be a set of nodes that are "surprisingly well connected" in some sense, but what one means by "surprising" and "well connected" is often application-dependent and subjective. In many cases, a precise definition of "community" depends on the method that one uses to detect communities. In particular, many popular community-detection approaches in single-layer and multilayer networks define communities as sets in a partition of a network that optimizes an objective function such as modularity [1,23,36]; stability [37][38][39]; InfoMap and its variants [9,40]; likelihood functions that are derived from stochastic block models (SBMs), which are models for partitioning a network into sets of nodes with statistically homogeneous connectivity patterns [21,[41][42][43][44]; and others [4,45,46]. In this paper, we refer to a partition of the set of nodes of a single-layer network as a single-layer partition and to a partition of the set of state nodes of a multilayer network as a multilayer partition. The two primary differences between a community in a multilayer partition and a community in a single-layer partition are that (1) the former can include state nodes from different layers; and (2) "induced communities" (see Sec. III) in one layer of a multilayer network may depend on connectivity patterns in other layers. For many notions of (single-layer or multilayer) community structure, including the most prominent methods, one cannot exactly solve a community-assignment problem in polynomial time (unless P = NP) [29,[47][48][49]; and popular scalable heuristics currently have few or no theoretical guarantees on how closely an identified partition resembles an optimal partition. These issues apply more generally to the field of cluster analysis, such as in graph partitioning [50], and many of the problems that plague community detection also apply to detecting other types of mesoscale structures.
Throughout our paper, to make it clear which results and observations apply to network community structure in particular and which apply to network mesoscale structure more generally, we use the term "community" when referring to a set in a partition that corresponds to community structure and the term "meso-set" (see Sec. III) when referring to a set in a partition that corresponds to any type of mesoscale structure. (In particular, a community is a type of meso-set.)

C. Generative models for mesoscale structure
The ubiquity and diversity of mesoscale structures in empirical networks make it crucial to develop generative models of mesoscale structure that can produce features that one encounters in empirical networks. Broadly speaking, the goal of such generative models is to construct synthetic networks that resemble real-world networks when one appropriately constrains and/or calibrates the parameters of a model using information about a scenario or application. Generative models of mesoscale structure can serve a variety of purposes, such as (1) generating benchmark network models for comparing meso-set-detection methods and algorithms [10,[51][52][53][54][55]; (2) undertaking statistical inference on empirical networks [10,42,56]; (3) generating synthetic networks with a desired set of properties [20,57]; (4) generating null models to take into account available information about an empirical network [58]; and (5) investigating "detectability limits" for mesoscale structure, as one can plant partitions that, under suitable conditions, cannot subsequently be detected algorithmically, despite the fact that they exist by construction [10,59,60].
One of the main challenges in constructing a realistic generative model (even for single-layer networks) is the breadth of possible empirical features in networks. The available generative models for mesoscale structure in single-layer networks usually focus on replicating a few empirical features at a time (rather than all of them at once): heterogeneous degree distributions and community-size distributions [42,52,61]; edge-weight distribution [51,56,62]; spatial embeddedness [12,63]; and so on. Multilayer networks inherit all of the empirical features of single-layer networks, and they also have a key additional one: dependencies between layers. These interlayer dependencies can be ordered (as in most models of temporal networks), unordered (as in multiplex networks), or combinations of these. However, despite this variety, existing generative models for mesoscale structure in multilayer networks allow only a restrictive set of interlayer dependencies. They assume either a temporal structure [10,12,53,[64][65][66], a simplified multiplex structure with the same planted partitions across all layers [20,43,[67][68][69], or independent groups of layers such that layers in the same group have identical planted partitions [9,21]. Using an alternative approach, a very recent model generates multilayer partitions that satisfy the constraint that nonempty "induced meso-sets" (see Sec. III) in different layers are identical if and only if they correspond to the same meso-set [70]. Recent work [13] on the link between multilayer modularity maximization and maximum-likelihood estimation of multilayer SBMs allows either temporal or multiplex interlayer dependencies with induced partitions that can vary across layers, but it makes restrictive assumptions on interlayer dependencies (e.g., all layers have the same set of nodes, interlayer dependencies are "diagonal" and "layer-coupled," and so on). Importantly, in all aforementioned generative models of mesoscale structure in multilayer networks, interlayer dependencies are either (1) not explicitly specifiable or (2) special cases of the framework that we discuss in this paper.

III. NOTATION
We now present a comprehensive set of notation for multilayer networks. Different subsets of notation are useful for different situations.
We consider a multilayer network M = (V M , E M , V, L) with n = |V| nodes (i.e., physical nodes) and l = |L| layers (i.e., states). We use d to denote the number of aspects and L a = {1, . . . , l a } to denote the labels of aspect a (where a ∈ {1, . . . , d}). We use O to denote the set of ordered aspects and U to denote the set of unordered aspects. For each ordered aspect a ∈ O, we assume that the labels L a reflect the ordering of the aspect. That is, for all α, β ∈ L a , we require that α < β if and only if α precedes β. We say that a multilayer network is fully ordered if U = ∅, unordered if O = ∅, and partially ordered otherwise. The set L = L 1 × · · · × L d of states is the Cartesian product of the aspects, where a state α ∈ L is an integer-valued vector of length d and each of its d entries specifies an element of the corresponding aspect. Note that l = |L| = d a=1 l a . We use (i, α) ∈ V M ⊆ V × L to denote the state node (i.e., "node-layer tuple" [4]) of physical node i ∈ V in state α ∈ L. We include a state node in V M if and only if the corresponding node exists in that state. The edges E M ⊆ V M × V M in a multilayer network connect state nodes to each other. We use ((i, α), ( j, β)) to denote a directed edge from (i, α) to ( j, β). For two state nodes, (i, α) and ( j, β), that are connected to each other via a directed edge ((i, α), ( j, β)) ∈ E M , we say that (i, α) is an in-neighbor of ( j, β) and that ( j, β) is an out-neighbor of (i, α). We categorize the edges into intralayer edges E L , which have the form ((i, α), ( j, α)) and link entities i and j in the same state α, and interlayer (i.e., coupling) edges E C , which have the form ((i, α), ( j, β)) for α = β. We thereby decompose the edge set as We define a weighted multilayer network by introducing a weight function w : E M → R (which is analogous to the weight function for weighted single-layer networks), which encodes the edge weights within and between layers. For an unweighted multilayer network, w(e) = 1 for all e ∈ E M . We encode the connectivity pattern of a multilayer network using an adjacency tensor A, which is analogous to an adjacency matrix for a single-layer network, with entries Note that G M = (V M , E M ) is a graph on the state nodes of the multilayer network M. We refer to G M as the flattened network of M. The adjacency matrix of a flattened network is the "supra-adjacency matrix" [4,15,71,72] of its associated multilayer network. One obtains a supra-adjacency matrix by flattening 3 an adjacency tensor [see Eq. (1)] of a multilayer network. The multilayer network and the corresponding flattened network encode the same information [4,73], provided one keeps track of the correspondence between the nodes in the flattened network and the physical nodes and layers of the multilayer network. In a similar vein, aspects provide a convenient way to keep track of the correspondence between state nodes and, for example, temporal and multiplex relationships.
We denote a multilayer partition with n set sets by S = {S 1 , . . . , S n set }, where n set s=1 S s = V M and S s ∩ S r = ∅ for s = r. We represent a partition S using a partition tensor S with entries S i,α , where S i,α = s if and only if the state node (i, α) is in the set S s . We refer to a partition of a temporal network (a common type of multilayer network) as a temporal partition and a partition of a multiplex network (another common type of multilayer network) as a multiplex partition. A multilayer partition induces a partition We refer to a set S s of a partition S as a meso-set, to S| α as the induced partition of S in layer α, and to S s | α as the induced meso-set of S s in layer α. A community is a set S s in a partition that corresponds to community structure. We call s ∈ {1, . . . , n set } the label of meso-set S s . One can examine overlapping meso-sets by identifying multiple state nodes from a single layer with the same physical node.

IV. GENERATING MULTILAYER PARTITIONS
The systematic analysis of dependencies between layers is a key motivation for analyzing a single multilayer network, instead of examining several single-layer networks independently. The goal of our partition-generation process is to model interlayer dependency in a way that can incorporate diverse types of dependencies. We now motivate our partitionmodeling approach; we describe it in detail in Secs. IV A-IV D.
The complexity of dependencies between layers can make it difficult to explicitly specify a joint probability distribution for meso-set assignments, especially for unordered or partially ordered multilayer networks. To address this issue, we require only the specification of conditional probabilities for a state node's meso-set assignment, given the assignments of all other state nodes. The idea of conditional models is old and follows naturally from Markov-chain theory [74,75]. Specifying conditional models (which capture different dependency features separately) rather than joint models (which try to capture many dependency features at once) is convenient for numerous situations. For example, conditionally specified distributions have been applied to areas such as spatial data modeling [76], imputation of missing data [75], secure disclosure of information [77], dependency networks for combining databases from different sources [78], and Gibbs sampling [79,80].
A problem with conditional models is that their conditional probability distributions are not necessarily compatible, in the sense that there may not exist a joint distribution that realizes all conditional distributions [78,[80][81][82][83]. Although several methods have been developed for checking the compatibility of discrete conditional distributions, these either make restrictive assumptions on the conditional distributions or have scalability constraints that significantly limit practical use [75,84]. Nevertheless, the employment of conditional models (even if potentially incompatible) is common [85], and they have been useful in many applications, provided that one cautiously handles any potential incompatibilities [75]. For our use case, potential incompatibilities arise in unordered or partially ordered multilayer networks. An issue that can result from incompatibility is nonuniqueness of a joint distribution. We carefully design our partition-generation process such that it is well-defined irrespective of whether conditional distributions are compatible. In particular, we show that convergence is guaranteed and we address nonuniqueness by appropriately sampling initial conditions (see Appendix B). We also suggest empirical checks to ensure that a generated partition reflects planted interlayer dependencies.
In our framework, we define conditional probabilities in two parts: we separately specify (1) independent layer-specific random components and (2) interlayer dependencies. For a choice of interlayer dependencies, there are several features that one may want to allow in a multilayer partition. These include variation in the numbers and sizes of meso-sets across layers (e.g., meso-sets can gain state nodes, lose state nodes, appear, and disappear) and the possibility for these meso-set variations to incorporate features of the application at hand. (For example, in a temporal network, one may not want a meso-set to reappear after it disappears.) We ensure that such features are possible for any choice of interlayer dependency via independent layer-specific null distributions that specify the set of possible meso-set assignments for each layer and determine the expected sizes and expected number of mesosets in the absence of interlayer dependencies.
Interlayer dependencies should reflect the type of multilayer network that one is investigating. For example, in a temporal network, dependencies tend to be stronger between contiguous layers and weaker for pairs of layers that are farther apart. It is common to assume that interlayer dependencies are uniform across state nodes for a given pair of layers (i.e., interlayer dependencies are "layer-coupled" [4]) and occur only between state nodes that correspond to the same physical node (i.e., interlayer dependencies are "diagonal" [4]) [12,16,21,32,60,71,72,86]. However, these assumptions are too restrictive in many cases (e.g., situations in which dependencies depend on the state nodes or in which dependencies exist between state nodes that do not correspond to the same physical node) [9,45,69,[87][88][89]. To ensure that one can relax these assumptions, we allow dependencies to be specified at the level of state nodes (or at the level of layers, when a user assumes that dependencies are layer-coupled). We encode these interlayer dependencies in a user-specified interlayerdependency tensor that determines the extent to which a state node's meso-set assignment in one layer "depends directly" on the assignments of state nodes in other layers (see Secs. IV B and IV C). Our independence assumption on the null dis-tributions allows us to encode all of the interlayer dependencies in a single object (namely, an interlayer-dependency tensor). The entries of an interlayer-dependency tensor specify an "interlayer-dependency network" and correspond to the causal links for the flow of information (in the information-theoretic sense) between different layers of a multilayer network. Therefore, they should reflect any constraints that one wishes to impose on the direct flow of information between layers. Longer paths in an interlayer-dependency network yield indirect dependencies between structures in different layers.
After specifying the interlayer-dependency structure, we define the conditional meso-set assignment probabilities such that we either sample a state node's meso-set assignment from the corresponding null distribution or obtain it by copying the assignment of another state node (based on the interlayerdependency tensor). Using these conditional probabilities, we define an iterative update process on the meso-set assignments of state nodes to generate multilayer partitions with dependencies between induced partitions in different layers. When updating the meso-set assignments of state nodes, we respect the ordering (e.g., temporal ordering) of an ordered aspect. For a fully ordered multilayer network, our update process reduces to sequentially sampling an induced partition for each layer based on the induced partitions in previous layers. For an unordered multilayer network, our update process defines a Markov chain on the space of multilayer partitions. We sample partitions from a stationary distribution of this Markov chain. This sampling strategy is known as (pseudo-)Gibbs sampling [78][79][80]85]. (We use the word "pseudo" because the conditional probabilities that we use to define the Markov chain are not necessarily compatible.) For a partially ordered multilayer network, our update process combines these two sampling strategies.
In Sec. IV A, we describe possible choices for the independent, layer-specific null distributions. In Sec. IV B, we explain our framework for generating a multilayer partition with general interlayer dependencies. In Sec. IV C, we focus on the specific situation in which interlayer dependencies are layer-coupled and diagonal, and we also assume that a physical node is present in all layers (i.e., the network is "fully interconnected" [4]). In Sec. IV D, we illustrate the properties of example temporal and multiplex partitions that are generated by models that we construct from our framework. Additionally, we take advantage of the tractability of our generative model for the special case of temporal partitions to highlight some of this special case's properties analytically.

A. Null distribution
We denote the null distribution of layer α by P α 0 and the set of all null distributions by P 0 = {P α 0 : α ∈ L}. A simple choice for the null distributions is a categorical distribution, where for each layer α and each meso-set label s, we fix the probability p α s that an arbitrary state node in layer α is assigned to a meso-set s in the absence of interlayer dependencies. That is, where n set is the total number of meso-sets in the multilayer partition and n set s=1 p α s = 1 for all α ∈ L. The set {1, . . . , n set } 023100-6 is the set of meso-set labels, and the support of a null distribution is the set of labels that have nonzero probability. In the absence of interlayer dependencies, a categorical null distribution corresponds to specifying independent multinomial distributions of the sizes of induced meso-sets for each layer and fixing the expected size np α s of each induced meso-set. (Recall that n is the number of physical nodes.) Therefore, by choosing the probabilities p α s , one has some control over the expected number and sizes of meso-sets in a sampled multilayer partition. A natural choice for p α is to sample it from a Dirichlet distribution, which is the conjugate prior of a categorical distribution [90,91]. One can think of the Dirichlet distribution, which is the multivariate form of the beta distribution, as a probability distribution over the space of all possible categorical distributions with a given number of categories. Any other (probabilistic or deterministic) choice for p α is also possible. We give further detail about the Dirichlet distribution in Appendix A, where we discuss how one can vary its parameters to control the expected number and sizes of meso-sets in the absence of interlayer dependencies. This can allow a user to generate, for example, a null distribution with equally-sized meso-sets or a null distribution with a few large meso-sets and many small meso-sets. Furthermore, irrespective of the particular choice of categorical distribution, it may be desirable to have meso-sets that are in the support of the null distribution in some, but not all, layers. In Appendix A, we give examples of how one can sample the support of the null distributions before sampling the probability vectors to incorporate this property when modeling the birth and/or death of meso-sets in temporal networks and the presence and/or absence of meso-sets in multiplex networks.
In general, the choice of the null distributions can have a large effect on the set of sampled multilayer partitions and ought to be guided by one's use case. For our numerical examples in Sec. VI, we fix a value of n set ∈ {1, . . . , nl}, where n is the number of physical nodes in each layer and l the number of layers (see Sec. III); and we use a symmetric Dirichlet distribution of order n set with concentration parameter θ = 1 (see Appendix A) to sample probability vectors p α of length n set . This produces multilayer partitions in which the expected meso-set labels are the same across layers (and are given by {1, . . . , n set }) in the absence of interlayer dependencies and for which the expected induced meso-set sizes (which are equal to np α s in the absence of interlayer dependencies) differ across layers.

B. General interlayer dependencies
We denote the user-specified interlayer-dependency tensor by P, where P j,β i,α is the probability that state node ( j, β) copies its meso-set assignment from state node (i, α) for any two state nodes (i, α), ( j, β) ∈ V M . Note that P is fixed throughout the copying process. The probability that state node ( j, β) copies its meso-set assignment from an arbitrary state node when state node ( j, β)'s meso-set assignment is updated iŝ where we require thatp j,β 1 for all state nodes ( j, β) ∈ V M . We also require that all intralayer probabilities are 0; that is, P j,α i,α = 0 for all i, j ∈ V and α ∈ L. We say that a state node ( j, β) depends directly on a state node (i, α) if and only if P j,β i,α is nonzero. By extension, we say that a layer β depends directly on a layer α if there exists at least one state node in layer β that depends directly on a state node in layer α. An interlayer-dependency tensor induces an associated interlayer-dependency network, whose edges are all interlayer and directed; these edges point in the direction of information flow between layers. The in-neighbors (see Sec. III) of a state node ( j, β) in this network consist of the state nodes from which ( j, β) can copy a meso-set assignment with nonzero probability. The support of the null distribution P β 0 corresponds to the set of possible meso-set assignments of a state node ( j, β) when the state node does not copy its assignment from one of its in-neighbors in the interlayer-dependency network.
For a given set P 0 of null distributions and a given interlayer-dependency tensor P, a multilayer partition that results from our sampling process depends on four choices: (1) the way in which we update a state-node assignment at a given step; (2) the order in which we update state-node assignments; (3) the initial multilayer partition; and (4) the criteria for convergence of the iterative update process. We discuss points (1), (2), and (3) in the remainder of this section. We describe the sampling process in more detail and discuss convergence in Appendix B.

Update equation
A single meso-set-assignment update depends only on the choice of state node to update and on the current multilayer partition. Let τ be an arbitrary update step of the assignment copying process. Suppose that we are updating the meso-set assignment of state node ( j, β) at step τ and that the current multilayer partition is S (τ ) [with associated partition tensor S(τ )]. We update the meso-set assignment of state node ( j, β) either by copying the meso-set assignment in S (τ ) from one of its in-neighbors in the interlayer-dependency network or by obtaining a new, random meso-set assignment from the null distribution P β 0 for layer β. In particular, with probabilityp j,β , a state node ( j, β) copies its meso-set assignment from one of its in-neighbors in the interlayer-dependency network; and with probability 1 −p j,β , it obtains its meso-set assignment from the null distribution P β 0 . This yields the following update equation at step τ of our copying process: The update equation (4) is at the heart of our framework. It is clear from Eq. (4) that the set P 0 of null distributions is responsible for the specification of meso-set assignments in the absence of interlayer dependencies [i.e., if P j,β i,α = 0 for all (i, α), ( j, β)]. In general,p j,β determines the relative importances of the interlayer dependencies and the null distribution on the meso-set assignment of state node ( j, β). Specifically, whenp j,β = 0, the meso-set assignment of ( j, β) depends only on the null distribution; and whenp j,β = 1, the meso-set assignment of ( j, β) depends only on the meso-set assignments of its in-neighbors in the associated interlayerdependency network.

Update order
The order in which we update meso-set assignments of state nodes via Eq. (4) can influence a generated multilayer partition. As we mentioned in Sec. II A, an aspect of a multilayer network can be either ordered or unordered, and the update order is particularly important when generating a multilayer partition with at least one ordered aspect. In particular, the structure of an interlayer-dependency tensor should reflect the causality that arises from the order of an aspect's elements. For an ordered aspect, structure in a given layer should depend directly only on structure in previous layers. Formally, for each ordered aspect a ∈ O of a multilayer network, we require that where α a denotes the element of state α that corresponds to aspect a and where (as we stated in Sec. III) we require that the labels of an aspect's elements reflect the order of those elements. 4 For example, if we think of the interlayer edges in Fig. 1 as interlayer edges of an interlayer-dependency network, then Eq. (5) states that all edges with nonzero edge weights must respect the arrow of time. That is, it is impossible to have edges of the form (( j, (2, α 2 )), (i, (1, α 2 ))), with nodes i, j ∈ {1, 2, 3} and aspect α 2 ∈ {1, 2, 3}. The update order for state nodes also needs to reflect the causality from the ordered aspects. In particular, if a ∈ O is an ordered aspect and we update a state node in a layer β, then the final meso-set assignments of state nodes in any layer α with α a < β a must already be fixed. For example, in Fig. 1, an edge direction respects the arrow of time; and it is also necessary that the assignments of all state nodes in the first time layer are fixed before one updates the assignments of state nodes in the second time layer. To satisfy this requirement for the assignments, we divide the layers into classes of "order-equivalent" layers, such that two layers α and β are order equivalent if and only if α a = β a for all a ∈ O. We also say that layer α precedes layer β if and only if α a < β a for all a ∈ O. (For example, the three layers in Fig. 1 that correspond to the first time point are order equivalent and precede the three layers that correspond to the second time point.) Based on this equivalence relation, we obtain an ordered set of equivalence classes by inheriting the ordering from the ordered aspects. As a consequence of Eq. (5), the meso-set assignments of state nodes in a given layer depend only on the assignments of state nodes in order-equivalent layers or preceding layers. Therefore, by sorting the classes of order-equivalent layers based on "lexicographic" ordering [92], we ensure that our partition-generation process reflects notions of causality that arise in multilayer networks with at least one ordered aspect. 4 For a fully ordered multilayer network (i.e., when U = ∅), condition (5) is equivalent to the existence of an upper-triangular matrix representation of its flattened interlayer-dependency tensor. [See Fig. 2(a) for an example.] For a partially ordered multilayer network, the same equivalence holds if we "aggregate" [15] its interlayerdependency tensor over its unordered aspects.

Sampling process
The general idea of our sampling process is to simultaneously sample the meso-set assignments of state nodes within each class of order-equivalent layers and to sequentially sample the meso-set assignments of state nodes in non-orderequivalent layers, conditional on the fixed assignments of state nodes in preceding layers.
More specifically, our sampling algorithm proceeds as follows. First, we sample an initial multilayer partition from the null distribution [i.e., S i,α (0) ∼ P α 0 ]. Second, we sample a partition for the first class of order-equivalent layers using (pseudo-)Gibbs sampling. In particular, we iteratively sample a layer uniformly at random from the first class of orderequivalent layers and update the meso-set assignments of all state nodes in that layer 5 based on Eq. (4). This defines a Markov chain on a subspace of multilayer partitions. We repeat the update process for sufficiently many iterations such that we approximately sample the meso-set assignments of the first class of order-equivalent layers from a stationary distribution of this Markov chain. Third, we sample state-node assignments for subsequent classes of order-equivalent layers in the same way, based on fixed state-node assignments from preceding layers. In particular, at each update step τ in Eq. (4), a state node can either copy a fixed meso-set assignment from an in-neighbor in a preceding layer, copy a current meso-set assignment from an in-neighbor in an order-equivalent layer, or obtain a meso-set assignment from the null distribution.
For a fully ordered multilayer network (i.e., a multilayer network with U = ∅), each class of order-equivalent layers consists of a single layer. Consequently, one needs only a single update for each state node for the sampling process to converge. (Subsequent updates would constitute independent samples from the same distribution.) For an unordered multilayer network, our sampling algorithm reduces to (pseudo-)Gibbs sampling, because all layers are order equivalent and there is thus only one class of order-equivalent layers. For a partially ordered multilayer network, there are multiple nonsingleton classes of order-equivalent layers. We use (pseudo-)Gibbs sampling for each class, conditional on the meso-set assignments of state nodes in preceding layers.
In Appendix B, we explain our sampling process in more detail. We describe our (pseudo-)Gibbs sampling procedure for each class of order-equivalent layers; and we show pseudocode that one can use to sample partitions in unordered, partially ordered, and fully ordered multilayer networks in Algorithm 1. We also discuss the convergence properties of our sampling procedure for unordered and partially ordered multilayer networks. As part of this discussion, we examine the effect of the potential incompatibility of the distributions that are defined by Eq. (4).

C. Layer-coupled and diagonal interlayer dependencies
As we mentioned at the beginning of Sec. IV, a particularly useful restriction of an interlayer-dependency 5 The order in which we update state nodes in the same layer has no effect on the sampling process because their conditional distributions are independent. One can even update them in parallel.
Layer-coupled interlayer-dependency tensors (which, in this case, are matrices) for different types of multilayer networks with a single aspect. (a) In a typical temporal network, an induced partition in a layer depends directly only on the induced partition in the previous layer. Therefore, the only nonzero elements of its associated layer-coupled interlayer-dependency tensor occur in the first superdiagonal. (b) In a typical multiplex network, an induced partition in a layer depends directly on the induced partitions in all other layers. We show a layer-coupled example. The copying probability in Eq. (6) isp β = p β for layer β in the temporal layercoupled interlayer-dependency matrix in (a) andp = (l − 1)p for the multiplex layer-coupled interlayer-dependency matrix in (b). We suppress the subscript β in (b), as the copying probability is the same for all layers.
tensor that still allows us to study many situations of interest is to assume that it is layer-coupled and diagonal. For simplicity, we also assume that we have a fully interconnected multilayer network. Under these assumptions, we can express the interlayer-dependency tensor P using a layer-coupled interlayer-dependency tensor P with elements P j,β i,α = δ(i, j) P β α . The probability that state node ( j, β) copies its meso-set assignment from an arbitrary state node when ( j, β)'s meso-set assignment is updated iŝ As before, we require thatp β 1 and P α α = 0 for all α ∈ L. The update equation [see Eq. (4)] then simplifies to which depends on the interlayer-dependency tensor P (which is now independent of state nodes) and the null distributions P 0 . Each term P β α quantifies the extent to which an induced partition in layer β depends directly on an induced partition in layer α. By considering different P, one can generate multilayer networks that correspond to several important scenarios, including temporal networks, multiplex networks, and multilayer networks with more than one aspect (e.g., with combinations of temporal and multiplex features). That is, the above restriction reduces the dimensionality of the interlayerdependency tensor significantly, while still allowing one to analyze several very important situations.
FIG. 3. Block-matrix representation of a layer-coupled interlayer-dependency tensor for a multilayer network that is both temporal and multiplex. This example has more than one aspect, and it combines the features of the examples in Fig. 2. It has two classes of order-equivalent layers (where each class corresponds to one time point), as indicated by the presence of two diagonal blocks. An induced partition in a layer depends directly on the induced partitions in its order-equivalent layers and on the induced partition in the corresponding layer in the preceding class of order-equivalent layers. If we think of the interlayer edges in Fig. 1 as edges in an interlayer-dependency network, then the tensor P (with blocks of size 3 × 3) is the corresponding layer-coupled interlayer-dependency tensor.
In Fig. 2, we show layer-coupled interlayer-dependency tensors for two types of single-aspect multilayer networks: a temporal network and a multiplex network. As we mentioned in Sec. IV, it is useful to think of interlayer dependencies as causal links for the flow of information between layers. In a temporal network, it is typical to assume that an induced partition in a layer depends directly only on induced partitions in the previous layer. There are thus l − 1 copying probabilities (one for each pair of contiguous layers), which we are free to choose. Common examples include choosing the same probability for each pair of contiguous layers [10,12] or making some of the probabilities significantly smaller than the others to introduce change points [43,93]. In a multiplex network, an induced partition in any layer can depend directly on induced partitions in all other layers. This yields l (l − 1) copying probabilities to choose. In Fig. 2(b), we illustrate the simplest case, in which each layer depends equally on every other layer. The layer-coupled interlayer-dependency tensors in Fig. 2 are matrices, so we sometimes refer to them and other examples of layer-coupled interlayer-dependency tensors with a single aspect as layer-coupled interlayer-dependency matrices.
We can also generate multilayer networks with more than one aspect and can thereby combine temporal and multiplex features. In Fig. 3, we illustrate how to construct a layer- We choose a node ordering for each visualization that (whenever possible) emphasizes "persistent" mesoscale structure [16]. We show only the first 15 layers of each multilayer partition, because (as one can see in the NMI heat maps) similarities between induced partitions for p < 1 decay steeply with the number of layers when there are dependencies only between contiguous layers. coupled interlayer-dependency tensor to generate such a multilayer network on a simple example with two aspects, one of which is multiplex and the other of which is temporal.

D. Temporal and multiplex partitions
In Figs. 4 and 5, we show example multilayer partitions that we obtain with the interlayer-dependency tensors from  99, and (f)p = 1, wherep = (l − 1)p is the probability that a state node copies its assignment from another state node, we show color-coded meso-set assignments ( , top) for a single example output partition and NMI values [94] (0 1, bottom) between induced partitions in different layers averaged over a sample of 10 output partitions. (For the temporal example in Fig. 4, note thatp = p.) The parameter values match those in the numerical examples in Sec.VI (with the exception ofp = 0, which we include for completeness). We choose a node ordering that (whenever possible) emphasizes "persistent" mesoscale structure [16].
partitions in different layers. NMI is a measure of similarity between partitions, so we expect the values of NMI to reflect the planted dependencies. As expected, for the temporal examples, partitions are most similar for contiguous layers. Similarities decay with the distance between layers and increase with the value of p. To provide a detailed illustration of the steps that result in a multilayer partition, we focus on the temporal examples in Fig. 4. For the important special case of temporal interlayer dependencies, we can significantly simplify our generative model for multilayer partitions. In particular, there is a single 023100-10 ordered aspect, so the layer index α ∈ N is a scalar and the order of the layers corresponds to temporal ordering. Furthermore, as we mentioned in Sec. IV B, for a fully ordered multilayer network, we require that the order of the mesoset-assignment update process in Eq. (4) respects the order of the layers. The update order of state nodes (1, α), . . . , (n, α) in any layer α is arbitrary (so we can update them simultaneously), but each update is conditional on the meso-set assignments of state nodes in layer α − 1.
The update process that we described in Sec. IV B reduces to three steps: (1) initialize (in an arbitrary order) the meso-set assignments in layer α = 1; (2) take the meso-set assignment of (i, α + 1) to be that of (i, α) with probability p, and sample the meso-set assignment of (i, α + 1) from P α+1 0 with complementary probability 1 − p; and (3) increment α by 1 and repeat steps (2) and (3) until we have updated the meso-set assignments in all layers. (See Algorithm 2 in Appendix C.) As we mentioned in Sec. IV B, convergence is not an issue for this case (or, more generally, for any fully ordered multilayer network), as we need only one iteration through the layers. This three-step generative model for temporal partitions was also suggested by Ghasemian et al. [49]. They used it to derive a detectability threshold when the null distributions are uniform across meso-sets (i.e., θ → ∞ in Appendix A) and intralayer edges are generated independently using a standard SBM. Specifically, one replaces the degree-corrected SBM that we discuss in Sec. V with the non-degree-corrected SBM of [41]. In Appendix C, we highlight properties of a sampled temporal partition that illustrate how the interplay between p and the null distributions affects the evolution of meso-sets across layers (e.g., growth, birth, and death of induced meso-sets). The properties that we highlight hold for any choice of null distribution. In the same appendix, we subsequently illustrate that the particular choice of null distributions can greatly influence resulting partitions. For example, a nonempty overlap between the supports of null distributions of contiguous layers is a necessary condition for meso-sets to gain new state nodes over time. The properties and observations in Appendix C are independent of one's choice of planted-partition network model.

V. GENERATING NETWORK EDGES
There are diverse types of multilayer networks [4,8]. One common type of multilayer network is one with only intralayer edges. There are also many multilayer networks with both intralayer and interlayer edges (e.g., multimodal transportation networks [95]), as well as ones with only interlayer edges (e.g., temporal networks with edge delays, such as departures and arrivals of flights between airports [5,96]). One can use our framework to generate all three types of examples, provided the underlying edge-generation model is a planted-partition network model.
Having generated a multilayer partition S with dependencies between induced partitions in different layers, the simplest way to generate edges is to use any single-layer plantedpartition network model (e.g., SBMs [41,42,48,56,62,97], models for spatially embedded networks [12,63], and so on) and to generate edges for each layer independently. This yields a multilayer network with only intralayer edges, such that any dependencies between different layers result only from dependencies between induced partitions in the different layers. For our numerical experiments in Sec. VI, we use a single-layer network model that is a slight variant (avoiding the creation of self-edges and multi-edges) of the degreecorrected SBM (DCSBM) benchmark from [42], where the term "DCSBM benchmark" designates the specific type of DCSBM that was used in the numerical experiments of [42].
One can also include dependencies between layers other than those that are induced by planted mesoscale structures. For example, one can introduce dependencies between the parameters of a single-layer planted-partition network model by (1) sampling them from a common probability distribution (e.g., to incorporate interlayer degree correlations [4,98] in a DCSBM) or by (2) introducing interlayer edge correlations, given a single-layer partition on each layer [99]. For temporal networks, one can also incorporate "burstiness" [2,3,100] in the interevent-time distribution of edges. In such a scenario, the probability for an edge to exist in a given layer depends not only on the induced partition in that layer, but also on the existence (and weights) of the edge in previous layers. For example, one can use a Hawkes process to specify the time points at which edges are active [101,102].
In Appendix D, we describe a multilayer generalization of the DCSBM of Ref. [42] that one can use to generate intralayer edges and/or interlayer edges. Our generalization gives a framework to formulate the parameters of a singlelayer DCSBM in a multilayer setting. One can use our multilayer DCSBM (M-DCSBM) framework to incorporate some of the features (e.g., degree correlations) that we described in the previous paragraph in a multilayer network with intralayer and/or interlayer edges.

VI. NUMERICAL EXAMPLES
In this section, we use our framework to construct benchmark models for multilayer community-detection methods and algorithms. We use the examples of interlayerdependency tensors from Figs. 2 and 3 to generate benchmark models. We also consider a variant of Fig. 2(b) in which we split the layers into groups, use uniform dependencies between layers in the same groups, and treat layers in different groups as independent of each other. These examples cover commonly studied temporal and multiplex dependencies, and they illustrate how one can generate benchmark multilayer networks with more than one aspect. In our numerical examples, we focus on a couple of popular objective functions for multilayer community detection, rather than investigating any given method or algorithm in detail.
We compare the behaviors of several variants of GenLouvain algorithms [103], which are similar to the the locally greedy Louvain computational heuristic [104], to optimize a multilayer modularity objective function [16,32] using the standard Newman-Girvan null model (which is a variant of a "configuration model" [105]). Modularity is an objective function that is popular for partitioning sets of nodes into communities that have a larger total internal edge weight than the expected total internal edge weight of the same sets in a "null network" [16], which is generated from some null model. Modularity maximization consists of finding a partition that maximizes this difference. For our numerical experiments, we use the generalization of modularity to multilayer networks of [32]. For multilayer modularity, the strengths of the interactions between different layers are governed by an interlayer-coupling tensor that controls the incentive for state nodes in different layers to be assigned to the same community. We use multilayer modularity with uniform interlayer coupling, so the strength of the interactions between different layers of a network depends on a layer-independent and node-independent interlayer coupling weight ω 0. We use diagonal and categorical (i.e., between all pairs of layers) interlayer coupling with weight ω for the multiplex examples in Sec. VI A. We use diagonal and ordinal interlayer coupling (so the coupling is between contiguous layers) with weight ω for the temporal examples in Sec. VI B. In Appendix E, we describe the Louvain algorithm and the variants of it (GENLOUVAIN and GENLOUVAINRAND) that we use in this paper.
We compare the results of multilayer modularity maximization with those of multilayer INFOMAP 6 [9], which uses an objective function called the "map equation" (which is not an equation) that is based on a discrete-time random walk and ideas from coding theory, to coarse-grain sets of nodes into communities [106]. In multilayer INFOMAP, one uses a probability r ∈ [0, 1] called the "relaxation rate" to control the relative frequency with which a random walker remains in the same layer or moves to other layers. (A random walker cannot change layers when r = 0.) The relaxation rate thus controls the interactions between different layers of a multilayer network. We allow the random walker to move to any other layer when r = 0 for the multiplex examples in Sec. VI A, and we allow the random walker to move only to contiguous layers for the temporal examples in Sec. VI B. In contrast to multilayer modularity (where ω = 0 yields single-layer modularity for each layer), computing singlelayer INFOMAP for each layer is not equivalent to choosing r = 0 in multilayer INFOMAP, because placing state nodes that correspond to the same physical node in the same community can contribute positively to the objective function even when r = 0 [9]. Consequently, we compute single-layer INFOMAP for each layer separately and reference it with the label "s" on the horizontal axis in our figures with results from INFOMAP.
In all experiments in this section, we generate a multilayer partition using our copying process in Sec. IV (see Algorithm 1 in Appendix B) and a multilayer network for a fixed planted partition using a slight variant of the DCSBM benchmark of Ref. [42]. We generate only intralayer edges, and our variant avoids creating self-edges and multi-edges (see Algorithm 3 in Appendix D). This multilayer plantedpartition network model is our M-DCSBM benchmark. This M-DCSBM benchmark produces multilayer networks that have only intralayer edges and whose connectivity patterns in different layers depend on each other. Following [42], we parametrize our M-DCSBM benchmark in terms of its distribution of expected degrees and a community-mixing parameter μ ∈ [0, 1] that controls the strength of the community structure in the sampled network edges. For μ = 0, all edges lie within communities; for μ = 1, edges are distributed independently of the communities, where the probability of observing an edge between two state nodes in the same layer depends only on the expected degrees of those two state nodes. We use a truncated power law for the distribution of expected degrees (see Appendix D).
The DCSBM benchmark has several features that make our M-DCSBM benchmark an interesting test for the ability of multilayer community-detection methods to aggregate information from multiple layers. The DCSBM benchmark imposes community structure as an expected feature of an ensemble of networks that it generates. Furthermore, the definition of the mixing parameter μ of the DCSBM benchmark ensures that the planted partition from the benchmark remains community-like for any μ < 1, as one is more likely to observe intracommunity edges and less likely to observe intercommunity edges than in a network that is generated from a single-block DCSBM with the same expected degree for each node. Consequently, given sufficiently many samples from the same DCSBM benchmark (i.e., all samples have the same planted partition and the same expected degrees), one should be able to identify the planted community structure for any μ < 1 (where the necessary number of samples goes to infinity as μ → 1).
We use NMI [94] to compare the performance of different community-detection algorithms. For each of our partitions, we compute the mean of the NMI between the partition induced on each layer by the output partition and that induced by the planted partition. That is, The quantity NMI is invariant under permutations of the meso-set labels within a layer. Consequently, NMI is wellsuited to comparing multilayer community-detection methods with single-layer community-detection methods. In particular, it allows us to test whether multilayer community-detection methods can exploit dependencies between layers of a multilayer network whenp 0 [see Eq. (6)]. In Appendix F, we show numerical experiments in which we compute NMI between multilayer partitions. We denote the NMI between two multilayer partitions by mNMI.
In all of our numerical experiments in Secs. VI A and VI B, we sample the benchmark networks in the following way. For each value ofp, we generate 10 sample partitions. For each sample partition and each value of μ, we generate 10 sample multilayer networks. This yields 100 benchmark instantiations for each pair (p, μ). We run each communitydetection algorithm 10 times on each instantiation. In Figs. 6-9, we show NMI between planted and recovered multilayer partitions averaged over sample partitions, sample networks, and algorithmic runs for each value of μ. Our results for different planted partitions tend to be similar to each other. The only exceptions are the results from using INFOMAP on temporal networks with p = 0.99 and p = 1, where we FIG. 6. Multiplex networks with uniform interlayer dependencies. We illustrate the effects of the interlayer coupling strength ω and relaxation rate r on the ability of different community-detection algorithms to recover planted partitions as a function of the community-mixing parameter μ in benchmark networks with uniform multiplex dependencies [see Fig. 2(b)]. Each multilayer network has 1000 nodes and 15 layers, and each node is present in all layers. All NMI values are means over 10 runs of the algorithms and 100 benchmark instantiations. (See the introduction of Sec. VI.) Each curve corresponds to the mean NMI values that we obtain for a given value of μ, and the shaded area around a curve indicates the minimum and maximum NMI values that we obtain with the 10 sample partitions for a given value ofp. observe large differences between our results for different partitions for certain values of μ. Figures 6-10 illustrate results for different choices of the interlayer-dependency tensor. In Figs. 6-9, rows correspond to different choices of community-detection algorithm, and columns correspond to different values of the copying probabilities (with the strength of interlayer dependencies increasing from left to right). All benchmark instantiations that we use for Figs. 6-10 are available online. 7

A. Multiplex examples
In this section, we consider two stylized examples of multiplex networks. Multiplex networks arise in a variety of different applications, including international relations and trade [107][108][109], social networks [110], and ecological networks [111]. In our multiplex examples, we consider simple dependency structures in which we expect multilayer 7 See https://dx.doi.org/10.5281/zenodo.3304059. community-detection methods to outperform single-layer methods by exploiting interlayer dependencies.
In Fig. 6, we consider multiplex networks with uniform dependencies between community structure in different layers. In Fig. 7, we consider multiplex networks with nonuniform dependencies between community structure in different layers. In both figures, we parametrize the amount of interlayer dependency in a network by the probabilityp [see Eq. (6)] that a state node copies its community assignment from an inneighbor in the interlayer-dependency network. All multilayer networks in these examples have n = 1000 physical nodes and l = 15 layers. Each node is present in every layer, so there are a total of 15 000 state nodes. In Fig. 2(b), we show the layer-coupled interlayer-dependency matrix that we use to generate the uniform multiplex networks. For the nonuniform multiplex networks, we split the layers into 3 groups of 5 layers each. We use uniform dependencies between layers in the same group, and layers in different groups are independent of each other. The resulting layer-coupled interlayerdependency matrix is block diagonal with diagonal blocks as in Fig. 2(b) and 0 entries in the off-diagonal blocks. We FIG. 7. Multiplex networks with nonuniform interlayer dependencies. We illustrate the effects of the interlayer coupling strength ω and relaxation rate r on the ability of different communitydetection algorithms to recover planted partitions as a function of the community-mixing parameter μ in a multiplex benchmark with nonuniform interlayer dependencies. Each multilayer network has 1000 nodes and 15 layers, and each node is present in all layers. The layer-coupled interlayer-dependency matrix is a block-diagonal matrix with diagonal blocks of size 5 × 5. Each diagonal block corresponds to a copy of the matrix in Fig. 2(b) with p =p/4. In each off-diagonal block, we set the value of p to p c = 0, thereby incorporating an abrupt change in community structure. All NMI values are means over 10 runs of the algorithms and 100 benchmark instantiations. Each curve corresponds to the mean NMI values that we obtain for a given value of μ, and the shaded area around a curve indicates the minimum and maximum NMI values that we obtain with the 10 sample partitions for a given value ofp. use a Dirichlet null distribution with n set = 10 and θ = 1 (see Appendix A) to specify expected community sizes and the M-DCSBM benchmark (see Appendix D) with η k = 2, k min = 3, and k max = 150 to generate intralayer edges. We perform 200 iterations of our update process (see Appendix B).
Comparing our results for GENLOUVAIN [see Figs. 6(a)-6(e) and 7(a)-7(c)] and GENLOUVAINRAND with "reiteration" and "postprocessing" 8 [see Figs. 6(f)-6(j) and 7(d)-7(f)], we see that the choice of optimization heuristic for a given objective function can significantly affect the quality of the 8 We describe reiteration and postprocessing in Appendix E. resulting output partitions. In particular, for GENLOUVAIN, we observe two distinct regimes and what appears to be a sharp transition between them. 9 For ω 1, we obtain partitions that are of similar quality to what we obtain by maximizing single-layer modularity for each layer; however, for ω 1, each partition that we obtain has identical induced partitions in each layer. We call the former the single-layer regime, and we call the latter the aggregate regime. Although partitions from the aggregate regime may be more similar to a planted partition than partitions from the single-layer regime whenp is sufficiently large [e.g., this occurs in Figs. 6(d) and 6(e)], we do not observe an interval of ω values between the two regimes in which GENLOUVAIN is better at recovering the planted partition than it is in the aggregate and single-layer regimes. By contrast, in most cases, GENLOUVAINRAND with reiteration and postprocessing identifies partitions that match the planted partition more closely than it does in either the single-layer or aggregate regimes. The exceptions are uniform multiplex networks withp = 0.5 [see Fig. 6 where GENLOUVAINRAND is unable to exploit the weak interlayer dependencies to outperform the partitions from the single-layer regime, andp = 1 [see Fig. 6(j)], where the partitions from the aggregate regime are always best. Most of the improvement in the results for GENLOUVAINRAND over those for GENLOUVAIN comes from reiteration. The additional randomization helps smooth out the transition (which tends to occur when ω ≈ 1) between the single-layer and aggregate regimes. Postprocessing only yields a minor improvement in the value of NMI . However, the effect of postprocessing is more pronounced when we compute mNMI between multilayer partitions (see Appendix F).
Multilayer INFOMAP has some problematic behavior in these benchmark experiments. Our results for multilayer INFOMAP are noticeably worse than those that we obtain with single-layer INFOMAP (corresponding to the label "s" on the horizontal axis) for networks with relatively weak community structure (specifically, for our computations with μ 0.4), unless the planted partition for different layers is very similar (specifically, unlessp 0.99). 10 Furthermore, the methods that are based on multilayer modularity outperform multilayer INFOMAP in our experiments for networks with weak community structure and similar layers (specifically, when both μ andp are close to 1). INFOMAP does not identify meaningful community structure in networks with μ 0.8 in any of our numerical examples, whereas GENLOUVAIN and GENLOUVAINRAND identify meaningful structure even for networks with very weak community structure (e.g., with μ = 0.9) for sufficiently large values of ω. 9 We observe some erratic behavior as the interlayer coupling ω approaches 1 from below. For values of ω near 1 but smaller than 1, the GENLOUVAIN algorithm has a tendency to place all state nodes into a single community. This observation is related to the transition behavior that was described in Ref. [16]. For values of ω that are above a certain threshold, only interlayer merges occur in the first phase of GENLOUVAIN (because of the locally greedy nature of the algorithm). 10 When we give inequalities like these for parameter values, we refer specifically to values at which we do numerical computations.

023100-14
Our results for nonuniform multiplex networks in Fig. 7 are similar to those that we showed for the uniform multiplex networks in Fig. 6. In particular, we see that GENLOUVAINRAND can exploit dependencies between community structure in different layers (although the improvement over single-layer modularity is less pronounced than in the uniform examples), whereas GENLOUVAIN cannot. The main difference is that for p = 1, the partitions from the aggregate regime for GENLOU-VAIN and GENLOUVAINRAND do not recover the planted partition, because the partitions induced by the planted partition in different layers are not identical. In principle, multilayer INFOMAP should have an advantage over methods that are based on maximizing multilayer modularity in this case study, as the former's objective function is designed to detect breaks in community structure, whereas multilayer modularity forces some persistence of community labels between any pair of layers. In practice, however, multilayer INFOMAP correctly identifies the planted community structure only when it is particularly strong (as we show in Appendix F, it correctly identifies the different groups of layers for our computations on networks with μ 0.3); and it is outperformed by singlelayer INFOMAP for our computations on networks with μ 0.4.
We suspect that at least some of the shortcomings of multilayer INFOMAP in these experiments are due to the use of a Louvain-type optimization heuristic, rather than from flaws in INFOMAP's objective function. As we have seen in our results with the heuristics GENLOUVAIN and GENLOU-VAINRAND, seemingly minor adjustments of an optimization heuristic can have large effects on the quality of the results. We are thus hopeful that similar adjustments can also improve the results for multilayer INFOMAP.

B. Temporal examples
Temporal networks arise in many different applications [2,3], such as the study of brain dynamics [17], financialasset correlations [16], and scientific citations [18]. When representing a temporal network as a multilayer network, one orders the layers in a causal way, such that structure in a particular layer depends directly only on structure in previous layers (and not on structure in future layers). To generate multilayer networks that have such structure, we use our sampling process for fully ordered multilayer networks and the interlayer-dependency tensor of Fig. 2(a).
We consider two stylized examples of temporal networks. In Fig. 8, we show results for temporal networks that have uniform dependencies between contiguous layers (and hence tend to evolve gradually). To generate these networks, we set p β = p for all layers in Fig. 2(a). In Fig. 9, we show results for temporal networks with change points. To generate these networks, we set p β = p for all layers except layers 25, 50, and 75, for which p β = p c = 0, resulting in abrupt changes in community structure. Each multilayer network in these two examples has n = 150 physical nodes and l = 100 layers. Each node is present in every layer, so there are a total of 15 000 state nodes. We use a Dirichlet null distribution to specify expected community sizes and set n set = 5 and θ = 1 (see Appendix A). For the temporal networks with change points in Fig. 9, we choose the supports of the null distributions so that communities after a change point have new labels. We use the M-DCSBM benchmark (see Appendix D) with η k = 2, k min = 3, and k max = 30 to generate intralayer edges.
Both multilayer-modularity-based algorithms (i.e., GEN-LOUVAIN and GENLOUVAINRAND) can exploit interlayer dependencies for these temporal benchmark networks and identify partitions with significantly larger NMI values than those that we obtain with single-layer modularity (i.e., with ω = 0). Typically, the peak of NMI seems to occur when 1 ω 4. When p < 1, one expects NMI to decrease for sufficiently large values of ω, as increasing ω further favors "persistence" [16] in the output partition that is not present in the multilayer planted partition. , we see that the difference in results between the two optimization heuristics for multilayer modularity is even more pronounced in these temporal examples than in the multiplex examples in Sec. VI A. In particular, GENLOUVAIN has an abrupt change in behavior near ω = 1; this is related to the transition behavior that was described in Ref. [16]. As we explained in Sec. VI A (where we also observed such a phenomenon), this transition occurs for the following reason: for values of ω that are above a certain threshold, only interlayer merges occur in the first phase of GENLOUVAIN. This abrupt transition no longer occurs with the additional randomization in GENLOUVAINRAND, which exhibits smooth behavior as a function of ω. However, unlike in our multiplex experiments, we observe benefits from the behavior of GENLOUVAIN in some situations, especially for networks with weak community structure (specifically, for μ 0.8). This phenomenon first becomes noticeable for networks with p = 0.95, and it becomes more pronounced for progressively larger values of p. Our results for temporal benchmark networks with change points (see Fig. 9) are mostly similar to those for temporal benchmark networks with uniform dependencies when we calculate NMI to compare planted and recovered partitions.  Fig. 2(a)]. Each multilayer network has 150 nodes and 100 layers, and each node is present in all layers. All NMI values are means over 10 runs of the algorithms and 100 benchmark instantiations. (See the introduction of Sec. VI.) Each curve corresponds to the mean NMI values that we obtain for a given value of μ, and the shaded area around a curve indicates the minimum and maximum NMI values that we obtain with the 10 sample partitions for a given value of p.
for each layer are not identical, even when p = 1. Consequently, unlike in Fig. 8, the results for p = 1 in Fig. 9 are qualitatively similar to those for smaller values of p. In principle, one might expect multilayer INFOMAP to have an advantage over multilayer-modularity-based methods in this case study, as the former's objective function is designed to detect abrupt changes in community structure, whereas maximizing multilayer modularity favors some persistence of community labels across layers. However, this theoretical advantage is not borne out in our computations. When comparing the planted and output partitions by calculating mNMI, we see that multilayer INFOMAP can correctly recover the change points only when both p = 1 and μ ∈ {0, 0.1} (see Appendix F).

C. Multiaspect examples
In this section, we illustrate the ability of our framework to generate multilayer networks with more than one aspect. A common type of multilayer network with two aspects is a multiplex network that changes over time (e.g., citations patterns that change over time [112]). In our framework, one can model such a situation by using an interlayer-dependency tensor like the one in Fig. 3. For the illustrative examples in this section, we use uniform multiplex and temporal dependencies. We parametrize the interlayer-dependency tensor as follows: where a ∈ [0, 1] controls the balance between multiplex and temporal dependencies andp controls the overall strength of the interlayer dependencies.
To recover communities, we optimize a multiaspect generalization, 023100-16 FIG. 9. Temporal networks with nonuniform interlayer dependencies. We illustrate the effects of the interlayer coupling strength ω and relaxation rate r on the ability of different communitydetection algorithms to recover planted partitions as a function of the community-mixing parameter μ in a temporal benchmark with nonuniform interlayer dependencies. Each multilayer network has 150 nodes and 100 layers, and each node is present in all layers. The layer-coupled interlayer-dependency matrix is the one in Fig. 2(a). We choose the values of p β in Fig. 2(a) such that we introduce three abrupt changes in community structure. Every 25th layer (specifically, for layers 25, 50, and 75), we set the value of p β to p c = 0, thereby introducing an abrupt change in community structure; we set all other values of p β to p. All NMI values are means over 10 runs of the algorithms and 100 benchmark instantiations. Each curve corresponds to the mean NMI values that we obtain for a given value of μ, and the shaded area around a curve indicates the minimum and maximum NMI values that we obtain with the 10 sample partitions for a given value of p.
of the multilayer modularity function that was introduced in Ref. [32].
We specify the interlayer-coupling tensor C (a multiaspect generalization of the "interslice coupling" of [32]), such that it reflects the planted interlayer dependencies, by writing where ω t denotes the coupling parameter for the temporal dependencies and ω m denotes the coupling parameter for the multiplex dependencies.
FIG. 10. Multiaspect networks with uniform interlayer dependencies. We illustrate the effects of the strengths of the temporal interlayer coupling ω t and the multiplex interlayer coupling ω m on the ability of multilayer modularity maximization to recover a planted partition in two-aspect multilayer networks that we generate using our framework. The first aspect is temporal (and thus ordered) and the second aspect is multiplex (and thus unordered). The heat maps show the mean ( NMI , 0 1) over all layers of the NMI between induced planted partitions and recovered community structure. A black cross indicates the maximum value of NMI . The dependence between structure in different layers increases from top to bottom, and the importance of the temporal aspect versus the multiplex aspect increases from left to right. We average our results for each pair of parameters over 10 runs of GENLOUVAINRAND with reiteration. The generated networks have n = 100 nodes and l = 100 layers with |L T | = |L M | = 10, so there are a total of 10 000 state nodes. We fix the community-mixing parameter of the M-DCSBM benchmark to be μ = 0.5.
In Fig. 10, we illustrate the effects of the two interlayer coupling parameters, ω t and ω m , on the extent to which multilayer modularity maximization can recover a planted partition. We use GENLOUVAINRAND with reiteration (but without postprocessing) to find a local optimum of Eq. (8). For all of the values ofp and a that we consider in Fig. 10, we see some improvement in the performance of multilayer modularity maximization with nonzero interlayer coupling compared with the case in which both ω t = 0 and ω m = 0. As the latter case corresponds to independently maximizing modularity on each layer of a network, this demonstrates that multilayer modularity is able to exploit the interlayer dependencies to better recover a planted partition. As expected, increasing a 023100-17 (which gives more weight to temporal dependencies) leads to a shift of the region of the parameter space with good planted-partition recovery to larger values of ω t and smaller values of ω m . For progressively largerp, we observe a small overall increase in the value of NMI , but its dependence on ω t and ω m remains similar.

VII. CONCLUSIONS AND DISCUSSION
We introduced a unifying and flexible framework for the construction of generative models for mesoscale structures in multilayer networks. The three most important features of our framework are the following: (1) it includes an explicitly parametrizable tensor P that controls interlayer-dependency structure; (2) it can generate an extremely general, diverse set of multilayer networks (including, e.g., temporal, multiplex, "multilevel" [113], and multiaspect multilayer networks with uniform or nonuniform dependencies between state nodes); and (3) it is modular, as the process of generating a partition is separate from the process of generating edges, enabling a user to first generate a partition and then use any planted-partition network model. Along with our paper, we provide publicly available code [14] that users can modify to readily incorporate different types of null distributions (see Sec. IV A), interlayer-dependency structures (see Sec. IV B), and plantedpartition network models (see Sec. V).
The ability to explicitly specify interlayer-dependency structure makes it possible for a user to control which layers depend directly on each other (by deciding which entries in the interlayer-dependency tensor are nonzero) and the strengths of such dependencies (by varying the values of the entries in the interlayer-dependency tensor). One can thereby generate multilayer networks with either a single aspect or multiple aspects (e.g., temporal and/or multiplex networks) and vary dependencies between layers from the extreme case in which induced partitions of a planted multilayer partition are the same across layers to the opposite extreme, in which induced partitions of a planted multilayer partition are generated independently for each layer from a null distribution for that layer. To the best of our knowledge, this level of generality is absent from existing generative models for mesoscale structures in multilayer networks, as those models tend to consider only networks with a single aspect (e.g., temporal or multiplex) or networks with limited interlayer-dependency structures (e.g., in which a planted multilayer partition has identical induced partitions in all layers).
We illustrated several examples of generative models that one can construct from our framework, with a focus on a few special cases of interest, rather than on trying to discuss as many situations as possible. We focused on community structure in our numerical experiments because it is a commonly studied mesoscale structure, but one can also use our framework to construct generative models of mesoscale structures other than community structure (e.g., core-periphery structure, bipartite structure, and so on) by taking advantage of our framework's ability to use any planted-partition network model. For our single-aspect examples, we assumed that interlayer dependencies exist either between all contiguous layers (a special case of temporal networks), between all layers (a special case of multiplex networks), or between all contiguous groups of layers (another special case of multiplex networks). For both our temporal and multiplex examples, we considered both uniform and nonuniform interlayer dependencies. We also combined some of these scenarios in an example with two aspects (a temporal aspect and a multiplex aspect). Our framework's flexibility allows us to construct generative models of multilayer networks with more realistic features. For example, for temporal networks, one can introduce dependencies between a layer and all layers that follow it (such that the interlayer-dependency tensor in Fig. 2(a) is an upper-triangular matrix with nonzero entries above the diagonal) to incorporate memory effects [114]. One can also consider interlayer dependencies that are not layer-coupled. For example, dependencies can be diagonal but nonuniform, or they can be nonuniform and exist only between sets of related nodes.
In Sec. VI, we considered the commonly studied case of a multilayer network with only intralayer edges and connectivity patterns in the different layers that depend on each other. We used a slight variant of the DCSBM benchmark from [42] to generate edges for each layer. However, other types of multilayer networks are also important [4], and one can readily combine our approach for generating multilayer partitions with different network generative models that capture various important features. For example, one can use an SBM to generate interlayer edges (e.g., using our M-DCSBM framework, which we discuss in Appendix D), or one can replace the degree-corrected SBM in Sec. V with any other planted-partition network model or other interesting models (e.g., other variants of SBMs [41,48] and models for networks whose structure is affected by space [12] or latent features [63]). In all of these examples, dependencies between connectivity patterns in different layers arise only from a planted multilayer partition. It is also possible to modify our network-generation process (see Sec. V) to incorporate additional dependencies between layers beyond those that are induced by a planted mesoscale structure. For example, one can introduce dependencies between a node's degree in different layers [98,115] or burstiness [100] in the intereventtime distributions of edges.
Our work has the potential for many useful and interesting extensions, and we highlight three of these. First, although we have given some illustrative numerical examples in Sec. VI, the area of benchmarking community-detection methods in multilayer networks is far from fully developed. Generative models are useful tools for understanding the behavior of community-detection methods in detail and thus for suggesting ways of improving heuristic algorithms without losing scalability. One can use our framework to construct benchmark models that provide a test bed for gaining insight into the advantages and shortcomings of communitydetection methods and algorithms (and, more generally, of meso-set-detection methods and algorithms). Importantly, we expect these benchmark models to be very informative for detecting potential artifacts of algorithms that can sometimes be masked in real-world applications. Second, a wellunderstood generative model can be a powerful tool for statistical inference (specifically, for inferring the structure of a multilayer network, rather than generating a multilayer network with a planted structure) [23]. For temporal networks,  closed forms for the joint distribution of meso-set assignments have been derived for models that are special cases of our framework [10,13]. These results may be useful for statistical inference. Additionally, it seems likely that it is possible to adapt Bayesian inference techniques (such as Gibbs sampling [116][117][118] and variational methods [119]) that have been developed for SBMs to derive methods both for inferring a multilayer partition and for inferring an interlayer-dependency tensor. A key advantage of the generality of the framework that we have developed in this paper is that it may be possible to frame many model-selection questions in terms of posterior estimation of an interlayer-dependency tensor. Finally, it is important to model interdependent data streams and not just fixed data sets. For example, for any fully ordered multilayer network, our generative model respects the causality of layers (e.g., temporal causality), and one can thus update a multilayer network with a new layer without the need to update any previous layers. It is critical to develop generative models that are readily adaptable to such situations, and our work in this paper is a step in this direction. Computational resources that we used for this research were supported in part by the National Science Foundation under Grant No. CNS-0521433; in part by Lilly Endowment, Inc., through its support for the Indiana University Pervasive Technology Institute; and in part by the Indiana METACyt Initiative. The Indiana METACyt Initiative at Indiana University Bloomington was also supported in part by Lilly Endowment, Inc. We thank Sang Hoon Lee, Peter Mucha (and his research group), Brooks Paige, and Roxana Pamfil for helpful comments.

APPENDIX A: EXAMPLE NULL DISTRIBUTIONS
In this appendix, we discuss parameter choices for categorical null distributions and give concrete examples that can be useful for modeling mesoscale structures in temporal and multiplex networks.
In Sec. IV A, we described the general form of a categorical null distribution: where p α s is the probability for a state node in layer α to be assigned to a meso-set s in the absence of interlayer dependencies, n set is the number of meso-sets in the multilayer partition, and n set s=1 p α s = 1 for each α ∈ L. The support G α of the null distribution P α 0 is G α = {s : P α 0 [s] = 0}. We say that a label s is active in a layer α if it is in the support of the null distribution P α 0 (i.e., if P α 0 [s] = 0), and we say that a label is inactive in layer α if it is in the complement of the support of P α 0 (i.e., if P α 0 [s] = 0). A natural choice is to sample p α from a Dirichlet distribution, which is the conjugate prior of a categorical distribution [90,91]. The Dirichlet distribution of order q has q parameters θ 1 , . . . , θ q . Its probability density function is where denotes the gamma function, x i ∈ (0, 1), and θ i > 0 for each i ∈ {1, . . . , q}. The case in which all θ i are equal is called a "symmetric Dirichlet distribution," which we parametrize by the common value θ (the so-called "concentration parameter") of the parameters and the number q of variables.
The concentration parameter θ determines the types of discrete probability distributions that one is likely to obtain from the symmetric Dirichlet distribution of order q = n set . For θ = 1, the symmetric Dirichlet distribution is the continuous uniform distribution over the space of all discrete probability distributions with n set states. As θ → ∞, the symmetric Dirichlet distribution becomes increasingly concentrated near the discrete uniform distribution, such that all entries in p α are approximately equal. As θ → 0, it becomes increasingly concentrated away from the uniform distribution, such that p α tends to have 1 (or a few) large entries and all other entries are close to 0. Consequently, to have very heterogeneous meso-set sizes, one wants θ ≈ 1. To have all meso-sets be of similar sizes, one uses a large value of θ . To have a few large meso-sets and many small meso-sets, one wants θ to be sufficiently smaller than 1. The value of n set also affects the amount of meso-set label overlap across layers in the absence of interlayer dependencies. For example, if p α is the same for all layers, then larger values of n set incentivize less label overlap across layers (because there are more possible labels for each layer) and smaller values of n set incentivize more label overlap across layers (because there are fewer possible labels for each layer).
In some situations-for example, when modeling the birth and death of communities in temporal networks or the appearance and disappearance of communities in multiplex networks-it is desirable to have meso-set labels that have a nonzero probability in Eq. (A1) only in some layers. For these situations, we suggest sampling the support of the distributions before sampling the probabilities p α . Given the supports for each layer, one then samples the corresponding probabilities from a symmetric Dirichlet distribution (or any other distribution over categorical distributions). That is, where G α = {s ∈ {1, . . . , n set } : P α 0 [s] = 0} is the complement of the support, with n set = max α∈L max(G α ).
A simple example of a process for modeling the birth and death of meso-sets is the following one. First, fix a number of meso-sets and a support (i.e., active community labels) for the first layer. One then sequentially initializes the supports for the other layers by removing each meso-set in the support of the previous layer with probability r d ∈ [0, 1] and adding a number, sampled from a Poisson distribution with rate r b ∈ [0, ∞), of new meso-sets (with new labels that are not active in any previous layer). In temporal networks, for example, this ensures that if a meso-set label is not in the support of a given layer after being in the support of previous layers, then the label is also not in the support of any subsequent layers. For this process, the expected number |G α | of mesosets in a layer approaches r b /r d as one iterates through the layers. Therefore, for the first layer, one should initialize the support of the null distribution with about r b /r d elements to avoid transients in the number of meso-sets. For this process, the lifetime of meso-sets follows a geometric distribution. The nature of the copying process that we use to introduce dependencies between induced partitions in different layers typically implies that meso-sets that have been removed from the support of the null distribution do not lose all of their members instantly, but instead shrink at a rate that depends on the values of the copying probabilities in the interlayerdependency tensor.
One can also allow labels to appear and disappear when examining multiplex partitions. For example, given a value for n set , one can generate the support for each layer by allowing each label s ∈ {1, . . . , n set } to be present with some probabilityq and absent with complementary probability 1 −q. This yields a sets of active and inactive meso-set labels for each layer. One can then sample the nonzero probabilities in p α that correspond to active labels from a Dirichlet distribution and set p α s to 0 for each inactive label s. Because multiplex partitions are unordered, there is no notion of one layer occurring after another one, so we do not need to ensure that an inactive label in a given layer is also inactive in "subsequent" layers.

APPENDIX B: PARTITION SAMPLING PROCESS
In this appendix, we provide a detailed description of the way in which we sample multilayer partitions, including a discussion of the convergence properties of our sampling process. As we mentioned in Sec. IV B, we assume that the multilayer partitions are generated by a copying process on the meso-set assignments of state nodes. Recall that the conditional probability distribution that governs the new meso-set assignment of a state node, given the current state of the copying process [see Eq. (4) in Sec. IV B], is where P 0 is a given set of independent layer-specific null distributions and P is a specified interlayer-dependency tensor. Recall as well that a multilayer network can have both ordered aspects and unordered aspects, where we denote the set of ordered aspects by O and the set of unordered aspects ALGORITHM 1. Pseudocode for generating multilayer partitions with an interlayer-dependency tensor P, a set P 0 of null distributions, and ordered aspects O. The parameter n upd is the expected number of updates for each state node.
Initialize partition tensor using null distributions end for for Loop over ordered aspects in lexicographic order Sample uniformly from order-equivalent layers Loop over nodes in some order Update (i, α) according to Eq. (B1) end for end for end for return S end function by U . Our goal is to sample multilayer partitions that are consistent with generation by Eq. (B1) while respecting the update order from its ordered aspects (if any are present). Recall that two layers, α and β, are order equivalent (which we denote by α ∼ O β) if α a = β a for all a ∈ O. Based on this equivalence relation, we obtain an ordered set L/∼ O of equivalence classes by inheriting the order from the ordered aspects (where, for definiteness, we consider lexicographic ordering [92] over ordered aspects). We want to simultaneously sample the meso-set assignments of state nodes in each class of order-equivalent layers; and we want to sequentially sample the meso-set assignments of state nodes in non-orderequivalent layers, conditional on the fixed assignments of state nodes in preceding layers (see Algorithm 1).

Scan order and compatibility
To sample the induced partitions for each class of orderequivalent layers, we use a technique that resembles Gibbs sampling [79]. We define a Markov chain, as in Gibbs sampling, on the space of multilayer partitions by updating the meso-set assignments of state nodes by sampling from Eq. (B1). In Gibbs sampling, one assumes that the conditional distributions that one uses to define such a Markov chain are compatible, in the sense that there exists a joint distribution that realizes these conditional distributions [80][81][82][83]. However, for a general interlayer-dependency tensor, there is no guarantee that the distributions that are defined by Eq. (B1) are compatible. Following [78], we use the term "pseudo-Gibbs sampling" to refer to Gibbs sampling from a set of potentially incompatible conditional distributions. We can define different Markov chains by changing the way in which we select which conditional distribution to apply at each step (i.e., which state node to update at each step). The order in which one cycles through conditional distributions is known as the scan order [85]. We use the term sampling Markov chain for the Markov chain that is defined by pseudo-Gibbs sampling with a specified scan order. Common scan orders are cycling through conditional distributions in a fixed order and sampling the update distribution uniformly at random from the set of conditional distributions. In pseudo-Gibbs sampling, the choice of scan order deserves careful attention. In such sampling, different scan orders correspond to sampling from different distributions [85]. (By contrast, in Gibbs sampling, the choice of scan order influences only the speed of convergence.) In particular, when cycling through conditional distributions in a fixed order, later conditional distributions have more influence than earlier ones on the stationary distributions of the sampling Markov chain [85]. Sampling the update distribution uniformly at random mitigates this problem, at the expense of introducing computational overhead.
We can improve on a purely random sampling strategy by exploiting the structure of the interlayer-dependency tensor. In particular, note that the conditional distributions for state nodes in the same layer are independent (as we assume that the interlayer-dependency tensor has no intralayer contributions) and hence compatible. When updating a set of state nodes from a given layer, we can update them in any order or even update them concurrently. To take advantage of this, we first sample a layer uniformly at random from the current class of order-equivalent layers and then update all state nodes in that layer. The sequence of layer updates defines a Markov chain on the space of multilayer partitions, and its convergence properties determine the level of effectiveness of our sampling algorithm.

Convergence guarantees
We now discuss two key properties, aperiodicity and ergodicity, that determine the convergence behavior of finite Markov chains.
First, note that sampling layers uniformly at random guarantees that we have an aperiodic sampling Markov chain. To see this, note that a sufficient condition for aperiodicity is that, for all states, there is a nonzero probability to transition from the state to itself. This clearly holds for our sampling Markov chain, as it is possible for two consecutive updates to update the same layer and [because the probability distributions in Eq. (B1) remain unchanged after the first update] there is a positive probability that the second update does not change the partition.
Second, note that ifp i,α < 1 for all state nodes (i, α) ∈ V M with α ∈ [α], where [α] is a class of order-equivalent layers, then the sampling Markov chain for that class is ergodic over a subspace of multilayer partitions that includes the support G [α] = α∈[α] G α of the null distributions. We have already shown that the sampling Markov chain is aperiodic, so we only need to verify that there exists a sample path from any arbitrary initial partition S 0 to any arbitrary multilayer partition S g ∈ G [α] in the support of the null distributions. Such a sample path clearly exists, as Therefore, one can achieve the desired transition using one update for each state node. Sample paths of a finite, aperiodic Markov chain converge in distribution to a stationary distribution of the Markov chain. However, different sample paths can converge to different distributions, and the eventual stationary distribution can depend both on the initial condition and on the transient behavior of the sample path. If a Markov chain is ergodic, its stationary distribution is unique, so it is independent of the initial condition and transient behavior.
After sufficiently many updates (i.e., by specifying a sufficiently large value for the expected number n upd of updates for each state node in Algorithm 1), we approximately sample the state of the Markov chains for each class of order-equivalent layers from their stationary distributions. For fully ordered multilayer networks (i.e., for multilayer networks with U = ∅), setting n upd = 1 is sufficient for convergence, because the equivalence classes include only a single layer in this case and subsequent updates are thus independent samples from the same distribution.
Sampling initial conditions in an unbiased way can be important for ensuring that one successfully explores the space of possible partitions, as the updating process is not necessarily ergodic over the support of the null distributions when (i,α)∈V M P j,β i,α = 1 for some state nodes. For example, ifp i,α = 1 for all (i, α) ∈ V M , then any partition in which all state nodes in each weakly connected component of the interlayer-dependency network have the same meso-set assignment is an absorbing state of the Markov chain.
However, provided the updating process has converged, any partition that one generates should still reflect the desired dependencies between induced partitions in different layers. Therefore, to circumvent the problem of nonunique stationary distributions when the updating process is not ergodic, we reinitialize the updating process by sampling from the null distribution for each partition that we want to generate. Reinitializing the initial partition from a fixed distribution for each sample always defines a unique joint distribution, which is a mixture of all possible stationary distributions when the updating process is not ergodic. When the updating process is ergodic, it is equivalent to sample partitions by reinitializing or to sample multiple partitions from a single long chain. Using a single long chain is usually more efficient when one needs many samples and it is not problematic for there to be some dependence between samples [120]. In our case, however, we usually need only a few samples with the same parameters, and ensuring independence tends to be more important than ensuring perfect convergence (provided the generated partitions have the desired interlayer-dependency structure). Using multiple chains thus has clear advantages even when the updating process is ergodic, and it is necessary to use multiple chains when it is not ergodic.

Convergence tests
It is difficult to determine whether a Markov chain has converged to a stationary distribution, mostly because it is difficult to distinguish the case of a slow-mixing chain becoming stuck in a particular part of state space from the case in which the chain has converged to a stationary distribution. FIG. 11. Autocorrelation of the sampling Markov chain, as measured by mNMI. We use the layer-coupled interlayer-dependency matrix in Fig. 2(b) to generate partitions and show the corresponding value ofp above each panel. On the horizontal axis, we show the value of n upd that we use to generate the first partition in the comparison. The lag is the number of additional updates per state node that we use to generate the second partition in the comparison. We average results over a sample of 100 independent chains with different initial partitions. The shaded area indicates one standard deviation above and below the mean. Note that the large amount of noise for small values of the lag is not a result of differences between chains in the sample; instead, it is an inherent feature of all chains in the sample.
There has been much work on trying to define convergence criteria for Markov chains [121], but none of the available approaches are entirely successful. In practice, one usually runs a Markov chain (or chains) for a predetermined number of steps. One checks manually for a few examples that the resulting chains exhibit behavior that is consistent with convergence by examining autocorrelations between samples of the same chain and cross correlations between samples of independent chains with the same initial state. When feasible, one can also check whether parts of different, independent chains or different parts of the same chain are consistent with sampling from the same distribution.
To estimate the number of updates that we need in the numerical experiments of Sec. VI A, we examine autocorrela- . We use the layer-coupled interlayer-dependency matrix in Fig. 2(b) to generate partitions and show the corresponding value ofp above each panel. On the horizontal axis, we show the value of n upd that we use to generate the first partition in the comparison. The lag is the number of additional updates per state node that we use to generate the second partition in the comparison. We average results over a sample of 100 independent chains with different initial partitions. The shaded area indicates one standard deviation above and below the mean.
tions of the Markov chain in two different ways. In Fig. 11, we calculate multilayer NMI (mNMI) between sampled partitions at different steps of the Markov chain. We are interested in the behavior of the mNMI both as a function of n upd (the expected number of updates per state node that we use to generate the first partition in the comparison) and the lag (the number of additional updates per state node that we use to generate the second partition in the comparison). Because the Markov chain converges to a stationary distribution, the value of mNMI for a given lag should become independent of n upd (so the curves in Fig. 11 should become flat). The results in Fig. 11 suggest that the number of updates that we need for convergence increases moderately withp, where n upd ≈ 50 is sufficient for convergence even when p = 1. In Fig. 12, we examine whether the dependency pattern  between induced partitions in different layers (we characterize such a pattern by calculating NMI between induced partitions in different layers) has converged using the mean absolute distance between the NMI matrices between induced partitions in a multilayer partition. (We showed examples of such NMI matrices in Figs. 4 and 5.) As with mNMI, the mean absolute distance d should become independent of n upd once the Markov chain has converged to a stationary distribution. This yields estimates for the minimum number of updates for convergence that are consistent with those from the results in Fig. 11. It is always safe to use a larger number of updates to sample partitions, and the results from Figs. 11 and 12 suggest that our choice of n upd = 200 is conservative. Figures 11 and 12 reveal additional important information about the behavior of our sampling process. We can estimate the mixing time of the sampling Markov chain at stationarity from Fig. 11 because the mixing time corresponds to the lag that is necessary for the mNMI to converge to 0. The mixing time at stationarity increases withp and becomes infinite aŝ p → 1 because the Markov chain converges to an absorbing state. From Fig. 12, we observe that, for all values ofp, the mean absolute distance d between NMI matrices converges to a value near 0 and becomes essentially independent of the lag. This corroborates our assumption that different partitions that we generate from our sampling process have similar interlayer-dependency structures.

APPENDIX C: A GENERATIVE MODEL FOR TEMPORAL NETWORKS
In this appendix, we examine the generative model for temporal networks that we considered in Sec. IV D. In this model, we assume uniform dependencies between contiguous layers [i.e., p β = p ∈ [0, 1] for all β ∈ {2, . . . , l} in Fig. 2(a)]. For this choice of interlayer-dependency tensor, Algorithm 1 simplifies to Algorithm 2. In this case (and, more generally, for any fully ordered multilayer network), convergence is automatic because we need only one iteration through the layers.
A first important feature of the generative model in Algorithm 2 is that it respects the arrow of time. In particular, meso-set assignments in a given layer depend only on meso-set assignments in the previous layer (e.g., the previous temporal snapshot) and on the null distributions P 0 . That is, for all s ∈ {1, . . . , n set } and all α ∈ {2, . . . , l}, we have where we recall that S| α is the n-node partition that is induced in layer α by the multilayer partition S. The value of p determines the relative importances of the previous layer and the null distribution. When p = 0, meso-set assignments in a given layer depend only on the null distribution of that layer [i.e., on the second term of the right-hand side of ALGORITHM 2. Pseudocode for generating partitions of a temporal network with uniform interlayer dependencies p between successive layers [i.e., p β = p for all β ∈ {2, . . . , l} in Fig. 2(a) Eq. (C1)]. When p = 1, meso-set assignments in a given layer are identical to those in the previous layer (and, by induction, to meso-set assignments in all previous layers).
Using Eq. (C1), we see that the marginal probability that a given state node has a specific meso-set assignment is Computing marginal probabilities can be useful for computing expected meso-set sizes.
We now highlight how the copying update (i.e., Step C) and the reallocation update (i.e., Step R) in Algorithm 2 govern the evolution of meso-set assignments between contiguous layers. Steps C and R deal with the movement of nodes by first removing some nodes ("subtraction") and then reallocating them ("addition"). In Step C, a meso-set with label s in layer α can lose a number of nodes that ranges from 0 to all of them. It can keep all of its nodes in layer α + 1 (i.e., S s | α+1 = S s | α ), lose some of its nodes (i.e., S s | α+1 ⊂ S s | α ), or disappear entirely (i.e., S s | α+1 = ∅ and S s | α = ∅). The null distribution in Step R is responsible for a meso-set gaining new nodes (i.e., S s | α+1 ⊂ S s | α ) or for the appearance of a new meso-set label (i.e., S s | α+1 = ∅ and S s | α = ∅). When defining the null distributions P 0 , it is necessary to consider the interplay between Step C and Step R.
To illustrate how the meso-set-assignment copying process and the null distribution in Algorithm 2 can interact with each other, we give the conditional probability that a label disappears in layer α and the conditional probability that a label appears in layer α. For all s ∈ {1, . . . , n set } and all α ∈ {2, . . . , l}, the conditional probability that a label s is not present in layer α is The expression (C2) depends only on our copying process and simplifies to (1 − p) |S s | α−1 | when P α 0 [S i,α = s] = 0 (i.e., when the probability of being assigned to label s is 0 using the null distribution of layer α). Furthermore, for progressively larger P α 0 [S i,α = s], the conditional probability that a label s is not present is progressively smaller.
For all s ∈ {1, . . . , n set } and all α ∈ {2, . . . , l}, the conditional probability that at least one new label appears in layer α is , the conditional probability that a new label appears depends only on our meso-setassignment copying process and is given by 1 − p n . Furthermore, larger values of S r | α−1 ∈S| α−1 P α 0 [S i,α = r] reduce the conditional probability that a new label appears in layer α.
All discussions thus far in this appendix hold for any choice of P 0 . We now give two examples to illustrate some features of the categorical null distribution from Sec. IV A. In particular, we focus on the effect of the support of a categorical null distribution on a sampled multilayer partition. (See Sec. IV A for definitions of some of the notation that we use below.) Recall that the support G α of a categorical null distribution P α 0 is G α = {s : p α s = 0}, where s ∈ {1, . . . , n set }. An important property of the supports of the null distributions for the interlayer dependencies in Fig. 2(a) is that overlap between G α and G α+1 (i.e., G α+1 ∩ G α = ∅) is a necessary condition for a meso-set in layer α to gain new members in layer α + 1.
Let c α denote the vector of expected sizes of the induced meso-sets in layer α (i.e., c α s = np α s ), and suppose that the probabilities p α are the same in each layer (i.e., p α = p for all α). The expected number of meso-set labels is then the same for each layer, and the expected number of state nodes with meso-set label s is also the same in each layer and is given by c α s . This choice produces a temporal network in which state nodes change meso-set labels across layers in a way that preserves both the expected number of induced meso-sets in a layer and the expected sizes of the induced meso-sets in a layer. Now suppose that we choose the p α values such that their supports do not overlap (i.e., G α ∩ G β = ∅ for all α = β). At each Step C in Algorithm 2, an existing meso-set label can then only lose members; and, in each layer except for the first, with probability 1 − p n , at least one new label will appear. For this case, one expects pc α s members of meso-set s in layer α to remain in meso-set s in layer α + 1 and (1 − p)c α s members of meso-set s in layer α to be assigned to new meso-sets (because labels are nonoverlapping) in layer α + 1. This choice thus produces multilayer partitions in which the expected number of new meso-set labels per layer is nonzero (unless p = 1) and the expected size of a given induced meso-set decreases in time.

APPENDIX D: M-DCSBM SAMPLING PROCEDURE
In this appendix, we discuss a generalization to multilayer networks of the degree-corrected SBM (DCSBM) from Ref. [42]. In other words, it is a multilayer DCSBM (M-DCSBM). The parameters of a general, directed M-DCSBM are a multilayer partition S (which determines the assignments of state nodes to meso-sets), a block tensor W (which determines the expected number of edges between meso-sets in different layers), and a set σ of state-node parameters (which determine the allocation of edges to state nodes in meso-sets). The probability of observing an edge (or the expected number of edges, if we allow multi-edges) from state node (i, α) to state node ( j, β) with meso-set assignments r = S i,α and s = S j,β in an M-DCSBM is where W s,β r,α is the expected number of edges from state nodes in layer α and meso-set S r to state nodes in layer β and mesoset S s , the quantity σ β i,α is the probability for an edge starting in meso-set S r in layer α and ending in layer β to be attached to state node (i, α) (note that the dependence on S r is implicit in σ β i,α ), and σ j,β α is the probability for an edge starting in layer α and ending in meso-set S s in layer β to be attached to state node ( j, β). (Note that the dependence on S s is implicit in σ j,β α .) For an undirected M-DCSBM, both the block tensor W and the state-node parameters σ are symmetric. That is, W r,α s,β = W s,β r,α and σ i,α β = σ β i,α . The above M-DCSBM can generate multilayer networks with arbitrary expected layer-specific in-degrees and outdegrees for each state node. (The DCSBM of [42] can generate single-layer networks with arbitrary expected degrees.) Given a multilayer network with adjacency tensor A, the layer-α-specific in-degree of state node ( j, β) is and the layer-β-specific out-degree of state node (i, α) is The layer-α-specific in-degree and out-degree of state node (i, α) are the "intralayer in-degree" and "intralayer outdegree" of (i, α). For an undirected multilayer network, layerβ-specific in-degrees and out-degrees are equal (i.e., k i,α β = k β i,α ). We refer to their common value as the "layer-β-specific degree" of a state node. For an ensemble of networks that are generated by an M-DCSBM, the associated means are and As we mentioned in Sec. V, we consider undirected multilayer networks with only intralayer edges for our experiments in Sec. VI. The block tensor W thus does not have any interlayer contributions (i.e., W s,β r,α = 0 if α = β). Furthermore, we can reduce the number of node parameters that we need to specify the M-DCSBM to a single parameter σ i,α = σ α i,α = σ i,α α for each state node (i, α).
We construct the M-DCSBM benchmark that we use in our numerical experiments by sampling the expected intralayer degrees e i,α = k α i,α of the state nodes from a truncated power law 11 with exponent η k , minimum cutoff k min , and maximum cutoff k max . We then construct the block tensor W and statenode parameters σ for the M-DCSBM from the sampled expected intralayer degrees e and the meso-set assignments S. Let be the expected degree of meso-set s in layer α; and let be the expected number of edges in layer α. Consequently, is the probability for an intralayer edge in layer α to be attached to state node (i, α), given that the edge is attached to a state node that is in layer α and is part of meso-set S i,α . The elements of the block tensor give, for r = s, the expected number of edges between state nodes in meso-set s in layer β and state nodes in meso-set r in layer α. For s = r and β = α, the block-tensor element W s,β r,α gives twice the expected number of edges. One way to think of the M-DCSBM benchmark is that we categorize each edge that we want to sample as an intra-meso-set edge with probability 1 − μ or as a "random edge" (i.e., an edge that can be either an intra-meso-set edge or an inter-meso-set edge) with probability μ. To sample an edge, we sample two state nodes (which we then join by an edge). We call these two state nodes the "end points" of the 11 A power law with cutoffs x min and x max and exponent τ is a continuous probability distribution with probability density function edge. We sample the two end points of an intra-meso-set edge with a frequency that is proportional to the expected degree of their associated state nodes, conditional on the end points being in the same meso-set. By contrast, we sample each end point of a random edge with a frequency that is proportional to the expected degree of its associated state node (without conditioning on meso-set assignments). We assume that the total number of edges in layer α is sampled from a Poisson distribution 12 with mean w α . Although our procedure for sampling edges describes a potential algorithm for sampling 12 The choice of a Poisson distribution for the number of edges ensures that this approach for sampling edges is approximately consistent with sampling edges independently from Bernoulli distributions with success probabilities (D1). The exact distribution of the number of edges is a Poisson-binomial distribution, from which it is difficult to sample. A Poisson-binomial distribution is well-approximated by a Poisson distribution if none of the individual edge probabilities are too large. networks from the M-DCSBM benchmark, it usually is more efficient to sample edges separately for each pair of meso-sets.
In Algorithm 3, we show pseudocode for the specific instance of the M-DCSBM benchmark that we use to sample intralayer network edges (independently for each layer) for a given multilayer partition in our numerical experiments of Sec. VI. The only difference between Algorithm 3 and the sampling algorithm of [42] is that we use rejection sampling to avoid creating self-edges and multi-edges. (In other words, if we sample an edge that has already been sampled or that is a self-edge, we do not include it in the multilayer network; instead, we resample.) Rejection sampling is efficient as long as all blocks of a network's adjacency tensor are sufficiently sparse, such that the probability of generating multi-edges is small. For dense blocks, we instead sample edges from independent Bernoulli distributions with success probabilities (D1). This algorithm for sampling networks from an M-DCSBM is very efficient; it scales linearly with the number of edges in a network.
In the above discussion, we defined the parameters in Eq. (D1) to generate intralayer edges in a multilayer network by specifying expected layer-specific in-degrees and expected layer-specific out-degrees. That is, the expected layerβ-specific in-degree and expected layer-β-specific out-degree of state node (i, α) are both 0 if β = α. We can also extend the M-DCSBM benchmark to generate interlayer edges. For a directed multilayer network with interlayer edges, we sample expected layer-β-specific in-degrees e i,α β and expected layer-β-specific out-degrees e β i,α for each state node (i, α) and layer β from appropriate distributions. Given expected layer-specific degrees e and a multilayer partition S, we then construct the block tensor W and state-node parameters σ of the M-DCSBM analogously to the special case on which we have been focusing. Let For directed networks, the expected layer-specific in-degrees and expected layer-specific out-degrees that are generated by this model do not correspond exactly to the values from the input parameters e, except when μ = 1.

APPENDIX E: THE LOUVAIN ALGORITHM AND ITS VARIANTS
In this appendix, we describe the Louvain communitydetection algorithm and the variants of the Louvain algorithm that we used in our numerical experiments in Sec. VI.
The Louvain algorithm [104] for maximizing (single-layer or multilayer) modularity proceeds in two phases, which are repeated iteratively. Starting from an initial partition of a network, one considers the state nodes one by one (in some order) and places each state node in a set that results in the largest increase of modularity. (If there is no move that improves modularity, a state node keeps the same assignment.) One repeats this first phase of the algorithm until reaching a local maximum. In the second phase of the Louvain algorithm, one obtains a new modularity matrix by aggregating the sets of state nodes that one obtains after the convergence of the first phase. One then applies the algorithm's first phase to the new modularity matrix and iterates both phases until one converges to a local maximum. The two Louvain-type algorithms that we use in Sec. VI differ in how they select which moves to make. The first is GENLOUVAIN, which always chooses a move that maximally increases modularity; the second is GENLOUVAINRAND, which chooses modularityincreasing moves at random, such that the probability of a particular move is proportional to the resulting increase in the modularity objective function. (The latter is a variant of the algorithm "LouvainRand" in Ref. [16]; in that algorithm, one chooses modularity-increasing moves uniformly at random.) We use "reiteration" and "postprocessing" to improve the output of GENLOUVAINRAND. Reiteration entails taking an output partition of GENLOUVAINRAND as an initial partition and reapplying the algorithm until the output partition no longer changes. Postprocessing entails using the Hungarian algorithm [122] to optimize the "persistence" [16] (i.e., the number of interlayer edges for which both end points have the same community assignment) of each output partition without changing the induced partitions in each layer. All Louvain variants and related functionality that we use in this paper are available at [103].

APPENDIX F: MULTILAYER NMI EXPERIMENTS
In this appendix, we calculate multilayer NMI (mNMI) between output multilayer partitions and corresponding planted multilayer partitions for the numerical experiments in Figs. 6-9. We include only a subset of the values of p andp that we considered in Figs. 6 and 8, as these suffice for the purpose of this appendix. One needs to be cautious when interpreting values of mNMI; they depend on how one determines correspondences between communities in different layers, and these correspondences can be ambiguous. For example, in a temporal network, if a community in one layer splits into multiple communities in the next layer, there are multiple plausible ways to define labels across layers (e.g., all new communities get new labels or the largest new community keeps the original label). Moreover, different objective functions reward community-assignment labeling across layers in different ways, and mNMI conflates this issue with how well an objective function recovers structure within layers. In 023100-26 FIG. 13. Effects of the interlayer coupling strength ω and relaxation rate r on the ability of different community-detection algorithms to recover planted partitions as a function of the community-mixing parameter μ in a uniform multiplex benchmark [see Fig. 2(b)]. Each multilayer network has 1000 nodes and 15 layers, and each node is present in all layers. (We use the shorthand "NO PP" for experiments with GENLOUVAINRAND that do not include postprocessing.) particular, a small value of mNMI can hide an optimal value of NMI (as NMI = 1 is necessary, but not sufficient, to guarantee that mNMI = 1).
In the numerical experiments of Sec. VI, we observed that postprocessing does not have a major effect on NMI . This suggests that any misassignment across layers that underemphasizes the "persistence" of a multilayer partition (which postprocessing tries to mitigate [16]) does not affect an algorithm's ability to identify structure within layers. (Recall that postprocessing relabels an assignment without changing the induced partitions.) We include figures without postprocessing (which we indicate with the shorthand "NO PP") in FIG. 14. Effects of the interlayer coupling strength ω and relaxation rate r on the ability of different community-detection algorithms to recover planted partitions as a function of the communitymixing parameter μ in a multiplex benchmark with nonuniform interlayer dependencies. Each multilayer network has 1000 nodes and 15 layers, and each node is present in all layers. The layercoupled interlayer-dependency matrix is a block-diagonal matrix with diagonal blocks of size 5 × 5. Each diagonal block corresponds to a copy of the matrix in Fig. 2(b) with p =p/4. In each offdiagonal block, we set the value of p to p c = 0, thereby incorporating an abrupt change in community structure. We define the null distributions in these experiments so that there is no overlap in community labels between different groups of layers. (We use the shorthand "NO PP" for experiments with GENLOUVAINRAND that do not include postprocessing.) the experiments in this appendix, as we expect the effect of postprocessing to be more noticeable for mNMI than for NMI (because mNMI accounts for meso-set assignments both within and across layers).
In Figs. 13 and 14, we use the same benchmark networks as in the multiplex examples in Figs. 6 and 7, respectively. FIG. 15. Effects of the interlayer coupling strength ω and relaxation rate r on the ability of different community-detection algorithms to recover planted partitions as a function of the communitymixing parameter μ in a temporal benchmark with uniform interlayer dependencies [i.e., p β = p ∈ [0, 1] for all β ∈ {2, . . . , l} in Fig. 2(a)]. Each multilayer network has 150 nodes and 100 layers, and each node is present in all layers. (We use the shorthand "NO PP" for experiments with GENLOUVAINRAND that do not include postprocessing.) In Figs. 15 and 16, we use the same benchmark networks as in the temporal examples of Figs. 8 and 9, respectively. All mNMI values are means over 10 runs of the algorithms and 100 benchmark instantiations. (See the introduction of Sec. VI for more details on how we generate these instantiations.) Each curve in each figure corresponds to the mean mNMI values that we obtain for a given value of the community-mixing parameter μ, and the shaded area around a curve corresponds to the minimum and maximum mNMI values that we obtain with the 10 sample partitions for a given value of p orp.
For experiments with nonuniform interlayer dependencies in Figs. 14 and 16, it seems that INFOMAP is better at detecting FIG. 16. Effects of the interlayer coupling strength ω and relaxation rate r on the ability of different community-detection algorithms to recover planted partitions as a function of the communitymixing parameter μ in a temporal benchmark with nonuniform interlayer dependencies. Each multilayer network has 150 nodes and 100 layers, and each node is present in all layers. The layer-coupled interlayer-dependency matrix is given by Fig. 2(a). We choose the values of p β in Fig. 2(a) such that we introduce three abrupt changes in community structure. Every 25th layer (specifically, for layers 25, 50, and 75), we set the value of p β to p c = 0 (thereby incorporating an abrupt change in community structure); we set all other values of p β to p. We define the null distributions in these experiments so that there is no overlap in community labels between different groups of layers. (We use the shorthand "NO PP" for experiments with GENLOUVAINRAND that do not include postprocessing.) abrupt differences (so-called "breaks," such as in the form of change points in temporal examples) in community structure than both GENLOUVAIN and GENLOUVAINRAND, especially in the multiplex examples (where this manifests for all examined values ofp). However, multilayer INFOMAP correctly identifies the planted community structure only when it is 023100-28 particularly strong (specifically, for our computations with μ 0.3). Recall from our experiments in Figs. 7 and 9 that both GENLOUVAIN and GENLOUVAINRAND are better than INFOMAP at detecting an induced partition in a layer (especially when the planted community structure is somewhat weak). The reason that multilayer modularity maximization does not perfectly recover breaks in community structure of a planted multilayer partition is because multilayer modularity maximization with uniform ordinal or categorical coupling (see the introduction of Sec. VI) incentivizes "persistence" of community labels between layers (and nonzero persistence is achievable even when one independently samples induced partitions in different layers) [16]. However, by design, our benchmark networks have no persistence whenever there is a break in community structure because of how we defined the support of the null distributions. Therefore, we do not expect multilayer modularity maximization to recover the planted labels between layers when p c = 0; this manifests as a de facto maximum in the mNMI between an output partition and its corresponding planted partition.