diff --git a/docs/source/api/generators.rst b/docs/source/api/generators.rst
index 2af796e3..2db05123 100644
--- a/docs/source/api/generators.rst
+++ b/docs/source/api/generators.rst
@@ -20,4 +20,5 @@ generators package
    ~xgi.generators.uniform
    ~xgi.generators.simplicial_complexes
    ~xgi.generators.randomizing
+   ~xgi.generators.growth
    
\ No newline at end of file
diff --git a/docs/source/api/generators/xgi.generators.growth.rst b/docs/source/api/generators/xgi.generators.growth.rst
new file mode 100644
index 00000000..e8743f11
--- /dev/null
+++ b/docs/source/api/generators/xgi.generators.growth.rst
@@ -0,0 +1,12 @@
+xgi.generators.growth
+=====================
+
+.. currentmodule:: xgi.generators.growth
+
+.. automodule:: xgi.generators.growth
+   
+   .. rubric:: Functions
+   
+   .. autofunction:: apollonian_complex
+   .. autofunction:: network_geometry_flavor
+   .. autofunction:: pseudofractal_simplicial_complex
\ No newline at end of file
diff --git a/tests/generators/test_growth.py b/tests/generators/test_growth.py
new file mode 100644
index 00000000..c904fbc0
--- /dev/null
+++ b/tests/generators/test_growth.py
@@ -0,0 +1,184 @@
+import pytest
+
+import xgi
+from xgi.exception import XGIError
+
+
+def test_pseudofractal():
+
+    # initial
+    S0 = xgi.pseudofractal_simplicial_complex(order=2, n_iter=0)
+    triangles = set(S0.edges.filterby("order", 2).members(dtype=list))
+
+    assert isinstance(S0, xgi.SimplicialComplex)
+    assert S0.num_nodes == 3
+    assert S0.num_edges == 4
+    assert triangles == {frozenset({0, 1, 2})}
+
+    # first iteration
+    S1 = xgi.pseudofractal_simplicial_complex(order=2, n_iter=1)
+    triangles = set(S1.edges.filterby("order", 2).members(dtype=list))
+
+    assert isinstance(S1, xgi.SimplicialComplex)
+    assert S1.num_nodes == 6
+    assert xgi.num_edges_order(S1, d=2) == 4
+    assert triangles == {
+        frozenset({1, 2, 5}),
+        frozenset({0, 2, 4}),
+        frozenset({0, 1, 3}),
+        frozenset({0, 1, 2}),
+    }
+
+    # second iteration
+    S2 = xgi.pseudofractal_simplicial_complex(order=2, n_iter=2)
+    triangles = set(S2.edges.filterby("order", 2).members(dtype=list))
+
+    assert isinstance(S2, xgi.SimplicialComplex)
+    assert S2.num_nodes == 15
+    assert xgi.num_edges_order(S2, d=2) == 13
+    assert triangles == {
+        frozenset({0, 4, 10}),
+        frozenset({0, 2, 7}),
+        frozenset({0, 1, 2}),
+        frozenset({1, 5, 11}),
+        frozenset({1, 2, 5}),
+        frozenset({0, 2, 4}),
+        frozenset({0, 3, 12}),
+        frozenset({1, 3, 14}),
+        frozenset({2, 4, 9}),
+        frozenset({1, 2, 8}),
+        frozenset({0, 1, 6}),
+        frozenset({2, 5, 13}),
+        frozenset({0, 1, 3}),
+    }
+
+
+def test_apollonian():
+
+    # initial
+    S0 = xgi.apollonian_complex(order=2, n_iter=0)
+    triangles = set(S0.edges.filterby("order", 2).members(dtype=list))
+
+    assert isinstance(S0, xgi.SimplicialComplex)
+    assert S0.num_nodes == 3
+    assert S0.num_edges == 4
+    assert triangles == {frozenset({0, 1, 2})}
+
+    # first iteration
+    S1 = xgi.apollonian_complex(order=2, n_iter=1)
+    triangles = set(S1.edges.filterby("order", 2).members(dtype=list))
+
+    assert isinstance(S1, xgi.SimplicialComplex)
+    assert S1.num_nodes == 6
+    assert xgi.num_edges_order(S1, d=2) == 4
+    assert triangles == {
+        frozenset({1, 2, 5}),
+        frozenset({0, 2, 4}),
+        frozenset({0, 1, 3}),
+        frozenset({0, 1, 2}),
+    }
+
+    # second iteration
+    S2 = xgi.apollonian_complex(order=2, n_iter=2)
+    triangles = set(S2.edges.filterby("order", 2).members(dtype=list))
+
+    assert isinstance(S2, xgi.SimplicialComplex)
+    assert S2.num_nodes == 12
+    assert xgi.num_edges_order(S2, d=2) == 10
+    assert triangles == {
+        frozenset({2, 4, 6}),
+        frozenset({1, 5, 8}),
+        frozenset({0, 3, 9}),
+        frozenset({0, 1, 2}),
+        frozenset({1, 2, 5}),
+        frozenset({0, 2, 4}),
+        frozenset({2, 5, 10}),
+        frozenset({1, 3, 11}),
+        frozenset({0, 1, 3}),
+        frozenset({0, 4, 7}),
+    }
+
+
+@pytest.mark.parametrize(
+    "generator",
+    [xgi.pseudofractal_simplicial_complex, xgi.apollonian_complex],
+)
+def test_order_too_small(generator):
+    for order in [0, 1]:
+        with pytest.raises(ValueError):
+            generator(order=order, n_iter=1)
+
+
+def test_network_geometry_flavor():
+    order, n_iter = 2, 20
+    for s in [-1, 0, 1]:
+        S = xgi.network_geometry_flavor(order, s, beta=1, n_iter=n_iter, seed=1)
+        assert isinstance(S, xgi.SimplicialComplex)
+        # one new node and one new d-simplex per iteration
+        assert set(S.nodes) == set(range(order + 1 + n_iter))
+        assert len(S.edges.filterby("order", order)) == n_iter + 1
+
+    # s = -1: each (d-1)-simplex is shared by at most two d-simplices
+    S = xgi.network_geometry_flavor(order, -1, beta=1, n_iter=50, seed=1)
+    triangles = S.edges.filterby("order", order).members(dtype=list)
+    face_counts = {}
+    for face in xgi.subfaces(triangles, order=order - 1):
+        face = frozenset(face)
+        face_counts[face] = face_counts.get(face, 0) + 1
+    assert max(face_counts.values()) <= 2
+
+    # order 1 grows a tree
+    S = xgi.network_geometry_flavor(1, 0, beta=1, n_iter=10, seed=1)
+    assert S.num_nodes == 12
+    assert len(S.edges.filterby("order", 1)) == 11
+
+    # large beta does not underflow
+    S = xgi.network_geometry_flavor(order, 0, beta=1000, n_iter=10, seed=1)
+    assert S.num_nodes == order + 1 + 10
+
+
+def test_network_geometry_flavor_seed():
+    S1 = xgi.network_geometry_flavor(2, 1, beta=1, n_iter=30, seed=42)
+    S2 = xgi.network_geometry_flavor(2, 1, beta=1, n_iter=30, seed=42)
+    S3 = xgi.network_geometry_flavor(2, 1, beta=1, n_iter=30, seed=43)
+    assert set(S1.edges.members(dtype=list)) == set(S2.edges.members(dtype=list))
+    assert set(S1.edges.members(dtype=list)) != set(S3.edges.members(dtype=list))
+
+    # custom energy distribution
+    S = xgi.network_geometry_flavor(
+        2,
+        0,
+        beta=1,
+        n_iter=10,
+        energy_distribution=lambda rng: rng.exponential(1.0),
+        seed=1,
+    )
+    assert S.num_nodes == 13
+
+    # energy distribution with an rng argument is reproducible with seed
+    def energy(rng):
+        return rng.exponential(1.0)
+
+    S1 = xgi.network_geometry_flavor(2, 1, 1, 30, energy_distribution=energy, seed=7)
+    S2 = xgi.network_geometry_flavor(2, 1, 1, 30, energy_distribution=energy, seed=7)
+    assert set(S1.edges.members(dtype=list)) == set(S2.edges.members(dtype=list))
+
+    # energy distribution without arguments also works
+    calls = []
+
+    def energy_no_arg():
+        calls.append(1)
+        return 1.0
+
+    S = xgi.network_geometry_flavor(
+        2, 0, beta=1, n_iter=10, energy_distribution=energy_no_arg, seed=1
+    )
+    assert S.num_nodes == 13
+    assert len(calls) == 13  # one energy per node
+
+
+def test_network_geometry_flavor_errors():
+    with pytest.raises(ValueError):
+        xgi.network_geometry_flavor(0, 0, beta=1, n_iter=5)
+    with pytest.raises(ValueError):
+        xgi.network_geometry_flavor(2, 2, beta=1, n_iter=5)
diff --git a/xgi/generators/__init__.py b/xgi/generators/__init__.py
index 77df5860..7fbf7ef0 100644
--- a/xgi/generators/__init__.py
+++ b/xgi/generators/__init__.py
@@ -1,5 +1,6 @@
 from . import (
     classic,
+    growth,
     lattice,
     random,
     randomizing,
@@ -8,6 +9,7 @@
     uniform,
 )
 from .classic import *
+from .growth import *
 from .lattice import *
 from .random import *
 from .randomizing import *
@@ -43,4 +45,7 @@
     "uniform_erdos_renyi_hypergraph",
     "uniform_HSBM",
     "uniform_HPPM",
+    "pseudofractal_simplicial_complex",
+    "apollonian_complex",
+    "network_geometry_flavor",
 ]
diff --git a/xgi/generators/growth.py b/xgi/generators/growth.py
new file mode 100644
index 00000000..227b1291
--- /dev/null
+++ b/xgi/generators/growth.py
@@ -0,0 +1,333 @@
+"""Generate growing simplicial complexes."""
+
+from inspect import Parameter, signature
+
+import numpy as np
+
+from ..core import SimplicialComplex
+from ..utils.utilities import subfaces
+
+__all__ = [
+    "pseudofractal_simplicial_complex",
+    "apollonian_complex",
+    "network_geometry_flavor",
+]
+
+
+def pseudofractal_simplicial_complex(order, n_iter):
+    r"""
+    Generate the pseudofractal simplicial complex of order `order`.
+
+    Starting with a single d-simplex (d = `order`), at each iteration a new
+    node is attached to every existing (d-1)-simplex, forming a new d-simplex
+    (with all its subfaces). This process is deterministic.
+
+    Parameters
+    ----------
+        order : int
+            The order of the simplices to add (e.g., 2 for triangles, 3 for tetrahedra, etc.).
+        n_iter : int
+            The number of iterations to generate simplices.
+
+    Returns
+    -------
+        S : xgi.SimplicialComplex
+            Generated simplicial complex
+
+    Raises
+    ------
+    ValueError
+        If `order` is less than 2.
+
+    Notes
+    -----
+    The size grows exponentially with `n_iter`. After :math:`t` iterations, the
+    complex has :math:`((d+1)^{t+1} - 1)/d` simplices of order :math:`d`, i.e.,
+    :math:`\mathcal{O}((d+1)^t)`, and :math:`d` more nodes than that.
+    For example, ``order=2, n_iter=10`` gives 88,573 triangles and
+    ``order=3, n_iter=10`` gives 1,398,101 tetrahedra.
+
+    See also
+    --------
+    apollonian_complex
+    network_geometry_flavor
+
+    References
+    ----------
+    Nurisso, M., Morandini, M., Lucas, M., Vaccarino, F., Gili, T., & Petri, G. (2024).
+    "Higher-order Laplacian Renormalization."
+    arXiv preprint arXiv:2401.11298.
+    https://arxiv.org/abs/2401.11298
+    """
+
+    if order < 2:
+        raise ValueError("order must be at least 2.")
+
+    S = SimplicialComplex()
+
+    # initialize the first d-simplex
+    first_simplex = tuple(range(order + 1))
+    S.add_simplex(first_simplex)
+
+    # generate simplices iteratively
+    for it in range(1, n_iter + 1):
+        # Find all (order - 1)-simplices present in the complex
+        nodes = S.nodes
+        subfaces = S.edges.filterby("order", order - 1).members(dtype=list)
+        max_index = max(nodes)
+        new_simplices = []
+
+        for subface in subfaces:
+            # create a new simplex by adding the new vertex to the existing d-simplex
+            max_index += 1  # new vertex index
+
+            new_simplex = (*subface, max_index)
+            new_simplices.append(new_simplex)
+
+        S.add_simplices_from(new_simplices)
+
+    return S
+
+
+def apollonian_complex(order, n_iter):
+    r"""
+    Generate the apollonian complex of order `order`.
+
+    Starting with a single d-simplex (d = `order`), at each iteration a new
+    node is attached to every (d-1)-simplex that contains a node added in the
+    previous iteration, forming a new d-simplex (with all its subfaces).
+    This process is deterministic.
+
+    Parameters
+    ----------
+        order : int
+            The order of the simplices to add (e.g., 2 for triangles, 3 for tetrahedra, etc.).
+        n_iter : int
+            The maximum iteration to generate simplices.
+
+    Returns
+    -------
+        S : xgi.SimplicialComplex
+            Generated simplicial complex
+
+    Raises
+    ------
+    ValueError
+        If `order` is less than 2.
+
+    Notes
+    -----
+    The size grows exponentially with `n_iter`, but more slowly than for
+    :func:`pseudofractal_simplicial_complex`. After :math:`t` iterations, the
+    complex has :math:`1 + (d+1)(d^t - 1)/(d-1)` simplices of order :math:`d`,
+    i.e., :math:`\mathcal{O}(d^t)`, and :math:`d` more nodes than that.
+    For example, ``order=2, n_iter=10`` gives 3,070 triangles and
+    ``order=3, n_iter=10`` gives 118,097 tetrahedra.
+
+    See also
+    --------
+    pseudofractal_simplicial_complex
+    network_geometry_flavor
+
+    References
+    ----------
+    Nurisso, M., Morandini, M., Lucas, M., Vaccarino, F., Gili, T., & Petri, G. (2024).
+    "Higher-order Laplacian Renormalization."
+    arXiv preprint arXiv:2401.11298.
+    https://arxiv.org/abs/2401.11298
+    """
+
+    if order < 2:
+        raise ValueError("order must be at least 2.")
+
+    S = SimplicialComplex()
+
+    # initialize the first d-simplex
+    first_simplex = tuple(range(order + 1))
+    S.add_simplex(first_simplex)
+
+    new_simplices = [first_simplex]
+    new_indices = list(first_simplex)
+
+    # generate simplices iteratively
+    for it in range(1, n_iter + 1):
+        # find all (order - 1)-simplices present in the complex
+        nodes = S.nodes
+        subfaces_previous_iter = S.edges.filterby("order", order - 1).members(
+            dtype=list
+        )
+
+        # keep only those attached to new nodes
+        subfaces_previous_iter = [
+            subface
+            for subface in subfaces_previous_iter
+            if any(new_index in subface for new_index in new_indices)
+        ]
+
+        max_index = max(nodes)
+        new_simplices = []
+        new_indices = []
+
+        for subface in subfaces_previous_iter:
+            # create a new simplex by adding the new vertex to the existing d-simplex
+            max_index += 1  # New vertex index
+
+            new_simplex = (*subface, max_index)
+            new_simplices.append(new_simplex)
+            new_indices.append(max_index)
+
+        S.add_simplices_from(new_simplices)
+
+    return S
+
+
+def network_geometry_flavor(
+    order, s, beta, n_iter, energy_distribution=None, seed=None
+):
+    r"""
+    Generate a Network Geometry with Flavor (NGF) simplicial complex.
+
+    The model grows a d-dimensional simplicial complex (d = `order`), starting
+    from a single d-simplex. At each iteration, a new node is attached to an
+    existing (d-1)-simplex :math:`\alpha`, forming a new d-simplex (with all
+    its subfaces). Which (d-1)-simplex is chosen depends on the energies of
+    its nodes and on the flavor `s`, which sets whether (d-1)-simplices
+    already shared by several d-simplices are forbidden (``s=-1``), neutral
+    (``s=0``), or favored (``s=1``).
+
+    Each node :math:`i` gets a fixed energy :math:`\epsilon_i`, drawn from
+    `energy_distribution` when it is added. The energy of a (d-1)-simplex
+    is the sum of the energies of its nodes,
+    :math:`\epsilon_\alpha = \sum_{i \in \alpha} \epsilon_i`. At each
+    iteration, :math:`\alpha` is chosen with probability
+
+    .. math::
+
+        \Pi_\alpha \propto e^{-\beta \epsilon_\alpha} (1 + s\, n_\alpha),
+
+    where :math:`n_\alpha` is the number of d-simplices containing
+    :math:`\alpha`, minus one. There is no energy threshold: lower-energy
+    (d-1)-simplices are more likely to be chosen, and `beta` sets how strongly.
+    For ``beta=0``, energies play no role.
+
+    The flavor `s` sets how the choice depends on :math:`n_\alpha`:
+
+    * ``s=-1``: only (d-1)-simplices contained in a single d-simplex can be
+      chosen, so each one ends up in at most two d-simplices. The result is a
+      discrete manifold.
+    * ``s=0``: the choice depends on the energies only.
+    * ``s=1``: (d-1)-simplices contained in many d-simplices are more likely
+      to be chosen (preferential attachment).
+
+    Parameters
+    ----------
+    order : int
+        The order d of the simplices to add. Must be at least 1.
+    s : {-1, 0, 1}
+        The flavor parameter.
+    beta : float
+        The inverse temperature. Higher values favor low-energy
+        (d-1)-simplices more strongly.
+    n_iter : int
+        The number of new nodes (and d-simplices) to add.
+    energy_distribution : callable, optional
+        A function that returns one node energy. It can either take a
+        :class:`numpy.random.Generator` as its only argument, e.g.
+        ``lambda rng: rng.exponential(1.0)``, in which case energies are
+        reproducible with `seed`, or take no argument, e.g.
+        ``lambda: np.random.exponential(1.0)``, in which case it handles its
+        own randomness. By default, energies are uniform in [0, 10).
+    seed : int, numpy.random.Generator, or None, optional
+        The seed for the random number generator. By default, None.
+
+    Returns
+    -------
+    S : xgi.SimplicialComplex
+        The generated NGF simplicial complex.
+
+    Raises
+    ------
+    ValueError
+        If `order` is less than 1 or `s` is not -1, 0, or 1.
+
+    See also
+    --------
+    apollonian_complex
+    pseudofractal_simplicial_complex
+
+    References
+    ----------
+    Bianconi, G., & Rahmede, C. (2016).
+    "Network geometry with flavor: from complexity to quantum geometry."
+    Physical Review E, 93(3), 032315.
+    https://arxiv.org/abs/1511.04539
+    """
+    if order < 1:
+        raise ValueError("order must be at least 1.")
+    if s not in {-1, 0, 1}:
+        raise ValueError("s must be -1, 0, or 1.")
+
+    rng = np.random.default_rng(seed)
+
+    if energy_distribution is None:
+        sample_energy = lambda: rng.uniform(0, 10)
+    elif _n_required_args(energy_distribution) == 0:
+        sample_energy = energy_distribution
+    else:
+        sample_energy = lambda: energy_distribution(rng)
+
+    S = SimplicialComplex()
+
+    # create initial d-simplex
+    initial_simplex = list(range(order + 1))
+    energy_nodes = {node: sample_energy() for node in initial_simplex}
+    S.add_simplex(initial_simplex)
+
+    # number of d-simplices attached to each (d-1)-simplex
+    face_counts = {}
+
+    def count_faces(simplex):
+        for face in subfaces([simplex], order=order - 1):
+            face = tuple(sorted(face))
+            face_counts[face] = face_counts.get(face, 0) + 1
+
+    count_faces(initial_simplex)
+
+    for new_node in range(order + 1, order + 1 + n_iter):
+        # attachment weights exp(-beta * energy) * (1 + s * (count - 1)),
+        # computed in log space to avoid underflow at large beta
+        faces = []
+        log_weights = []
+        for face, count in face_counts.items():
+            multiplicity = 1 + s * (count - 1)
+            if multiplicity <= 0:
+                continue  # s = -1: face already shared by two simplices
+            energy = sum(energy_nodes[node] for node in face)
+            faces.append(face)
+            log_weights.append(-beta * energy + np.log(multiplicity))
+
+        log_weights = np.array(log_weights)
+        probs = np.exp(log_weights - log_weights.max())
+        probs /= probs.sum()
+
+        # attach a new node to the chosen (d-1)-simplex
+        chosen_face = faces[rng.choice(len(faces), p=probs)]
+        energy_nodes[new_node] = sample_energy()
+        new_simplex = list(chosen_face) + [new_node]
+        S.add_simplex(new_simplex)
+        count_faces(new_simplex)
+
+    return S
+
+
+def _n_required_args(func):
+    """Number of positional arguments `func` needs (0 if it cannot be inspected)."""
+    try:
+        params = signature(func).parameters.values()
+    except (TypeError, ValueError):
+        return 0
+    return sum(
+        p.default is Parameter.empty
+        and p.kind in (Parameter.POSITIONAL_ONLY, Parameter.POSITIONAL_OR_KEYWORD)
+        for p in params
+    )