\section{Introduction}
This manuscript establishes a general approach to formulate partial differential equations (PDEs) on networks of (hyper)surfaces, referred to as hypergraphs. Such PDEs consist of differential expressions with respect to all hyperedges (surfaces) and compatibility conditions on the hypernodes (joints, intersections of surfaces). These compatibility conditions ensure conservation properties (in case of continuity equations) or incorporate other properties---motivated by physical or mathematical modeling.  We illuminate how to discretize such equations numerically using hybrid discontinuous Galerkin (HDG) methods, which appear to be a natural choice, since they consist of local solvers (encoding the differential expressions on hyperedges) and a global compatibility condition (related to our hypernode conditions). We complement the physically motivated compatibility conditions by a derivation through a singular limit analysis of thinning structures yielding the same results.

Albeit many physical, sociological, engineering, and economic processes have been described by partial differential equations posed on domains which cannot be described as subsets of linear space or smooth manifolds, there is still a lack of mathematical tools and general purpose software specifically addressing the challenges arising from the discretization of these models.

Fractured porous media (see \cite{BerreDK2019} for a comprehensive review) have gained substantial attention and have become an active field of research due to their critical role with respect to flow patterns in several applications in the subsurface, in material science, and in biology. Most commonly, a fracture is described as a very thin, not necessarily planar object in which, for example Darcy's equation holds. This motivates the singular limit approximation in which a fracture is assumed to be a two dimensional \emph{surface} within the three dimensional \emph{space}. When several of these fractures meet, they form a fracture network of two-dimensional surfaces. Thus, fracture networks illustrate a physical application of the type of problem we investigate in this publication.
Moreover, a model in which the joints of two (or more) fractures are assigned additional physical properties can be found in \cite{ReichenbergerJBH2006}. Beyond this, fracture networks have been simulated using hybrid high order (HHO) methods \cite{HedinPE2019}.

Graph based models for porous media (without fractures) consist of simulating preferential flow paths within the porous matrix. One of the first publications implementing this idea is \cite{Fatt1956} who observed that a network of tubes might approximate the flow of porous media better than the classical model of tube bundles, which has also been used in the most common upscaling techniques---see \cite{SchulzRZRK2019,RayRSK2018} and the references therein for a discussion of those tube models in upscaling procedures. The tube network approach \cite{Fatt1956} has been successfully applied to couple porous media flow to free (Navier--) Stokes flow \cite{WeishauptJH2019}.

PDEs on hypergraphs are especially suitable to be used in the description of elastic networks \cite{Eremeyev2019}: Here, we discriminate between one dimensional elastic beam (rod) networks, trusses, etc.\ and two dimensional elastic plate (shell) networks \cite{LagneseL1993}, respectively. Beam networks have been used to model truss bridges and towers (most prominently the Eiffel tower) and other mechanical structures, originating the field of elastic beam theory (for instance \cite{BauchauC2009} for an introduction which also covers elastic plate models). Elastic plate models describe the stability of houses and have several engineering applications such as the description of the stability of (bend) plates (used in automobile industries and several others). They have even been used to understand interseismic surface deformation at subduction zones \cite{KandaS2010}.

Elastic beam networks have been used to evaluate elastic constants in amorphous materials. That is, the elastic properties of stiff, beam like polymers have been investigated. Such polymers are key to understanding the cytoskeleton which is an important part of biological cells \cite{Heussinger2007,Lieleg2007}, but they are also important for the healing of wounds (fibrin), for skin stability (collagen), and for the properties of paper. Moreover, such models can be used for modelling rubber \cite{PhysRevE.76.031906}, foams, and fiber networks \cite{PhysRevLett.96.017802}. 

Conservation laws in the form of PDEs on hypergraphs have been used in the simulation and optimization of gas networks \cite{RuefflerMH2018} and other networks of pipelines. They have been extended to networks of traffic (streets and data), (tele-)communication, and blood flow. For an overview of the main ideas that are related to these applications, the reader may consult \cite{Garavello2010,BressanCGHP2014}. Additionally, rigorous mathematical analysis of such problems is developing to a field of current research \cite{SkreMR2021}.

We conclude the overview over some applications by stating that regular surfaces and volumes can also be interpreted as hypergraphs. Thus, PDEs on surfaces \cite{DziukE2013} and standard ``volume" problems (in which the hyperedges have the same dimension as the surrounding space and at most two hyperedges meet in a common hypernode) are also covered by our approach.

Hypergraph models usually are approximations of problems in higher dimensional networks of thin structures, for example a network of thin pipes or thin plates in 3D. As a model example we give a rigorous  derivation of a diffusion equation on a hypergraph. More precisely, we consider a network of thin plates in three dimensions, where the thickness of the plates is small compared to their length. We denote the ratio between the thickness and the length by the small parameter $0 < \epsilon \ll 1$. Due to the different scales the computational effort for numerical simulations is very high. To overcome this problem the idea is to replace the thin-structure by a hypergraph. For this we give a rigorous mathematical justification using asymptotic analysis. We pass to the limit $\epsilon \to 0$ in the weak formulation of the problem, and derive a limit problem stated on the hypergraph. The solution of this limit-problem is an approximation of the model in the higher-dimensional thin domain. Singular limits for thin plates and shells (leading to lower-dimensional manifolds in the limit $\epsilon \to 0$) in elesticity can be found in \cite{ciarlet1997mathematical,ciarlet2000theory}. Dimension reduction for a folded elastic plate is treated in \cite{le1989folded}.
Singular limits leading to hypergraphs for fluid equations can be found in \cite{maruvsic2003rigorous}, where a Kirchhoff law in a junction of thin pipes is derived, and \cite{maruvsic2019mathematical} where junctions of thin pipes and plates are treated using the method of two-scale convergence.

The remainder of this manuscript is structured as follows: First, we discuss conservation equations on hypergraphs. Second, we rigorously formulate an elliptic model equation and investigate some of its properties in Section \ref{SEC:model_eq}. Third, we discuss its discretization by means of the HDG method in Section \ref{SEC:hdg_graph}. Fourth, we discuss how PDEs on hypergpahs can be obtained by a model reduction approach, in particular, by considering singular limits. The publication is wrapped up, by a section on possible conclusions.

\section{Conservation equations on geometric hypergraphs}
\subsection{Hypergraphs}
A hypergraph $\ensuremath{{\mathcal G}}=(\ensuremath{{\mathcal N}},\ensuremath{{\mathcal E}})$ consists of a finite set $\ensuremath{{\mathcal E}}$ of hyperedges and a finite set $\ensuremath{{\mathcal N}}$ of hypernodes. We refer to it as a \emph{geometric hypergraph} if the hyperedges are smooth, open manifolds of dimension $\dimEdge$ with piecewise smooth, Lipschitz boundary and the hypernodes can be identified with smooth subsets of the boundaries of these hyperedges. More specifically, the boundary of each hyperedge $\ensuremath{E}_e\in \ensuremath{{\mathcal E}}$ is subdivided into $k_e$ nonoverlapping subsets $\Gamma^e_{i}$ such that $\partial\ensuremath{E}_e = \bigcup \overline{\Gamma^e_{i}}$. We associate to $\ensuremath{E}_e$ an index vector $\ensuremath{\eta}^e_1,\dots,\ensuremath{\eta}^e_{k_e}$ and isometries
\begin{gather}
 \iota^e_i\colon \Gamma^e_{i} \to \ensuremath{N}_{\ensuremath{\eta}^e_i}, \qquad i=1,\dots,k_e.
\end{gather}
The hypernodes are thus identified with the closures of the subsets of the boundaries of one or more hyperedges. Their dimension is $\dimEdge-1$.

$\ensuremath{{\mathcal G}}$ has the structure of a hypergraph in the classical sense as each edge $\ensuremath{E}_e\in\ensuremath{{\mathcal E}}$ connects a set of nodes $\{\ensuremath{N}_{\ensuremath{\eta}^e_1},\dots,\ensuremath{N}_{\ensuremath{\eta}^e_{k_e}}\}\subset \ensuremath{{\mathcal N}}$. The dual hypergraph $\ensuremath{{\mathcal G}}^*=(\ensuremath{{\mathcal E}},\ensuremath{{\mathcal N}})$ describes the situation where each hypernode $\ensuremath{N}_n \in \ensuremath{{\mathcal N}}$ connects $\ell_n$ hyperedges with indices $\ensuremath{\nu}^n_1,\dots,\ensuremath{\nu}^n_{\ell_n}$.

We call a hypernode $\ensuremath{N}_n \in \ensuremath{{\mathcal N}}$ a boundary hypernode if $\ell_n=1$, i.e., it is part of the boundary of only a single hyperedge.
Accordingly, we define the set of boundary hypernodes $\ensuremath{\SetNode_\textup B}$ and the set of interior hypernodes $\ensuremath{\SetNode_\textup I} = \ensuremath{{\mathcal N}} \setminus \ensuremath{\SetNode_\textup B}$.
\begin{figure}
 \includegraphics[width=.8\textwidth]{pictures/hypernode}
 \caption{A hypergraph with three hyperedges of dimension 2, a hypernode (red) connecting them, and 9 boundary hypernodes (blue). Embedded hypergraph (left) in $\ensuremath{\mathbb R}^3$ and without embedding (right). Isometries $\iota$ only shown for the interior hypernode.}\label{fig:hypergraph1}
\end{figure}

As special cases: a geometric graph is a geometric hypergraph where the edges are smooth curves and the nodes are their end points. If every hypernode is either at the boundary or connects exactly two hyperedges, the hypergraph represents a piecewise smooth manifold.

The structure might become more evident if we consider an embedded geometric hypergraph in some ambient space $\ensuremath{\mathbb R}^\ensuremath{\mathfrak D}$, as in Figure~\ref{fig:hypergraph1} on the left. In this case, the isometries $\iota^e_i$ are identical mappings and the hypernodes are identified with the boundary pieces $\Gamma^e_i$. On the right of this figure, the same hypergraph is displayed without embedding. In this case, the hyperedges are objects in $\ensuremath{\mathbb R}^2$, possibly with a non-flat metric. Hypernodes are intervals in $\ensuremath{\mathbb R}$, inheriting their metric through the isometries $\iota$.

Due to the isometries $\iota$, every point of a hypernode $\ensuremath{N}$ is uniquely identified with a point on the boundary of each of the hyperedges it connects. Thus, convergence of a point sequence in the union of these hyperedges to a point on the hypernode is well-defined, for instance by considering the (finitely many) subsequences on each hyperedge. Also, a distance between two points on different hyperedges sharing a hypernode is defined locally by these isometries and triangle inequality.

The domain $\Omega$ of the hypergraph, its closure, and its boundary are 
\begin{gather}\begin{split}
 \overline \Omega &= \bigcup_{\ensuremath{E}\in\ensuremath{{\mathcal E}}} \ensuremath{E} \cup \bigcup_{\ensuremath{N}\in\ensuremath{{\mathcal N}}} \ensuremath{N},
 \qquad \quad
 \partial\Omega = \bigcup_{\ensuremath{N}\in\ensuremath{\SetNode_\textup B}} \ensuremath{N},\\
  \Omega &= \overline \Omega \setminus \partial \Omega.
 \label{EQ:def_omega}
\end{split}
\end{gather}
In this definition, the hyperedges are considered open with respect to their topology and do not contain their boundaries.
The hypernodes are closed. We introduce the skeletal domain
\begin{gather}
    \Sigma = \bigcup_{\ensuremath{N}\in\ensuremath{{\mathcal N}}} \ensuremath{N}.
\end{gather}
We make the assumption that $\Omega$ is connected. Note that this implies that any two hyperedges are either connected by a common node or not connected, since $\Omega$ is open, see \eqref{EQ:def_omega}.
Without such an assumption, the problems of partial differential equations below separate into subproblems, which then can be analyzed and solved independently.

In \autoref{fig:hypergraph1}, $\partial \Omega$ comprises all blue hypernodes, which also include the end points of the red hypernode. The union of the red and blue hypernodes is $\Sigma$. The domain $\Omega$ consists of the interior of the red hypernode and the the three hyperedges.

Many concepts of standard domains in $\ensuremath{\mathbb R}^\dimEdge$ transfer to $\Omega$, even if it is not a manifold. In particular, the notion of a small open ball $B_r(x)$ with radius $r > 0$ around $x \in \Omega$, see \autoref{fig:hg_ball}, in $\Omega$ is maintained by construction and thus the notion of open subsets. A subset is called compactly embedded in $\Omega$ if its closure is contained in $\Omega$ and thus has a positive distance to $\partial\Omega$.
\begin{figure}
 \includegraphics[width=.8\textwidth]{pictures/ball}
 \caption{The hypergraph of \autoref{fig:hypergraph1} with the illustration of an open ball. The ball's center is located on the shared hypernode, and its radius is $r$.}\label{fig:hg_ball}
\end{figure}

A function is continuous on $\Omega$, if it is continuous inside each hyperedge and its limits on a hypernode are consistent between all hyperedges connected by this hypernode. Analogously, a function is in $L^2(\Omega)$ if it is in $L^2(\ensuremath{E})$ for all $\ensuremath{E} \in \ensuremath{{\mathcal E}}$ and it is in $L^2(\Sigma)$ if it is in $L^2(\ensuremath{N})$ for all $\ensuremath{N} \in \ensuremath{{\mathcal N}}$.
\begin{remark}[Comparison to standard nomenclatures]
 In this article, we mix concepts from graph theory, partial differential equations, and finite elements. Thus, a clash of names was unavoidable. What is referred to as a hypernode here, is a face ---an edge in two dimensions--- in finite element literature, while the hyperedges here correspond to mesh cells or elements. In order to reduce ensuing confusion, we consistently use the term ``hyperedge''. Another difference to finite element literature is established by the fact that we consider the hypergraph fixed and are not concerned with refinement limits. Finally, we would like to point out that there has been a concept of geometric hypergraphs in the literature; it is nevertheless very limited, such that we coin this term in a new way here, meaning a hypergraph whose elements are geometric shapes themselves.
\end{remark}
\subsection{Continuity equations on hypergraphs}\label{SEC:cont_eq}
Next, we conduct a heuristic derivation, employing control volumes $V$ in the shape of infinitesimal, open hyperballs.
Let $\rho$ be a conserved quantity and $\vec J$ be its flux. Then, the conservation property of $\rho$ is usually stated in integral form such that for any such control volume $V$ there holds
\begin{gather}
    \label{EQ:conservation-1}
    \frac{d}{dt} \int_V \rho\ensuremath{\, \textup d x} = - \int_{\partial V} \vec J\cdot\mathbf n \ensuremath{\, \textup d \sigma}.
\end{gather}
When $V$ is a subset of a Lipschitz manifold $\ensuremath{E}$, the meaning of this statement is clear if $\vec J$ is a smooth tangential vector field in $\bar \ensuremath{E}$ and $\mathbf n$ is the outer normal vector to $\partial V$ in the tangential plane of $\ensuremath{E}$. The term $\ensuremath{\, \textup d x}$ denotes the volume element of the manifold, and $\ensuremath{\, \textup d \sigma}$ is the induced surface element.

If the hyperball $V$ intersects a hypernode $\ensuremath{N}$ in which several edges meet, meaning can be given to equation~\eqref{EQ:conservation-1} by the following observation: if $\ensuremath{E}_1,\dots,\ensuremath{E}_\ell$ are the hyperedges which meet in $\ensuremath{N}$ inside $V$, then for $i\in\{1,\dots,\ell\}$ the intersection $V_i = V \cap \ensuremath{E}_i$ has a piecewise smooth boundary $\partial V_i$. We observe that $\ensuremath{N}_V = V\cap\ensuremath{N}$ is in the interior of $V$ (see Figure~\ref{fig:hg_ball} for an illustration) and the boundary of $V$ is nowhere tangential to $\ensuremath{N}$. Thus, with the assumption that no mass is created or destroyed in the hypernode $\ensuremath{N}$, the conservation property~\eqref{EQ:conservation-1} can be restated as
\begin{gather}
 \frac{d}{dt} \int\limits_V \rho\ensuremath{\, \textup d x}
 = \frac{d}{dt}\sum_{i=1}^\ell \int\limits_{V_i} \rho\ensuremath{\, \textup d x}
 = - \sum_{i=1}^\ell \int\limits_{\partial V_i\setminus \ensuremath{N}_V} \vec J\cdot\mathbf n \ensuremath{\, \textup d \sigma}.\label{EQ:conservation-2}
\end{gather}
Again, the flux $\vec J$ and the outer normal vector $\mathbf n$ to $V$ are well defined along $\partial V_i$ in the tangential plane of $\ensuremath{E}_i$.

As a generalization of \eqref{EQ:conservation-2}, we allow for sinks and sources $f$ living in the hyperedges and $g$ living within hypernode $\ensuremath{N}$: This can be implemented by setting
\begin{equation}
 \frac{d}{dt} \int\limits_V \rho \ensuremath{\, \textup d x} = - \sum_{i=1}^\ell \int\limits_{\partial V_i \setminus \ensuremath{N}_V} \vec J\cdot\mathbf n \ensuremath{\, \textup d \sigma} 
 + \int\limits_V f \ensuremath{\, \textup d x} + \int\limits_{\ensuremath{N}_V} g \ensuremath{\, \textup d \sigma},\label{EQ:balance_cond}
\end{equation}
where positive $f$ and $g$ describe sources, while negative $f$ and $g$ describes sinks.

Before we convert~\eqref{EQ:conservation-1} into a problem of partial differential equations, we make the simplifying assumption that the hyperedges and hypernodes are planar and that $\ensuremath{\, \textup d x}$ is the standard Lebesgue measure.
This way, we avoid delving into the complexities of surface partial differential equations. This simplification is purely for the ease of presentation and we refer the readers to \cite{DziukE2013} and \cite{BENARTZI2007989} for more general surfaces in the elliptic and hyperbolic settings, respectively.

Thus, in the interior of each hyperedge $\ensuremath{E}$, we can apply Gauss' divergence theorem in standard form to obtain
\begin{gather}
    \int_{\partial V} \vec J \cdot\mathbf n\ensuremath{\, \textup d \sigma}
    = \int_V \div \vec J \ensuremath{\, \textup d x}.
\end{gather}
If on the other hand $V$ overlaps a hypernode $\ensuremath{N}$ which connects hyperedges $\ensuremath{E}_1\dots,\ensuremath{E}_\ell$, we can still apply the divergence theorem in each hyperedge to obtain
\begin{gather}
    \int_{\partial V} \vec J \cdot\mathbf n\ensuremath{\, \textup d \sigma}
    = \sum_{i=1}^{\ell} \left[
      \int_{V_i}\div \vec J \ensuremath{\, \textup d x}
      - \int_{\ensuremath{N}_V}
      \vec J\cdot\mathbf n\ensuremath{\, \textup d \sigma}
      \right].
\end{gather}
This notion also extends to control volumes $V$ intersecting with several hypernodes in a natural way. Then, rearranging the sum over boundary integrals yields
\begin{gather}
    \label{EQ:conservation-3}
    \int_{\partial V} \vec J \cdot\mathbf n\ensuremath{\, \textup d \sigma}
    = \sum_{\ensuremath{E}\in\ensuremath{{\mathcal E}}} \int\limits_{V\cap \ensuremath{E}}\div \vec J \ensuremath{\, \textup d x}
    - \sum_{\ensuremath{N}\in\ensuremath{\SetNode_\textup I}}
    \int\limits_{V \cap \ensuremath{N}} \jump{\vec J\cdot \mathbf n} \ensuremath{\, \textup d \sigma},
\end{gather}
for any control volume $V\subset\Omega$.
Here, $\div \vec J$ is the standard divergence of the differentiable vector field and $\jump{\cdot}$ is the summation operator such that on a hypernode $\ensuremath{N}$ with hyperedges $\ensuremath{E}_1,\dots,\ensuremath{E}_\ell$ there holds
\begin{gather}
    \label{EQ:divergence_node}
    \jump{\vec J\cdot \mathbf n}
    = \sum_{i=1}^{\ell} \vec J_{|\ensuremath{E}_i}\cdot\mathbf n_{\ensuremath{E}_i}.
\end{gather}
A vector field $\vec J$ is usually called solenoidal, if the left side of equation~\eqref{EQ:conservation-3} vanishes for any control volume $V$. The right hand side of this equation generalizes this notion from standard domains to hypergraphs.
Therefore, we call a piecewise smooth vector field $\vec J$ solenoidal, if
\begin{subequations}
\begin{align}
    \div \vec J(x) &= 0 && \text{ for all } x\in \ensuremath{E} \text{ and } \ensuremath{E} \in \ensuremath{{\mathcal E}},\\
    \label{EQ:Kirchhoff-node}
    \jump{\vec J\cdot \mathbf n}(x) &= 0& & \text{ for all } x\in \ensuremath{N} \text{ and } \ensuremath{N} \in \ensuremath{\SetNode_\textup I}.
\end{align}
\end{subequations}
Note that the second condition is an extension of Kirchhoff's junction rule from points to higher dimensional hypernodes.

To put it in a nutshell, assuming there are no leaks and sources in hypernodes and hyperedges, the conservation condition~\eqref{EQ:conservation-2} induces the PDE--interface problem to find $\rho$ and $\vec J$ such that
\begin{subequations}\label{EQ:continuity}
\begin{align}
 \partial_t \rho + \div \vec J &= 0 && \text{ in all } \ensuremath{E} \in \ensuremath{{\mathcal E}},\\
 \jump{\vec J\cdot \mathbf n} &= 0& & \text{ on all } \ensuremath{N} \in \ensuremath{\SetNode_\textup I}.
\end{align}
\end{subequations}

Analogously, the continuity condition \eqref{EQ:balance_cond} induces the PDE interface problem to find $\rho$ and $\vec J$ such that
\begin{subequations}\label{EQ:balance}
\begin{align}
 \partial_t \rho + \div \vec J &= f && \text{ in all } \ensuremath{E} \in \ensuremath{{\mathcal E}},\\
 \jump{\vec J\cdot \mathbf n} &= g & & \text{ on all } \ensuremath{N} \in \ensuremath{\SetNode_\textup I}.
\end{align}
\end{subequations}

In \eqref{EQ:continuity} and \eqref{EQ:balance}, $\rho$ and $\vec J$ might be linked by some phenomenological description, i.e., $\vec J = \vec J(\rho)$ (depending on the specific application). Both equations are complemented by appropriate initial and boundary conditions. Beyond this, additional continuity constraints might be formulated, such as $\rho \in C(\Omega)$, $\rho \in C^\infty(\bigcup \ensuremath{E})$, \ldots.
\begin{remark}\label{REM:hybrid_formulation}
 The interface problems in \eqref{EQ:continuity} and \eqref{EQ:balance} resemble the \emph{hybrid} or \emph{hybridized} formulation of a PDE, which was introduced for instance in context of the mixed elements of Raviart--Thomas and Brezzi--Douglas--Marini in \cite{RaviartT1977a,RaviartT1977b,BrezziDM1985}.
\end{remark}
We are, exemplary, going to discuss the continuity equation~\eqref{EQ:continuity} in the context of diffusion
problems in section \ref{SEC:model_eq}.
\section{Elliptic model equation}\label{SEC:model_eq}
The standard diffusion equation in mixed form defined on a hypergraph $\ensuremath{{\mathcal G}} = (\ensuremath{{\mathcal E}}, \ensuremath{{\mathcal N}})$ is a conservation equation of type~\eqref{EQ:continuity} for the flux $\vec J = -\kappa\nabla u$ of a scalar function $u$. This is for instance known as Fourier's law of thermal conduction, where $u$ is the temperature and $\kappa$ is the dimensionless heat conductivity of the material. It is also Fick's law of diffusion where $u$ is a concentration and $\kappa$ is the diffusion coefficient.

Like in the previous section, we simplify the presentation by assuming that all hyperedges are flat and thus can be identified with a domain in $\ensuremath{\mathbb R}^\dimEdge$. In the more general case, the differential operators must be replaced by their differential geometric counterparts as in \cite{DziukE2013}.

We focus on the stationary case and set the time derivative in~\eqref{EQ:continuity} to zero. Thus, the discussion of the previous section leads to the following problem: find $u$ satisfying
\begin{subequations}\label{EQ:diffusion_primal}
\begin{align}
 - \div (\kappa \nabla u ) & = f && \text{ in all } \ensuremath{E}\in\ensuremath{{\mathcal E}}, \label{EQ:diffusion_eq}\\
 u & = u_\textup D && \text{ on all } \ensuremath{N} \in \ensuremath{\SetNode_\textup D}, \label{EQ:Dirichlet}\\
 u_{|\ensuremath{E}_1} & = u_{|\ensuremath{E}_2} && \text{ on all } \ensuremath{N} \in \ensuremath{{\mathcal N}},\; \ensuremath{N} \subset \partial \ensuremath{E}_1 \cap \partial \ensuremath{E}_2,\label{EQ:diffusion_cont}\\
 - \jump{\kappa \nabla u \cdot \mathbf n} & = g && \text{ on all } \ensuremath{N} \in \ensuremath{{\mathcal N}} \setminus \ensuremath{\SetNode_\textup D}, \label{EQ:interface}
\end{align}
\end{subequations}
for all $\ensuremath{E}_1, \ensuremath{E}_2 \in \ensuremath{{\mathcal E}}$, right hand sides $f$ and $g$, and a diffusion coefficient $\kappa \ge \kappa_0 > 0$. A justification by taking the limit of thin domains can be found in Section~\ref{SEC:limit} below.

We observe that in~\eqref{EQ:diffusion_primal} the diffusion equation \eqref{EQ:diffusion_eq} is complemented by three  boundary and interface conditions.
First, it is closed by a ``Dirichlet'' boundary condition~\eqref{EQ:Dirichlet}: We choose a non-empty set $\ensuremath{\SetNode_\textup D} \subset \ensuremath{{\mathcal N}}$ of ``Dirichlet'' hypernodes, on which we impose $u = u_\textup D$ for a prescribed boundary value $u_D$.
In~\eqref{EQ:diffusion_cont}, we employ a continuity constraint. This constraint prohibits jumps in the primary unknown across interior nodes, and therefore, loosely speaking, imitates the standard constraint that $u \in H^1$ of the domain.

On interior nodes $\ensuremath{N} \in \ensuremath{\SetNode_\textup I} \subset \ensuremath{{\mathcal N}} \setminus \ensuremath{\SetNode_\textup D}$, we set out with Kirchhoff's junction law, but with the option of a concentrated source $g$ in~\eqref{EQ:interface}. This equation also incorporates the Neumann condition $\vec -\kappa \nabla u \cdot \mathbf n = g$, since on a boundary hypernode the sum in the definition~\eqref{EQ:divergence_node} of the operator $\jump{\cdot}$ reduces to a single hyperedge. Note that \eqref{EQ:interface} for $g=0$ on interior nodes serves as compatibility condition for mimicking $-\kappa \nabla u \in H^\textup{div}$.
\begin{definition}[Function spaces on hypergraphs]\label{DEF:h}
 For each $\ensuremath{E}\in\ensuremath{{\mathcal E}}$ let $H^1(\ensuremath{E})$ be the standard Sobolev space on $\ensuremath{E}$ and $\gamma\colon H^1(\ensuremath{E})\to H^{1/2}(\partial\ensuremath{E})$ be the standard trace operator.
 
 Then, we define
 \begin{equation}
  \ensuremath{\mathcal H} = \left\{ u \in \bigoplus_{\ensuremath{E} \in \ensuremath{{\mathcal E}}} H^1(\ensuremath{E}) \;\middle|\; \begin{array}{@{\,}c@{\, }} \gamma_1 u = \gamma_2u \\ \text{on } \ensuremath{N} = \partial \ensuremath{E}_1 \cap \partial \ensuremath{E}_2, \; \ensuremath{N} \in \ensuremath{{\mathcal N}} \end{array} \right\},
 \end{equation}
 where $\gamma_1u$ and $\gamma_2u$ are the traces of $u$ from the hyperedges $\ensuremath{E}_1$ and $\ensuremath{E}_2$ on $\ensuremath{N}$, respectively.
 Due to the equality of traces in the definition of $\ensuremath{\mathcal H}$, we can define the trace operator to the skeleton
 \begin{equation}
  \ensuremath{\gamma} \colon \ensuremath{\mathcal H} \to \ensuremath{\mathcal M} := \left\{ \mu \in L^2(\Sigma) \;\middle|\; \begin{array}{@{\,}c@{\, }} \mu_{|\partial \ensuremath{E}} \in H^{1/2}(\partial \ensuremath{E}) \\ \text{for all } \ensuremath{E} \in \ensuremath{\mathcal E} \end{array} \right\}.
 \end{equation}
 
 Additionally, the spaces $\ensuremath{\mathcal H}_0$ and $\ensuremath{\mathcal M}_0$ are defined as
 \begin{align}
  \ensuremath{\mathcal M}_0 := & \{ \mu \in \ensuremath{\mathcal M} \colon \mu_{|\ensuremath{N}} = 0 \text{ for all } \ensuremath{N} \in \ensuremath{\SetNode_\textup D} \},\\
  \ensuremath{\mathcal H}_0 := & \{ u \in \ensuremath{\mathcal H} \colon \ensuremath{\gamma} u \in \ensuremath{\mathcal M}_0 \},
 \end{align}
 and we denote the  dual spaces of $\ensuremath{\mathcal H}_0$ by $\ensuremath{\mathcal H}_0^\star$ and of $\ensuremath{\mathcal M}_0$ by $\ensuremath{\mathcal M}_0^\star$.
 
 Norms ($\|\cdot\|_\ensuremath{\mathcal H}$ and $\|\cdot\|_\ensuremath{\mathcal M}$) on the respective spaces ($\ensuremath{\mathcal H}$ and $\ensuremath{\mathcal M}$) are induced by summed versions of the local scalar-products:
 \begin{align}
  ( u, v )_\ensuremath{\mathcal H} := & \sum_{\ensuremath{E} \in \ensuremath{{\mathcal E}}} (u_{|\ensuremath{E}}, v_{|\ensuremath{E}})_{H^1(\ensuremath{E})}, &
  \| u \|^2_\ensuremath{\mathcal H} := & (u, u)_\ensuremath{\mathcal H}, \\
 \langle \lambda, \mu \rangle_\ensuremath{\mathcal M} := & \sum_{\ensuremath{E} \in \ensuremath{{\mathcal E}}} \langle \lambda_{|\ensuremath{E}}, \mu_{|\ensuremath{E}}\rangle_{H^{1/2}(\partial \ensuremath{E})}, &
 \| \mu \|^2_\ensuremath{\mathcal M} := & \langle \mu, \mu \rangle_\ensuremath{\mathcal M},
 \end{align}
 such that
 \begin{gather}
     \| u \|^2_\ensuremath{\mathcal H} = \sum_{\ensuremath{E} \in \ensuremath{{\mathcal E}}} \| u_{|\ensuremath{E}} \|^2_{H^1(\ensuremath{E})}
     \quad\text{and}\quad
     \| \mu \|^2_\ensuremath{\mathcal M} = \sum_{\ensuremath{E} \in \ensuremath{{\mathcal E}}} \| \mu_{|\partial \ensuremath{E}} \|^2_{H^{1/2}(\partial \ensuremath{E})}.
 \end{gather}
\end{definition}
These definitions have a few immediate consequences:
\begin{enumerate}
 \item $\ensuremath{\gamma}\colon \ensuremath{\mathcal H} \to \ensuremath{\mathcal M}$ is a well-defined and surjective, linear, and continuous operator.
 \item We have the Gelfand triple relations
 \begin{gather}
  \ensuremath{\mathcal H}_0 \hookrightarrow L^2(\Omega) \cong [L^2(\Omega)]^\star \hookrightarrow \ensuremath{\mathcal H}_0^\star,\\
  \ensuremath{\mathcal M}_0 \hookrightarrow L^2(\Sigma) \cong [L^2(\Sigma)]^\star \hookrightarrow \ensuremath{\mathcal M}_0^\star.
 \end{gather}
\end{enumerate}
 Note that $\ensuremath{\mathcal M}^\star$ is analogous to space $M$ in of Raviart and Thomas \cite{RaviartT1977a}.
 
\begin{lemma}
 The space $\ensuremath{\mathcal H}$ with inner product $(\cdot, \cdot)_\ensuremath{\mathcal H}$ is a Hilbert space.
\end{lemma}
\begin{proof}
 Obviously, $\ensuremath{\mathcal H}$ is a subspace of the Hilbert space $\bigoplus H^1(\ensuremath{E})$, and the function
 \begin{equation*}
  h \colon \bigoplus_{\ensuremath{E} \in \ensuremath{{\mathcal E}}} H^1(\ensuremath{E}) \ni u \mapsto \sum_{\ensuremath{N} \in \ensuremath{\SetNode_\textup I}} \sum_{\bar \ensuremath{E}_1,\bar \ensuremath{E}_2 \supset \ensuremath{N}} \| u|_{\ensuremath{E}_1} - u|_{\ensuremath{E}_2} \|^2_{L^2(\ensuremath{N})} \in \ensuremath{\mathbb R}
 \end{equation*}
 is continuous and $\ensuremath{\mathcal H}$ is its kernel. Thus, $\ensuremath{\mathcal H}$ is closed.
\end{proof}
\begin{definition}
\label{DEF:weak_sol}
 A weak solution to the primal formulation of \eqref{EQ:diffusion_primal} with $\kappa \in L^\infty(\Omega)$, $f \in \ensuremath{\mathcal H}_0^\star$, and $g \in \ensuremath{\mathcal M}_0^\star$ is a function $u \in \ensuremath{\mathcal H}$ with $\gamma u = u_\textup D$ on all $\ensuremath{N} \in \ensuremath{\SetNode_\textup D}$, and
 \begin{equation}\label{EQ:weak_diff}
  \sum_{\ensuremath{E} \in \ensuremath{{\mathcal E}}} \int_\ensuremath{E} \kappa \nabla u \cdot \nabla v \ensuremath{\, \textup d x} = \langle f, v\rangle_{\ensuremath{\mathcal H}_0^\star, \ensuremath{\mathcal H}_0}  - \langle g, v \rangle_{\ensuremath{\mathcal M}_0^\star, \ensuremath{\mathcal M}_0} \qquad \forall v \in \ensuremath{\mathcal H}_0.
 \end{equation}
\end{definition}
In particular, if $f \in L^2(\Omega)$ and $g \in L^2(\Sigma)$, we can rewrite \eqref{EQ:weak_diff} as
\begin{equation}
  \sum_{\ensuremath{E} \in \ensuremath{{\mathcal E}}} \int_\ensuremath{E} \kappa \nabla u \cdot \nabla v \ensuremath{\, \textup d x} = \int_\Omega f v \ensuremath{\, \textup d x}  - \sum_{\ensuremath{N} \in \ensuremath{{\mathcal N}}} \int_\ensuremath{N} g v \ensuremath{\, \textup d \sigma} \qquad \forall v \in \ensuremath{\mathcal H}_0.
\end{equation}
\subsection{Existence and uniqueness of solutions}
\label{SEC:ExUnique}
\begin{theorem}\label{TH:ex_uni}
 Assume for $u_\textup D$ that there is a lifting $\bar u_\textup D \in \ensuremath{\mathcal H}$ with $\bar u_\textup D = u_\textup D$ on all $\ensuremath{N} \in \ensuremath{\SetNode_\textup D}$. If $\kappa \in L^{\infty}(\Omega)$ with $\kappa \ge \kappa_0 > 0$ a.~e., $f \in \ensuremath{\mathcal H}_0^\star$, $g \in \ensuremath{\mathcal M}_0^\star$, and all $\ensuremath{E} \in \ensuremath{{\mathcal E}}$ are Lipschitz domains, there is an unique weak solution $u$ according to Definition~\ref{DEF:weak_sol}, which continuously depends on the data.
\end{theorem}
\begin{proof}
Due to the existence of $\bar u_\textup D$, we can reduce the problem to the one with homogeneous Dirichlet values if we replace $u$ by $u-\bar u_\textup D$ and modifying the right hand side accordingly. Since the right hand side is bounded and $\ensuremath{\mathcal H}_0$ is a Hilbert space, it suffices to show ellipticity of the weak form to conclude the proof by the Lax--Milgram lemma. We note that for $v\in \ensuremath{\mathcal H}_0$ there holds
\begin{gather}
    \sum_{\ensuremath{E} \in \ensuremath{{\mathcal E}}} \int_\ensuremath{E} \kappa \nabla v \cdot \nabla v \ensuremath{\, \textup d x} 
    \ge \kappa_0 \sum_{\ensuremath{E} \in\ensuremath{{\mathcal E}}} \| \nabla v \|^2_{L^2(\ensuremath{E})}
\end{gather}
Thus, the following Poincaré--Friedrichs inequality implies ellipticity and concludes the proof.
\end{proof}

\begin{lemma}[Poincaré--Friedrichs inequality for $\ensuremath{\mathcal H}_0$]\label{LEM:Poincare1}
For all $v \in \ensuremath{\mathcal H}_0$ it holds that
\begin{align*}
 \Vert v \Vert_{L^2(\Omega)} \le C  \sum_{\ensuremath{E} \in \ensuremath{{\mathcal E}}} \| \nabla v \|_{L^2(\ensuremath{E})}.
\end{align*}
 \end{lemma}
 \begin{proof}
 Similar to the standard case of subdomains in $\ensuremath{\mathbb R}^{\ensuremath{\mathfrak d}}$,
 this inequality follows easily by contradiction:
 To this end, we assume that there is a sequence $(v_n)_{n=1,\ldots} \subset \ensuremath{\mathcal H}_0$ with
 \begin{equation}
  \lVert v_n \rVert_{L^2(\Omega)} =1
  \qquad\text{and}\qquad
  \sum_{\ensuremath{E} \in \ensuremath{{\mathcal E}}} \lVert \nabla v_n \rVert_{L^2(\ensuremath{E})} \le \frac1n.
 \end{equation}
 Thus, $v_{n|\ensuremath{E}}$ is bounded in $H^1(\ensuremath{E})$ for all $\ensuremath{E} \in \ensuremath{{\mathcal E}}$. Hence, by the weak compactness of the unit ball in $H^1(E)$ and the Rellich-Kondrachov theorem, there exists a subsequence (also denoted $v_n$) such that
 \begin{equation}
  v_{n|\ensuremath{E}} \rightarrow \tilde v_{|\ensuremath{E}} \text{ in } {L^2(\ensuremath{E})},
  \qquad \text{and} \qquad
  v_{n|\ensuremath{E}} \rightharpoonup \tilde v_{|\ensuremath{E}} \text{ in } H^1(\ensuremath{E}).
 \end{equation}
 We have that $\tilde v \in \ensuremath{\mathcal H}_0$ (due to its completeness), and that the seminorm $\sum_{\ensuremath{E} \in \ensuremath{{\mathcal E}}} \|\nabla \tilde v\|_{L^2(\ensuremath{E})} = 0$. Thus, $\tilde v$ is constant in all $\ensuremath{E} \in \ensuremath{{\mathcal E}}$ and overall continuous. Therefore it is overall constant and has to be zero, due to the zero boundary condition on Dirichlet nodes and the connectedness of $\Omega$. Hence, the strong convergence of $v_{n|E}$ in $L^2(E)$ implies $\lVert v_{n|E} \rVert_{L^2(E)} \rightarrow 0$, which contradicts $\| \tilde v \|_{L^2(\Omega)} = 1$. Therefore, the  Poincaré--Friedrichs inequality is valid.
\end{proof}
\section{HDG method for elliptic model equation}\label{SEC:hdg_graph}
When we derived PDE problems on hypergraphs, we were led to a formulation local on each hyperedge with coupling conditions on hypernodes. This is a structure which is nicely reflected in hybridized methods. Indeed, there the separation goes one step further. By putting degrees of freedom on the hypernode, values on hyperedges are not coupling anymore to other hyperedges across these hypernodes, but only to the values on the hypernodes constituting their boundary. Thus, differing from standard or discontinuous finite element methods, the number of hyperedges attached to a hypernode does not affect the solution process on a single hyperedge. Therefore, we consider hybridized methods ideally suited to PDEs on hypergraphs.

Hybridized discontinuous Galerkin (HDG) methods break the continuity condition~\eqref{EQ:diffusion_cont} by introducing Lagrange multipliers on each hypernode which enforce the continuity of fluxes~\eqref{EQ:interface} weakly.
It turns out though, that the Lagrange multiplier is an approximation to the solution $u$ of~\eqref{EQ:diffusion_primal} on the skeleton itself.

With such methods, the actual PDE \eqref{EQ:diffusion_eq} is represented locally on each hyperedge by Steklov-Poincaré operators on the hyperedges, which transform function values to flux values on the boundary of the hyperedges, a process called ``local solver'' in HDG terminology.
The global problem is posed in terms of the degrees of freedom on the hypernodes only, yielding a square, linear system of equations.

In this respect, HDG methods have a similar structure as the family of HHO methods. These are based on hybridizing the primal formulation and lead to a rather simple error analysis on polytopic meshes where only $L^2$ projections are used (as opposed to the rather complicated projections used for HDG). This is achieved by a novel stabilization design \cite{DiPietro2015}. For recent developments in hybrid high-order and HDG methods, the reader may consult \cite{Burman2018,Qiu2016}.

The separation of the local solution of bulk problems from the global coupling of interface variables is also achieved by the virtual element method \cite{VEM0,VEM1}. Thus, it fits into our view of coupled differential equations on connected hyperedges. Different to the methods discussed so far, it does not rely on polynomial shape functions inside mesh cells but rather on forms of fundamental solutions of any shape \cite{VEM2}. Accordingly, when applied to hypergraphs, the actual type of local solvers and of the specific boundary trace operators will differ from our approach, but remain within the same principal concept.
\subsection{The hybridized dual mixed formulation}\label{SEC:hdg_mixed}
In physical applications, there often is a need to receive reasonable approximations for both the primal unknown $u$ and the dual unknown $\vec q = - \kappa \nabla u$. In other words, considering diffusion, we would like to know both the distribution of some species' concentration and the species ``movement" (flux).  This becomes particularly important if we interpret $u$ as pressure and $\vec q$ as fluid flow through a porous medium (Darcy's equation). In this situation, the flow field $\vec q$ will govern the movement of chemical species dissolved within the fluid. It is often the main quantity of interest and conservativity is crucial. Therefore, we turn to the mixed formulation.

The mixed HDG methods use the weak, dual, mixed, hybrid formulation of \eqref{EQ:diffusion_primal}, i.e., find $(u,\vec q, \lambda) \in L^2(\Omega) \times \bigoplus \mathbf H^{\textup{div}}(\ensuremath{E}) \times \ensuremath{\mathcal M}$ with $\lambda = u_\textup D$ on all $\ensuremath{N} \in \ensuremath{\SetNode_\textup D}$ such that
\begin{subequations}\label{EQ:diffusion_mixed}
\begin{align}
 \int_\ensuremath{E} \left[ u (\div \vec p ) - \kappa^{-1} \vec q \cdot \vec p \right] \ensuremath{\, \textup d x} & = \int_{\partial \ensuremath{E}} \lambda \vec p \cdot \mathbf n \ensuremath{\, \textup d \sigma} && \forall \vec p \in \bigoplus_{\ensuremath{E}\in\ensuremath{{\mathcal E}}} \mathbf H^{\textup{div}}(\ensuremath{E}),
 \label{EQ:primal_dual}\\
  \int_\ensuremath{E} v\div \vec q \ensuremath{\, \textup d x} & = \int_\ensuremath{E} f v \ensuremath{\, \textup d x} && \forall v \in L^2(\Omega),
  \label{EQ:div_is_f}\\
  \sum_{\ensuremath{E} \in \ensuremath{{\mathcal E}}} \int_{\partial\ensuremath{E}} (\vec q \cdot \mathbf n) \mu \ensuremath{\, \textup d \sigma} &= \langle g, \mu \rangle_{\ensuremath{\mathcal M}^\star_0, \ensuremath{\mathcal M}_0} && \forall \mu \in \ensuremath{\mathcal M}_0.\label{EQ:interface_mixed}
\end{align}
\end{subequations}

Well-posedness of this formulation can be deduced from Theorem \ref{TH:ex_uni} if $f \in L^2(\Omega)$. Indeed, on the one hand, this implies that the (uniquely existing) solution $\overline{u}$ of Definition \ref{DEF:weak_sol} solves \eqref{EQ:diffusion_mixed} with $u = \overline{u}$, $\vec q = - \kappa \nabla \overline{u}$, and $\lambda = \gamma \overline{u}$. On the other hand, for any solution $(u,\vec q, \lambda) \in L^2(\Omega) \times \bigoplus \mathbf H^{\textup{div}}(\ensuremath{E}) \times \ensuremath{\mathcal M}$ of \eqref{EQ:diffusion_mixed}, we have $u \in \ensuremath{\mathcal H}$ (by the space's definition), and $\vec q = -\kappa \nabla u$ in the weak sense. Therefore, any solution to \eqref{EQ:diffusion_mixed} satisfies Definition \ref{DEF:weak_sol}.

Equations \eqref{EQ:primal_dual} \& \eqref{EQ:div_is_f} are local equations on the hyperedge, like in the standard case of a domain. They only couple to the Lagrange multipliers on the boundary of the hyperedge. Thus, we can eliminate them locally in the fashion of the Schur complement method.
To this effect, we introduce the local solution operator $S_\ensuremath{E}\colon \ensuremath{\mathcal M}\to\ensuremath{\mathcal M}^*$ for the right hand side $f=0$. It is in fact a Steklov-Poincaré operator on $\ensuremath{E}$ mapping the Dirichlet data $\lambda$ to the normal trace of the flux in~\eqref{EQ:interface_mixed}.
  
Then, the solution $\lambda$ of~\eqref{EQ:diffusion_mixed} can be characterized as
\begin{gather}
  \label{EQ:equation-steklov}
      \sum_\ensuremath{E} \langle S_\ensuremath{E}\lambda, \mu \rangle_{\ensuremath{\mathcal M}^\star_0, \ensuremath{\mathcal M}_0}
      =
      \langle g, \mu \rangle_{\ensuremath{\mathcal M}^\star_0, \ensuremath{\mathcal M}_0}.
\end{gather}

\begin{remark}
  The Steklov-Poincaré operators $S_\ensuremath{E}$ in this equation are the same ones as in the case of a manifold. They do not depend on the connectivity of a hypernode to other hyperedges. Thus, their implementation does not differ from that of a standard finite element method. The  only difference lies in the structure of the sum on the left, and is thus almost purely of algebraic nature.
\end{remark}

For inhomogeneous right hand side $f\neq 0$, we can define the operators
\begin{gather}
    \begin{split}
    (\tilde \ensuremath{\mathcal U}, \tilde \ensuremath{\vec{\mathcal Q}}) \colon L^2(\Omega) &\to L^2(\Omega) \times \bigoplus_{\ensuremath{E} \in \ensuremath{{\mathcal E}}} \mathbf H^{\textup{div}}(\ensuremath{E})\\
    f &\mapsto (u,\vec q)
    \end{split}
\end{gather}
by an edge-wise solution of \eqref{EQ:primal_dual} and \eqref{EQ:div_is_f} with $\lambda \equiv 0$.

In order to obtain a better understanding of the Steklov-Poincaré operators, we follow the route laid out in \cite{CockburnGL2009} for the discrete version and define the solution operators
\begin{gather}
    \begin{split}
    (\ensuremath{\mathcal U}, \ensuremath{\vec{\mathcal Q}}) \colon \ensuremath{\mathcal M} &\to L^2(\Omega) \times \bigoplus_{\ensuremath{E} \in \ensuremath{{\mathcal E}}} \mathbf H^{\textup{div}}(\ensuremath{E})\\
    \lambda &\mapsto (u,\vec q)
    \end{split}
\end{gather}
which map a given $\lambda$ to the element-wise solution of \eqref{EQ:primal_dual} and \eqref{EQ:div_is_f} with $f \equiv 0$.

The well-posedness and linearity of all local solution operators follow directly from the fact (see for instance \cite{BoffiBrezziFortin13}) that the mixed formulation on a single hyperedge is well-posed for any given $\lambda \in H^{1/2}(\partial \ensuremath{E})$ and its solution depends continuously on $\lambda$. 
Entering these solution operators into~\eqref{EQ:interface_mixed} yields
\begin{equation}
 - \sum_{\ensuremath{E} \in \ensuremath{{\mathcal E}}} \int_{\partial\ensuremath{E}} (\ensuremath{\vec{\mathcal Q}}\lambda \cdot \mathbf n) \mu \ensuremath{\, \textup d \sigma} = \sum_{\ensuremath{E} \in \ensuremath{{\mathcal E}}} \int_{\partial\ensuremath{E}} (\tilde \ensuremath{\vec{\mathcal Q}} f \cdot \mathbf n) \mu \ensuremath{\, \textup d \sigma} - \langle g, \mu \rangle_{\ensuremath{\mathcal M}^\star_0, \ensuremath{\mathcal M}_0},
\end{equation}
since $\vec q = \ensuremath{\vec{\mathcal Q}} \lambda + \tilde \ensuremath{\vec{\mathcal Q}} f$.
By some simple transformations of~\eqref{EQ:diffusion_mixed}, we can write~\eqref{EQ:equation-steklov} with inhomogeneous right hand sides in terms of bilinear and linear forms. This argument allows to reduce the problem to finding $\lambda \in \ensuremath{\mathcal M}$ with $\lambda = u_\textup D$ on all $\ensuremath{N} \in \ensuremath{\SetNode_\textup D}$, such that
\begin{subequations}
\begin{equation}
 a(\lambda, \mu) = b(\mu) \qquad \forall \mu \in \ensuremath{\mathcal M}_0,
\end{equation}
\begin{align}
 a(\lambda, \mu ) & = \sum_{\ensuremath{E} \in \ensuremath{{\mathcal E}}} \int_\ensuremath{E} \kappa^{-1} \ensuremath{\vec{\mathcal Q}} \lambda \cdot \ensuremath{\vec{\mathcal Q}} \mu \ensuremath{\, \textup d x},\label{EQ:cond_bil}\\
 b(\mu) & = \sum_{\ensuremath{E} \in \ensuremath{{\mathcal E}}} \int_{\partial\ensuremath{E}} (\tilde \ensuremath{\vec{\mathcal Q}} f \cdot \mathbf n) \mu \ensuremath{\, \textup d \sigma} - \langle g, \mu \rangle_{\ensuremath{\mathcal M}^\star_0, \ensuremath{\mathcal M}_0}.
\end{align}
\end{subequations}
Obviously, bilinear form $a$ and linear form $b$ are continuous due to the continuity of operators $\ensuremath{\vec{\mathcal Q}}$ and $\tilde \ensuremath{\vec{\mathcal Q}}$. Surprisingly, we have recovered a symmetric bilinear form. The following lemma is a key to the discrete well-posedness and adds the fact that this form is even $\ensuremath{\mathcal M}_0$-elliptic.
\begin{lemma}
 If $\kappa \ge \kappa_0 > 0$ and $\kappa \in L^\infty(\Omega)$, bilinear form $a$ from \eqref{EQ:cond_bil} is $\ensuremath{\mathcal M}_0$ elliptic.
\end{lemma}
\begin{proof}
Like in the proof of the Poincaré-Friedrichs inequality, we prove ellipticity of $a(\cdot,\cdot)$ by a contradiction argument. To this end, let $\lambda_n$ a sequence in $\ensuremath{\mathcal M}_0$ such that
\begin{gather}\label{EQ:contradiction-poincare2}
    \lVert\lambda_n\rVert_{\ensuremath{\mathcal M}} = 1
    \qquad\text{and}\qquad
    \lVert\ensuremath{\vec{\mathcal Q}}\lambda_n\rVert_{L^2(\Omega)} \to 0.
\end{gather}
Thus, there exists a subsequence $\lambda_n \rightharpoonup \tilde\lambda$ in $\ensuremath{\mathcal M}$, and by the compact embedding of $H^{1/2}(\partial\ensuremath{E})$ in $L^2(\partial\ensuremath{E})$ there holds again for a subsequence $\lambda_n \to \tilde \lambda$ in $L^2(\Sigma)$. Since $\ensuremath{\vec{\mathcal Q}}$ is continuous and $\lambda_n$ converges weakly in $H^{\frac12}(\partial E)$ we obtain $\ensuremath{\vec{\mathcal Q}} \lambda_n \rightharpoonup Q\tilde\lambda$ weakly in $\mathbf H^{\textup{div}}(\ensuremath{E})$ for each hyperdege and therefore also in $L^2(\Omega)$. Thus,~\eqref{EQ:contradiction-poincare2} implies that $\ensuremath{\vec{\mathcal Q}}\tilde\lambda = 0$. We denote by $(u,\vec q)$ the solution to~\eqref{EQ:primal_dual} and~\eqref{EQ:div_is_f} associated to $\tilde\lambda$, especially we have $\vec q = \ensuremath{\vec{\mathcal Q}} \tilde\lambda = 0$. Hence, for every 
$\ensuremath{E} \in \ensuremath{{\mathcal E}}$ we have
\begin{equation}
 \int_\ensuremath{E} u (\div \vec p) \ensuremath{\, \textup d x} = \int_{\partial \ensuremath{E}} \lambda \vec p \cdot \mathbf n \ensuremath{\, \textup d \sigma} \qquad \forall \vec p \in \mathbf H^{\textup{div}}(\ensuremath{E}).\label{EQ:const_u}
\end{equation}
Moreover, we have that the divergence operator
\begin{equation}
 \mathrm{div}\colon H^1_0(\ensuremath{E})^\ensuremath{\mathfrak d} \to L^2_0(\ensuremath{E}) := \{ u \in L^2(\ensuremath{E})\colon \int_\ensuremath{E} u = 0 \}
\end{equation} is surjective, and therefore
\begin{equation}
 \int_\ensuremath{E} u \varphi \ensuremath{\, \textup d x} = 0 \qquad \forall \varphi \in L^2_0(\ensuremath{E}).
\end{equation}
This, however, implies that $u$ is constant on $\ensuremath{E}$. Furthermore, since $u$ is constant on $\ensuremath{E}$, we can deduce by \eqref{EQ:const_u} and Gauss' divergence theorem that $\tilde\lambda = u$ is constant and
\begin{gather}
    \label{EQ:constant-poincare2}
    \tilde\lambda \equiv \operatorname{const} \quad\text{on}\;\partial \ensuremath{E}
    \qquad \forall \ensuremath{E} \in \ensuremath{{\mathcal E}}.
\end{gather}
The contradiction argument is concluded by the fact that some hyperedges are adjacent to the Dirichlet nodes and thus  $\tilde\lambda = 0$ on their boundary. For the other hyperedges, $\tilde\lambda = 0$ follows from connectedness of $\Omega$. Thus, $\tilde\lambda\equiv 0$ on $\Sigma$ in contradiction to $\lVert\lambda\rVert_{\ensuremath{\mathcal M}} = 1$. Altogether we showed that for all $\lambda \in \ensuremath{\mathcal M}_0$ it holds that
\begin{align*}
    \lVert \lambda \rVert_{\ensuremath{\mathcal M}} \lesssim \lVert \ensuremath{\vec{\mathcal Q}} \lambda \rVert_{L^2(\Omega)}.
\end{align*}
Together with \eqref{EQ:cond_bil} we obtain for a positive constant $\alpha$
\begin{align*}
    a(\lambda,\lambda) \geq \alpha \lVert \lambda \rVert_{\ensuremath{\mathcal M}}^2,
\end{align*}
i.e., the ellipticity of $a$ on $\ensuremath{\mathcal M}_0$.

Thus, the Lax-Milgram lemma guarantees a unique solution as soon as it is clear that $f$ and $g$ generate a right hand side in the dual of $\ensuremath{\mathcal M}_0$. But for $g$ this is obvious as $\ensuremath{\mathcal M}_0\subset L^2(\Sigma)$. For $f\in L^2(\Omega)$, there is a unique solution of~\eqref{EQ:primal_dual} and~\eqref{EQ:div_is_f} with $\lambda=0$. Its trace $\vec q\cdot\mathbf n$ is in $H^{-1/2}(\partial\ensuremath{E})$ for each hyperedge and thus bounded on $\ensuremath{\mathcal M}$.
\end{proof}
\subsection{HDG methods in dual mixed form}\label{SEC:HDG_form}
Let $\hat M$ be some finite dimensional, scalar function space. Then, we define the space of discrete functions on the skeleton $\Sigma$ by
\begin{gather}\label{EQ:skeletal_space}
 M := \left\{ \lambda \in L^2(\Sigma) \;\middle|\;
 \begin{array}{r@{\,}c@{\,}ll}
  \lambda_{|\ensuremath{N}} &\in& \hat M & \forall \ensuremath{N} \in \ensuremath{{\mathcal N}}\\
  \lambda_{|\ensuremath{N}} &=& 0 & \forall \ensuremath{N} \in \ensuremath{\SetNode_\textup D}    
 \end{array}
 \right\}.
\end{gather}
The mixed HDG methods involve a local solver on each hyperedge $\ensuremath{E} \in \ensuremath{{\mathcal E}}$, producing hyperedge-wise approximations $U_\ensuremath{E} \in V_\ensuremath{E}$ and and $\vec Q_\ensuremath{E} \in \vec W_\ensuremath{E}$ of the functions $u$ and $\vec q$ in equation \eqref{EQ:diffusion_mixed}, respectively. Here, $V_\ensuremath{E}$ is some finite dimensional, scalar function space, and $\vec W_\ensuremath{E}$ is some finite dimensional, vector valued function space. We will also use the concatenations of the spaces $V_\ensuremath{E}$ and $\vec W_\ensuremath{E}$, respectively, as a function space on $\Omega$, namely
\begin{gather}\label{EQ:dg_spaces}
 \begin{aligned}
  \ensuremath{V} &:=\bigl\{ v \in L^2(\Omega) & \big|\;v_{|\ensuremath{E}} &\in V_\ensuremath{E}, &\forall \ensuremath{E} &\in \ensuremath{{\mathcal E}} \bigr\},\\
  \vec W &:=\bigl\{ \vec q \in L^2(\Omega;\ensuremath{\mathbb R}^d) & \big|\;\vec q_{|\ensuremath{E}} &\in \vec W_\ensuremath{E}, &\forall \ensuremath{E} &\in \ensuremath{{\mathcal E}} \bigr\}.
 \end{aligned}
\end{gather}
The HDG scheme for \eqref{EQ:diffusion_mixed} on a hypergraph $\ensuremath{{\mathcal G}}$ consists of the local solver and a global coupling equation. The local solver is defined hyperedge-wise by a weak formulation of \eqref{EQ:diffusion_mixed} in the discrete spaces $V_\ensuremath{E} \times \vec W_\ensuremath{E}$ and defining suitable numerical traces and fluxes. Namely, given $\lambda \in M$ find $U_\ensuremath{E} \in V_\ensuremath{E}$ and $\vec Q_\ensuremath{E} \in \vec W_\ensuremath{E}$ , such that
\begin{subequations}\label{EQ:hdg_scheme}
 \begin{align}
  \int_\ensuremath{E} \frac{1}{\kappa} \vec Q_\ensuremath{E} \cdot \vec p \ensuremath{\, \textup d x} - \int_\ensuremath{E} U_\ensuremath{E} \ensuremath{\Nabla\!\cdot\!} \vec p \ensuremath{\, \textup d x} & = - \int_{\partial \ensuremath{E}} \lambda \vec p \cdot \mathbf n \ensuremath{\, \textup d \sigma} \label{EQ:hdg_primary}\\
  \int_{\partial \ensuremath{E}} ( \vec Q_\ensuremath{E} \cdot \mathbf n + \tau  U_\ensuremath{E} ) v \ensuremath{\, \textup d \sigma} - \int_\ensuremath{E} \vec Q_\ensuremath{E} \cdot \ensuremath{\nabla_\Edge} v \ensuremath{\, \textup d x}
  & = \tau \int_{\partial \ensuremath{E}} \lambda v \ensuremath{\, \textup d \sigma} \label{EQ:hdg_flux}
 \end{align}
\end{subequations}
hold for all $v \in V_\ensuremath{E}$, and all $\vec p \in \vec W_\ensuremath{E}$, and for all $\ensuremath{E} \in \ensuremath{{\mathcal E}}$. Here, $\tau \ge 0$ is the penalty coefficient. While the local solvers are implemented hyperedge by hyperedge, it is helpful for the analysis to combine them by concatenation. Thus, the local solvers define a mapping
\begin{gather}
 \begin{split}
  M & \to \ensuremath{V} \times \vec W\\
  \lambda &\mapsto (\ensuremath{\mathcal U} \lambda, \ensuremath{\vec{\mathcal Q}} \lambda),
 \end{split}
\end{gather}
where for each hyperedge $\ensuremath{E} \in \ensuremath{{\mathcal E}}$ holds $\ensuremath{\mathcal U} \lambda = U_\ensuremath{E}$ and $\ensuremath{\vec{\mathcal Q}} \lambda = \vec Q_\ensuremath{E}$. Analogously, we set  $\ensuremath{\mathcal U} (f, u_\textup D)$ and $ \ensuremath{\vec{\mathcal Q}} (f, u_\textup D)$, where now the local solutions are defined by the system
\begin{subequations}\label{EQ:hdg_f}
 \begin{align}
  \int_\ensuremath{E} \frac{1}{\kappa} \vec Q_\ensuremath{E} \cdot \vec p \ensuremath{\, \textup d x} - \int_\ensuremath{E} U_\ensuremath{E} \ensuremath{\Nabla\!\cdot\!} \vec p \ensuremath{\, \textup d x} & = - \int_{\partial \ensuremath{E}} u_\textup D \vec p \cdot \mathbf n \ensuremath{\, \textup d \sigma}\label{EQ:hdg_f_primary}\\
  \int_{\partial \ensuremath{E}} ( \vec Q_\ensuremath{E} \cdot \mathbf n + \tau U_\ensuremath{E} ) v \ensuremath{\, \textup d \sigma} - \int_\ensuremath{E} \vec Q_\ensuremath{E} \cdot \ensuremath{\nabla_\Edge} v & \ensuremath{\, \textup d x} \label{EQ:hdg_f_flux}\\
  = & \int_{\ensuremath{E}} f v \ensuremath{\, \textup d x} +  \tau \int_{\partial \ensuremath{E}} u_\textup D v \ensuremath{\, \textup d \sigma}\notag
 \end{align}
\end{subequations}

Once $\lambda$ has been computed, the HDG approximation to \eqref{EQ:diffusion_mixed} on $\ensuremath{{\mathcal G}}$ will be computed as
\begin{equation}\label{EQ:construct_unk}
 U = \ensuremath{\mathcal U} \lambda + \ensuremath{\mathcal U} (f, u_\textup D), \qquad \vec Q = \ensuremath{\vec{\mathcal Q}} \lambda + \ensuremath{\vec{\mathcal Q}} (f, u_\textup D)
\end{equation}

The global coupling condition is derived through a discontinuous
Galerkin version of mass balance and reads: Find
$\lambda \in M$, such that for all
$ \mu \in M$
\begin{equation}
 \sum_{\ensuremath{\eta} \in \ensuremath{\mathcal E}} \sum^{\ensuremath{\nu} \in \ensuremath{\mathcal N} \setminus \ensuremath{\setNode_\textup D}}_{\ensuremath{\nu} \in \ensuremath{\eta}} \int_\ensuremath{N} \left[ \vec Q \cdot \mathbf n + \tau (U - \lambda) \right] \mu \ensuremath{\, \textup d \sigma} = \sum_{\ensuremath{\nu} \in \ensuremath{\mathcal N} \setminus \ensuremath{\setNode_\textup D}} \int_\ensuremath{N} g \mu \ensuremath{\, \textup d \sigma}.\label{EQ:hdg_global}
\end{equation}
\begin{remark}\label{REM:spaces}
Hybridized DG methods in dual mixed form differ by the choice of local polynomial spaces and the stabilization parameter $\tau$. Defining $\ensuremath{\mathcal P}_p$ as the space of multivariate polynomials of degree at most $p$, Table~\ref{tab:hdg-methods} lists some well-known combinations on simplices.
\begin{table}[tp]
    \centering
    \begin{tabular}{l|ccc|c}
         Method &  $\hat M$ & $V_\ensuremath{E}$ & $\vec W_\ensuremath{E}$ & $\tau$ \\\hline
         LDG-H & $\ensuremath{\mathcal P}_p$ & $\ensuremath{\mathcal P}_p$ & $\ensuremath{\mathcal P}^d_p$ & $>0$\\
         RT-H & $\ensuremath{\mathcal P}_p$ & $\ensuremath{\mathcal P}_p$
            & $\ensuremath{\mathcal P}^d_p + \vec x \ensuremath{\mathcal P}_p$ & $=0$\\
         BDM-H & $\ensuremath{\mathcal P}_p$ & $\ensuremath{\mathcal P}_{p-1}$
            & $\ensuremath{\mathcal P}^d_p $ & $=0$\\
    \end{tabular}
    \caption{Combinations of local polynomial spaces and stabilization parameters for various hybridized methods.}
    \label{tab:hdg-methods}
\end{table}
Well-posedness of the local solvers for all of them is proven in \cite{CockburnGL2009} and the works cited there.
Analogous combinations based on tensor product polynomials exist for hypercubes.
\end{remark}
Existence and uniqueness of the discrete solution $\lambda$, $U$, and $\vec Q$ to the HDG method can be shown repeating the arguments mentioned in Section \ref{SEC:hdg_mixed} in the finite-dimensional setting. A natural assumption is the well-posedness of the local problems~\eqref{EQ:hdg_scheme}, see Remark~\ref{REM:spaces}.

Given the local solvers, the HDG method for elliptic diffusion problems is consistent with respect to the solution to~\eqref{EQ:diffusion_mixed}. Using consistency, we can immediately apply the analysis in \cite{CockburnGW2009}, as it proceeds locally for each hyperedge. Thus, we obtain optimal convergence rates for LDG-H (and also RT-H by slight adaptions) on simplicial hypergraphs. They also transfer to quadrilateral hypergraphs, since these allow for a Raviart--Thomas projection satisfying equation (2.7) in \cite{CockburnGW2009}.

\subsection{Numerical convergence tests for LDG-H}\label{SEC:diff_conv_ell}
\begin{table}[t]
 \begin{tabular}{cc|@{\,}lcc@{\,}lcc@{\,}lcc}
  \toprule
  \multicolumn{2}{c|@{\,}}{}      && \multicolumn{2}{c}{$\dimEdge=1$}  && \multicolumn{2}{c}{$\dimEdge=2$}   && \multicolumn{2}{c}{$\dimEdge=3$} \\
  \cmidrule{4-5} \cmidrule{7-8} \cmidrule{10-11}
  \multicolumn{2}{c|@{\,}}{mesh}  && err & eoc && err & eoc && err & eoc   \\
  \midrule
  \multirow{6}{*}{\rotatebox[origin=c]{90}{filling}}
  & $i = 0$ && 2.71e-1 & --- && 2.64e-1 & --- && 1.31e-1 & ---  \\
  & $i = 1$ && 9.82e-2 & 1.5 && 7.96e-2 & 1.7 && 3.24e-2 & 2.0  \\
  & $i = 2$ && 4.05e-2 & 1.3 && 2.55e-2 & 1.6 && 8.07e-3 & 2.0  \\
  & $i = 3$ && 1.81e-2 & 1.2 && 8.56e-3 & 1.6 && 2.01e-3 & 2.0  \\
  & $i = 4$ && 8.57e-3 & 1.1 && 2.94e-3 & 1.5 && 5.04e-4 & 2.0  \\
  & $i = 5$ && 4.16e-3 & 1.0 && 1.02e-3 & 1.5 && 1.26e-4 & 2.0  \\
  \midrule
  \multirow{4}{*}{\rotatebox[origin=c]{90}{refinement}}
  & $r = 0$ && 4.05e-2 & --- && 2.55e-2 & --- && 8.07e-3 & ---  \\
  & $r = 1$ && 1.00e-2 & 2.0 && 6.38e-3 & 2.0 && 2.01e-3 & 2.0  \\
  & $r = 2$ && 2.52e-3 & 2.0 && 1.59e-3 & 2.0 && 5.04e-4 & 2.0  \\
  & $r = 3$ && 6.30e-4 & 2.0 && 3.98e-4 & 2.0 && 1.26e-4 & 2.0  \\
  \bottomrule
 \end{tabular}\vspace{1ex}
 \caption{$L^2$ errors (err) and estimated orders of convergence (eoc) of linear approximation to the diffusion equation for hypergraphs with hyperedge dimension $\dimEdge$.}\label{TAB:diff_hg_el_conv}
\end{table}
\begin{figure}
 \begin{tabular}{ccc}
  \includegraphics[width=0.3\textwidth,draft=false]{pictures/i1r0.png} & \includegraphics[width=0.3\textwidth,draft=false]{pictures/i2r0.png} & \includegraphics[width=0.3\textwidth,draft=false]{pictures/i2r1.png} \\
  $i = 1$ & $i = 2$ & $i = 2$ \\
  $r = 0$ & $r = 0$ & $r = 1$ 
 \end{tabular}
 \caption{Pictures of the computational domains for $\ensuremath{\mathfrak d} = 2$, and different combinations of filling $i$ and refinement $r$. I all three illustrations, the same corner of the ``cubes" has been removed in the plots to illustrate the interior structure of the ``cubes".}\label{FIG:diff_hg_el_conv}
\end{figure}
Next, we consider a convergence example on a hypergraph. It is constructed to approximate
\begin{equation}\label{EQ:hypergraph_example}
 - \nabla \cdot (\kappa \nabla u) = f \text{ in } \ensuremath{E} \in \ensuremath{{\mathcal E}}, \qquad u = u_\text D \text{ on } \ensuremath{N} \in \ensuremath{\SetNode_\textup D},
\end{equation}
where the Dirichlet nodes are those that are located on the boundary of $[0,1]^\ensuremath{\mathfrak D}$ with $\ensuremath{\mathfrak D} = 3$.

The filling $i$ indicates that the cube has been $i$ times uniformly refined (in the standard three dimensional sense), and the calculation is conducted on the $\ensuremath{\mathfrak d}$ dimensional ``surfaces" of this filling. These surfaces themselves might be further refined $r$ times, and these refined surfaces are identified to be our standard hyperedges, see Figure \ref{FIG:diff_hg_el_conv} for an illustration.

The $\ensuremath{\mathfrak d}-1$ dimensional faces of this approach are interpreted as nodes and the nodes located on the boundary of the unit cube are considered Dirichlet nodes. All other nodes are supposed to be interior nodes. The solution is constructed to be $u = -x^2 - y^2 - z^2$, diffusion coefficient $d = 1$, and right-hand side $f = 2 \ensuremath{\mathfrak d}$. Of course, polynomial degrees $\ge 2$ are supposed to exactly reproduce the given solution, which is true in our numerical experiments. Thus, we only plot the errors for $p = 1$ in Table \ref{TAB:diff_hg_el_conv}.

Interestingly, the $L^2$ errors converge although with filling $i$, also the computational domain increases for $\ensuremath{\mathfrak d} = 1$ and $\ensuremath{\mathfrak d} = 2$. However, the rate of convergence deteriorates by $1$ if $\ensuremath{\mathfrak d} = 1$, and $\tfrac{1}{2}$ if $\ensuremath{\mathfrak d} = 2$. The optimal order is obtained for $\ensuremath{\mathfrak d} = 3$. 

Beyond this, the refinement indicated with $r$ uses filling level $i = 2$ and then uniformly refines the respective faces. This does not lead to an increase of the computational domain (even if $\ensuremath{\mathfrak d} < \ensuremath{\mathfrak D}$) and, therefore, gives the optimal convergence rate $p = 2$.

The aforementioned results have been obtained using our code HyperHDG \cite{HyperHDGgithub}.
\section{Hypergraph PDE as singular limit}
\label{SEC:limit}
\input{limit}
\section{Conclusions}\label{SEC:conclusions}
We have motivated the formulation of PDEs on geometric hypergraphs, which generalize the notions of ``domains'', ``graphs'', and ``network of surfaces''. Using a simple singular limit example in which several thin subdomains meet in a common inter-domain, we underlined that problems on hypergraphs might evolve from practical applications, and that the hybrid formulations of PDEs is particular useful for their expression. Thus, we used hybrid methods (in particular HDG methods), which intrinsically fit to this formulation, to approximate the solution of simple PDEs on hypergraphs. Doing so, we obtained the optimal convergence behavior, which is predicted by the theory of these methods applied to PDEs on standard domains.
\section*{Acknowledgements}
This work is supported by the German Research Foundation under Germany's Excellence Strategy EXC 2181/1 - 390900948 (the Heidelberg STRUCTURES Excellence Cluster).

The authors thank Dominic Kempf and the Scientific Software Center of Heidelberg University for their support in developing HyperHDG.
\bibliographystyle{ARalpha}

\subsection{Description of the 3D model problem}
We consider the simplified case of one hypernode $\ensuremath{N}$ with length $L$ connecting $m$ hyperedges $\ensuremath{E}_i$ for $i=1,\ldots,m$ which are rectangles with side lengths $L$ and $L_i$. Thus, Figure \ref{fig:hg_limit} shows the case with $m=3$. The opposite node of $\ensuremath{N}$ with respect to $\ensuremath{E}_i$ is denoted by $\ensuremath{N}_{i,e}$ (and is a boundary node).  Without loss of generality, we assume that $\ensuremath{N}$ lies in the $x_1$-axis and we have 
\begin{gather}
    \ensuremath{N}= \{s e_1 \, : \, s \in (0,L)\}.
\end{gather}
We denote by $\nu_i$ a unit normal vector to $E_i$ and define extruded hyperedges for $0 <\epsilon \ll 1$ and $0 < d_i$
\begin{equation*}
\tilde{E}_i^{\epsilon}:= \left\{ x \in \ensuremath{\mathbb R}^3 \,: \, x = y + s\nu_i \mbox{ for } s \in \left(-\dfrac{d_i \epsilon}{2},\dfrac{d_i \epsilon}{2}\right), \, y \in E_i \right\}.
\end{equation*}
Hence, $\tilde{E}_i^{\epsilon}$ is a hexahedron with side lengths $L$ and $L_i$, and with thickness $d_i \epsilon$. We construct now a domain $\Omega^\epsilon$ which contains the union of all these extruded hyperedges and a nonoverlapping decomposition of this domain.
To this end, let $\alpha>0$ be chosen such that the sets 
\begin{gather}
  \label{EQ:def-eei}
  E_i^{\epsilon}:= \left\{x \in \tilde{E}_i^{\epsilon} \, : \, \mathrm{dist}(x,\tilde{S}_i^{\epsilon}) > \alpha \epsilon \right\},
\end{gather}
do not overlap.
We denote the side of $\tilde{E}_i^{\epsilon}$ that contains $\ensuremath{N}$ by $\tilde{S}_i^{\epsilon}$, and define
\begin{align*}
S_i^{\epsilon}&:= \mathrm{int}\left\{x \in \partial E_i^{\epsilon} \, : \,  \mathrm{dist}(x,\tilde{S}_i^{\epsilon}) = \alpha \epsilon \right\}.
\end{align*}
The side lengths of $E_i^{\epsilon}$ are $L$ and $L_i - \alpha\epsilon$.
Additionally, we define the convex hull of the node $\ensuremath{N}$ and the sides $S_i^{\epsilon}$:
\begin{equation*}
N^{\epsilon}:= \mathrm{int}\left( \mathrm{conv}\{\overline{\ensuremath{N}},\overline{S_1^{\epsilon}},\ldots,\overline{S_m^{\epsilon}}\}\right).
\end{equation*}
By construction, we have
\begin{equation*}
N^{\epsilon} = (0,L)\times \epsilon \omega.
\end{equation*}
Then, we define the thin domain $\Omega^{\epsilon}$ as
\begin{equation*}
\Omega^{\epsilon} := N^{\epsilon} \cup \bigcup_{i=1}^m \left( E_i^{\epsilon} \cup S_i^{\epsilon}\right).
\end{equation*}
On $\Omega^{\epsilon}$ we define a  diffusion problem and then to pass to the limit $\epsilon \searrow 0$ in order to derive a problem on the hypergraph $\Omega = N \cup \bigcup_{i=1}^m E_i$.
To ensure uniqueness for our model we assume a zero-Dirichlet boundary condition on one face $S_D^{\epsilon}$.  We consider the following problem for the unknown $u^{\epsilon}: \Omega^{\epsilon} \rightarrow \ensuremath{\mathbb R}$
\begin{subequations}\label{EQ:diff_eps_model}
\begin{align}
- \nabla \cdot \big( \kappa^{\epsilon} \nabla u^{\epsilon} \big) &= f^{\epsilon} &\mbox{ in }& \Omega^{\epsilon},
\\
u^{\epsilon} &= 0 &\mbox{ on }& S_D^{\epsilon},
\\
-\kappa^{\epsilon} \nabla u^{\epsilon} \cdot \nu &= 0 &\text{ on }& \partial \Omega^{\epsilon} \setminus  S_D^{\epsilon}.
\end{align}
\end{subequations}
Here we assume that $\kappa^{\epsilon}$ is piecewise constant on $\ensuremath{E}^{\epsilon}_i$ and $\ensuremath{N}$, i.e. we have $\kappa^{\epsilon}\vert_{\ensuremath{E}^{\epsilon}_i} = \kappa_i  > 0$ and $\kappa^{\epsilon}\vert_{\ensuremath{N}} = \kappa_{\ensuremath{N}}  > 0$.  For $f^{\epsilon}$ we assume that $f^{\epsilon}\vert_{\ensuremath{E}^{\epsilon}_i} = f_i$ with $f_i \in C^0(\ensuremath{\mathbb R}^3)$ and 
  \begin{gather}
      f^\epsilon|_{\ensuremath{N}^\epsilon} (x)=
      \frac1{\epsilon\lvert\omega\rvert} g(x)
  \end{gather}
with $g \in C^0(\ensuremath{\mathbb R}^3)$.  
\begin{remark}\label{REM:scale_f}
 We emphasize that the function $f^{\epsilon}_N$ is of order $\epsilon^{-1}$. That is, a nonzero source term $g$ on a hypernode can only be caused by a large sink/source which converges to a measure on the hypernode.
\end{remark}
\begin{definition}
Let
\begin{gather}
    H^1_D(\Omega^{\epsilon}) =
    \bigl\{ u\in H^1(\Omega^{\epsilon}) \bigm\lvert u^{\epsilon}_{|S_D^{\epsilon}} = 0
    \bigr\}.
\end{gather}
We call $u^{\epsilon} \in H^1_D(\Omega^{\epsilon})$ a weak solution of $\eqref{EQ:diff_eps_model}$  if for all $\phi^{\epsilon} \in H^1_D(\Omega^{\epsilon})$ there holds
\begin{align}\label{EQ:Var_eps}
\int_{\Omega^{\epsilon}} \kappa^{\epsilon} \nabla u^{\epsilon} \cdot \nabla \phi^{\epsilon} \ensuremath{\, \textup d x}
= \int_{\Omega^{\epsilon}} f^{\epsilon} \phi^{\epsilon} \ensuremath{\, \textup d x}.
\end{align}
\end{definition}
As the standard solution theory for elliptic equations applies to $\Omega^\epsilon$, the Lax--Milgram Lemma immediately implies the existence of a unique weak solution of $\eqref{EQ:diff_eps_model}$ for all $\epsilon > 0$.
\subsection{Transformation to $\epsilon$ independent domains}
First, we map the domains $E_i^{\epsilon}$ and $N^{\epsilon}$ to fixed domains $E_{\mathrm{ref}}$ and $N_{\mathrm{ref}}$, respectively:
\begin{align*}
 E_{\mathrm{ref}} &:= (0,L) \times (0,L) \times \left(-\frac12,\frac12\right),\\
 N_{\mathrm{ref}} &:= (0,L) \times \omega.
\end{align*}
For $N^{\epsilon}$ we just use a simple scaling:
\begin{equation*}
\Phi_N^{\epsilon}: N^{\epsilon} \rightarrow N_{\mathrm{ref}},\qquad \Phi_N^{\epsilon}(x)= A_N^{\epsilon} x,
\end{equation*}
where the matrix $A_N^{\epsilon} \in \ensuremath{\mathbb R}^{3 \times 3} $ is given via
\begin{align*}
A_N^{\epsilon} = \begin{pmatrix}
1 & 0 & 0 \\
0  & \epsilon^{-1} & 0 \\
0 & 0 & \epsilon^{-1}
\end{pmatrix}.
\end{align*}
The transformation between $E_i^{\epsilon}$ and $E_{\mathrm{ref}}$ is defined by  ($i=1,\ldots,m$)
\begin{align}\label{TransformationReferenceElement}
\Phi_i^{\epsilon} : E_i^{\epsilon} \rightarrow E_{\mathrm{ref}}, \qquad \Phi_i^{\epsilon}(x) = A_i^{\epsilon} R_i x + a_i^{\epsilon},
\end{align}
where $R_i \in \ensuremath{\mathbb R}^{3\times 3}$ is a rotation around the $x_1$-axis, $A_i^{\epsilon} \in \ensuremath{\mathbb R}^{3\times 3}$ and $a_i^{\epsilon} \in \ensuremath{\mathbb R}^3$ are given by
\begin{align*}
R_i = \begin{pmatrix}
1 & 0 \\
0 & \bar{R}_i 
\end{pmatrix},
\quad
A_i^{\epsilon} = \begin{pmatrix}
1 & 0 & 0 \\
0 & \frac{L}{L_i - \alpha\epsilon} & 0 \\
0 & 0 & \frac{1}{d_i \epsilon} 
\end{pmatrix}, \quad a_i^{\epsilon} = \begin{pmatrix}
0 \\ - \frac{L}{L_i - \alpha\epsilon} \alpha\epsilon \\ 0
\end{pmatrix},
\end{align*}
where the matrix $\bar{R}_i \in \ensuremath{\mathbb R}^{2\times 2}$ is a rotation.
More precisely, $R_i$ can be determined by the equation $R_i\nu_i = e_3$, where $e_3$ is the unit vector in $x_3$-direction and $\nu_i$ is a unit normal on $E_i$.

Further, we write 
\begin{align*}
\Psi_N^{\epsilon}(x) &:= (\Phi_N^{\epsilon})^{-1}(x) = (A_N^{\epsilon})^{-1} x, 
\\
\Psi_i^{\epsilon} (x) &:= (\Phi_i^{\epsilon})^{-1}(x) = R_i^{-1} (A_i^{\epsilon})^{-1} \big(x - a_i^{\epsilon}\big).
\end{align*} 
The side of $E_{\mathrm{ref}}$ with the homogeneous Dirichlet boundary condition on $\ensuremath{E}_m$ is denoted by $S_D$.
It is characterized via
\begin{align*}
S_D := \Phi_i^{\epsilon}(S_{D}^{\epsilon})= (0,L)\times \{L\} \times \left(-\frac12,\frac12\right).
\end{align*} 
The sides of $N_{\mathrm{ref}}$ where $\ensuremath{N}^\epsilon$ interfaces to the hyperedges are
\begin{align*}
S_i = \Phi_N^{\epsilon}(S_i^{\epsilon}) \subset (0,L)\times \partial \omega.
\end{align*}


Now we make a change of variables to transform equation $\eqref{EQ:Var_eps}$ to the fixed domains $E_{\mathrm{ref}}$ and $N_{\mathrm{ref}}$. Then we define
\begin{align*}
u^{\epsilon}_i(x) &:= u^{\epsilon} \big(\Psi_i^{\epsilon}(x)\big) &&\mbox{ for almost every } x \in E_{\mathrm{ref}},
\\
u^{\epsilon}_N(x) &:= u^{\epsilon} \big(\Psi_N^{\epsilon} (x)\big) &&\mbox{ for almost every } x \in N_{\mathrm{ref}}.
\end{align*}
In other words, we identify the function $u^{\epsilon} \in H^1_D(\Omega^{\epsilon})$ with the tuple 
\begin{align*}
(u^{\epsilon}_1,\ldots,u^{\epsilon}_m, u^{\epsilon}_N) \in H^1(E_{\mathrm{ref}})^m \times H^1(N_{\mathrm{ref}}),
\end{align*}
together with the interface and boundary conditions
\begin{align}\label{EQ:inter_bound_eps}
u^{\epsilon}_N = u^{\epsilon}_i \circ \Phi_i^{\epsilon} \circ \Psi_N^{\epsilon} \,\, \mbox{ on } S_i, \qquad 
u^{\epsilon}_m = 0 \,\, \text{ on } S_D.
\end{align}

Next, we use the fact that
\begin{subequations}
\begin{equation}\label{EQ:def_s}
 (\Phi_i^{\epsilon} \circ \Psi_N^{\epsilon})|_{S_i}\colon S_i \to S := (0,L)\times \{0\} \times \left(-\frac12,\frac12\right)
\end{equation}
is an isomorphism between a face of $N_{\mathrm{ref}}$ and $E_{\mathrm{ref}}$. Thus, we can write it in the $\epsilon$-independent form
\begin{equation}\label{EQ:Composition}
 (\Phi_i^{\epsilon} \circ \Psi_N^{\epsilon})|_{S_i} = \underbrace{ \begin{pmatrix} 1  & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & d_i^{-1} \end{pmatrix} }_{ =: A_i } R_i x - de_2.
\end{equation}
\end{subequations}
Note that here, the first and second component of the mapping can be obtained by a simple concatenation of $\Phi^\epsilon_i$ and $\Psi^\epsilon_N$, in particular, the second component is zero ---cf.\ the definition of S in \eqref{EQ:def_s}, since $(R_i x)_2 = d$ for all $x \in S_i$.

By a change of coordinates and an elemental calculation we obtain
\begin{align}
\epsilon^2 \int_{N_{\mathrm{ref}}} & \kappa_{\ensuremath{N}} (A_N^{\epsilon})^2 \nabla u^{\epsilon}_N \cdot \nabla \phi_N \ensuremath{\, \textup d x} \notag\\
& + \sum_{i=1}^m \frac{(L_i - \alpha\epsilon) d_i \epsilon}{L}\int_{E_{\mathrm{ref}}} \kappa_i (A_i^{\epsilon})^2 \nabla u^{\epsilon}_i \cdot \nabla \phi_i \ensuremath{\, \textup d x} 
\label{EQ:var_eps_ref}\\
= &\frac{\epsilon}{\vert \omega\vert} \int_{N_{\mathrm{ref}}}   g \circ \Psi^{\epsilon}_N \phi_N \ensuremath{\, \textup d x} 
+ \sum_{i=1}^m \frac{(L_i - \alpha\epsilon) d_i \epsilon}{L} \int_{E_{\mathrm{ref}}} f_i \circ \Psi^{\epsilon}_i\phi_i \ensuremath{\, \textup d x}\notag
\end{align}
for all $\phi =(\phi_1,\ldots,\phi_m,\phi_N) \in H^1(E_{\mathrm{ref}})^m \times H^1(N_{\mathrm{ref}})$ which fulfills the boundary and interface conditions in $\eqref{EQ:inter_bound_eps}$.
\subsection{A priori estimates that are uniform in $\epsilon$}
In the following we use for $x \in \ensuremath{\mathbb R}^3$ the notation $x_\parallel = (x_1,x_2)$ and $x_{\perp} = (x_2,x_3)$, as well as $\nabla_{\parallel} = (\partial_1,\partial_2)^T$ and $\nabla_{\perp} = (\partial_2,\partial_3)^T$.

In a first step, we derive \textit{a priori} estimates for $u^{\epsilon}_N$ and $u^{\epsilon}_i$ uniformly with respect to $\epsilon$. We define the space
\begin{align*}
H^1(N_{\mathrm{ref}},\nabla_{\perp}) := \left\{\phi_N \in L^2(N_{\mathrm{ref}}) \, : \, \nabla_{\perp} \phi_N \in L^2(N_{\mathrm{ref}})^2 \right\},
\end{align*}
with the norm $\Vert \phi_N \Vert^2_{H^1(N_{\mathrm{ref}},\nabla_{\perp})} := \Vert \phi_N\Vert_{L^2(N_{\mathrm{ref}})}^2 + \Vert \nabla_{\perp} \phi_N\Vert_{L^2(N_{\mathrm{ref}})}^2$. The proof of the trace theorem in \cite[Chapter 5.5, Theorem 1]{EvansPartialDifferentialEquations} implies the existence of a bounded linear trace operator
\begin{align*}
    \gamma_{\nabla{\perp}} : H^1(N_{\mathrm{ref}},\nabla_{\perp}) \rightarrow L^2(S_N)
\end{align*}
for $S_N:= (0,L)\times \partial \omega$.
We use the abbreviations $\phi_{N\vert S_N} := \gamma_{\nabla_{\perp}}(\phi_N)$ and $\phi_{N\vert S_i}$ for the restriction to $S_i$, respectively. Additionally, we have the following Poincaré-inequality:
\begin{lemma}\label{LemmaPoincareInequality}
For all $v = (v_1,\ldots,v_m,v_N) \subset H^1(E_{\mathrm{ref}})^m \times H^1(N_{\mathrm{ref}})$ with the boundary conditions $\eqref{EQ:inter_bound_eps}$ it holds 
\begin{multline*}
\sum_{i=1}^m\Vert v_i \Vert_{L^2(E_{\mathrm{ref}})}+\Vert v_N \Vert_{L^2(N_{\mathrm{ref}})} \\
\le C \left(\Vert \nabla_{\perp} v_N \Vert_{L^2(N_{\mathrm{ref}})} + \sum_{i=1}^m \Vert \nabla v_i \Vert_{L^2(E_{\mathrm{ref}})} \right).
\end{multline*}
\end{lemma}
\begin{proof}
As for the  Poincaré-inequality in the proof of Theorem \ref{TH:ex_uni} we use a contradiction argument. We assume that there exists a sequence  $(v_1^n,\ldots,v_m^n,v_N^n)_n \subset H^1(E_{\mathrm{ref}})^m \times H^1(N_{\mathrm{ref}})$ with the boundary conditions $\eqref{EQ:inter_bound_eps}$, such that
\begin{align}
\begin{aligned}\label{EQ:Poincare_Contradiction}
1 = \sum_{i=1}^m\Vert v_i^n \Vert_{L^2(E_{\mathrm{ref}})}+\Vert v_N^n & \Vert_{L^2(N_{\mathrm{ref}})}  
\\
&\geq n \left(\Vert \nabla_{\perp} v_N^n \Vert_{L^2(N_{\mathrm{ref}})} + \sum_{i=1}^m \Vert \nabla v_i^n \Vert_{L^2(E_{\mathrm{ref}})} \right).
\end{aligned}
\end{align}
Since $v_i^n$ is bounded in $H^1(E_{\mathrm{ref}})$ and $v_N^n$ is bounded in $H^1(N_{\mathrm{ref}},\nabla_{\perp})$, there exist $v_i$ and $v_N$, such that up to a subsequence 
\begin{align*}
v_i^n &\rightharpoonup v_i &\mbox{ weakly in }& H^1(E_{\mathrm{ref}}),
\\
v_i^n &\rightarrow v_i &\mbox{ in }& L^2(E_{\mathrm{ref}}),
\\
v_N^n &\rightharpoonup v_N &\mbox{ weakly in }& H^1(N_{\mathrm{ref}},\nabla_{\perp}).
\end{align*}
Further, due to \eqref{EQ:Poincare_Contradiction} we have $\nabla v_i = 0 $ and $\nabla_{\perp} v_N = 0$. This implies that $v_i$ is constant on $E_{\mathrm{ref}}$ (and $v_m = 0$ due to the zero boundary condition on $S_D$) and there exists $\bar{v}_N \in L^2(N)$ such that $v_N(x) = \bar{v}_N(x_1)$ for almost every $x \in N_{\mathrm{ref}}$. The continuity of the usual trace operator on $H^1(E_{\mathrm{ref}})$ and the continuity of $\gamma_{\nabla_{\perp}}$ imply that $v_i(A_iR_ix - de_2)\vert_{S_i} = \bar{v}(x_1)$ on $S_i$.
Further, the continuity of the trace operator implies the weak convergence of the traces (from both sides). Additionally we used, that we have $v_N^n\vert_{S_i} = v_i^n\vert_{S_i}$, which can be shown by a density argument. This implies $v_i = 0$ and $v_N = 0$, which contradicts  $1 = \sum_{i=1}^m\Vert v_i \Vert_{L^2(E_{\mathrm{ref}})}+\Vert v_N  \Vert_{L^2(N_{\mathrm{ref}})}  $.
\end{proof}
Now, we obtain the following \textit{a priori} estimates:
\begin{lemma}\label{LE:apriori_estimates}
For $u^{\epsilon}_i$ and $u^{\epsilon}_N$ it holds that
\begin{align*}
\epsilon & \Vert \partial_1 u^{\epsilon}_N \Vert_{L^2(N_{\mathrm{ref}})} + \Vert \nabla_{\perp} u^{\epsilon}_N \Vert_{L^2(N_{\mathrm{ref}})} 
\\
&+ \sum_{i=1}^m \left\{ \sqrt{\epsilon} \Vert \nabla_{\parallel} u^{\epsilon}_i \Vert_{L^2(E_{\mathrm{ref}})} + \Vert \partial_3 u^{\epsilon}_i \Vert_{L^2(E_{\mathrm{ref}})} \right\} \le C\sqrt{\epsilon}
\end{align*}
for a constant $C > 0$ independent of $\epsilon$.
\end{lemma}
\begin{proof}
Choosing  $\phi_i = u^{\epsilon}_i$ and $\phi_N = u^{\epsilon}_N$ as a test-function in $\eqref{EQ:var_eps_ref}$ and using the positivity of $\kappa^{\epsilon}$ and the assumptions on $f^{\epsilon}$, as well as the trace-inequality and the Poincar\'e-inequality from Lemma \ref{LemmaPoincareInequality}, we obtain
\begin{align*}
\epsilon^2 &\Vert \partial_1 u_N^{\epsilon}\Vert_{L^2(N_{\mathrm{ref}})}^2 + \Vert \nabla_{\perp} u^{\epsilon}_N \Vert_{L^2(N_{\mathrm{ref}})}^2
\\
&+ \sum_{i=1}^m \left\{ \epsilon \Vert \nabla_{\parallel} u^{\epsilon}_i \Vert_{L^2(E_{\mathrm{ref}})}^2 + \Vert \partial_3 u^{\epsilon}_i \Vert_{L^2(E_{\mathrm{ref}})} \right\} 
\\
\le & C \epsilon \Vert g \circ \Psi_{\ensuremath{N}}^{\epsilon} \Vert_{L^2(N_{\mathrm{ref}})} \Vert u^{\epsilon}_N \Vert_{L^2(N_{\mathrm{ref}})} + C\epsilon \sum_{i=1}^m \Vert f_i \circ \Psi_i^{\epsilon} \Vert_{L^2(E_{\mathrm{ref}})} \Vert u^{\epsilon}_i \Vert_{L^2(E_{\mathrm{ref}})} 
\\
\le & C \epsilon \Vert u^{\epsilon}_N \Vert_{L^2(N_{\mathrm{ref}})} + C\epsilon \sum_{i=1}^m \Vert u^{\epsilon}_i \Vert_{L^2(E_{\mathrm{ref}})} 
\\
\le & C\epsilon + \frac12 \left(  \Vert \nabla_{\perp} u^{\epsilon}_N \Vert_{L^2(N_{\mathrm{ref}})}^2 + \epsilon \sum_{i=1}^m \Vert \nabla  u^{\epsilon}_i \Vert_{L^2(E_{\mathrm{ref}})}^2 \right)
\end{align*}
The second term on the right-hand side can be absorbed from the left-hand side and we obtain the desired result.
\end{proof}
\subsection{Convergence by compactness results and characterization of limit problem}
Using the weak compactness of the unit ball in $L^2$ and the Rellich-Kondrachov theorem we immediately obtain the following compactness result:

\begin{corollary}\label{CO:convergences_eps}
There exists $u_N^0 \in L^2(N_{\mathrm{ref}})$ with $\nabla_{\perp} u_N^0 = 0$ in the weak sense, and $u_i^0 \in H^1(E_{\mathrm{ref}})$ for $i=1,\ldots,m$, such that up to a subsequence
\begin{align*}
u^{\epsilon}_N &\rightharpoonup u_N^0 &\mbox{ weakly in }& L^2(N_{\mathrm{ref}}),
\\
u^{\epsilon}_i &\rightarrow u_i^0 &\mbox{ strongly in }& L^2(E_{\mathrm{ref}}),
\\
\nabla u^{\epsilon}_i &\rightharpoonup \nabla u_i^0 &\mbox{ weakly in }& L^2(E_{\mathrm{ref}}),
\\
\nabla_{\perp} u^{\epsilon}_N &\rightarrow 0 &\mbox{ strongly in }& L^2(N_{\mathrm{ref}}),
\\
\partial_3 u^{\epsilon}_i &\rightarrow 0 &\mbox{ strongly in }& L^2(E_{\mathrm{ref}}).
\end{align*}
Especially the last two convergences imply $\nabla_{\perp} \tilde{u}_N^0 = 0$ and $\partial_3 \tilde{u}_i^0 = 0$ in the weak sense. Hence, there exist $u_N^0 \in L^2(N)$ and $u_i^0 \in H^1((0,L)^2)$ such that $\tilde{u}_N^0(x) = u_N^0(x_1)$  for almost every $x \in N_{\mathrm{ref}}$ and $\tilde{u}_i^0(x) = u_i^0(x_\parallel)$ for almost every $x \in E_{\mathrm{ref}}$.
\end{corollary}
In the following, we drop the notation $\tilde{\cdot}$ and just use the notation $u_N^0 \in L^2(N_{\mathrm{ref}})$ and $u_i^0\in H^1(E_{\mathrm{ref}})$ for the limit functions. Let us consider the interface and boundary conditions for these limits. Obviously, the zero boundary condition $u^{\epsilon}_m = 0 $ on $S_D$ is inherited to  $u_0^m$ by the continuity of the trace operator on $H^1(E_{\mathrm{ref}})$. 

The weak convergence of $u^{\epsilon}_N$ in $H^1(N_{\mathrm{ref}},\nabla_{\perp})$ to $u_N^0$ implies the weak convergence $u^{\epsilon}_N\vert{S_i} \rightharpoonup u_N^0 $ in $L^2(S_i)$. Further, the weak convergence of $u^{\epsilon}_i $ in $H^1(E_{\mathrm{ref}})$ and the compactness of the embedding $H^1(E_{\mathrm{ref}}) \hookrightarrow L^2(\partial E_{\mathrm{ref}})$ implies the strong convergence   $u^{\epsilon}_i\vert_{S_i} \rightarrow u_i^0 $ in $L^2(S_i)$. Using $\eqref{EQ:inter_bound_eps}$ and $\eqref{EQ:Composition}$ we obtain for all $\phi \in C_0^{\infty}(S_i)$
\begin{align*}
\int_{S_i} u_N^0 \phi d\sigma &= \lim_{\epsilon\to 0} \int_{S_i} u^{\epsilon}_N \phi d\sigma 
= \lim_{\epsilon\to 0} \int_{S_i} u^{\epsilon}_i (A_iR_i x - d e_2) \phi d\sigma
\\
&= \int_{S_i} u_i^0(A_iR_i x - de_2) \phi d\sigma.
\end{align*}
This implies the interface condition
\begin{align}\label{EQ:inter_limit}
u_0^N\vert_{S_i}(x) = u_i^0\vert_{S}  (A_i R_i x - de_2) \quad \mbox{ for almost every } x \in S_i.
\end{align}
Now, let us pass to the limit in the variational equation $\eqref{EQ:var_eps_ref}$ for suitable test-functions. We choose $\phi_N \in C^{\infty}(\overline{N})$ (hence $\phi_N $ is constant on every section $\{x_1\}\times \omega$ of $N_{\mathrm{ref}}$) and $\phi_i \in C^{\infty}(\overline{(0,L)^2})$ with $\phi_i = \phi_N$ on $ (0,L)\times \{0\} $ ($\phi_i$ is constant in $x_3$-direction) and $\phi_m $ has compact support away from $S_D$. Obviously, $(\phi_1,\ldots,\phi_m,\phi_N)$ is an admissible test-function for $\eqref{EQ:var_eps_ref}$. We multiply $\eqref{EQ:var_eps_ref}$ with $1/\epsilon$ and obtain
\begin{align}
\begin{aligned}
\label{EQ:aux_var_eps_limit}
\epsilon &\int_{N_{\mathrm{ref}}} \kappa_N \partial_1 u^{\epsilon}_N \partial_1 \phi_N \ensuremath{\, \textup d x} 
+ \sum_{i=1}^m \frac{(L_i - \alpha\epsilon) d_i }{L}\int_{E_{\mathrm{ref}}} \kappa_i (A_i^{\epsilon})^2 \nabla u^{\epsilon}_i \cdot \nabla \phi_i \ensuremath{\, \textup d x} 
\\
=&  \int_{\ensuremath{N}}  g \circ \Psi^{\epsilon}_N \phi_N \ensuremath{\, \textup d x}_1
+ \sum_{i=1}^m \frac{(L_i - \alpha\epsilon) d_i }{L} \int_{E_{\mathrm{ref}}} f_i \circ \Psi^{\epsilon}_i \phi_i \ensuremath{\, \textup d x}.
\end{aligned}
\end{align}
Now, we pass to the limit $\epsilon \to 0$ in every single term. From Lemma \ref{LE:apriori_estimates}  we immediately obtain that the first term on the left-hand side is of order $\sqrt{\epsilon}$ and vanishes for $\epsilon \to 0$. 
The  convergence of $\nabla u^{\epsilon}_i$ from Corollary \ref{CO:convergences_eps} implies
\begin{align*}
 \frac{(L_i - \alpha\epsilon) d_i }{L}&\int_{E_{\mathrm{ref}}} \kappa_i (A_i^{\epsilon})^2 \nabla u^{\epsilon}_i \cdot \nabla \phi_i \ensuremath{\, \textup d x} 
\\
\overset{\epsilon \to 0}{\longrightarrow}& \frac{L_i d_i}{L} \int_{(0,L)^2} \kappa_i \left( \partial_1 u_i^0 \partial_1 \phi_i + \frac{L^2}{L_i^2}\partial_2 u_i^0 \partial_2 \phi_i \right) \ensuremath{\, \textup d x}_\parallel.
\end{align*}
Further, we have 
\begin{align*}
  \frac{(L_i - \alpha\epsilon) d_i }{L} \int_{E_{\mathrm{ref}}} f_i \circ \Psi^{\epsilon}_i \phi_i \ensuremath{\, \textup d x} \overset{\epsilon \to 0 }{\longrightarrow} \frac{L_i}{L}d_i \int_{(0,L)^2} f_i \circ \Psi_i^0 \phi_i \ensuremath{\, \textup d x}_\parallel,  
\end{align*}
with 
\begin{align*}
    \Psi_i^0(x):= R_i^{-1} \begin{pmatrix}
1 & 0 & 0 \\
0 & \frac{L}{L_i} & 0 \\
0 & 0 & 0
\end{pmatrix}x,
\end{align*}
and in a similar way we get
\begin{align*}
    \int_{\ensuremath{N}}  g \circ \Psi^{\epsilon}_N \phi_N \ensuremath{\, \textup d x}_1 \overset{\epsilon \to 0 }{\longrightarrow} \int_{\ensuremath{N}} g(x_1,0,0) \phi_N \ensuremath{\, \textup d x}_1.
\end{align*}
Now, we define the function $u^0$ on the hypergraph $\Omega$ in the following way: For almost every $x \in E_i $ we define (we emphasize that $(R_ix)_3 = 0$)
\begin{align*}
u^0(x) := u_i^0 (C_iR_i x) \quad \mbox{ with } C_i = \begin{pmatrix}
1 & 0 & 0 \\
0 & \frac{L}{L_i} & 0 \\
0 & 0 & 0
\end{pmatrix}.
\end{align*}
Hence, we have $u^0 \in \ensuremath{\mathcal H}$ with $u^0\vert E_m = 0$ on $\ensuremath{N}_{m,e}$ (see Definition \ref{DEF:weak_sol}), and altogether, we obtain for $\epsilon \to 0$ in $\eqref{EQ:aux_var_eps_limit}$ after a change of coordinates for all $\phi \in \ensuremath{\mathcal H}$ with $\phi_m = 0$ on $\ensuremath{N}_{m,e}$
\begin{multline}
\sum_{i=1}^m d_i \int_{E_i} \kappa_i \nabla_{E_i} u^0_i \cdot \nabla_{E_i} \phi \ensuremath{\, \textup d x}
= \\
 \int_{\ensuremath{N}} g(x_1,0,0)\phi_{\vert\ensuremath{N}} \ensuremath{\, \textup d x}_1
+  \sum_{i=1}^n d_i \int_{\ensuremath{E}_i}f_i \phi_i \ensuremath{\, \textup d x}.
\end{multline}
In other words, $u^0$ is the unique weak solution of ($\tilde\kappa_i=\tilde\kappa_{\vert E_i}= d_i \kappa_i$)
\begin{subequations}
\begin{align}
- \nabla_{E_i} \cdot (\tilde\kappa_i \nabla_{E_i} u^0) &= d_if_i &\mbox{ for }& i=1,\ldots,m,
\\
u^0 &= 0 &\mbox{ on }& \ensuremath{N}_{m,e},
\\
-  \jump{\tilde\kappa  \nabla_{E} u^0 \cdot \mathbf{n}} &=g &\mbox{ on }& N,
\\
-\tilde\kappa_i \nabla_{E_i} u^0 \cdot \mathbf{n} &= 0 &\mbox{ on }& \ensuremath{N}_{i,e} \mbox{ for } i=1,\ldots,m-1,
\\
u^0\vert_{E_i} &= u^0\vert_{E_j} &\mbox{ on }& \ensuremath{N} \mbox{ for } i,j=1,\ldots,m.
\end{align}
\end{subequations}
\subsection{Concluding remarks}
For the ease of presentation, we have just performed the limit analysis for a very simple model case. Therefore, we have to discuss the validity of our analysis in more general cases.
\begin{enumerate}[leftmargin=*]
    \item Different combinations of boundary conditions, Dirichlet, (inhomogeneous) Neumann, or Robin will yield the same result as long as the variation of boundary functions in direction normal to the hyperedges $\ensuremath{E}_i$ vanishes for $\epsilon \searrow 0$.
    \item A hyperedge with several ``interior'' hypernodes: after changing the definition of $\ensuremath{E}_i^\epsilon$ in~\eqref{EQ:def-eei} and the reference mapping $\Phi_i^\epsilon$ in~\eqref{TransformationReferenceElement}, two hypernode domains $\ensuremath{N}_i^\epsilon$ meeting in the same corner will have an intersection, which we have to treat separately in order to have a nonoverlapping decomposition of $\Omega^\epsilon$. But then, the integrals over these corner domains converge to zero by one order faster than the first term in~\eqref{EQ:aux_var_eps_limit}. Thus, they do not enter the limit equation. Obviously, not all nodes can be located in the $x$-axis, but an affine transformation to the reference node $N_{\mathrm{ref}}$ can always be found and the analysis works in the same way as shown.
    \item General hypergraphs with planar hyperedges: after the previous point, it is clear that we can construct and decompose $\Omega^\epsilon$ in the same fashion for any finite hypergraph, provided $\epsilon$ sufficiently small.
    \item Hypergraphs with smooth nonplanar hyperedges: in this case, the reference mappings become nonlinear mappings and many aspects become technically much more involved. The cross section $\epsilon\omega$ of the node $\ensuremath{N}^\epsilon$ may depend smoothly on the tangential coordinate, but it will always have positive diameter and will be bounded as long as the hyperedges are smooth manifolds. Again, there will be an upper bound for $\epsilon$, but the limit properties will not be affected.
    \item Higher dimensional hypergraphs and graphs: here we end up in a situation with new reference domains
    \begin{align*}
 E_{\mathrm{ref}} &:= (0,L)^{\ensuremath{\mathfrak d}}  \times \left(-\frac12,\frac12\right)^{\ensuremath{\mathfrak D}-\ensuremath{\mathfrak d}},\\
 N_{\mathrm{ref}} &:= (0,L)^{\ensuremath{\mathfrak d}-1} \times \omega
\end{align*}
with $|\epsilon \omega| \sim \epsilon^{\ensuremath{\mathfrak D}+1-\ensuremath{\mathfrak d}}$.
\item More general assumptions on $\kappa^{\epsilon}$ are possible. For example   continuity of $\kappa^{\epsilon}$ on $\ensuremath{E}_i^{\epsilon}$ and $\ensuremath{N}$ as in the definition of $f^{\epsilon}$  is enough (jumps between the different compartments are valid). Another possible choice is to construct $\kappa^{\epsilon}$ from a $L^{\infty}$-function on $\ensuremath{E}_i$ and $\ensuremath{N}$, constantly extended in normal direction with respect to $\ensuremath{E}_i$ resp. $\ensuremath{N}$.
\end{enumerate}

We conclude this section with some observations about the singular limit model and the 3D model problem:

First, we observe that in the limit problem the angles under which two or more hyperedges meet have become irrelevant, while these angles certainly have been relevant for the 3D model problem. Note that the definition of $\omega$ needs that there are no angles of zero degrees and if there are small angles, the shape of $\omega$ will compensate this drawback. Since we take $\epsilon \searrow 0$, this shape as well as the volume of $\omega$ is does not affect the limit problem. This effect implies that the $\epsilon$-limit of solutions will be approached slower when angles become small.

Second, nodal sources $g$ correspond to ``strong" sources $f^\epsilon_\ensuremath{N}$: Since the measures (and thereby the effects) of hypernodes are scaled by $\epsilon^2$ and vanish more quickly than those of hyperedges, sources on hypernodes are only relevant if they are of order of $\epsilon^{-1}$, cf.\ Remark \ref{REM:scale_f}. That is, analytically sources even become stronger during the limiting process.
