# Generators of $H^1(\Gamma, \partial\Gamma^c)$ with $\partial\Gamma^c \subset \partial\Gamma$ for Triangulated Surfaces $\Gamma$: Construction and Classification of Global Loops

Silvano Pitassi

► To cite this version:

Silvano Pitassi. Generators of $H^1(\Gamma, \partial\Gamma^c)$ with $\partial\Gamma^c \subset \partial\Gamma$ for Triangulated Surfaces $\Gamma$: Construction and Classification of Global Loops. 2025. hal-04876473v2

HAL Id: hal-04876473

https://hal.science/hal-04876473v2

Preprint submitted on 7 Apr 2025 (v2), last revised 23 Dec 2025 (v3)

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Generators of $H^1(\Gamma, \partial\Gamma^c)$ with $\partial\Gamma^c \subset \partial\Gamma$ for Triangulated Surfaces $\Gamma$: Construction and Classification of Global Loops

Silvano Pitassi *

## Abstract

Given a compact surface $\Gamma$ embedded in $\mathbb{R}^3$ with boundary $\partial\Gamma$, our goal is to construct a set of representatives for a basis of the relative cohomology group $H^1(\Gamma, \partial\Gamma^c)$, where $\Gamma^c$ is a specified subset of $\partial\Gamma$. To achieve this, we propose a novel graph-based algorithm with two key features: it is applicable to non-orientable surfaces, thereby generalizing previous approaches, and it has a worst-case time complexity that is linear in the number of edges of the mesh $\mathcal{K}$ triangulating $\Gamma$. Importantly, this algorithm serves as a critical pre-processing step to address the low-frequency breakdown encountered in boundary element discretizations of integral equation formulations.

**Keywords:** relative cohomology; loop-star decomposition; global loops; electrical field integral equation.

## 1 Introduction

This paper bridges the fields of algebraic topology and computational electromagnetics. Specifically, it addresses a challenge arising from boundary element method computations for an electromagnetic scattering problem: computing a set of representatives for a basis of the first relative cohomology group $H^1(\Gamma, \partial\Gamma^c)$, where $\Gamma$ is a triangulated surface embedded in $\mathbb{R}^3$ and $\partial\Gamma^c \subset \partial\Gamma$. To tackle this, we propose a worst-case linear-time algorithm.

**Background and motivation.** Throughout this paper, we consider the following prototypical electromagnetic scattering problem, which serves as the primary motivation for the algorithm discussed here. Specifically, an incident harmonic electric field $E^i$, propagating in a space with permittivity $\epsilon$ and permeability $\mu$, impinges upon a triangulated surface $\Gamma$ without boundary, which is a perfect conductor. This results in a scattered electric field $E^s$. To determine $E^s$, we solve the Electric Field Integral Equation (EFIE)

$$ \gamma_T \circ \mathfrak{U}_k(J) = -\gamma_T E^i, $$

which involves finding the induced surface current density $J$ on $\Gamma$ determining $E^s$ via the EFIE integral operator $\gamma_T \circ \mathfrak{U}_k$ [1]. Here, $\gamma_T$ denotes the tangential trace operator on $\Gamma$.

Crucially, when the EFIE is discretized using the Boundary Element Method (BEM), two numerical challenges arise as the frequency of the involved harmonic fields approaches zero—a well-documented phenomenon referred to as the low-frequency breakdown. The first challenge is the growth in the condition number of the resulting linear system after BEM discretization, while the second is the occurrence of round-off errors due to finite precision in the numerical integration of the EFIE’s source term. These two effects can lead to inefficiencies in standard numerical solvers, requiring significantly higher computational resources (e.g., more memory or more solver iterations), or even causing the solution to fail to converge [2].

To address the low-frequency breakdown of the EFIE, the key idea is to decompose the current density $J$ as the sum $J = J_\Lambda + J_\Sigma$, where $J_\Lambda$ is solenoidal (i.e., $\mathrm{div}_\Gamma J_\Lambda = 0$, with $\mathrm{div}_\Gamma$ denoting the surface divergence operator on $\Gamma$) and $J_\Sigma$ is non-solenoidal. This decomposition is essentially

<footnote>*IMAG, University of Montpellier, Montpellier, France (silvano.pitassi@umontpellier.fr)</footnote>
derived from the Hodge decomposition theorem applied to the L²-based *de Rham complex* on Γ, along with the definition of the first (absolute) cohomology group H¹(Γ) as a quotient space, so that we can write:

$$ \mathbf{J} = \underbrace{\mathrm{curl}_{\Gamma} \psi + \mathbf{g}}_{=:\mathbf{J}_{\Lambda}} + \underbrace{\mathrm{grad}_{\Gamma} \varphi}_{=:\mathbf{J}_{\Sigma}} \quad (1) $$

where grad_Γ is the L²-adjoint of div_Γ, and g represents a cohomology class in H¹(Γ). The vector field ψ is referred to as the *stream function* or *local loop*, and g is called a (*cohomology*) *generator* or *global loop*.

Two approaches for computing this decomposition have been proposed in the literature. Initially, the *loop-star* decomposition was introduced [3, 4]. This approach requires computing a set of generators for the image spaces Im curl_Γ and Im grad_Γ, as well as a set of generators for the cohomology space H¹(Γ). More recently, *quasi-Helmholtz projectors* have been proposed as an alternative to avoid explicitly computing cohomology generators [5]. This method computes the non-solenoidal part J_Σ via an orthogonal projection P_Σ onto Im grad_Γ, while the solenoidal part J_Λ is found using the complementary orthogonal projection P_Λ.

However, we emphasize that the computation of generators of H¹(Γ) is, in practice, strictly necessary to address the low-frequency breakdown of the EFIE, even when quasi-Helmholtz projectors are used. This point is supported by recent work [6], which introduces a technique for reducing round-off errors in the numerical evaluation of the EFIE's source term. Notably, the methodology outlined in Sections IV A and B of that paper clearly shows the need for a set of generators of H¹(Γ), as they appear explicitly in the definition of the proposed quasi-Helmholtz projectors (cf. [6, Equations (24)–(28)]).

**Review of cohomology computations in the numerical analysis community.** We begin by pointing out early works by Hiptmair and Ostrowski in [7] (which we refer to as the *HO-algorithm* hereafter), along with a non-exhaustive list of related works in [8–10]. These articles focus on a surface Γ without boundary and aim to compute generators of the first homology group H₁(Γ). However, in numerical analysis and engineering applications, what is typically needed are generators of H¹(Γ), rather than those of H₁(Γ).

Since Γ is orientable (i.e., the boundary is a closed subset of R³), it is possible to recover generators of H¹(Γ) from those of H₁(Γ) using the Poincaré-Lefschetz Duality Theorem [11]. This approach is mentioned in the introduction of [7], though the actual implementation details are missing. Instead, an algorithm that implements this construction was first developed by Dłotko in [12] and later reorganized by Dłotko and Specogna in the series of works [13, 14] (referred to as the *DS-algorithm*).

Now, consider the case where Γ is a surface with a non-empty boundary ∂Γ. In this case, note that the DS-algorithm can still be applied with a simple modification of the construction, as described in [12, Section 6]. However, in the context of our electromagnetic application, the cohomology space H¹(Γ) is no longer appropriate. We must instead account for the boundary condition J · n_∂Γ = 0 on ∂Γ, which arises from the EFIE model.

To satisfy this boundary condition, it can be shown that the vector field g in (1) must belong to the first *relative cohomology space* H¹(Γ, ∂Γ) [11]. Therefore, algorithms like the DS-algorithm or the HO-algorithm are unsuitable for computing generators of H¹(Γ, ∂Γ), since H¹(Γ) and H¹(Γ, ∂Γ) are generally quite different spaces; see Fig. 1. As far as we are aware, no combinatorial algorithm has been devised to address this case.

In this paper, we build upon this observation and consider the more general problem of computing generators of H¹(Γ, ∂Γᶜ), where ∂Γᶜ is an arbitrary subset of the boundary ∂Γ of Γ. This general version arises in the EFIE model, specifically for cases where the surface Γ is connected to external lumped-element circuits through suitable *contacts* or *ports*, such as in antennas and their connections with electronic circuits.

From a physical perspective, these connections result in net currents flowing through the body from one electrode to another, a phenomenon that would otherwise be absent. Crucially, this modifies the topology of the surface. For example, in the simplest case of a simply connected surface Γ, these connections make the surface non-simply connected. We discuss the details of the EFIE model with contact regions in Section 3.
<figure><img src="image_2.png" /><figcaption>**Figure 1:** A triangulation of a hollow disk $\Gamma$. The support of the representative of the $H^1(\Gamma)$ generator, as formally defined in Section 2.2, is shown in red, while the support of the representative of the $H^1(\Gamma, \partial\Gamma)$ generator is shown in blue. The green part shows the intersection of the two supports. The shaded grey triangles indicate the support of the corresponding generators as functions in $L^2$, obtained by interpolation using the standard Raviart-Thomas basis functions.</figcaption></figure>

**The main results and their novelty w.r.t. literature.** The primary objective of this paper is to provide an efficient algorithm for constructing and classifying the generators of the first relative cohomology group $H^1(\Gamma, \partial\Gamma^c)$, where $\Gamma$ is an arbitrary triangulated surface embedded in $\mathbb{R}^3$ and $\partial\Gamma^c$ is a subset of its boundary. Specifically, we present two key contributions: one for the engineering community and the other for the mathematical community.

The first contribution is a rigorous definition of *global loops* [15] in the context of *cohomology theory*. In Definition 3.1 at the end of Section 3, we demonstrate that when global loops are defined as generators of the first relative cohomology group $H^1(\Gamma, \partial\Gamma^c)$ of a triangulated surface that excludes contact regions, the decomposition of the current density $J$ in (1) holds. This formal definition has several important implications. First, it covers cases where $\Gamma$ has no contact regions and no boundary, in which the space $H^1(\Gamma, \partial\Gamma^c)$ reduces to both $H^1(\Gamma, \partial\Gamma)$ and $H^1(\Gamma)$, since we have $\partial\Gamma^c = \partial\Gamma$ and $\partial\Gamma = \emptyset$, respectively. Thus, our definition naturally encompasses a variety of modeling contexts, correctly identifying the appropriate cohomology space for each case, whether relative or absolute.

Second, our definition is based directly on the *cochain complex* (see [16, Section 2.6]) supported on a triangulation $\mathcal{K}$ of $\Gamma$, rather than the standard $L^2$ de Rham complex typically used in the engineering literature (see [2, Section IV]). In the latter framework, conforming basis functions are added to the set of standard *Raviart-Thomas basis functions* to span the de Rham first cohomology group. By contrast, our approach, grounded in the cochain complex, provides new insights into the construction of algorithms, which can be efficiently implemented using graph-based techniques. In particular, our algorithm has a worst-case time complexity that is linear in the number of edges of the triangulation, making it highly efficient for engineering applications, such as those involving the EFIE model. In fact, regarding the problem of computing global loops in engineering applications, it has even been suggested that “it is impossible to provide a complete description of this family of functions” [5].

Regarding the mathematical contribution, a significant novel aspect is that our algorithm is provably correct for non-orientable surfaces, such as the *Möbius strip*. Since we assume that $\Gamma$ is embedded in $\mathbb{R}^3$, it follows that $\Gamma$ must have a non-empty boundary. This is because a non-orientable surface without boundary cannot be embedded in $\mathbb{R}^3$ (e.g., the Klein bottle), while every surface with boundary can be embedded in $\mathbb{R}^3$. The key challenge for non-orientable surfaces is that a generator of the *torsion subgroup* appears when considering relative homology with integer coefficients (a fact evident from a canonical representation of the Möbius strip). The problem,
then, is how to distinguish the torsion generator from other homology generators.

Crucially, in the proof of correctness presented in Section 5, we explicitly identify the torsion
generator. In doing so, we also resolve an open problem mentioned in [12, Section 7]. The proof
relies on concepts from *Discrete Morse Theory*, which has several far-reaching implications. First,
it generalizes the construction used in the HO-algorithm, making our algorithm a natural extension
of the HO-algorithm to the case of non-orientable surfaces. Second, it provides an alternative proof
of the classical Classification Theorem for compact surfaces [17], which follows as a corollary of
our specific constructions.

Another important aspect is the treatment of contact regions. While the EFIE model in the presence of contacts has been studied in engineering works [18, 19], the focus has primarily been on the design of the coupling between the EFIE model and external circuit equations, as well as the quasi-static assumptions near contact regions. Mathematically, contributions have been made to establish links with relative cohomology, beginning with seminal ideas in [11] and revisited in works such as [20–22], among others. However, these contributions have largely focused on orientable domains, where the Poincaré-Lefschetz Duality Theorem can be applied.

In contrast, when torsion generators are present, additional complications arise, particularly the
introduction of an extra generator associated with contact regions. This stands in stark contrast
to the orientable case, where a strong physical interpretation is also available. Indeed, techniques
that rely on duality theorems do not necessarily work in non-orientable cases, where the hypotheses
underlying these theorems are not satisfied.

We believe that the combination of all these features opens up new possibilities for efficient
numerical procedures to solve the EFIE.

2 Preliminaries from algebraic topology

In this preliminary section, we introduce concepts from algebraic topology and provide specific
definitions tailored to our purposes. In particular, we define the concept of a cochain as a model
for a discrete vector field.

In this paper, we consider two types of cochains. The first are integer-valued cochains, such as
representatives of the generators of the first cohomology group over the integers. The second are
real-valued cochains, which are used to model physical variables.

**2.1 CW complex and triangulations**

A 2-manifold $\Gamma$ with boundary is a Hausdorff topological space in which every point has a neighborhood homeomorphic to either $\mathbb{R}^2$ or the half-plane $\{(a, b) \in \mathbb{R}^2 | a \ge 0\}$. The boundary of $\Gamma$, denoted by $\partial\Gamma$, consists of all points with neighborhoods homeomorphic to the half-plane. We say that $\Gamma$ is a 2-manifold without boundary, or simply a 2-manifold, if $\partial\Gamma$ is empty.

An (open) $m$-cell, for $m \in \{0, 1, 2\}$, is a topological space homeomorphic to the interior of an $m$-ball $\text{int } B^m$ (note that $B^0$ is a single point).

Let us now introduce the concept of a CW complex $\mathcal{K}$ defined on $\Gamma$.

**Definition 2.1 (CW complex).** A (2-dimensional) *CW complex* $\mathcal{K}$ (on $\Gamma$) is a partition of $\Gamma$ into $m$-cells $x$ such that the following holds: for each $m$-cell $x$, there exists a continuous surjection (the *gluing map*) $g_x : B^m \to x$ satisfying:

(i) The restriction of $g_x$ to the interior, $g_x : \text{int } B^m \to x$, is a homeomorphism.

(ii) The image of the boundary, $g_x(\partial B^m)$, is covered by the closure of a finite number of $k$-cells
with $0 \le k < m$.

For $m \in \{0, 1, 2\}$, we denote by $\mathcal{K}^{(m)}$ the set of all $m$-cells of $\mathcal{K}$.

A CW complex $\mathcal{K}$ is called *regular* if, for every $m$-cell $x$, the restriction of the gluing map $g_x : \partial B^m \to g_x(\partial B^m)$ is a homeomorphism.

A regular CW complex $\mathcal{K}$ where each $m$-cell is an $m$-simplex is also referred to as a *triangulation*.
In this case, we adopt a more familiar notation and terminology for the $m$-cells. Specifically: a
node n is a 0-cell in $V := K^{(0)}$, an edge e is a 1-cell in $E := K^{(1)}$, and a face f is a 2-cell in $F := K^{(2)}$. Accordingly, we denote a triangulation $K$ by $K = (V, E, F)$.

A CW subcomplex $S$ of $K$ is a union of $m$-cells $x \in K$ such that, if $x \subseteq S$, then the closure of $x$ is also a subset of $S$.

If $\Gamma$ has a non-empty boundary $\partial\Gamma$, we denote by $\partial K$ the subcomplex of $K$ induced on $\partial\Gamma$. Note that $\partial\Gamma$ is itself a manifold because $\Gamma$ is.

The 1-skeleton of $K$ is the union $K^{(0)} \cup K^{(1)}$, i.e., the set of all 0- and 1-cells. It is straightforward to observe that the 1-skeleton of $K$ can be regarded as a graph whose vertices and edges are $K^{(0)}$ and $K^{(1)}$, respectively.

Using this identification, a spanning tree $T = (V_T, E_T)$ of $K$ refers to a spanning tree of the 1-skeleton of $K$, treated as a graph.

## 2.2 (Co)chains and (co)homology

Let $K$ be a CW complex, and let $G$ be either the module of integers $\mathbb{Z}$ or reals $\mathbb{R}$. We start by fixing an orientation of all its $m$-cells (see, e.g., [28, Chapter 2]). The $m$-th chain group $C_m(K; G)$ is the group of all formal sums $c = \sum_{x \in K^{(m)}} a_x x$ with $a_x \in G$. We call each element of $C_m(K; G)$ a $m$-chain. For two $m$-chains $c = \sum_{x \in K^{(m)}} a_x x$ and $c' = \sum_{x \in K^{(m)}} b_x x$, we define their scalar product by

$$ \langle c, c' \rangle := \sum_{x \in K^{(m)}} a_x b_x. $$

The $m$-th cochain group $C^m(K; G)$ is the group of all (linear) maps $c : C_m(K; G) \to G$ such that $c(c) = c(\sum_{x \in K^{(m)}} a_x x) = \sum_{x \in K^{(m)}} a_x c(x)$ for every $c \in C_m(K; G)$. We call each element of $C^m(K; G)$ a $m$-cochain. For a $m$-cochain $c \in C^m(K; G)$, its support is the set of all $m$-cells $x \in K^{(m)}$ such that $c(x) \neq 0$. Notice that to compute the value $c(c)$ it is sufficient to know the values $c(x)$ on every $m$-cell. Therefore, we can associate to the $m$-cochain $c$ the $m$-chain $\phi(c) := \sum_{x \in K^{(m)}} a_x x$ with $a_x = c(x)$ for $x \in K^{(m)}$ in such a way that, for each $m$-chain $c' \in C_m(K; G)$, it holds

$$ c(c') = \langle \phi(c), c' \rangle $$

In what follows, we will simply write $\langle c, c' \rangle$ in place of $\langle \phi(c), c' \rangle$ so that, by the above equality, this expression denotes also the value $c(c')$.

For $m \in \{1, 2\}$, let $\iota : K^{(m)} \times K^{(m-1)} \to \mathbb{Z}$ be the incidence number between a $m$-cell and a $(m-1)$-cell (see [29, Section 2.2], and, in particular, the cellular boundary formula on p. 140; equivalently, see [28]). If $\iota(x, x') \neq 0$, then we say that $x$ is incident on $x'$, or simply that $x$ and $x'$ are incident. We use incidence numbers to define the boundary map $\partial_m : C_m(K; G) \to C_{m-1}(K; G)$ for $m \in \{1, 2\}$. For each $m$-cell $x \in K^{(m)}$ we define $\partial_m x := \sum_{x' \in K^{(m-1)}} \iota(x, x')x'$, and for a $m$-chain $c = \sum_{x \in K^{(m)}} a_x x \in C_m(K; G)$, by linearity, we define $\partial_m c := \sum_{x \in K^{(m)}} a_x \partial_m x$. It can be checked that $\partial_{m+1} \circ \partial_m = 0$ for $m \in \{0, 1\}$. The coboundary map $\delta^m : C^{m-1}(K; G) \to C^m(K; G)$, for $m \in \{0, 1\}$, is defined as follows. For a $m$-cochain $c \in C^m(K; G)$, $\delta^m c$ is the (unique) $(m+1)$-cochain such that $\langle \delta^m c, c' \rangle = \langle c, \partial_{m+1} c' \rangle$ for all $c' \in C_{m+1}(K; G)$. Using $\partial_{m+1} \circ \partial_m = 0$, it can be checked that $\delta^m \circ \delta^{m+1}$ for $m \in \{0, 1\}$, so that the (graded) family $(C^m(K; G), \delta^m)$ with $m \in \{0, 1, 2\}$ is a (cochain) complex once we also set $\delta^2 := 0$ [16].

The boundary operator gives rise to a classification of $m$-chains in $K$. We define the subgroup $Z_m(K; G) = \{c \in C_m(K; G) \mid \partial_m c = 0\}$ and the group $B_m(K; G) = \{\partial_{m+1} c \mid c \in C_{m+1}(K; G)\}$. The $m$-th homology group is then the quotient group $H_m(K; G) := Z_m(K; G)/B_m(K; G)$.

In a similar way, we define the group $Z^m(K; G) = \{c \in C^m(K; G) \mid \delta^m c = 0\}$ and the group $B^m(K; G) = \{\delta^{m-1} c \mid c \in C^{m-1}(K; G)\}$. We call each element of $Z^m(K; G)$ a $m$-cocycle. The $m$-th cohomology group is then the quotient group $H^m(K; G) := Z^m(K; G)/B^m(K; G)$. When the group $G$ is the group of reals $\mathbb{R}$, in the following we will simply write $H^m(K)$ in place of $H^m(K; \mathbb{R})$. A cohomology basis of $H^m(K; G)$ is a set of equivalence classes of $m$-cocycles $\{[g_1], \dots, [g_N]\} \subset H^m(K; G)$ such that every equivalence class in $H^m(K; G)$ can be obtained in a unique way as a linear combination $\sum_{i=1}^N \alpha_i [g_i]$ with $\alpha_i \in G$. In this case, we say that the set of $m$-cocycles $\{g_1, \dots, g_N\}$ are cohomology generators of $H^m(K; G)$, and we call the number $N$ the rank of $H^m(K; G)$ and we denote it by $\mathrm{rank}(H^m(K; G))$.
In the following, we will also need the standard concept of relative cohomology. In relative cohomology, some parts of the complex may be considered irrelevant. In our case, we will focus on the CW subcomplex $\partial K^c$ induced by $K$ on a closed subset of $\partial \Gamma^c \subset \partial \Gamma$. We define the m-th relative cochain group by $C^m(K, \partial K^c; G) := C^m(K; G)/C^m(\partial K^c; G)$ (notice that $C_2(\partial K^c; G)$ is the zero group since $\partial K^c$ is 1-dimensional). Accordingly, the definition of the m-th relative cohomology group $H^m(K, \partial K^c; G)$ remains unchanged with respect to the absolute version once relative 1-chains are used in the corresponding definitions. The following characterization of elements of $C^1(K, \partial K^c; G)$ will be useful later: a 1-cochain $c \in C^1(K; G)$ is in $C^1(K, \partial K^c; G)$ if and only if the support of $c$ does not contain any 1-cell in $\partial K^c$.

The m-th relative homology group $H^m(K, \partial K^c; \mathbb{Z})$ is an Abelian group by construction. Therefore, using the Fundamental Theorem for Finitely Generated Abelian Groups (see e.g. [11]), it is isomorphic to the direct sum of a free Abelian group made of $\beta_1 \in \mathbb{N}$ cohomology generators and a torsion group made of $\tau_1 \in \mathbb{N}$ torsion generators. In our target application, the relationship between $H^m(K, \partial K^c; \mathbb{Z})$ and $H^m(K, \partial K^c; \mathbb{R})$ will be important. In fact, for the numerical approximation of physical fields, we employ real coefficients, and hence generators of $H^m(K, \partial K^c; \mathbb{R})$, while computational topology algorithms employ integer coefficients and hence generators of $H^m(K, \partial K^c; \mathbb{Z})$. Despite this distinction, we notice that when $H^m(K, \partial K^c; \mathbb{Z})$ is found, the cohomology generators $H^m(K, \partial K^c; \mathbb{R})$ can be deduced from those of $H^m(K, \partial K^c; \mathbb{Z})$ by not taking into account the torsion ones; see Section 5.5 and, in particular, Corollary 1 for the formal statement of this fact.

# 3 Definition of global loops

Let $\Gamma$ be a 2-manifold with boundary, embedded in $\mathbb{R}^3$, and satisfying the conditions described in [1, Section 1.1]. As explained in the introduction, an induced current density $J$ on $\Gamma$ arises either from an incident electric field or from an external circuit. In this section, we focus on the latter case and specifically review how the EFIE model on $\Gamma$ can be coupled with lumped-element circuit models.

This coupling requires defining contact regions as the intersections between $\Gamma$ and the physical connections of the circuits. To rigorously define these contact regions, we first assume that $\Gamma$ is the boundary of a domain $\Omega$ in $\mathbb{R}^3$. We then extend the discussion to the general case, where $\Gamma$ represents a screen. This is modeled by assuming $\Omega$ to be a thin conductor, where the thickness is small compared to the wavelength and other characteristic lengths.

This distinction is necessary to justify the inclusion of contact regions in the definition of a suitable triangulation of $\Gamma$. It also lays the foundation for our general definition of global loops, presented in Definition 3.1.

## 3.1 EFIE modeling in presence of contacts

Assume that $\Gamma = \partial \Omega$. The electric field $E$ must satisfy the following condition:

$$ \gamma_T \mathbf{E} = \mathbf{0} \quad \text{on } \Gamma. \qquad (2) $$

Let $E^i$ and $E^s$ denote the incident and scattered electric fields from $\Gamma$, respectively, so that the total electric field $E$ is given by

$$ \mathbf{E} = \mathbf{E}^i + \mathbf{E}^s. \qquad (3) $$

In terms of electromagnetic potentials $A$ and $\varphi$, the scattered electric field can be written as

$$ \mathbf{E}_s = -i\omega\mathbf{A} - \text{grad } \varphi.s \qquad (4) $$

The potentials $A$ and $\varphi$ are expressed in terms of the equivalent current density $J$ and charge density $\sigma$ on $\Gamma$ as follows:

$$ \mathbf{A}(\mathbf{x}) = \mu \int_{\Gamma} G_k(\mathbf{x}, \mathbf{x}') \mathbf{J}(\mathbf{x}') d\mathbf{x}', \quad \varphi(\mathbf{x}) = \frac{1}{\epsilon} \int_{\Gamma} G_k(\mathbf{x}, \mathbf{x}') \sigma(\mathbf{x}') d\mathbf{x}', $$

where

$$ G_k(\mathbf{x}, \mathbf{x}') = \frac{e^{ik|\mathbf{x} - \mathbf{x}'|}}{4\pi|\mathbf{x} - \mathbf{x}'|} $$
is the free-space Green's function, $k = \omega/c$, and $c = 1/\sqrt{\epsilon\mu}$.

We also need to consider the standard continuity equation, which relates the charge density $\sigma$ to the current density $J$:

$$i\omega\sigma = -\operatorname{div}_{\Gamma} J \quad \text{on } \Gamma.$$

Next, we introduce the concept of contact. Let $\Gamma^c \subset \Gamma$ be the subset where the physical connections of lumped circuits intersect $\Omega$. Thus, $\Gamma$ is partitioned as

$$\Gamma = \Gamma^c \cup \Gamma^{nc},$$

where $\Gamma^{nc} := \Gamma \setminus \Gamma^c$. We assume that $\Gamma^c$ decomposes into disjoint open sets $\{\Gamma_i^c\}_{i=1}^{N^c}$, where each $\Gamma_i^c$ is homeomorphic to an (open) disk. Each such $\Gamma_i^c$ is referred to as a contact.

Each contact $\Gamma_i^c$ represents an electrically small region, meaning that its diameter in $\mathbb{R}^3$ is small compared to the smallest characteristic wavelength. Consequently, for all $i \in \{1, \dots, N^c\}$, we have

$$\varphi = \operatorname{const} = V_i \quad \text{and} \quad \gamma_T \mathbf{A} = \mathbf{0} \quad \text{on} \quad \Gamma_i^c, \qquad (5)$$

where $V_i$ is the constant scalar potential enforced on $\Gamma_i^c$, which corresponds to the voltage of the corresponding circuit node. The potential difference is calculated with respect to the potential at infinity, assumed to be zero.

Now, consider the case of a conducting body $\Omega$ whose thickness is much smaller than the smallest characteristic wavelength. In this case, $\Omega$ can be modeled as a screen $\Gamma$ with possibly a non-empty boundary $\partial\Gamma$.

In the presence of contacts, the screen $\Gamma$ is again partitioned as $\Gamma = \Gamma^c \cup \Gamma^{nc}$, as described above. Furthermore, the boundary $\partial\Gamma$ can be also partitioned as $\partial\Gamma = \partial\Gamma^c \cup \partial\Gamma^{nc}$, where $\partial\Gamma^c$ is the disjoint union of the contact regions. This partitioning arises when a contact $\Gamma_i^c$ intersects $\Omega$ in a subset of $\Gamma$, and we still want to include its effects in the model.

In this case, the following boundary condition holds on $\partial\Gamma$:

$$J \cdot n_{\partial\Gamma} = 0 \quad \text{on} \quad \partial\Gamma^{nc}, \qquad (6)$$

where $n_{\partial\Gamma}$ is the unit outer normal to $\partial\Gamma$.

## 3.2 EFIE discretization and interpretation of discrete unknowns as real-valued cochains

To motivate the general definition of global loops, we first recall the initial steps of the Galerkin method for the EFIE, considering the case where $\Gamma$ is a piecewise-smooth 2-manifold.

Let $\mathcal{K}$ be a triangulation of $\Gamma$, and let $W_h^1$ denote the Sobolev space spanned by the standard Raviart-Thomas basis functions, supported on $\mathcal{K}$ and satisfying the boundary condition on $\partial\Gamma$ as given in (6). The boundary condition (2) may be enforced on the Galerkin space $W_h^1$ by requiring that

$$\int_{\Gamma} E \cdot j = 0, \quad \forall j \in W_h^1. \qquad (7)$$

Substituting the expressions for $E$ from (3) and $E^s$ from (4) into (7), and using integration by parts with the two conditions (5) on $\Gamma^c$, we obtain the following equation

$$i\omega\mu \int_{\Gamma^{nc}} \mathbf{A} \cdot j \, dS + \int_{\Gamma^{nc}} \varphi \operatorname{div}_{\Gamma} j \, dS = \sum_{i=1}^{N} \int_{\partial\Gamma_i^c} \varphi j \cdot \mathbf{n}_{\partial\Gamma_i^c} \, dl + \int_{\Gamma^{nc}} E^i \cdot j \, dS, \qquad (8)$$

where we have used the conditions (5) on each contact $\Gamma_i^c$. Notice that if the electric potential $\varphi$ is assigned on $\partial\Gamma_i^c$, then the right-hand side of the equation is known. On the other hand, specifying the normal components of the current density, $J \cdot n_{\partial\Gamma_i^c}$, is related to our formal definition of global loops. Also, note that $\Gamma^{nc}$ is generally the sub-manifold of $\Gamma$ over which the field integrals are evaluated in the weak formulation (8). Indeed, in the absence of contacts, we have $\Gamma^c = \emptyset$, so $\Gamma^{nc} = \Gamma$. Therefore, we can focus on a triangulation $\mathcal{K}$ of $\Gamma^{nc}$.

Now, let $F_h^1: W_h^1 \to C^1(\mathcal{K})$ be the De Rham map [16, 30], restricted to the space $W_h^1$. Crucially, this map is a cochain isomorphism between the cochain complex and the finite element complex
supported on the triangulation $\mathcal{K}$. More precisely, taking into account the boundary condition (6) on $\partial\mathcal{K}$, the map $F_h^1$ defines an isomorphism between $W_h^1$ and the set of relative 1-cochains $C^1(\mathcal{K}, \partial\mathcal{K}^{nc})$ [31].

We are now ready to present the general definition of global loops.

**Definition 3.1** (Definition of global loops). The set $GL = \{g_i\}_{i=1}^{\beta_1(\Gamma^{nc})}$ of global loops on $\Gamma$ is defined as a set of 1-cocycles, representing a basis of the first relative cohomology group $H^1(\mathcal{K}, \partial\mathcal{K}^{nc})$ of the triangulation $\mathcal{K}$ of the closed sub-manifold $\Gamma^{nc}$, excluding the contact regions.

A crucial aspect of this definition is that the global loops are defined as 1-cochains supported on a triangulation $\mathcal{K}$, rather than as vector fields on $\Gamma^{nc}$. As we will show in the subsequent sections, this combinatorial definition provides a more effective approach for handling topologically complex surfaces, as it suggests efficient graph-based algorithms.

# 4 General algorithm to construct global loops

As formalized in Definition 3.1, global loops should be understood as generators of the first cohomology group $H^1(\mathcal{K}, \partial\mathcal{K}^{nc})$, where $\mathcal{K}$ is a triangulation of the sub-manifold $\Gamma^{nc}$ of $\Gamma$ that excludes contact regions. The types of generators required depend on the topology of $\mathcal{K}$ and the presence or absence of contact regions.

In this section, we present a general algorithm to compute these generators, which is detailed in Algorithm 5. The key idea behind the main procedure in Algorithm 5 is that the set $GL$ of all global loops of $\mathcal{K}$ can be partitioned into three subsets that are non-empty according to the presence of, informally, “handles”, “holes”, and contacts on $\mathcal{K}$, and that they can be computed independently using three distinct procedures.

We assume that $\mathcal{K}$ is connected for the purposes of this discussion. This assumption is not restrictive in general; indeed, if $\mathcal{K}$ consists of $p$ connected components $K_i$ for $i \in \{1, \dots, p\}$, we have the following decomposition

$$H^1(\mathcal{K}) \cong \bigoplus_{i=1}^{p} H^1(K_i),$$

which implies that it suffices to consider each connected component $K_i$ individually.

## 4.1 Intepreting cocyles as cycles on the dual graph $\tilde{G}$

In the constructions used in Algorithm 5, we require the definition of a special graph associated with the triangulation $\mathcal{K}$. This graph, known as the dual graph, is denoted by $\tilde{G} = (\tilde{V}, \tilde{E})$, where $\tilde{V}$ and $\tilde{E}$ are defined using the (bijective) duality map $D: E \cup F \to \tilde{V} \cup \tilde{E}$ as follows:

* For each face $f \in F$, the corresponding dual node is $\tilde{n} = D(f) \in \tilde{V}$.

* For each edge $e \in E$, which is a common edge of two faces $f_e^{(1)}, f_e^{(2)} \in F$, the corresponding (non-oriented) dual edge is $\tilde{e} = \{D(f_e^{(1)}), D(f_e^{(2)})\} = D(e) \in \tilde{E}$.

Additionally, we need to describe how the duality map $D$ acts on edges that lie on the boundary $\partial\mathcal{K}$ of $\mathcal{K}$ by defining the set $\tilde{V}_\partial$:

* For each edge $e \in \partial\mathcal{K}$, the corresponding dual node is $\tilde{n} = D_\partial(e) \in \tilde{V}_\partial$.

* For each edge $e \in \partial\mathcal{K}$, the corresponding dual edge is $\tilde{e} = \{D_\partial(e), D(f)\} = D(e) \in \tilde{E}$, where $f$ is the unique face in $F$ incident to $e$.

A path $\tilde{p}$ of $\tilde{G}$ is a set of $l \in \mathbb{N}, l \ge 1$ distinct dual nodes $\tilde{n}_1, \dots, \tilde{n}_l$ and $l-1$ distinct dual edges $\tilde{e}_1, \dots, \tilde{e}_{l-1}$ such that $\tilde{e}_i = \{\tilde{n}_i, \tilde{n}_{i+1}\}$ for $i \in \{1, \dots, l-1\}$ ($l=1$ means the path reduces to the single dual node $\tilde{n}_1$). By duality via the map $D$, there is a corresponding set of distinct $l$-faces $f_1 = D^{-1}(\tilde{n}_1), \dots, f_l = D^{-1}(\tilde{n}_l)$ and $l-1$ distinct edges $e_1 = D^{-1}(\tilde{e}_1), \dots, e_{l-1} = D^{-1}(\tilde{e}_{l-1})$ of $\mathcal{K}$.

We use paths on the dual graph $\tilde{G}$ to construct 1-cochains of $\mathcal{K}$ according to the following procedure.
**Algorithm 1** Construct a 1-cocycle g on K, taking as input a path $\tilde{p}_{e,e'}$ in $\tilde{G}$ and a pair of edges $(e, e')$ of K, where $e \subset f_1$ and $e' \subset f_l$.

1: initialize g = 0

2: set < g, e > to 1

3: **for** each dual edge $\tilde{e}_i$ of $\tilde{p}_{e,e'}$ **do**

4: compute

▷ initialize the value of g on the edge e to 1

▷ each $\tilde{e}_i$ is selected by walking along $\tilde{p}$ from $\tilde{e}_1$ to $\tilde{e}_l$

$$< g, e_i > := - \frac{\iota(f_i, e_i)}{\iota(f_i, e_{i-1})} < g, e_{i-1} >$$

▷ where $e_i = D^{-1}(\tilde{e}_i)$, and the conventions $e_0 := e$ and $e_l := e'$ are applied

5: **end for**

6: **if** < g, e' > = 0 **then**

7: **return** 0

8: **else**

9: **return** g

10: **end if**

## 4.2 Global loops corresponding to “handles”: the set $GL^{ha}$

The generators arising from the presence of “handles” are the more difficult to find, and the novel algorithm to compute them reads as follows.
<caption>**Algorithm 2** Construct the set $GL^{ha}$</caption>

1: construct a spanning tree $T$ on $\partial K$
2: extend $T$ to a spanning tree on the whole $K$
3: construct a spanning tree $\tilde{T}$ on the subgraph $\tilde{G}' = (\tilde{V}, \tilde{E} \setminus \bigcup_{e \in E_T} \{D(e)\}) \triangleright$ i.e., the subgraph of $\tilde{G}$ made of all dual edges in $\tilde{G}$ that are not dual to an edge in $T$
4: **for** each connected component $\partial K_k$ of $\partial K$ **do**
5: &nbsp;&nbsp;&nbsp;&nbsp;select an edge $e$ of $\partial K_k$
6: &nbsp;&nbsp;&nbsp;&nbsp;add the dual edge $D(e)$ to $E_{\tilde{T}}$ and the dual node $D_{\partial}(e)$ to $V_{\tilde{T}}$
7: **end for**
8: let $E_c := \{e \in E \mid e \notin E_T, D(e) \notin E_{\tilde{T}}\} \triangleright$ i.e., the set of edges not in $T$ whose dual is not in $\tilde{T}$
9: **for** each $e \in E_c$ **do**
10: &nbsp;&nbsp;&nbsp;&nbsp;let $f_e^{(1)}, f_e^{(2)}$ be the unique two faces of $K$ incident on $e$
11: &nbsp;&nbsp;&nbsp;&nbsp;construct the unique path $\tilde{p}_e$ in $\tilde{T}$ from $D(f_e^{(1)})$ to $D(f_e^{(2)})$
12: &nbsp;&nbsp;&nbsp;&nbsp;apply Algorithm 1 with input $\tilde{p}_e$ and the pair $(e, e)$ to get $g_e^{ha} \triangleright$ the pair $(e, e)$ consists of the same edge $e$
13: **end for**
14: let $E_c^{\text{II}} := \{e \in E_c \mid g_e^{ha} = 0\} \triangleright$ i.e., the set of edges for which the 1-cochain $g_e^{ha}$ computed from Algorithm 1 is the zero 1-cochain
15: **if** $|E_c^{\text{II}}| = 0$ **then**
16: &nbsp;&nbsp;&nbsp;&nbsp;**return** $GL^{ha} := \{g_1^{ha}, \dots, g_{N^{ha}}^{ha}\} = \{g_e^{ha}\}_{e \in E_c}$
17: **else**
18: &nbsp;&nbsp;&nbsp;&nbsp;select an edge $e^* \in E_c^{\text{II}}$
19: &nbsp;&nbsp;&nbsp;&nbsp;**for** each $e \in E_c^{\text{II}} \setminus \{e^*\}$ **do**
20: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;let $f_e^{(1)}, f_e^{(2)}$ be the unique two faces of $K$ incident on $e$
21: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;let $f_{e^*}^{(1)}, f_{e^*}^{(2)}$ be the unique two faces of $K$ incident on $e^*$
22: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;construct the unique path $\tilde{p}_{e,e^*}^{(1)}$ in $\tilde{T}$ from $D(f_e^{(1)})$ to $D(f_{e^*}^{(1)})$
23: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;construct the unique path $\tilde{p}_{e,e^*}^{(2)}$ in $\tilde{T}$ from $D(f_e^{(2)})$ to $D(f_{e^*}^{(2)})$
24: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;apply Algorithm 1 with input $\tilde{p}_{e,e^*}^{(1)}$ and $(e, e^*)$ to get $g_e^{(1)}$
25: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;apply Algorithm 1 with input $\tilde{p}_{e,e^*}^{(2)}$ and $(e, e^*)$ to get $g_e^{(2)}$
26: &nbsp;&nbsp;&nbsp;&nbsp;define $g_e^{ha}$ by
$$ \langle g_e^{ha}, e' \rangle := \begin{cases} \langle g_e^{(1)}, e' \rangle & \text{if } e' \in \{e, e^*\}, \\ \langle g_e^{(1)}, e' \rangle + \langle g_e^{(2)}, e' \rangle & \text{if } e' \notin \{e, e^*\} \end{cases} $$
27: **end for**
28: **return** $GL^{ha} := \{g_1^{ha}, \dots, g_{N^{ha}}^{ha}\} = \{g_e^{ha}\}_{e \in E_c \setminus \{e^*\}} \triangleright$ $N^{ha}$ is equal to $|E_c| - 1$
29: **end if**

## 4.3 Global loops corresponding to “holes”: the set $GL^h$

The global loops that arise from the presence of “holes” are straightforward to construct. Let $\partial K_k$ for $k \in \{1, \dots, N^h\}$ represent the connected components of $\partial K$, and fix one of them, say $\partial K_{N^h}$.

For each connected component $\partial K_k$ with $k \in \{1, \dots, N^h - 1\}$, there corresponds a global loop $g_k^h$. The support of this global loop consists of all edges that have exactly one of their nodes in $\partial K_k$.
**Algorithm 3** Construct the set $GL^h$

1: **for** each connected component $\partial K_k$ of $\partial K$ except $\partial K_{N^h}$ **do**

2: define $\mathbf{c}_k^h$ by

$$ < \mathbf{c}_k^h, n > := \begin{cases} 1 & \text{if } n \in \partial K_k, \\ 0 & \text{if } n \notin \partial K_k. \end{cases} $$

▷ $\mathbf{c}_k^h$ is equal to 1 on each boundary node $n$ of $\partial K_k$ and 0 elsewhere.

3: compute $\mathbf{g}_k^h := \delta^1 \mathbf{c}_k^h$

▷ Apply the coboundary map $\delta^1$ to $\mathbf{c}_k^h$.

4: **end for**

5: **return** $GL^h := \{\mathbf{g}_1^h, \dots, \mathbf{g}_{N^h-1}^h\}$

## 4.4 Global loops corresponding to contacts: the set $GL^c$

The global loops arising from the presence of contacts are also straightforward to construct. Let $\partial K_j^c$ for $j \in \{1, \dots, N^c\}$ denote the contacts of $\partial K^c$, and fix one of them, say $\partial K_{N^c}^c$.

Then, for each contact $\partial K_j^c$ with $j \in \{1, \dots, N^c\}$, there corresponds a global loop $\mathbf{g}_j^c$, whose support consists of all edges whose duals belong to a specific path in the dual graph $\tilde{G}$, associated with the pair of contacts $\partial K_j^c$ and $\partial K_{N^c}^c$ of $\partial K^c$. Moreover, an additional global loop needs to be added depending on the cardinality of a special set defined in Algorithm 2.

**Algorithm 4** Construct the set $GL^c$

1: let $\tilde{T}$ be a spanning tree of $\tilde{G}$

2: select an edge $e$ of $\partial K_{N^c}^c$ and let $f_e$ be the unique face of $K$ incident on it

3: **for** each connected component $\partial K_j^c$ of $\partial K^c$ except $\partial K_{N^c}^c$ **do**

4: select an edge $e_j$ of $\partial K_j^c$

5: let $f_{e_j}$ be the unique face of $K$ incident on $e_j$

6: construct the unique path $\tilde{p}_{e_j,e}$ in $\tilde{T}$ from $D(f_{e_j})$ to $D(f_e)$

7: apply Algorithm 1 with input $\tilde{p}_{e_j,e}$ and $(e_j, e)$ to get $\mathbf{g}_j^c$

8: **end for**

9: let $E_c^{\mathrm{II}}$ be the set computed in Algorithm 2 at line 14 ▷ the set $E_c$ and its subset $E_c^{\mathrm{II}}$ have already been computed in Algorithm 2 (see the order in which the algorithms are called in Algorithm 5)

10: **if** $|E_c^{\mathrm{II}}| = 0$ **then**

11: **return** $GL^c := \{\mathbf{g}_1^c, \dots, \mathbf{g}_{N^c-1}^c\}$

12: **else**

▷ If $|E_c^{\mathrm{II}}| > 0$

13: select an edge $e^* \in E_c^{\mathrm{II}}$

14: let $f_{e^*}^{(1)}, f_{e^*}^{(2)}$ be the unique two faces of $K$ incident on $e^*$

15: construct the unique path $\tilde{p}_{e^*,e^*}^{(1)}$ in $\tilde{T}$ from $D(f_e)$ to $D(f_{e^*}^{(1)})$

16: construct the unique path $\tilde{p}_{e^*,e^*}^{(2)}$ in $\tilde{T}$ from $D(f_e)$ to $D(f_{e^*}^{(2)})$

17: apply Algorithm 1 with input $\tilde{p}_{e^*,e^*}^{(1)}$ and $(e, e^*)$ to get $\mathbf{g}_e^{(1)}$

18: apply Algorithm 1 with input $\tilde{p}_{e^*,e^*}^{(2)}$ and $(e, e^*)$ to get $\mathbf{g}_e^{(2)}$

19: define $\mathbf{g}_{N^c}^c$ by

$$ < \mathbf{g}_{N^c}^c, e' > := \begin{cases} < \mathbf{g}_e^{(1)}, e' > & \text{if } e' = e^*, \\ < \mathbf{g}_e^{(1)}, e' > + < \mathbf{g}_e^{(2)}, e' > & \text{if } e' \neq e^*. \end{cases} $$

20: **return** $GL^c := \{\mathbf{g}_1^c, \dots, \mathbf{g}_{N^c}^c\}$

21: **end if**
## 4.5 The main algorithm to construct the set GL of global loops

The algorithm to compute global loops combines Algorithms 2, 3, and 4 and can be described as a three-step procedure.

<caption>**Algorithm 5** Construct of the set GL</caption>

1: apply Algorithm 2 to get global loops $GL^{ha}$ corresponding to handles
2: **if** $\partial K \neq \emptyset$ **then** ▷ if $K$ has a boundary
3: apply Algorithm 3 to get global loops $GL^h$ corresponding to holes
4: **if** $\partial K^c \neq \emptyset$ **then** ▷ if $K$ has contacts
5: apply Algorithm 4 to get global loops $GL^c$ corresponding to contacts
6: **end if**
7: **end if**
8: **return** $GL := GL^{ha} \sqcup GL^h \sqcup GL^c$ ▷ return the set of all global loops

In the next section we prove the correctness of Algorithm 5.

# 5 Proof of correctness of Algorithm 5

We divide the proof of in nine steps.

## 5.1 Tools from Discrete Morse Theory

To start with, we recall a specific version of Discrete Morse Theory introduced by Kozlov in his book [32].

Let $K$ be a CW complex, and for $m \in \{1, 2\}$, let $x$ be a $m$-cell and $x'$ be a $(m-1)$-cell of $K$. If $\iota(x, x') \neq 0$, we write $x' \prec x$.

**Definition 5.1** (Matching, acyclic matching). A matching $M$ of $m$-cells on $K$ is a set of pairs $(x, x')$ made of an $m$-cell and a $(m-1)$-cell of $K$ such that:

(i) If $(x, x') \in M$, then $x' \prec x$.

(ii) If $(x_1, x'), (x_2, x') \in M$, then $x_1 = x_2$.

When $(x, x') \in M$, we write $x' = d(x)$. A matching $M$ is *acyclic* if there does not exist a cycle of matched elements $x_i$ for $i \in \{1, \dots, l\}$ of the form

$$x_1 \succ d(x_1) \prec x_2 \succ d(x_2) \prec \cdots \prec x_l \succ d(x_l) \prec x_1,$$

with $l \ge 2$.

We say that an $m$-cell or $(m-1)$-cell of $K$ is *matched* in $M$ if it appears as the first or second member of a pair in $M$. The set of *critical cells* of $K$ is set of all cells of $K$ that are not matched in $M$.

The next lemma says that to every spanning tree on $K$ we can associate an acyclic matching.

**Lemma 1** (Acyclic matching from spanning tree). Let $T$ be a spanning tree of $K$. Then, there exist an acyclic matching $M_T$ such that every vertex of $T$, except one (denoted $n$), is matched in $M_T$. In particular, $n$ can be chosen arbitrarily. Similarly, let $\tilde{G} = (\tilde{V}, \tilde{E})$ be the dual graph of $K$ and let $\tilde{T}$ be a spanning tree on $\tilde{G}$. Then, there exist an acyclic matching $M_{\tilde{T}}$ such that every dual vertex of $\tilde{T}$, except one (denoted $\tilde{n}$), is matched in $M_{\tilde{T}}$. In particular, $\tilde{n}$ can be chosen arbitrarily.

*Proof.* The proof of these two results follows directly from the reasoning in [33, Section 4.2] for the cases $k=0$ and $k=2$ (notice that the construction for the case $k=2$ is stated for 3-dimensional domains but it can be easily adapted to 2-dimensional ones). □

The next theorem can be regarded as the main theorem of Discrete Morse Theory for CW complexes, as it highlights the topological role of acyclic matchings.
**Theorem 5.1.** [32, Theorem 11.13 (b) and (c)] Let M be an acyclic matching on K, and denote by $d_m$ the number of critical m-cells of K with respect to M. For every critical m-cell x and (m-1)-cell x' of K, let $S(x, x')$ be the set of all sequences s of matched cells of the form

$$s := (x \succ d(x_1) \prec x_1 \succ d(x_2) \prec x_2 \succ \cdots \succ d(x_l) \prec x_l \succ x').$$

For each $s \in S(x, x')$, define $w(s)$ as the m-chain computed from the sequence s as follows:

$$w(s) := w_l,$$

where, for all $i \in \{0, \dots, l\}$, $w_i$ is computed recursively as

$$\left\{ \begin{array}{l} w_0 := x, \\ w_i := w_{i-1} + \xi_i x_i, \end{array} \right. \qquad (9)$$

with

$$\xi_i := - \frac{\langle \partial_{m+1} w_{i-1}, d(x_i) \rangle}{\langle \partial_{m+1} x_i, d(x_i) \rangle}, \qquad (10)$$

for all $i \in \{1, \dots, l\}$. Then, $K$ is homotopy equivalent to a CW complex $K_c$ with exactly $d_m$ cells of dimension $m$. In particular, there is a natural identification of the cells of $K_c$ with the critical cells of $K$ such that, for any $x, x'$ as above, the incidence number $\iota^M(x, x')$ in $K_c$ is given by

$$\iota^M(x, x') := \sum_{s \in S(x, x')} \langle \partial_m w(s), x' \rangle. \qquad (11)$$

## 5.2 Properties of the 1-cochain computed by Algorithm 1

Let $\tilde{p}$ be a path in $\tilde{G}$, and denote by $f_1, \dots, f_l$ the corresponding set of 2-cells and by $e_1, \dots, e_{l-1}$ the associated set of 1-cells in $K$. Define $M_{\tilde{p}}$ as the set of all pairs $(f_i, e_i)$ for $i \in \{1, \dots, l-1\}$. It is straightforward to verify that $M_{\tilde{p}}$ forms an acyclic matching on $K$, and that $f_l$ is a critical 2-cell in $K$.

Accordingly, for each critical 1-cell $e$ such that $e \subset f_1$, we define the sequence $s_{\tilde{p}}$ associated with the path $\tilde{p}$ in $M_{\tilde{p}}$ as

$$s_{\tilde{p}} := (f_l \succ e_{l-1} \prec f_{l-1} \succ e_{l-2} \prec f_{l-2} \succ \cdots \succ e_1 \prec f_1 \succ e).$$

**Lemma 2** (Properties of 1-cochains of Algorithm 1). Let $\tilde{p}$ be a path in $\tilde{G}$, and let $(e, e')$ be a pair of 1-cells of $K$ such that $e \subset f_1$ and $e' \subset f_l$. Let $g$ be the 1-cochain of $K$ computed by Algorithm 1. If $g \neq 0$, then $\langle g, e \rangle = 1$, and the following holds

$$\langle g, e' \rangle = - \frac{\iota(\xi_l f_l, e')}{\iota(\xi_1 f_1, e)}, \qquad (12)$$

where the coefficients $\xi_1, \xi_l$ are defined as in (10) for the sequence $s_{\tilde{p}}$.

Furthermore, for all $i \in \{1, \dots, l-1\}$, the following equation holds

$$\langle g_e, \partial_2(\xi_i f_i) \rangle = 0. \qquad (13)$$

*Proof.* By the definition of the $g$ in Algorithm 1 at lines 1 and 4, we have that $\langle g, e \rangle = 1$, and for each $i \in \{1, \dots, l\}$, the following holds:

$$\langle g, e_i \rangle = - \frac{\iota(f_i, e_i)}{\iota(f_i, e_{i-1})} \langle g, e_{i-1} \rangle \qquad (14)$$

where we define $e_0 := e$, $e_l := e'$ as the initial and final edges, respectively. This formula expresses the transformation of the 1-cochain as we move through the sequence of 1-cells in $\tilde{p}$.

Note that (14) is independent of the specific orientation of each face $f_i$. In fact, for any coefficient $\chi_i \in \{-1, +1\}$ corresponding to an oriented face $\chi_i f_i$, we have $\iota(\chi_i f_i, e_i) = \chi_i \iota(f_i, e_i)$.
This shows that the expression in (14) remains unchanged under the orientation change of the faces.

Next, we multiply each face $f_i$ by the coefficient $\xi_i$, which is computed according to (10) using the sequence $s_{\tilde{p}}$. By applying (14) recursively, and using the fact that $\iota(\xi_i f_i, e_i) = -\iota(\xi_i f_i, e_{i-1})$ for all $i \in \{1, \dots, l\}$ (this follows from the specific choice of coefficients $\xi_i$ in (10), as well as from the recursion formula (9)), we obtain the desired result in (12).

To prove (13), we rewrite (14) as the following equation:

$$< g, e_i > \iota(\xi_i f_i, e_i) + < g, e_{i-1} > \iota(\xi_i f_i, e_{i-1}) = 0,$$

for each $i \in \{1, \dots, l-1\}$. This is simply a restatement of the transformation rule for the 1-cochain.

Finally, recalling that $g = 0$ as initialized in Algorithm 1 at line 1, and noting that each $f_i$ is distinct since $\tilde{p}$ is a path, we conclude that this leads to the desired result in (13). □

## 5.3 Definition of a new CW complex $\bar{\mathcal{K}}$ with empty boundary and its reduction to a homotopy-equivalent CW complex $\bar{\mathcal{K}}_c$

We first define the CW complex $\bar{\mathcal{K}}$, constructed from $\mathcal{K}$ as follows. For each connected component of the boundary $\partial \mathcal{K}_k$ with $k \in \{1, \dots, N^h\}$, we do the following operation: we consider a gluing map (in the CW sense as in Definition 2.1) attaching a closed disk $D_k$ (i.e. a 2-ball) along the (polygonal) boundary $\partial \mathcal{K}_k$ and we contract $D_k$ it to a single 0-cell, say $n_k$. Accordingly, the CW complex structure of $\bar{\mathcal{K}}$ is obtained from that of $\mathcal{K}$ by removing every 1-cell of $\partial \mathcal{K}$ and by identifying, for each $k \in \{1, \dots, N^h\}$, every 0-cell in $\partial \mathcal{K}_k$ with the 0-cell $n_k$ in $\bar{\mathcal{K}}$. Notice that a 1-cocycle $g$ of $\mathcal{K}$ is also a 1-cocycle of $\bar{\mathcal{K}}$ if and only if its support does not contain any 1-cell in $\partial \mathcal{K}$.

Let $T$ and $\tilde{T}$ be the spanning trees of $\mathcal{K}$ and its dual graph $\tilde{G}$, respectively, as defined in Algorithm 2, lines 1-7. These trees also induce corresponding spanning trees of $\bar{\mathcal{K}}$ and its dual graph $\tilde{G}$, which we denote by the same symbols for consistency.

Let $M_T$ and $M_{\tilde{T}}$ be the acyclic matchings constructed from the spanning trees $T$ and $\tilde{T}$, respectively, by applying Lemma 1 to $\bar{\mathcal{K}}$. Applying Theorem 5.1 to $\bar{\mathcal{K}}$ and the acyclic matching $M = M_{\tilde{T}} \cup M_T$, we obtain a CW complex $\bar{\mathcal{K}}_c$ that is homotopy equivalent to $\bar{\mathcal{K}}$, with its cells naturally indexed by the set of critical cells of $\bar{\mathcal{K}}$ with respect to $M$.

Since $T$ is a spanning tree of $\bar{\mathcal{K}}$, each edge of $T$ is matched in $M_T$ to a vertex (a 0-cell) of $\bar{\mathcal{K}}$. As a result, $M_T$ consists of $|\bar{\mathcal{K}}^{(0)}| - 1$ pairs. Similarly, since $\tilde{T}$ is a spanning tree of $\tilde{G}$, each edge of $\tilde{T}$ is matched in $M_{\tilde{T}}$ to a vertex (a 2-cell, by duality) of $\tilde{G}$. Therefore, $M_{\tilde{T}}$ consists of $|\tilde{V}| - 1 = |\bar{\mathcal{K}}^{(2)}| - 1$ pairs.

Consequently, $\bar{\mathcal{K}}_c$ consists of:

* One 0-cell, which is the unique non-matched vertex of $\bar{\mathcal{K}}$ in $M_T$.

* $|\bar{\mathcal{K}}^{(1)}| - (|\bar{\mathcal{K}}^{(0)}| - 1) - (|\bar{\mathcal{K}}^{(2)}| - 1) = 2 - \chi$ 1-cells, where $\chi$ is the Euler characteristic of $\bar{\mathcal{K}}$.

* One 2-cell, which is the unique non-matched vertex (again, by duality) of $\tilde{G}$ in $M_{\tilde{T}}$.

It is important to note that the set of all 1-cells of $\bar{\mathcal{K}}_c$ is in one-to-one correspondence with the set $E_c$, as defined in Algorithm 2 at line 8. For the remainder of this paper, we will identify these two sets.

Now, let $f$ denote the unique 2-cell of $\bar{\mathcal{K}}_c$, and consider the expression for the boundary map $\partial_2^M$ on the CW complex $\bar{\mathcal{K}}_c$ with incidence numbers $\iota^M(f, e)$, as given by (11). Specifically, we have

$$\partial_2^M f = \sum_{e \in E_c} \iota^M(f, e) e, \qquad (15)$$

where $\iota^M(f, e)$ is the incidence number associated with the boundary of $f$ along the 1-cell $e$.

Since $\bar{\mathcal{K}}_c$ is homotopy equivalent to $\bar{\mathcal{K}}$ (by Theorem 5.1) and the boundary of $\bar{\mathcal{K}}$ is empty, we conclude that each coefficient $\iota^M(f, e)$ in the above sum can only take one of two possible values:

$$1. \iota^M(f, e) = 0.$$
2. $\iota^M(f, e) = 2\eta_e$, where $\eta_e \in \{-1, +1\}$.

Therefore, we can partition $E_c$ into two subsets: $E_c = E_c^I \sqcup E_c^{II}$, where $E_c^I$ and $E_c^{II}$ consist of the edges $e \in E_c$ that satisfy conditions (1) and (2), respectively.

## 5.4 Differentiation between homology and torsion generators of $H_1(\bar{\mathcal{K}})$ passing through $\bar{\mathcal{K}}_c$

The characterization of the homology of $\bar{\mathcal{K}}$ is an immediate consequence of the Classification Theorem of Compact Surfaces [17]. In the next Lemma we provide an equivalent characterization which will be useful for our specific purposes.

**Lemma 3 (Homology characterization of $\bar{\mathcal{K}}$).** For each $e \in E_c$, let $c_e$ be the unique cycle in the graph $T \cup \{e\}$, naturally oriented according to the orientation induced by $e$. Note that $c_e$ can considered as a 1-cycle in $\bar{\mathcal{K}}$. Let $E_c = E_c^I \sqcup E_c^{II}$, and denote by $g := |E_c| = 2 - \chi$ the genus of $\bar{\mathcal{K}}_c$. We distinguish between the following two cases:

(i) If $|E_c^{II}| = 0$, then $H_1(\bar{\mathcal{K}}; \mathbb{Z}) = \mathbb{Z}^g$.

(ii) If $|E_c^{II}| > 0$, then $H_1(\bar{\mathcal{K}}; \mathbb{Z}) = \mathbb{Z}^{g-1} \oplus \mathbb{Z}/2\mathbb{Z}$. In this case, the torsion generator $\bar{c}^*$ is given by

$$ \bar{c}^* := \sum_{e \in E_c^{II}} \eta_e c_e, \qquad (16) $$

where each $\eta_e \in \{-1, +1\}$ is the sign of the coefficient $\iota^M(f, e)$ in the expression of $\partial_2^M f$ as a linear combination of 1-cells $e \in E_c^{II}$ in (15).

*Proof.* Since $\bar{\mathcal{K}}_c$ has only one 0-cell, every 1-cell $e \in E_c$ forms a 1-cycle.

In the first case, where $|E_c^{II}| = 0$, we have $E_c = E_c^I$. Therefore, each 1-cell $e \in E_c$ is a homology generator of $H_1(\bar{\mathcal{K}}_c; \mathbb{Z})$, as $\partial_2^M f \neq e$. It follows that each 1-cycle $c_e$ is a homology generator of $H_1(\bar{\mathcal{K}}; \mathbb{Z})$.

In the second case, where $|E_c^{II}| > 0$, fix a 1-cell $e^* \in E_c^{II}$. The boundary operator $\partial_2^M f$ satisfies

$$ \partial_2^M f = \sum_{e \in E_c^{II}} 2\eta_e e := 2c^*, $$

where $c^*$ is the (unique) torsion generator of $\bar{\mathcal{K}}_c$. For each $e \in E_c \setminus \{e^*\}$, the 1-cell $e$ is a homology generator of $H_1(\bar{\mathcal{K}}_c; \mathbb{Z})$, since $\partial_2^M f \neq e$. This means there are $g-1$ independent homology generators.

Next, let $c_e$ the unique 1-cycle of $\bar{\mathcal{K}}$ associated with each $e \in E_c$, as in the case where $|E_c^{II}| = 0$, and define the 1-cycle $\bar{c}^*$ as in (16). Using the fact that $\bar{\mathcal{K}}$ is homotopy equivalent to $\bar{\mathcal{K}}_c$, it follows that each $c_e$ for $e \in E_c \setminus \{e^*\}$ is a homology generator of $H_1(\bar{\mathcal{K}}; \mathbb{Z})$. Furthermore, $\bar{c}^*$ is the torsion generator of $H_1(\bar{\mathcal{K}}; \mathbb{Z})$, since $\bar{c}^*$ is homotopy equivalent to the torsion 1-cycle $c^*$ of $\bar{\mathcal{K}}_c$. $\square$

## 5.5 Characterization of generators of $H^1(\bar{\mathcal{K}})$

Let us remind the Universal Coefficient Theorem for cohomology.

**Theorem 5.2 ([29, Theorem 3.2]).** If a chain complex $C$ of free Abelian groups has homology groups $H_n(C; G)$, then the cohomology groups $H^n(C; G)$ of the cochain complex $\text{Hom}(C; G)$ are determined by the exact sequence

$$ 0 \to \text{Ext}(H_{n-1}(C); G) \to H^n(C; G) \xrightarrow{h} \text{Hom}(H_n(C); G) \to 0. $$

For our purposes, we do not need to delve into the definition of the Ext functor. It suffices to note the following property: $\text{Ext}(Q; G) = 0$ when $Q$ is a free group. For further details and a proof of this property, consult [29].

We now apply Theorem 5.2 to $\bar{\mathcal{K}}$ with $n = 1$ and $G = \mathbb{Z}$. In Theorem 5.2, the map $h: H^1(\bar{\mathcal{K}}; \mathbb{Z}) \to \text{Hom}(H_1(\bar{\mathcal{K}}); \mathbb{Z})$ assigns to each generator $g \in H^1(\bar{\mathcal{K}}; \mathbb{Z})$ the function $\langle g, \cdot \rangle$, which maps an element $c \in H_1(\bar{\mathcal{K}}; \mathbb{Z})$ to the value $\langle g, c \rangle$.
<figure><img src="image_3.png" /><figcaption>Figure 2: Illustration of the notation used in the proof of Point (ii) in Lemma 4 and in Algorithm 2.</figcaption></figure>

Note that if $\mathrm{Ext}(H_0(\bar{K}); \mathbb{Z}) = 0$, then $h$ is an isomorphism. Since we assume $\mathcal{K}$ is connected, $\bar{K}$ is also connected, implying that $H_0(\bar{K}; \mathbb{Z}) = \mathbb{Z}$. Consequently, $\mathrm{Ext}(H_0(\bar{K}); \mathbb{R}) = 0$, and from the exactness of the sequence in Theorem 5.2, we conclude that $h: H^1(\bar{K}; \mathbb{Z}) \to \mathrm{Hom}(H_1(\bar{K}); \mathbb{Z})$ is an isomorphism.

As a result, we obtain the following characterization of generators of $H^1(\bar{K}; \mathbb{R})$ with real coefficients in terms of a set of generators of $H_1(\bar{K}; \mathbb{Z})$.

**Corollary 1 (Characterization of cohomology generators).** Let $\{c_i\}_{i=1}^{\beta_1}$ be a set of homology generators of $H_1(\bar{K}; \mathbb{Z})$, and let $\{g_i\}_{i=1}^{\beta_1}$ be a set of 1-cycles of $\bar{K}$ such that

$$<g_i, c_j> = \delta_{i,j}, \qquad (17)$$

for all $i, j \in \{1, \dots, \beta_1\}$, where $\delta_{i,j}$ is the Kronecker delta. Then, the set $\{g_i\}_{i=1}^{\beta_1}$ forms a basis of $H^1(\bar{K}; \mathbb{Z})$. In particular, the same set yields a basis for $H^1(\bar{K}; \mathbb{R})$.

## 5.6 A first decomposition of global loops and proof of correctness of Algorithm 2

Recall the definition of the spanning trees $T$ and $\tilde{T}$ in Algorithm 2, at lines 1–7. It follows that every path $\tilde{p}$ in $\tilde{G}$ does not contain any boundary edge of $\mathcal{K}$. Accordingly, the support of each 1-cochain $g$ of $\mathcal{K}$ computed via Algorithm 1 does not contain any boundary edge of $\mathcal{K}$ either (recall the formula in Algorithm 1 at line 4). Hence, if $g$ is a 1-cocycle, then it is also 1-cocycle of $\bar{K}$.

**Lemma 4 (Properties of 1-cochains of Algorithm 2).** Let $e \in E_c$. The following two assertions hold:

(i) Let $\tilde{p}_e$ be the path defined in Algorithm 2 at lines 10 and 11, and let $g_e^{ha}$ be the corresponding 1-cochain of $\bar{K}$ computed at line 12 using Algorithm 1. Then, $e \in E_c^1$ if and only if $g_e^{ha}$ is a 1-cocycle of $\bar{K}$ whose restriction to $\bar{K}_c$ satisfies the condition

$$<g_e, e' > = \delta_{e,e'} \qquad (18)$$

for every $e' \in E_c$.
(ii) Assume $|E_c^{II}| > 0$, and let $e^*$ the 1-cell selected in Algorithm 2 at line 18. For each $e \in E_c^{II} \setminus \{e^*\}$, let $\tilde{p}_{e,e^*}$, $\tilde{p}_{e,e^*}^{(1)}$, $\tilde{p}_{e,e^*}^{(2)}$ be the paths defined in Algorithm 2, at lines 20-23. Let $g_e^{(1)}$ and $g_e^{(2)}$ denote the corresponding 1-cochains of $\bar{K}$ computed at lines 24 and 25 using Algorithm 1, and let $g_e^{ha}$ be the 1-cochain defined at line 26. Then, $g_e^{ha}$ is a 1-cocycle of $\bar{K}$ whose restriction to $\bar{K}_c$ satisfies the condition

$$< g_e, e' > = \delta_{e,e'} \quad (19)$$

for every $e' \in E_c^{II} \setminus \{e^*\}$.

*Proof.* (i) Let $e \in E_c$. Denote by $f_1, \dots, f_l$ and $e_1, \dots, e_{l-1}$ the set of 2-cells and 1-cells of $\tilde{p}_e$, respectively. Notice that we can identify the unique 2-cell $f$ of $\bar{K}_c$ by setting $f := f_l$. This follows from the fact that every 2-cell of $\bar{K}$, except one, is matched in the acyclic matching $M$, and that the unique non-matched 2-cell in $M$ can be chosen arbitrarily (as discussed in the property of $M_{\tilde{T}}$ in Lemma 1). This also matches with the natural identification of cells of $\bar{K}_c$ with the critical cells of $\bar{K}$ as stated in Theorem 5.1.

Let $s_{\tilde{p}_e}$ be the sequence associated with the path $\tilde{p}_e$ in $M$. Applying Point (ii) in Theorem 5.1 to the pair $f, e$, we observe that $s_{\tilde{p}_e}$ is the unique sequence in the set $S(f, e)$. This is because $f_1, f_l$ are the only 2-cells of $\bar{K}$ incident on $e$, with $f := f_l$.

Next, applying formula (11), we conclude that $e \in E_c^I$ if and only if

$$\iota^M(f, e) = \iota(\xi_1 f_1, e) + \iota(\xi_l f_l, e) = 0.$$

Thus, we obtain the following condition

$$\iota(\xi_1 f_1, e) = -\iota(\xi_l f_l, e). \quad (20)$$

By Lemma 2, we know that the 1-cochain $g_e^{ha}$ satisfies $< g_e^{ha}, e > = 1$ and $< g_e^{ha}, \partial_2(\xi_l f_l) > = 0$ for all $i \in \{1, \dots, l-1\}$. Since $e' = e$ in Algorithm 2 at line 12, it follows that $g_e^{ha}$ is a 1-cocycle if and only if $< g_e^{ha}, \partial_2(\xi_l f_l) > = 0$, which leads to the compatibility condition

$$1 = < g, e > = < g, e' >.$$

Using (12) from Lemma 2, we conclude that this condition is exactly (20).

(ii) Let $f_1^{(1)}, \dots, f_l^{(1)}$ and $e_1^{(1)}, \dots, e_{l-1}^{(1)}$ be the sets of 2-cells and 1-cells of $\tilde{p}_{e,e^*}^{(1)}$, and let $f_1^{(2)}, \dots, f_q^{(2)}$ and $e_1^{(2)}, \dots, e_{q-1}^{(2)}$ be the corresponding sets for $\tilde{p}_{e,e^*}^{(2)}$. See Fig. 2 for reference.

Let $\tilde{p}$ be the unique path in $\tilde{T}$ from $D(f_1^{(2)})$ to $D(f_l^{(1)})$ and let $s_{\tilde{p}}$ be the corresponding sequence in $M$. As in Point (i), we can identify $f := f_l^{(1)}$.

Apply Point (ii) of Theorem 5.1 to the pair $f, e$, and observe that $s_{\tilde{p}}$ and $s_{\tilde{p}_{e,e^*}^{(1)}}$ are the only two sequences in the set $S(f, e)$, since $f_1^{(1)}, f_1^{(2)}$ are the unique 2-cells of $\bar{K}$ incident on $e$. Then, applying formula (11) to the sequences $s_{\tilde{p}_{e,e^*}^{(1)}}$ and $s_{\tilde{p}}$, we have

$$\iota^M(f, e) = \iota(\xi_1^{(1)} f_1^{(1)}, e) + \iota(\xi_1^{(2)} f_1^{(2)}, e).$$

Since $e \in E_c^{II}$, by the definition of $E_c^{II}$ in Section 5.4, we have $\iota^M(f, e) = 2\eta_e$ with $\eta_e \in \{-1, 1\}$. Moreover, since both 2-cells $f_1^{(1)}, f_1^{(2)}$ are also 2-cells of $\bar{K}$, we deduce that $\iota(\xi_1^{(1)} f_1^{(1)}, e)$, $\iota(\xi_1^{(2)} f_1^{(2)}, e) \in \{-1, 1\}$. Therefore, we conclude that

$$\iota(\xi_1^{(1)} f_1^{(1)}, e) = \iota(\xi_1^{(2)} f_1^{(2)}, e) = \eta_e.$$

Similarly, let $\tilde{p}$ be the unique path in $\tilde{T}$ from $D(f_q^{(2)})$ to $D(f_l^{(1)})$. Apply again Point (ii) in Theorem 5.1 to the pair $f, e^*$, and observe that $s_{\tilde{p}}$ is the unique sequence of the set $S(f, e)$, since $f_l^{(1)}$ and $f_q^{(2)}$ are the only 2-cells of $\bar{K}$ incident on $e$.

Then, applying formula (11) to the sequence $s_{\tilde{p}}$, we have

$$\iota^M(f, e^*) = \iota(\xi_l^{(1)} f_l^{(1)}, e^*) + \iota(\xi_q^{(2)} f_q^{(2)}, e^*).$$
Since $e^* \in E_c^{II}$, we conclude, using similar reasoning as above, that

$$ \iota(\xi_l^{(1)} f_l^{(1)}, e^*) = \iota(\xi_q^{(2)} f_q^{(2)}, e^*) = \eta_{e^*} $$

with $\eta_{e^*} \in \{-1, 1\}$.

By (13) in Lemma 2, we have that $\langle \mathbf{g}_e^{(1)}, \partial_2(\xi_i^{(1)} f_i^{(1)}) \rangle = 0$ for all $i \in \{1, \dots, l-1\}$, and $\langle \mathbf{g}_e^{(2)}, \partial_2(\xi_i^{(2)} f_i^{(2)}) \rangle = 0$ for all $i \in \{1, \dots, q-1\}$. Recalling the definition of $\mathbf{g}_e^{ha}$ in Algorithm 2 at line 26, it follows that $\mathbf{g}_e^{ha}$ is a 1-cocycle if and only if $\langle \mathbf{g}_e^{(1)}, \partial_2(\xi_l^{(1)} f_l^{(1)}) \rangle = 0$ and $\langle \mathbf{g}_e^{(2)}, \partial_2(\xi_q^{(2)} f_q^{(2)}) \rangle = 0$. These conditions lead to the compatibility condition

$$ \langle \mathbf{g}_e^{(1)}, e^* \rangle = \langle \mathbf{g}_e^{(2)}, e^* \rangle, $$

along with $\langle \mathbf{g}_e^{(1)}, e \rangle = \langle \mathbf{g}_e^{(2)}, e \rangle$, which holds by the initialization in Algorithm 1 at line 2.

Applying formula (12) in Lemma 2 to the sequences $s_{\tilde{P}_{e,e^*}^{(1)}}$, $s_{\tilde{P}_{e,e^*}^{(2)}}$, we deduce that

$$ \langle \mathbf{g}_e^{(1)}, e^* \rangle = - \frac{\iota(\xi_l^{(1)} f_l^{(1)}, e^*)}{\iota(\xi_1^{(1)} f_1^{(1)}, e)} = - \frac{\eta_{e^*}}{\eta_e} = - \frac{\iota(\xi_q^{(2)} f_q^{(2)}, e^*)}{\iota(\xi_1^{(2)} f_1^{(2)}, e)} = \langle \mathbf{g}_e^{(2)}, e^* \rangle. $$

Thus, $\mathbf{g}_e^{ha}$ is a 1-cocycle, and its restriction to $\bar{\mathcal{K}}_c$ satisfies the desired condition, as stated in (19). □

**Lemma 5** (First decomposition of global loops). Let $GL^{ha} = \{\mathbf{g}_i^{ha}\}_{i=1}^{|GL^{ha}|}$ be the set of 1-cochains computed by Algorithm 2, and let $\langle [GL^{ha}] \rangle := \{\sum_{i=1}^{|GL^{ha}|} \gamma_i[\mathbf{g}_i^{ha}] | \gamma_i \in \mathbb{R}\}$ be the linear span of the equivalence classes associated with each generator $\mathbf{g}_i^{ha} \in GL^{ha}$. Then, the following isomorphism holds

$$ H^1(\bar{\mathcal{K}}) \cong \langle [GL^{ha}] \rangle. $$

*Proof.* For each $e \in E_c$, let $\mathbf{g}_e^{ha}$ be the 1-cocycle of $\bar{\mathcal{K}}$ computed in Algorithm 2 at line 12. Applying Point (i) in Lemma 4, it follows that, if $\mathbf{g}_e^{ha} \neq \mathbf{0}$, then $\mathbf{g}_e^{ha}$ is such that its restriction to $\bar{\mathcal{K}}_c$ satisfies (18), and $e \in E_c^I$. On the other hand, if $\mathbf{g}_e^{ha} = \mathbf{0}$, then $e \in E_c^{II}$. This defines the partition $E_c = E_c^I \cup E_c^{II}$, and we now analyze the two possible cases for this partition.

In the case where $|E_c^{II}| = 0$, observe that by construction, the support of $\mathbf{g}_e^{ha}$ does not include any edge of $T$. Consequently, for each $e' \in E_c$, the equality

$$ \langle \mathbf{g}_e^{ha}, c_{e'} \rangle = \langle \mathbf{g}_e^{ha}, e' \rangle = \delta_{e,e'} $$

holds, thanks to (18). Applying Corollary 1 in combination with Lemma 3, we deduce that $GL^{ha} = \{\mathbf{g}_e^{ha}\}_{e \in E_c}$ forms a set of generators of $H^1(\bar{\mathcal{K}})$.

In the case where $|E_c^{II}| > 0$, let $e^* \in E_c^{II}$ denote the 1-cell fixed in Algorithm 2 at line 18. For each $e \in E_c^{II} \setminus \{e^*\}$, let $\mathbf{g}_e^{ha}$ be the corresponding 1-cocycle of $\bar{\mathcal{K}}$ computed in Algorithm 1 at line 26. By Point (ii) of Lemma 4, each $\mathbf{g}_e^{ha}$ for $e \in E_c^{II} \setminus \{e^*\}$ is 1-cocycle of $\bar{\mathcal{K}}$ whose restriction to $\bar{\mathcal{K}}_c$ satisfies (19). Since the support of $\mathbf{g}_e^{ha}$ does not contain any edge of $T$, it follows that, for all $e' \in E_c \setminus \{e^*\}$, we have

$$ \langle \mathbf{g}_e^{ha}, c_{e'} \rangle = \langle \mathbf{g}_e^{ha}, e' \rangle = \delta_{e,e'}, $$

thanks to (19), as well as $\langle \mathbf{g}_e^{ha}, \bar{c}^* \rangle = \langle \mathbf{g}_e^{ha}, c^* \rangle$.

Using Corollary 1 and Lemma 3, we conclude that $GL^{ha} = \{\mathbf{g}_e^{ha}\}_{e \in E_c \setminus \{e^*\}}$ is a set of generators of $H^1(\bar{\mathcal{K}})$.

Thus, in both cases, the set $GL^{ha}$ provides a valid set of generators of $H^1(\bar{\mathcal{K}})$, determined by analyzing the contributions of 1-cocycles based on the partition $E_c = E_c^I \cup E_c^{II}$. □

## 5.7 A second decomposition of global loops and proof of correctness of Algorithm 3

We now focus on $H^1(\bar{\mathcal{K}}, \partial\bar{\mathcal{K}})$ and present a preparatory result for the second decomposition of global loops, which is associated with the connected components of the boundary $\partial\bar{\mathcal{K}}$.
**Lemma 6** (Dimension calculation via long exact sequences). Let $\mathcal{K}$ a 2-dimensional CW complex with a non-empty boundary $\partial\mathcal{K}$. Then, the relationship between the ranks of $H^1(\mathcal{K}, \partial\mathcal{K})$ and $H^1(\bar{\mathcal{K}})$ is given by the following expression

$$ \mathrm{rank}(H^1(\mathcal{K}, \partial\mathcal{K})) - \mathrm{rank}(H^1(\bar{\mathcal{K}})) = N^h - 1. \qquad (21) $$

*Proof.* The proof consists in expressing the ranks of $H^1(\mathcal{K}, \partial\mathcal{K})$ and $H^1(\bar{\mathcal{K}})$ in terms of that $H^1(\mathcal{K})$ and in computing the ranks of some elementary cohomology spaces.

We start by considering the Mayer-Victoris long exact sequence [34, Theorem 3, Chapter 11] for the triad $(\bar{\mathcal{K}}, \mathcal{K}, D)$, where $D$ is defined as $D := \bigcup_{k=1}^{N^h} D_k$. The exact sequence is as follows:

$$ \begin{aligned} 0 &\to \underbrace{H^0(\bar{\mathcal{K}})}_{1} \to \underbrace{H^0(\mathcal{K})}_{1} \oplus \underbrace{H^0(D)}_{N^h} \to \underbrace{H^0(\mathcal{K} \cap D)}_{N^h} \to \\ &\to H^1(\bar{\mathcal{K}}) \to H^1(\mathcal{K}) \oplus \underbrace{H^1(D)}_{0} \to \underbrace{H^1(\mathcal{K} \cap D)}_{N^h} \to \\ &\to H^2(\bar{\mathcal{K}}) \to \underbrace{H^2(\mathcal{K})}_{0} \oplus \underbrace{H^2(D)}_{0} \to \underbrace{H^2(\mathcal{K} \cap D)}_{0} \to 0, \end{aligned} \qquad (22) $$

where the ranks of the corresponding cohomology groups are indicated under the braces.

Next, we apply the long exact sequence in relative cohomology [34, Theorem 13] for the pair $(\partial\mathcal{K}, \mathcal{K})$

$$ \begin{aligned} 0 &\to \underbrace{H^0(\mathcal{K}, \partial\mathcal{K})}_{0} \to \underbrace{H^0(\mathcal{K})}_{1} \to \underbrace{H^0(\partial\mathcal{K})}_{N^h} \to \\ &\to H^1(\mathcal{K}, \partial\mathcal{K}) \to H^1(\mathcal{K}) \to \underbrace{H^1(\partial\mathcal{K})}_{N^h} \to \\ &\to H^2(\mathcal{K}, \partial\mathcal{K}) \to \underbrace{H^2(\mathcal{K})}_{0} \to \underbrace{H^2(\partial\mathcal{K})}_{0} \to 0, \end{aligned} \qquad (23) $$

where the ranks of the cohomology groups are indicated under the braces. In particular, $N^h$ denotes the number of connected components of $\partial\mathcal{K}$.

Finally, applying the Euler characteristic formula [29, Theorem 2.44] to the exact sequences (22) and (23), and noting that in (23), the rank of $H^2(\mathcal{K}, \partial\mathcal{K})$ equals that of $H^2(\bar{\mathcal{K}})$ in (22), regardless of whether $\mathcal{K}$ is orientable, we obtain the desired formula (21). $\square$

We are now ready to prove the second decomposition of global loops, which demonstrates the correctness of Algorithm 3.

**Lemma 7** (Second decomposition of global loops). Let $GL^h = \{g_k^h\}_{k=1}^{[GL^h]}$ be the set of 1-cochains computed by Algorithm 3, and define $\langle[GL^h]\rangle := \{\sum_{k=1}^{[GL^h]} \beta_k[g_k^h] | \beta_k \in \mathbb{R}\}$ as the linear span of the equivalence classes associated with each generator $g_k^h \in GL^h$. Then, the following group isomorphism holds

$$ H^1(\mathcal{K}, \partial\mathcal{K}) \cong H^1(\bar{\mathcal{K}}) \oplus \langle[GL^h]\rangle. \qquad (24) $$

*Proof.* We begin by recalling the definition of $g_k^h$ in Algorithm 3, specifically at lines 2 and 3. Each $g_k^h \in GL^h$ is a 1-cocycle, and it follows that $g_k^h \notin B^0(\mathcal{K}, \partial\mathcal{K})$. Therefore, we have that $\langle[GL^h]\rangle \subset H^1(\mathcal{K}, \partial\mathcal{K})$.

Additionally, since $g_k^h \in B^1(\bar{\mathcal{K}})$, we also deduce that $\langle[GL^h]\rangle \cap H^1(\bar{\mathcal{K}}) = \{0\}$.

Now, let $Z \subset H^1(\mathcal{K}, \partial\mathcal{K})$ be a subspace such that $H^1(\mathcal{K}, \partial\mathcal{K}) \cong H^1(\bar{\mathcal{K}}) \oplus Z$. By Lemma 6, we conclude that the rank of $Z$ is given by

$$ \mathrm{rank}(Z) = \mathrm{rank}(H^1(\mathcal{K}, \partial\mathcal{K})) - \mathrm{rank}(H^1(\bar{\mathcal{K}})) = N^h - 1. \qquad (25) $$

Note that $\langle[GL^h]\rangle \subset Z$. Since the set $GL^h$ consists of 1-cocycles $g_k^h$ corresponding to the $k$-th connected component $\partial\mathcal{K}_k$ of $\partial\mathcal{K}$, it is linearly independent. Therefore, the set $GL^h$ contains $N^h - 1$ independent 1-cocycles, and thus it forms a basis for $Z$.

Consequently, we conclude that $\langle[GL^h]\rangle = Z$, which gives us the decomposition Lemma 5. $\square$
## 5.8 A third decomposition of global loops and proof of correctness of Algorithm 4

We now discuss a preparatory result similarly to the one in Lemma 6.

**Lemma 8** (Dimension computation via long exact sequences). Let $\mathcal{K}$ a 2-dimensional CW complex with a non-empty boundary $\partial\mathcal{K}$. Recalling that $\partial\mathcal{K}^{nc} := \partial\mathcal{K} \cup \bigcup_{j=1}^{N^c} \partial\mathcal{K}_j^c$, where each $\partial\mathcal{K}_j^c$ is a contact region, the following expression relates the ranks of $H^1(\mathcal{K}, \partial\mathcal{K}^{nc})$ and $H^1(\mathcal{K}, \partial\mathcal{K})$:

$$ \mathrm{rank}(H^1(\mathcal{K}, \partial\mathcal{K}^{nc})) - \mathrm{rank}(H^1(\mathcal{K}, \partial\mathcal{K})) = \begin{cases} N^c - 1 & \text{if } \mathcal{K} \text{ is orientable,} \\ N^c & \text{otherwise.} \end{cases} \qquad (26) $$

*Proof.* We proceed in a similar manner as in the proof of Lemma 6, by expressing the ranks of $H^1(\mathcal{K}, \partial\mathcal{K}^{nc})$ and $H^1(\mathcal{K}, \partial\mathcal{K})$ in terms of $H^1(\mathcal{K})$.

First, let us consider the long exact sequence in relative cohomology [34, Theorem 12] for the pair $(\partial\mathcal{K}^{nc}, \mathcal{K})$:

$$ \begin{aligned} 0 &\to \underbrace{H^0(\mathcal{K}, \partial\mathcal{K}^{nc})}_{0} \to \underbrace{H^0(\mathcal{K})}_{1} \to \underbrace{H^0(\partial\mathcal{K}^{nc})}_{N^{nc}} \to \\ &\to H^1(\mathcal{K}, \partial\mathcal{K}^{nc}) \to H^1(\mathcal{K}) \to H^1(\partial\mathcal{K}^{nc}) \to \\ &\to \underbrace{H^2(\mathcal{K}, \partial\mathcal{K}^{nc})}_{0} \to \underbrace{H^2(\mathcal{K})}_{0} \to \underbrace{H^2(\partial\mathcal{K}^{nc})}_{0} \to 0, \end{aligned} \qquad (27) $$

where $N^{nc}$ represents the number of connected components of $\partial\mathcal{K}^{nc}$.

Next, let us note the following two facts:

1. In the exact sequence for $(\partial\mathcal{K}, \mathcal{K})$ given in (22) of Lemma 6, the rank of $H^2(\mathcal{K}, \partial\mathcal{K})$ is 1 if $\mathcal{K}$ is orientable and 0 if $\mathcal{K}$ is non-orientable.

2. In the exact sequence for $(\partial\mathcal{K}^{nc}, \mathcal{K})$ in (23), the rank of $H^1(\partial\mathcal{K}^{nc})$ satisfies the formula

$$ \mathrm{rank}(H^1(\partial\mathcal{K}^{nc})) = N^{nc} - N^c, \qquad (28) $$

which will now we prove.

Consider the set of all connected components $\partial\mathcal{K}_k$ of $\partial\mathcal{K}$, where $k \in \{1, \dots, N^h\}$. After potentially rearranging the indices, we assume that there is a number $p \in \{1, \dots, N^h\}$ such that the first $p$ connected components, $\{\partial\mathcal{K}_1, \dots, \partial\mathcal{K}_p\}$, contain some contact regions $\partial\mathcal{K}_j^c$ for $j \in \{1, \dots, N^c\}$, while the remaining components $\{\partial\mathcal{K}_{p+1}, \dots, \partial\mathcal{K}_{N^h}\}$ do not. If $p = N^h$, then the second set is empty.

For each $k \in \{1, \dots, N^h\}$, define $J_k^c$ as the set of indices $j \in \{1, \dots, N^c\}$ such that $\partial\mathcal{K}_j^c \subset \partial\mathcal{K}_k$. Let $|J_k^c|$ denote the cardinality of $J_k^c$, and define $N_k^{nc}$ as the number of connected components of the set $\partial\mathcal{K}_k \cup \bigcup_{j \in J_k^c} \partial\mathcal{K}_j^c$.

It is easy to check that $|J_k^c| = N_k^{nc}$ for all $k \in \{1, \dots, p\}$, and that $|J_k^c| = 0$ and $N_k^{nc} = 1$ for all $k \in \{p+1, \dots, N^h\}$. Hence, we have the relation

$$ N^h - p = \sum_{k=p+1}^{N^h} N_k^{nc}. $$

Also, since $\mathrm{rank}(H^1(\partial\mathcal{K}^{nc})) = N^h - p$, we obtain the formula:

$$ N^{nc} = \sum_{k=1}^{N^h} N_k^{nc} = \sum_{k=1}^{p} N_k^{nc} + \sum_{k=p+1}^{N^h} N_k^{nc} = N^c + \mathrm{rank}(H^1(\partial\mathcal{K}^{nc})), $$

which is the desired result (28).

Finally, applying the Euler characteristic formula [29, Theorem 2.44] to the exact sequences (27) and (22) in Lemma 6, and using Facts 1 and 2, we derive the formula (26). □
We are now ready to state the third decomposition of global loops, which also demonstrates the correctness of Algorithm 4.

**Lemma 9** (Third decomposition of global loops). Let $GL^c = \{g_j^c\}_{j=1}^{|GL^c|}$ be the set of 1-cochains computed by Algorithm 4, and define $\langle[GL^c]\rangle := \{\sum_{j=1}^{|GL^c|} \alpha_j[g_j^c] | \alpha_j \in \mathbb{R}\}$ as the linear span of the equivalence classes associated with each generator $g_j^c \in GL^c$. Then, the following group isomorphism holds

$$H^1(K, \partial K^{nc}) \cong H^1(K, \partial K) \oplus \langle[GL^c]\rangle. \qquad (29)$$

*Proof.* Let $E_c^{II}$ be the set defined in Algorithm 2 at line 14, and used in Algorithm 4 at line 9. We distinguish between the two cases $|E_c^{II}| = 0$ and $|E_c^{II}| > 0$.

If $|E_c^{II}| = 0$, then in Algorithm 4 each 1-cochain $g_j^c \in GL^c$ is computed using Algorithm 1 with input pairs $(e, e')$ where $e, e' \in \partial K$. Applying formulas (12) and (13) from Lemma 2, we conclude that $g_j^c$ is 1-cocycle that is nonzero on exactly one 1-cell in the boundary $\partial K_j^c \subset \partial K$.

This ensures that $g_j^c \notin B^0(K, \partial K^{nc})$, since otherwise it would be nonzero on more than 1-cell in $\partial K_j^c$, contradicting the earlier conclusion. Thus, we obtain the inclusion $\langle[GL^c]\rangle \subset H^1(K, \partial K^{nc})$.

Furthermore, since each $g_j^c$ is nonzero on exactly one 1-cell $\partial K_j^c$, it follows that $H^1(K, \partial K) \cap \langle[GL^c]\rangle = \{0\}$.

Next, consider a subspace $Y \subset H^1(K, \partial K^{nc})$ such that

$$H^1(K, \partial K^{nc}) \cong H^1(K, \partial K) \oplus Y.$$

Using Point (i) of Lemma 3, we deduce that $\bar{K}$, and consequently $K$, is orientable. This follows from the fact that $\bar{K}$ is obtained by gluing disks along the boundary of $K$. See [17] for a discussion of this property. Then, applying Lemma 8, we obtain

$$\mathrm{rank}(Y) = \mathrm{rank}(H^1(K, \partial K^{nc})) - \mathrm{rank}(H^1(K, \partial K)) = N^c - 1. \qquad (30)$$

Since $\langle[GL^c]\rangle \subset Y$ and the set $GL^c$ is linearly independent, where each $g_j^c$ is nonzero on exactly one 1-cell of the $j$-th connected component of $\partial K^{nc}$, it follows that $GL^c$ consists of $N^c - 1$ independent 1-cocycles, forming a basis of $Y$. Thus, we conclude that $\langle[GL^c]\rangle = Y$, leading to the decomposition Lemma 5.

On the other hand, if $|E_c^{II}| > 0$, we consider the 1-cochains $g_e^{(1)}$ and $g_e^{(2)}$ defined in Algorithm 4 at lines 17 and 18. Applying the same reasoning as in the proof of Point (ii) of Lemma 4, using the paths $\tilde{p}_{e,e^*}^{(1)}$ and $\tilde{p}_{e,e^*}^{(2)}$ defined at lines 15 and 16, we deduce that both 1-cochains $g_e^{(1)}$ and $g_e^{(2)}$ attain the same value on the 1-cell $e^*$, chosen at line 13. This implies that the 1-cochain $g_{N_c}^c$ defined at line 19, is a 1-cocycle. Using Point (ii) of Lemma 3, we deduce that in this case, $\bar{K}$, and hence $K$, must be non-orientable. Repeating the same reasoning as above for the case $|E_c^{II}| = 0$ and applying Lemma 8, we obtain the desired conclusion. $\square$

## 5.9 Final step

To conclude, it is sufficient to combine Lemmas 5, 7, and 9.

**Lemma 10** (Correctness of Algorithm 5). Let $GL^{ha}, GL^h, GL^c$ denote the sets of global loops computed by Algorithms 2, 3, and 4, respectively. Then, we have the following isomorphism

$$H^1(K, \partial K^{nc}) \cong \langle[GL^{ha}]\rangle \oplus \langle[GL^h]\rangle \oplus \langle[GL^c]\rangle.$$

## 6 Conclusions

In this final section, we provide an overview of the algorithm's key features, its impact on practical applications, and potential extensions that could guide future work.

The first aspect to highlight is the algorithm's ability to automatically compute global loops for surfaces, both with and without boundaries, and for orientable as well as non-orientable surfaces. The primary algorithm, presented in Section 4.5, utilizes a three-step procedure to partition global
loops into distinct disjoint classes. Informally, these classes correspond to the presence of “handles,” “holes,” or contact regions. The algorithm handles all these cases automatically, adapting to each situation to produce the correct output. This is a crucial feature for applications where manually verifying the explicit form and count of global loops is challenging. Furthermore, this reliability is underpinned by the robust theoretical analysis provided earlier in the paper.

For example, consider the case where $\Gamma$ is a triangulated Möbius strip. Our algorithm outputs no global loop, understood as a generator of $H^1(\Gamma, \partial\Gamma)$. This result is easily validated theoretically via Lemma 7. Since the Möbius strip has only one connected component, $N^h = 1$, it follows that $H^1(\mathcal{K}, \partial\mathcal{K}) \cong H^1(\bar{\mathcal{K}})$. However, $\bar{\mathcal{K}}$, obtained by gluing a disk to the Möbius strip, is homeomorphic to a projective plane (cf. [17, pg. 19]). From the homological characterization of the projective plane, we know that $H^1(\bar{\mathcal{K}}) \cong 0$. This observation corrects an error in [6], where it is claimed that the Möbius strip supports a global loop, as stated explicitly in the caption of Fig. 2 in Section V. This figure refers to their proposed procedure for computing global loops of the surface under consideration. To avoid such pitfalls and ensure provably correct results, the algorithm proposed in our work is therefore explicitly required.

The second crucial aspect of the algorithm is its linear time worst-case complexity. Specifically, the algorithms described in Section 4 require only a list-based representation of a triangulation $\mathcal{K}$ of the surface $\Gamma$. In this representation, each cell in $\mathcal{K}$ has an associated data structure that stores its neighboring cells. In practical applications, it is generally assumed that the number of neighboring cells for each cell in $\mathcal{K}$ is bounded by a constant, reflecting common mesh regularity assumptions in boundary element methods. Considering this constant, the use of standard graph algorithms for spanning trees results in a worst-case linear time complexity [23]. Additionally, all algorithms in Section 4 are entirely graph-based, which makes them straightforward to implement without requiring external libraries or specialized procedures.

Another consideration is that, while the algorithm performs efficiently in practice, it does not aim to produce a minimal or shortest set of cohomology generators. This means that the generated set of cohomology generators may not have the most compact support (e.g., minimal length or the smallest number of edges). For further insights, readers may consult Erickson et al. [24] for homology generators and Dłotko et al. [14] for cohomology generators, which are particularly relevant to our applications. Achieving a truly minimal basis is often computationally prohibitive, necessitating a trade-off between runtime and the quality of the basis obtained.

Finally, we emphasize that the algorithm's foundation in Discrete Morse Theory ensures its broad applicability. The algorithms presented here are effective as long as an *optimal discrete Morse function* can be computed. For example, such functions can be determined even for pathological cases, such as surfaces connected by a single vertex or edge [25]. Additionally, these ideas can extend to handle *non-manifold surfaces* (e.g., screens with junction points, as discussed in [26, 27]). We propose that slight modifications to our algorithms could extend their applicability to these cases while maintaining linear worst-case complexity. For instance, identifying manifold parts can be achieved in linear time by analyzing the coboundary elements of each edge. By orienting a random polygon in each manifold part and iteratively orienting others using a *breadth-first search* [23], we expect the algorithm's linear time complexity to remain intact.

# Acknowledgements

The author acknowledges the funding provided by the European Union through the MSCA EffECT, project number 101146324.

# Conflict of interest

The author declares no competing interests.

# References

[1] A. Buffa and S. H. Christiansen, “The electric field integral equation on lipschitz screens: definitions and numerical approximation,” *Numerische Mathematik*, vol. 94, pp. 229–267, 2003.
[2] S. B. Adrian, A. Dely, D. Consoli, A. Merlini, and F. P. Andriulli, "Electromagnetic integral equations: Insights in conditioning and preconditioning," *IEEE Open Journal of Antennas and Propagation*, vol. 2, pp. 1143-1174, 2021.

[3] G. Vecchi, "Loop-star decomposition of basis functions in the discretization of the efie," *IEEE Transactions on Antennas and Propagation*, vol. 47, no. 2, pp. 339-346, 1999.

[4] F. P. Andriulli, "Loop-star and loop-tree decompositions: Analysis and efficient algorithms," *IEEE Transactions on Antennas and Propagation*, vol. 60, no. 5, pp. 2347-2356, 2012.

[5] F. P. Andriulli, K. Cools, I. Bogaert, and E. Michielssen, "On a well-conditioned electric field integral operator for multiply connected geometries," *IEEE transactions on antennas and propagation*, vol. 61, no. 4, pp. 2077-2087, 2012.

[6] B. Hofmann, T. F. Eibert, F. P. Andriulli, and S. B. Adrian, "A low-frequency stable, excitation agnostic discretization of the right-hand side for the electric field integral equation on multiply-connected geometries," *IEEE Transactions on Antennas and Propagation*, vol. 71, no. 12, pp. 9277-9288, 2023.

[7] R. Hiptmair and J. Ostrowski, "Generators of $H_1(\Gamma_h, \mathbb{Z})$ for triangulated surfaces: Construction and classification," *SIAM Journal on Computing*, vol. 31, no. 5, pp. 1405-1423, 2002.

[8] D. Eppstein, "Dynamic generators of topologically embedded graphs," *ArXiv*, vol. cs.DS/0207082, 2002.

[9] F. Lazarus, M. Pocchiola, G. Vegter, and A. Verroust, "Computing a canonical polygonal schema of an orientable triangulated surface," in *Proceedings of the seventeenth annual symposium on Computational geometry*, pp. 80-89, 2001.

[10] T. K. Dey, F. Fan, and Y. Wang, "An efficient computation of handle and tunnel loops via reeb graphs," *ACM Transactions on Graphics (TOG)*, vol. 32, no. 4, pp. 1-10, 2013.

[11] P. W. Gross and P. R. Kotiuga, *Electromagnetic theory and computation: a topological approach*, vol. 48. Cambridge University Press, 2004.

[12] P. Dlotko, "A fast algorithm to compute cohomology group generators of orientable 2-manifolds," *Pattern Recognition Letters*, vol. 33, no. 11, pp. 1468-1476, 2012.

[13] P. Dlotko and R. Specogna, "Physics inspired algorithms for (co) homology computations of three-dimensional combinatorial manifolds with boundary," *Computer Physics Communications*, vol. 184, no. 10, pp. 2257-2266, 2013.

[14] P. Dlotko and R. Specogna, "Lazy cohomology generators: A breakthrough in (co) homology computations for cem," *IEEE transactions on magnetics*, vol. 50, no. 2, pp. 577-580, 2014.

[15] D. Wilton and A. Glisson, "On improving the stability of the electric field integral equation at low frequency," in *Proc. IEEE Antennas and Propagation Soc. National Symp*, pp. 124-133, 1981.

[16] D. N. Arnold, *Finite element exterior calculus*. SIAM, 2018.

[17] J. H. Gallier, D. Xu, et al., *A guide to the classification theorem for compact surfaces*. Springer, 2013.

[18] G. Miano and F. Villone, "A surface integral formulation of maxwell equations for topologically complex conducting domains," *IEEE transactions on antennas and propagation*, vol. 53, no. 12, pp. 4001-4014, 2005.

[19] Y. Wang, D. Gope, V. Jandhyala, and C.-J. Shi, "Generalized kirchoff's current and voltage law formulation for coupled circuit-electromagnetic simulation with surface integral equations," *IEEE transactions on microwave theory and techniques*, vol. 52, no. 7, pp. 1673-1682, 2004.

[20] S. Suuriniemi, J. Kangas, L. Kettunen, and T. Tarhasaari, "State variables for coupled circuit-field problems," *IEEE transactions on magnetics*, vol. 40, no. 2, pp. 949-952, 2004.

[21] S. Suuriniemi, J. Kangas, and L. Kettunen, "Driving a coupled field-circuit problem," *COMPEL-The international journal for computation and mathematics in electrical and electronic engineering*, vol. 26, no. 3, pp. 899-909, 2007.

[22] R. Hiptmair and J. Ostrowski, "Electromagnetic port boundary conditions: Topological and variational perspective," *International Journal of Numerical Modelling: Electronic Networks, Devices and Fields*, vol. 34, no. 3, p. e2839, 2021.

[23] T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein, *Introduction to algorithms*. MIT press, 2022.

[24] J. Erickson and K. Whittlesey, "Greedy optimal homotopy and homology generators," in *SODA*, vol. 5, pp. 1038-1046, 2005.

[25] T. Lewiner, H. Lopes, and G. Tavares, "Optimal discrete morse functions for 2-manifolds," *Comput. Geom.*, vol. 26, pp. 221-233, 2003.

[26] X. Claeys and R. Hiptmair, "Integral equations on multi-screens," *Integral equations and operator theory*, vol. 77, no. 2, pp. 167-197, 2013.

[27] X. Claeys and R. Hiptmair, "Integral equations for electromagnetic scattering at multi-screens," *Integral Equations and Operator Theory*, vol. 84, pp. 33-68, 2016.

[28] G. E. Cooke and R. L. Finney, *Homology of cell complexes*, vol. 2239. Princeton University Press, 2015.

[29] A. Hatcher, "Algebraic topology," 2001.

[30] J. Dodziuk, "Finite-difference approach to the hodge theory of harmonic forms," *American Journal of Mathematics*, vol. 98, no. 1, pp. 79-104, 1976.
[31] G. F. Duff, "Differential forms in manifolds with boundary," *Annals of Mathematics*, vol. 56, no. 1, pp. 115-127, 1952.

[32] D. Kozlov, "Combinatorial algebraic topology," in *Algorithms and computation in mathematics*, 2008.

[33] S. Pitassi, R. Ghiloni, and R. Specogna, "Inverting the discrete curl operator: a novel graph algorithm to find a vector potential of a given vector field," *Journal of Computational Physics*, vol. 466, p. 111404, 2022.

[34] M. Spivak, "A comprehensive introduction to differential geometry, publish or perish," *Inc., Berkeley*, vol. 2, 1979.