# Volume enclosed by a flux surface

R.S.MacKay  Mathematics Institute, University of Warwick, Coventry CV4 7AL, UK
<EMAIL_ADDRESS>
###### Abstract.

The paper describes ways that the computation of the volume enclosed by an
invariant torus (flux surface) for a magnetic field can be reduced from a 3D
integral to a 2D integral.

##  1\. Introduction

In generality, a magnetic field is a vector field $B$ on an orientable 3D
manifold $M$ that preserves a volume-form $\Omega$ (for the standard Euclidean
volume, one writes $\textrm{div}\,B=0$). Defining the associated flux 2-form
$\beta=i_{B}\Omega$, the volume-preservation condition can be written as
$d\beta=0$ ($\beta$ is closed). For a tutorial on use of differential forms in
plasma physics, see [M1].

A stronger condition (for general $M$) is that $\beta$ be exact, i.e.
$\beta=d\alpha$ for some 1-form $\alpha$. It is usual to write this in terms
of the associated vector potential $A$, related by $\alpha=A^{\flat}$ (where
for a Riemannian metric $g$, $A^{\flat}$ is defined by
$A^{\flat}(\xi)=g(A,\xi)$ for all tangent vectors $\xi$). Exactness is
equivalent to assuming $\int_{S}\beta=0$ for any closed surface $S$ (not just
those that bound a volume). Most discussions of magnetic fields assume this
stronger condition and we shall do so.

An invariant torus for $B$ is called a flux surface. If $B$ is smooth enough
(e.g. $C^{3}$ and the third derivative is Hölder continuous), then various
simple conditions imply a set of positive volume of flux surfaces, for
example, existence of a generic elliptic closed fieldline or existence of one
flux surface with smooth enough conjugacy to a Diophantine rotation (KAM
theory, for a semi-popular introduction, see [D]).

In the design of magnetic confinement devices for plasma, it is often desired
to quantify the volume $V$ enclosed by a flux surface $S$. In particular one
might want to know the volume enclosed by the outermost flux surface of given
class contained within the vacuum vessel, but we pose the question more
generally. Once a flux surface has been calculated, this can be done by 3D
integration, but the question arises whether there is a more efficient way and
that uses the invariance of $S$ under $B$.

One simple improvement is to reduce the problem to a surface integral: if $M$
is contractible then $\Omega$ can be written as $d\nu$ for a 2-form $\nu$ (in
many ways, e.g. for $\Omega=dx\wedge dy\wedge dz$ one can take
$\nu=x\,dy\wedge dz$ or cyclic permutations) and then
$\int_{V}\Omega=\int_{S}\nu$, by Stokes’ theorem. But this still does not use
the invariance under $B$.

One way to use invariance under $B$ is the formula

(1) $V=\int_{D}T\beta,$

where $D$ is any disk transverse to $B$ whose boundary is on the flux surface
$S$, and $T$ is the first-return time to $D$ along the fieldline flow
$\dot{x}=B(x)$. The formula can be obtained from the identity
$\Omega=\frac{B^{\flat}}{|B|^{2}}\wedge\beta$ by decomposing the domain into
the bundle of fieldline segments from $D$ to $D$. The same can be done for the
volume between two nested tori by replacing $D$ by any annulus transverse to
$B$ whose boundary components are on the respective tori; this includes the
limiting case where one of the tori is just an elliptic closed fieldline,
called a magnetic axis. Although (1) is a 2D integral, it involves computing
the return time function, which effectively makes the integration 3D, so there
is no real saving.

Here we give various formulae for the volume enclosed by a flux surface that
use the invariance under $B$ and promise to be more efficient. We build up
from the simplest situation to the general one.

Such an important topic has, of course, been addressed before, e.g. [He] in
the context of magnetohydrostatic (MHS) fields. We will mention connections as
we go along.

##  2\. Quasisymmetric fields

We start by addressing the special case of a quasisymmetric magnetic field
$B$, i.e. there is a vector field $u$ independent from $B$ almost everywhere
(in the region of interest), with $L_{u}\beta=0$, $L_{u}\Omega=0$ and
$L_{u}B^{\flat}=0$. Under mild additional conditions [BKM1], which we assume,
every orbit of $u$ is closed and has the same period $\tau>0$.

The simplest examples are the axisymmetric fields, those for which
$u=\partial_{\phi}$ in cylindrical coordinates (then $\tau=2\pi$). It is an
open question whether there are any other exactly quasisymmetric fields
(compare Kovalevskaya’s tops [K]), but fields with non-axisymmetric
quasisymmetry to a high degree of accuracy can be made [LP].

The conditions $L_{u}\beta=0$, $L_{u}\Omega=0$ imply that $i_{u}\beta$ is a
closed 1-form. Assuming there is no homological obstruction then it is exact,
i.e. there exists a function $\psi$ such that $i_{u}\beta=d\psi$. It follows
that $i_{u}d\psi$ and $i_{B}d\psi=0$, so $u$ and $B$ are tangent to the level
sets of $\psi$. The bounded regular components of level sets of $\psi$ are
tori, because they are orientable surfaces supporting a nowhere-zero vector
field.

Because both $B$ and $u$ are volume-preserving, $L_{u}\beta=0$ is equivalent
to the commutator $[u,B]$ being $0$. Then $(u,B)$ generate an action $\phi$ of
$\mathbb{R}^{2}$ on the flux surfaces; for $t=(t_{u},t_{B})\in\mathbb{R}^{2}$
and $x$ on a flux surface $S$, $\phi_{t}(x)$ is the point reached by flowing
for time $t_{u}$ with $u$ and $t_{B}$ with $B$ (the order does not matter
because the fields commute). We already have that $\phi_{(\tau,0)}$ is the
identity. Choose a $u$-line $\gamma$ and a point $x\in\gamma$. There is a
first $T>0$ such that the flow of $B$ starting at $x$ returns to $\gamma$. $T$
is independent of $x\in\gamma$ because if $y\in\gamma$ then
$y=\phi_{(t_{u},0)}(x)$ for some $t_{u}$, so
$\phi_{(0,T)}(y)=\phi_{(t_{u},T)}(x)=\phi_{(t_{u},0)}(\phi_{(0,T)}(x))$ is in
$\gamma$ and is not for smaller $T$. Furthermore, $T$ does not depend on the
choice of $u$-line $\gamma$ on $S$ (but it does in general depend on the flux
surface $S$).

Define $V$ to be the volume enclosed by a flux surface. It is locally a
function of $\psi$.

######  Theorem 1.

$dV=\tau T(\psi)d\psi.$

###### Proof.

The magnetic flux across an annulus with boundary components $\gamma$ and
$\gamma^{\prime}$ on two neighbouring flux surfaces is $\tau d\psi$. The
volume it sweeps out when following $B$ is the integral of the return time
with respect to the flux. ∎

So one can compute the volume enclosed by a flux surface by integrating
$\frac{dV}{d\psi}$ from some reference $\psi_{0}$ where the volume is known
(e.g. magnetic axis where it is zero), if one computes the return time
function $T(\psi)$. The result is a 2D integration, which is a saving over the
original 3D integration.

##  3\. Fields with weak quasisymmetry

A magnetic field has weak quasisymmetry $u$ if $L_{u}\beta=0$, $L_{u}\Omega=0$
and $L_{u}|B|=0$; this is a translation to differential forms of the
conditions defined in [RHB]. The latter is equivalent to
$i_{B}L_{u}B^{\flat}=0$, so weak quasisymmetry is a generalisation of
quasisymmetry. The status of weak quasisymmetry as the zeroth order conditions
for a velocity-dependent Hamiltonian symmetry of guiding-centre motion is
discussed in [BKM2].

The preceding arguments go through with the exception that the period $\tau$
of the $u$-lines is now in general a function of $\psi$, so it needs computing
too.

######  Theorem 2.

$dV=\tau(\psi)T(\psi)d\psi.$

Thus, computing $V$ is now a 1D integration of a pair of 1D integrations,
still equivalent to 2D.

##  4\. Fields with flux-form symmetry

Next we generalise to magnetic fields $B$ with a volume-preserving field $u$
independent from $B$ almost everywhere (in the region of interest), such that
$L_{u}\beta=0$ (dropping the requirement of $i_{B}L_{u}B^{\flat}=0$). We call
this a flux-form symmetry. Equivalently, $u$ is a volume-preserving field that
commutes with $B$.

This enlarges the previous context to include, for example, any non-degenerate
magnetohydrostatic (MHS) field $B$: simply take $u$ to be the current density
$J=\textrm{curl}\,B$. The MHS property $J\times B=\nabla p$ for pressure field
$p$ implies that $L_{J}\beta=0$. Non-degeneracy means $dp\neq 0$ almost
everywhere, which thereby makes $J$ and $B$ independent almost everywhere (it
should be noted, however, that no non-axisymmetric non-degenerate MHS fields
are known).

As before, assuming there is no homological obstruction then
$i_{u}\beta=d\psi$ for some function $\psi$ (in the MHS case with $u=J$ one
can take $\psi=p$), and the bounded regular components of level sets of $\psi$
are tori invariant under both $B$ and $u$.

Again, $(u,B)$ generate an action $\phi$ of $\mathbb{R}^{2}$ on the flux
surfaces. This implies that for each flux surface there is a lattice $\Gamma$
in $\mathbb{R}^{2}$ (discrete subset closed under subtraction) such that
$\phi_{t}=\textrm{Id}$ iff $t\in\Gamma$. It can no longer be assumed that
there is a $\tau>0$ with $(\tau,0)\in\Gamma$. Let
$T_{1},T_{2}\in\mathbb{R}^{2}$ be generators for $\Gamma$. They can be
computed by integration of ODEs for $u$ and $B$ and root finding. Readers may
recognise this as a strategy to prove the Liouville-Arnol’d theorem (e.g.
[AK]), which gives a parametrisation of $S$ by $\mathbb{R}^{2}/\mathbb{Z}^{2}$
with respect to which both $u$ and $B$ are constant. In the magnetic field
context, it was rediscovered by Hamada [Ha].

Make a matrix $\mathcal{T}$ with the vectors $T_{1},T_{2}$ as columns and let
$\Delta=\det\mathcal{T}$, which is a function of $\psi$.

######  Theorem 3.

$dV=\Delta(\psi)d\psi.$

###### Proof.

In a Liouville-Arnol’d parametrisation of a flux surface by
$(\theta^{1},\theta^{2})\in\mathbb{R}^{2}/\mathbb{Z}^{2}$, $u$ and $B$ have
constant components $(u^{1},u^{2}),(B^{1},B^{2})$, so make a matrix
$\mathcal{V}$ with the components of $u$ and $B$ as columns. Then
$d\psi=i_{u}i_{B}\Omega=\det\mathcal{V}\,i_{\partial_{\theta_{1}}}i_{\partial_{\theta_{2}}}\Omega=\det\mathcal{V}\,dV$.
But the lattice generators satisfy $\mathcal{T}\mathcal{V}=I$, the identity
matrix, so $\det\mathcal{V}\det\mathcal{T}=1$ and hence
$dV=\det\mathcal{T}d\psi$. ∎

From this one can again deduce the volume enclosed by a flux surface, by 1D
integration with respect to $\psi$ from a reference case. As the method
requires a 1D set of 1D integrations (to compute $\Delta(\psi)$), it can again
be considered a 2D integration overall, but that is still a saving compared to
the 3D integration of $\Omega$.

Quasi-symmetric and weak quasisymmetric fields are special cases of this, with
$[T_{1}\,T_{2}]=\left[\begin{array}[]{cc}\tau&c\\\ 0&T\end{array}\right]$ for
some $c$ (the value such that $\phi_{(0,T)}(x)=\phi_{(-c,0)}(x)$).

##  5\. Fields with a foliation by flux surfaces

If there is a foliation by flux surfaces (without knowing if there is a
symmetry $u$) then one can still compute the enclosed volume by an essentially
2D integration, as follows.

For any invariant surface $S$ under $B$, the 2-form $\beta=dA^{\flat}$
vanishes on its tangent spaces. It follows that

$\Phi=\int_{\gamma}A^{\flat}$

is the same for all closed curves $\gamma$ on $S$ in the same homology class.
It can be interpreted as a magnetic flux. Specifically, if $\gamma$ is the
boundary of a disk $D$ then $\int_{\gamma}A^{\flat}=\int_{D}\beta$; if
$\gamma-\gamma_{0}$ is the boundary of an annulus $A$ then
$\int_{\gamma}A^{\flat}-\int_{\gamma_{0}}A^{\flat}=\int_{A}\beta$. Think of
the cases of a poloidal loop and of a toroidal loop relative to a magnetic
axis $\gamma_{0}$.

We restrict attention to the case that $B$ is nowhere zero in the region of
interest. Using volume-preservation and the foliation by invariant tori, one
can deduce that the flow on each torus is of Poincaré type (in the terminology
of [BGKM]), i.e. it has a transverse section such that every $B$-line crosses
it in both directions. Specifically, label the flux surfaces as level sets of
a smooth function $\psi$ with $d\psi\neq 0$, choose any vector field $n$ such
that $i_{n}d\psi=1$ (e.g. choose a Riemannian metric and let
$n=\nabla\psi/|\nabla\psi|^{2}$), and let $\mathcal{A}=i_{n}\Omega$. Then
$\mathcal{A}$ is non-degenerate on each flux surface and preserved by $B$.
This rules out the case of Reeb components (an annulus bounded by periodic
orbits in opposite directions), because area would be contracted on
approaching a bounding periodic orbit, so for a nowhere-zero field, the only
case left is Poincaré type.

Given a closed curve $\gamma$ on a flux surface $S$, transverse to $B$, define
$\bar{T}$ to be the limiting average return time to $\gamma$ of an infinitely
long $B$-line starting at some point of $\gamma$. By the theory of Poincaré
flows on a torus, this average exists.

Furthermore, if $B$ winds irrationally on $S$ it is the same for all initial
conditions on $\gamma$. In this case, $\gamma$ can be deformed to one for
which the return time is constant, but for computational purposes it is enough
to estimate the average for a given $\gamma$. The weighted Birkhoff average
method [SM] is a promising route to do this accurately. Alternatively, we
propose approximating the integral of the return time from its values along
one orbit segment (see section 5.2).

In the case of a rational surface, all trajectories are closed (this follows
from preservation of $\mathcal{A}$ on flux surfaces, specifically
$\int_{\eta}i_{B}\mathcal{A}$ from a reference point to an arbitrary point $x$
of $S$ along any curve $\eta$ on $S$ is preserved by $B$ if $\eta$ is chosen
to depend continuously on $x$ [N]) and we define $\bar{T}$ to be the average
of the return time with respect to the area $\mathcal{A}$. But if the set of
rational surfaces have measure zero then we will not really need this.

If the space is locally foliated by flux surfaces and $\gamma$ is chosen
continuously then we obtain in the same way again

######  Theorem 4.

$dV=\bar{T}d\Phi.$

Note that in the case of a field with a flux-form symmetry, this result is
consistent with that of Theorem 3, because given $u$, construct Liouville-
Arnol’d coordinates $(\theta^{1},\theta^{2})$, suppose the second component
$B^{2}\neq 0$ (else interchange the coordinates), choose $\gamma$ to be
$\theta^{2}=0$, then $d\Phi=i_{\partial_{\theta_{1}}}i_{B}\Omega=B^{2}dV$ and
the return time $T=1/B^{2}$ (constant).

The work required to use Theorem 4 is slightly larger than 2D integration, in
the sense that estimating the mean return time for a torus requires more than
one revolution. But we think it will still be less than a full 3D integration.

If the vector potential $A$ is not provided then computing it at a point from
$B$ is a 1D integral (see the treatment of Poincaré’s lemma in section 6), so
$\Phi$ is a 2D integration and then $V$ is a 3D integration. But one can
replace the computation of $\Phi$ by $d\Phi=\int_{\gamma}i_{Y}\beta$ where $Y$
is a vector field on $\gamma$ pointing to its homologue on an infinitesimally
nearby flux surface, so use of $A$ was not really necessary, and this remains
a 2D integration.

###  5.1. Connections to established results

Let us connect to results surveyed in [He]. Eq. (19) of [He] is

$V^{\prime}(\psi)=\int_{0}^{2\pi}d\alpha\int\frac{d\ell}{|B|},$

where the $\ell$-integral is taken for one toroidal revolution, $\psi$ has
been chosen to be the toroidal flux enclosed by the surface divided by $2\pi$,
and $\alpha$ is a fieldline label such that $\beta=d\psi\wedge d\alpha$
(making $(\alpha,\psi)$ into “Clebsch coordinates” for $B$). Now,
$\int\frac{d\ell}{|B|}$ is the return time along fieldline flow to a poloidal
section, the average with respect to $\mathcal{A}$ is the average with respect
to $\alpha$, and $\Phi=2\pi\psi$. So the formula agrees with that of Theorem
4.

Eq. (32) of [He] says that in non-degenerate MHS, the period of all fieldlines
on a rational surface is the same (Hamada condition), so “all flux tubes
carrying the same magnetic flux on a given surface must have the same volume”.
Thus, $\bar{T}$ can be replaced by the period divided by the number of
toroidal revolutions, which recovers the result in the equation before (22) in
[He]. The Hamada condition is a particular consequence of the Liouville-
Arnol’d theorem used in section 4. A more general result for non-degenerate
MHS, cited just after (30) in [He], is that the Jacobian of the Hamada
coordinate system is a flux function. Indeed, Hamada coordinates are ones in
which $J$ and $B$ are constant on each flux surface and the Jacobian is
$\Delta(\psi)$ from Theorem 3.

###  5.2. Average return time from one orbit segment

Here is the proposed method to estimate the mean return time from one orbit
segment (more accurately than just taking the average of the return time). It
applies to tori with irrational winding ratio $\iota$. Then there is a
parametrisation of $\gamma$ by $\theta\in\mathbb{R}/\mathbb{Z}$ such that the
first return map to $\gamma$ is $\theta^{\prime}=\theta+\iota$. Express the
return time as a function $T(\theta)$. We desire its average $\bar{T}$ with
respect to $\theta$, i.e. $\bar{T}=\int T(\theta)\,d\theta$. Given an orbit
segment with successive return times $T_{n},n=1,\ldots,N$, we can label the
initial point by $\theta=0$ and then $T_{n}=\tau((n-1)\iota)$. The values
$(n-1)\iota$ can be reduced to the interval $[0,1)$ and
$\int\tau(\theta)\,d\theta$ approximated by the trapezoidal rule applied to
these points.

The method requires to first estimate $\iota$. This can be done by choosing
some coordinate $\phi:\gamma\to\mathbb{R}/\mathbb{Z}$, with initial point
$\phi=0$, and looking at the order in which the successive returns $\phi_{n}$
come. Let $n_{k}$, $k\geq 1$, be the indices for the successively closest
returns to the initial point (they will be on alternate sides of $0$) and
$m_{k}$ be the closest integer to the number of revolutions of $\gamma$ made
(this requires choosing the same convention as for $\iota$). Then $\iota$ can
be estimated from the linear approximation to $\phi(\theta)$ near $\theta=0$:

$\frac{n_{k}\iota-m_{k}}{n_{k-1}\iota-
m_{k-1}}\sim\frac{\phi_{n_{k}}}{\phi_{n_{k-1}}}.$

So

$\iota\sim\frac{m_{k}\phi_{n_{k-1}}-m_{k-1}\phi_{n_{k}}}{n_{k}\phi_{n_{k-1}}-n_{k-1}\phi_{n_{k}}}.$

Note that there will be a sequence $a_{k}$ of positive integers such that
$n_{k+1}=a_{k}n_{k}+n_{k-1}$ and the same for $m_{k+1}$ (defining
$n_{0}=1,m_{0}=0$). This sequence gives the continued fraction expansion of
$\iota$.

##  6\. General fields

For a general magnetic field there may be a set of flux surfaces of given
class, interspersed with island chains and “chaos” (the class can be defined
relative to a magnetic axis, for example).

Can one compute the integral $\int_{D^{\prime}}T\beta$ of the return time $T$
over the part $D^{\prime}$ of a transverse section $D$ that is inside a flux
surface in a more efficient way than the obvious 2D integral of the result of
1D integrations?

Yes. $T\beta$ is a top-form on $D$. If $D$ is contractible (e.g. a poloidal
section) then choose a vector field $X$ on $D$ whose flow $\phi$ contracts $D$
to a point as $t\to+\infty$, and for any tangent $\xi$ to $D$, let
$\eta(\xi)=-\int_{0}^{\infty}T\beta(\phi_{t*}X,\phi_{t*}\xi)\,dt$ (this is one
way to prove Poincaré’s lemma). Then $T\beta=d\eta$ on $D$ and so for any
subset $D^{\prime}\subset D$, $\int_{D^{\prime}}T\beta=\int_{\partial
D^{\prime}}\eta$.

To make it concrete, suppose coordinates $(u,v)$ on $D$ such that the forward
flow of $X(u,v)=(-u,-v)$ is defined on $D$. It is just
$\phi_{t}(u,v)=e^{-t}(u,v)$. Then for a tangent vector $\xi$ at $r=(u,v)$,
$\eta(\xi)=\int_{0}^{\infty}e^{-2t}T\beta(r,\xi)\,dt$, where $T$ and $\beta$
are evaluated at $e^{-t}r$. Recall that $\beta=i_{B}\Omega$, so
$\beta(r,\xi)=\Omega(B,r,\xi)$, the triple product in the Euclidean case.

Thus, we have reduced the calculation of the volume enclosed by a flux surface
to a 1D integral (around $\partial D^{\prime}$) of a 1D integral $\eta(\xi)$
with $\xi$ the tangent to $\partial D^{\prime}$. This is essentially a 2D
integral.

If $D$ is not contractible, e.g. an annulus between a closed curve on $S$ and
a magnetic axis $\gamma_{0}$, then can instead contract to $\gamma_{0}$ and
obtain a similar reduction to integration of a 1-form, that we do not spell
out here.

But the question remains whether this is any better than integrating a
primitive of $\Omega$ over the surface.

##  7\. Fluxes

It is interesting to compare computation of enclosed volume with the question
of computing the fluxes of a flux surface.

Given a “poloidal” closed curve $\gamma$ on a flux surface $S$ bounding a disk
$D$ in $M$, its flux is $\Phi=\int_{D}\beta$. This a 2D integral. It can be
written as $\Phi=\int_{\gamma}A^{\flat}$, reducing the computation to a 1D
integral (assuming that a vector potential $A$ is given, else its computation
is already a 1D integration and no saving results). By exactness of $\beta$,
$\Phi$ is independent of the disk $D$ chosen to span $\gamma$. By invariance
of the flux surface, it is also independent of the representative $\gamma$ of
its homology class in $H_{1}(S)$. This is because $\beta=0$ on all pairs of
tangent vectors to $S$.

To define the flux for a general homology class on a flux surface requires
some more discussion. For a flux surface $S$ bounding a solid torus with $B$
nowhere-zero inside, there is a closed fieldline inside. Choosing one,
designated the magnetic axis $\gamma_{0}$, one can define the flux for a non-
poloidal closed curve $\gamma$ on $S$ to be $\int_{A}\beta$ where $A$ is the
annulus with boundary $\gamma-\gamma_{0}$. Again, this can be reduced to 1D
integrals $\int_{\gamma}A^{\flat}-\int_{\gamma_{0}}A^{\flat}$.

As already mentioned, for a field with flux-form symmetry $u$ with closed
$u$-lines, $d\Phi=\tau d\psi$, where $\tau(\psi)$ is the period of the
$u$-lines.

###  7.1. Percival’s variational principle

The fluxes are associated with one strategy for computing flux surfaces,
namely an adaptation of Percival’s variational principle for invariant tori of
Hamiltonian systems [P1] to the magnetic field context (see [GKN] for another
extension of Percival’s variational principle).

Let $A$ be a vector potential for $B$. Given
$\omega\in\mathbb{R}^{2}\setminus\\{0\\}$ and a differentiable mapping
$x:\mathbb{T}^{2}=\mathbb{R}^{2}/\mathbb{Z}^{2}\to M$, define $P$ (for
Percival) by

$P_{\omega}(x)=\int A_{j}(x(\theta))\frac{\partial
x^{j}}{\partial\theta^{i}}(\theta)\omega^{i}\,d^{2}\theta.$

The first variation of $P_{\omega}(x)$ with respect to a variation $\delta x$
is

$\delta P_{\omega}=\int A_{j,k}\delta x^{k}x^{j}_{,i}\omega^{i}+A_{j}\delta
x^{j}_{,i}\omega^{i}\,d^{2}\theta,$

where subscript ,i denotes partial derivative in direction $i$. The second
term can be integrated by parts to produce

$\delta P_{\omega}=\int(A_{j,k}-A_{k,j})x^{j}_{,i}\omega^{i}\delta
x^{k}\,d^{2}\theta=\int\varepsilon_{jkl}B^{l}x^{j}_{,i}\omega^{i}\delta
x^{k}\,d^{2}\theta,$

where $\varepsilon$ is the Levi-Civita symbol. This is zero for all variations
$\delta x$ iff $x^{j}_{,i}\omega^{i}$ is parallel to $B^{j}$, i.e. there is a
function $c:\mathbb{T}^{2}\to\mathbb{R}$ such that

$x^{j}_{,i}(\theta)\omega^{i}=c(\theta)B^{j}(x(\theta)).$

In particular, if $x$ is a diffeomorphism then the image of $x$ is a flux
surface (furthermore, $c$ is nowhere zero).

In more detail, this says that the function $x$ conjugates
$\dot{\theta}=\omega/c(\theta)$ to $\dot{x}=B(x)$ on the flux surface.
Furthermore, assuming $B$ is nowhere zero on the surface, the fieldlines on
the surface have winding ratio $\iota=\omega^{2}/\omega^{1}$ with respect to
$(\theta^{1},\theta^{2})$. The $\theta$ motion has invariant probability
density $c(\theta)/\bar{c}$, where $\bar{c}$ is the average of $c$ with
respect to $d^{2}\theta$.

For irrational winding ratio, we can use ergodicity. $P_{\omega}(x)$ is the
$\theta$-average of $A\cdot\frac{dx}{d\tau}$ for
$\frac{d\theta}{d\tau}=\omega$. Transforming time by $dt=c(\theta)d\tau$, we
see that $P/\bar{c}$ is the average of the helicity $A\cdot B$ with respect to
the invariant probability on $S$. By unique ergodicity of the $\theta$ motion,
this is the $t$-average of $A\cdot B$ along any fieldline on $S$. Similarly,
$\bar{c}$ has the interpretation that $\omega_{j}/\bar{c}$ is the time-
averaged numbers of revolutions $\Delta\theta_{j}$ along a fieldline.

Note that both $P$ and $c$ are homogeneous of first degree in $\omega$. So one
can consider $P/\bar{c}$ to be a function of $\iota$.

In the differentiable context, there is no reason for $P$ to have a critical
point. A main point of Percival’s variational principle is that it allows
extension to functions $x$ that are differentiable in direction $\omega$ but
not necessarily continuous. Indeed, Mather used the principle to prove
existence of a set $M_{\rho}$ of quasiperiodic orbits for each irrational
rotation number $\rho$ for area-preserving twist maps [Ma1]. $M_{\rho}$ is
either a circle or a Cantor set, thus providing an answer in this context to
the question of what happens to invariant tori beyond where KAM theory
applies. The case of a Cantor set was christened “cantorus” by Percival [P2].

One can expect the same conclusion to apply to magnetic fields as long as they
have shear, meaning that in some coordinate system $(r,\phi_{1},\phi_{2})$,
$\frac{\partial}{\partial r}\frac{B^{\phi_{2}}}{B^{\phi_{1}}}\neq 0$.

One consequence of the variational principle is that if one has computed the
stationary value of $P_{\omega}$ for each $\omega$ (it is enough to do for
each $\iota$ by homogeneity in $\omega$) then one can read off the toroidal
and poloidal fluxes $\Phi_{j}$ as follows. By invariance of the flux surface,
the integral of $A^{\flat}$ along a long fieldline can be deformed into the
weighted sum of its integrals around a basis for cycles on the flux surface,
$P_{\omega}=\Phi_{1}\omega_{1}+\Phi_{2}\omega_{2}$. $P_{\omega}$ is stationary
with respect to variations in $x$, so assuming differentiability of the
critical value of $P_{\omega}$ with respect to $\omega$ (which could probably
be proved along the lines of [Ma2]) we deduce that $\frac{\partial
P}{\partial\omega_{j}}=\Phi_{j}$.

This formula is an analogue of that for the area under an invariant circle or
cantorus for area-preserving twist maps as derivative of mean action with
respect to rotation number [Ch], and the equivalent one for chemical potential
of an incommensurate structure for Frenkel-Kontorova chains as derivative of
the mean energy with respect to the mean spacing [Au].

##  8\. Computing flux surfaces

The paper begs the question of how to compute flux surfaces. This might depend
on the context, but one way is Percival’s variational principle. It can be
formulated as looking for a critical point of a functional in a space of
double periodic functions, which could be represented by Fourier series. One
may not know in advance, however, whether there exists a smooth solution. If
the field has shear then one can do a partial stationarity to replace a minor
radius coordinate by first derivatives with respect to the others, rendering a
functional that is bounded below and hence amenable to minimisation methods.
But a minimiser might still not be smooth.

Another method is to look for flux surfaces as limits of sequences of periodic
orbits. Again, this might not give a flux surface, but if the “residues” of
the periodic orbits go to zero, they are quite likely to accumulate on a flux
surface. A formal result is available in the other direction (if there is a
smooth flux surface with shear, then there are sequences of periodic orbits
that converge to it with residues going to zero [M2]). Outermost flux surfaces
can be found as boundaries of limits of sequences of periodic orbits whose
residues are bounded [GMS].

##  9\. Conclusion

The paper has presented methods to reduce the computation of the volume
enclosed by a flux surface of a magnetic field from 3D to 2D integration.

A question remains: are any of the methods more efficient than the first one
that does not use invariance of the flux surface under the field? The answer
is likely to depend on how the field is presented (e.g. directly or via a
vector potential) and how the flux surfaces are presented (e.g. as a graph of
minor radius against two angle variables, as the closure of the orbit of one
point…).

It would be interesting to test the methods numerically.

## Acknowledgements

This work was supported by a grant from the Simons Foundation (601970, RSM). I
am grateful to David Martinez for detailed comments.

## References

  * [AK] Arnol’d VI, Khesin BA, Topological methods in hydrodynamics (Springer, 1998).
  * [Au] Aubry S, The devil’s staircase transformation in incommensurate lattices, in: The Riemann problem, complete integrability and arithmetic applications, eds Chudnovsky, Chudnovsky, Springer Lect Notes Math 925 (1982) 221–245.
  * [BGKM] Baesens C, Guckenheimer J, Kim S, MacKay RS, Three coupled oscillators: Mode-locking, global bifurcations and toroidal chaos, Physica D 49 (1991) 387–475
  * [BKM1] Burby JW, Kallinikos N, MacKay RS, Some mathematics for quasisymmetry, J Math Phys 61 (2020) 093503
  * [BKM2] Burby JW, Kallinikos N, MacKay RS, Approximate symmetries of guiding centre motion, J Phys A 54 (2021) 125202
  * [Ch] Chen Q, Area as a devil’s staircase in twist maps, Phys Lett A 123 (1987) 444–450.
  * [D] Dumas HS, The KAM Story: A Friendly Introduction to the Content, History and Significance of Classical Kolmogorov-Arnold-Moser Theory (World Scientific Publishing Company, 2014)
  * [GKN] Gray CG, Karl G, Novikov VA, The four variational principles of mechanics, Ann Phys 251 (1996) 1–25.
  * [GMS] Greene JM, MacKay RS, Stark J, Boundary circles for area-preserving maps, Physica D 21 (1986) 267–295.
  * [Ha] Hamada S, Hydromagnetic equilibria and their proper coordinates, Nucl Fusion 2 (1962) 23–37.
  * [He] Helander P, Theory of plasma confinement in non-axisymmetric magnetic fields, Rep Prog Phys 77 (2014) 087001.
  * [K] Kowalevski S, Sur le problème de la rotation d’un corps solide autour d’un point fixe, Acta Math 12 (1889) 177–232.
  * [LP] Landreman M, Paul E, Magnetic fields with precise quasisymmetry for plasma confinement, Phys Rev Lett 128 (2022) 035001
  * [M1] MacKay RS, Differential forms for plasma physics, J Plasma Phys 86 (2020) 925860101
  * [M2] MacKay RS, Greene’s residue criterion, Nonlinearity 5 (1992) 161–187.
  * [Ma1] Mather JN, Existence of quasi-periodic orbits for twist homeomorphisms of the annulus, Topology 21 (1982) 457–467.
  * [Ma2] Mather JN, Differentiability of the minimal average action as a function of the rotation number, Bol Soc Bras Mat 21 (1990) 59–70.
  * [N] Newcomb WA, Magnetic differential equations, Phys Fluids 2 (1959) 362–5.
  * [P1] Percival IC, A variational principle for invariant tori of fixed frequency, J Phys A 12 (1979) L57–60.
  * [P2] Percival IC, Variational principles for invariant tori and cantori, in: Nonlinear dynamics and the beam-beam interaction, eds Month M, Herrera JC, AIP conf proc 57 (1980) 302–310.
  * [RHB] Rodriguez E, Helander P, Bhattacharjee A, Necessary and sufficient conditions for quasisymmetry, Phys Plasmas 27 (2020) 062501
  * [SM] Sander E, Meiss JD, Birkhoff averages and rotational invariant circles for area-preserving maps, Phys D 411 (2020) 132569


# Comparative Analysis of Sub-band Allocation Algorithms in In-body Sub-
networks Supporting XR Applications
††thanks: The work by Thomas Jacobsen and Ramoni Adeogun was supported by the
HORIZON-JU-SNS-2022-STREAM-B-01-03 6G-SHINE project (grant agreement No.
101095738). Ramoni Adeogun’s work was also partly supported by HORIZON-JU-
SNS-2022-STREAM-B-01-02 project - CENTRIC (grant agreement No. 101096379)

Saeed Bagherinejad1, Thomas Jacobsen2, Nuno K. Pratas2 and Ramoni O. Adeogun1
1Department of Electronic Systems, Aalborg University, 9220 Aalborg, Denmark
2Nokia Standards, 9220 Aalborg, Denmark  Email: 1{sbag<EMAIL_ADDRESS>2{nuno.kiilerich_pratas<EMAIL_ADDRESS>
###### Abstract

In-body subnetworks (IBS) are envisioned to support reliable wireless
connectivity for emerging applications including extended reality (XR) in the
human body. As the deployment of in-body sub-networks is uncontrollable by
nature, the dynamic radio resource allocation scheme in place becomes of the
uttermost importance for the performance of the in-body sub-networks. This
paper provides a comparative study on the performance of the state-of-the-art
interference-aware sub-band allocation algorithms in in-body sub-networks
supporting the XR applications. The study identified suitable models for
characterizing in-body sub-networks which are used in a snapshot-based
simulation framework to perform a comprehensive evaluation of the performance
of state-of-art sub-band allocation algorithms, including greedy selection,
sequential greedy selection (SG), centralized graph coloring (CGC), and
sequential iterative sub-band allocation (SISA). The study shows that for XR
requirements, the SISA and SG algorithms can support IBS densities up to 75%
higher than CGC.

###### Index Terms:

In-body sub-networks, In-X subnetworks, Sixth Generation (6G), Extended
Reality (XR), Sub-band Allocation

##  I Introduction

With the commercialization of the fifth generation (5G) of radio technologies,
the definition of the next generation has already begun. As before, the
motivation behind the development of sixth-generation (6G) wireless networks
is to enable the proliferation of more demanding services and technologies
that previous generations, including 5G, do not fully support. Fully-immersive
and multi-sensory eXtended Reality (XR) is among the main services that 6G is
expected to deliver and it includes various use cases such as Augmented
Reality (AR), Virtual Reality (VR), and Cloud Gaming (CG) [1, 2]. In these use
cases, the stringent requirements on the wireless networks originate from the
wireless transmission of high quality and delay sensitive video frames to the
XR display device (XRDD) from external computation resources such as
smartphones [3]. We envision that in-body sub-networks (IBS) could be a
potential platform for supporting the XR requirements.

IBSs are a part of the umbrella paradigm called 6G in-X sub-networks [4, 5]
which includes networks of specialized short-range low-power cells that serve
different types of use-cases inside entities such as vehicles, robots,
classrooms or the human body. Among them, IBS deals with wireless
communications in the vicinity of the human body. It should be noted that,
although there are existing standards such as IEEE 802.16.5 for wireless body
area networks (WBAN) [6], their supported quality of service (QoS) are nowhere
near the requirements of XR applications, e.g., average data rate of
$30-45$Mbps, packet delay budget (PDB) of 10ms, and 99% reliability.
Supporting such demanding QoSs compared to state-of-the-art is the main
motivation for developing IBS.

A key challenge for in-body sub-networks is interference, particularly in
situations where a high number of people (IBSs) are assembled in a limited
area, e.g., people gathering for a collective immersive social experience
using XR services. There has been some works in the past on interference
mitigation in 6G in-X sub-networks with a focus on in-factory sub-networks,
see e.g. [7, 8, 9, 10, 11, 12, 13]. In [7], a distributed and centralized
dynamic sub-band group selection algorithm has been proposed for the in-X sub-
networks based on the sensed interference level. A deep-learning-based
algorithm has been proposed in [8] for distributed dynamic sub-band allocation
which relies on offline training in a controller. In [10], a multi-agent
Q-learning approach is presented to joint sub-band and power allocation which
relies on limited sensing information in each subnetwork. Multi-agent
reinforcement learning-based algorithm in [12] utilizes received signal
strength indicator (RSSI) for joint sub-band and power allocation. A
centralized algorithm was proposed in [9] for sub-band allocation in sub-
networks which depends on the interference to signal ratio (ISR) level
reported from sub-networks. The authors of [13] introduced a hybrid framework
for resource management in in-X sub-networks in which the resource allocation
is done centrally or in a distributed manner based on the link state to the
central controller and the battery level of the AP.

As mentioned earlier, the development of interference mitigation algorithms
has been majorly focused on in-factory sub-networks. The use cases that IBS is
envisioned to support, such as XR, have different data traffic patterns and
also QoS levels. In addition, the signal propagation around the human body is
different from other environments. These necessitate a study on the
effectiveness of the resource allocation algorithms in the context of the in-
body sub-networks. To enable the comparative evaluation, a methodology must be
developed for modeling the IBS systems as a platform for supporting XR
applications. Thus, the main contributions of this work include: i) Developing
a method to model IBS networks that captures the most important aspects of
IBSs and their use cases, and ii) conducting a comparative evaluation on the
efficacy of different centralized and distributed sub-band allocation
algorithms in IBS networks. In particular, at first, the appropriate modeling
is defined for IBSs, each including a smartphone (acting as an AP and the
application server for video rendering) and an XRDD. Then, the reliability of
IBS networks is studied in supporting the video frame streams transmitted on
the downlink (DL), from smartphones to the XRDDs. The paper is organized as
follows: we start with the system model in Section II, followed by problem
formulation and a brief description of sub-band allocation algorithms in
Section III. In Section IV, the numerical evaluation methodology and
simulation results are presented and the conclusion is included in Section V.

Figure 1: Deployment and shape model of $N$ In-Body sub-networks.

##  II System Model

###  II-A Deployment and Channel Models

We consider a deployment comprising of $N$ in-body sub-networks (IBS) in an
$L\times L$ indoor environment. Each IBS contains one smartphone acting as the
AP which can coordinate with the other APs and/or the central controller. The
smartphone has also the computational power to run the XR application server,
as suggested in [3]. Furthermore, we assume that each IBS comprises of one
head-mounted XRDD that displays the received frames. Note that in general, an
IBS contains more devices such as haptic sensors, optical sensors,
smartwatches, and so forth that the AP would utilize for rendering the
appropriate frames. However, we focus our study on the video frame
transmission in the DL from AP to XRDD which is the most important and
demanding data flow of XR use cases compared to sensing feedback from other
devices. That is why we only consider only an AP and an XRDD in each IBS.
Furthermore, We assume that IBSs are randomly located following the Thomas
Cluster Process [14] where IBSs are located around different attraction
points. These attraction points are drawn from a Homogeneous Poisson Point
Process $\Phi_{p}$ with the density, $\lambda_{p}$ in the $\mathbb{R}^{2}$
plane. Then, the central points of IBSs are deployed as the offspring points
of these attraction points where the distance from an IBS’s center point to
its cluster center (attraction points) follows a normal distribution

$f_{n}(x,y|(x_{pn},y_{pn}))=\frac{1}{\sigma_{n}\sqrt{2\pi}}e^{-\frac{(x-x_{pn})^{2}+(y-y_{pn})^{2}}{2\sigma_{n}^{2}}},$
(1)

where the $(x_{pn},y_{pn})$ denotes the coordinates of the attraction point
and $\sigma_{n}$ is the standard deviation of IBSs relative distance from
parent points. It is worthwhile to mention that we enforce constraints on the
minimum distance between the clusters and the IBSs, i.e., $d_{c}^{min}$ and
$d_{IBS}^{min}$, respectively.

We consider the shape of an IBS as cylinders with a radius $R$ and height of
$H$ as depicted in Figure 1. Furthermore, IBSs use a dedicated spectrum with a
total bandwidth of $B_{t}$ which is divided into $K$ equally-sized sub-bands
as depicted in Figure 2. The channel from the $n$-th XRDD and the $l$-th AP
over the $k$-th sub-band is modeled as follows:

$\displaystyle h_{nl}^{k}=\sqrt{\beta_{nl}}g_{nl}^{k},$ (2)

where $\beta_{nl}$ is the large-scale fading coefficient and
$g_{nl}^{k}\sim\mathcal{CN}(0,1)$ denotes the small-scale fading coefficient
over the different sub-bands.

Figure 2: Frame structure and data traffic model for $T_{\nu}=3T_{f}$.

###  II-B XR application and data traffic model

The focus of our analysis is on the transmission of video frames on the DL
which constitutes the most important data flow in XR applications. We have
adopted the single eye buffer model where the video frames for both eyes
arrive at the same time as a single packet [15, 2]. According to this periodic
traffic model, packets arrive with the arrival rate equal to the frame rate,
$f_{\nu}$ which corresponds to a packet inter-arrival time of
$T_{\nu}=1/f_{\nu}$. Note that, we have considered the traffic model without
jitter in the arrival of the packets. This is due to the fact that the jitter
caused by changes in the routing of the packets over the network is not
present for IBS networks. Moreover, the PDB determines the maximum over-the-
air latency in which the packets must be transmitted successfully otherwise
they are considered as dropped. The PDB for XR applications has been defined
to be in the range of $5-15\text{ms}$ depending on the application, e.g.,
$10\text{ms}$ for AR/VR and $15$ms for CG. The packet size, $S_{\nu}$, is
determined by the average data rate $R_{\nu}$ and the packet arrival rate
$f_{\nu}$, i.e., $S_{\nu}=R_{\nu}/f_{\nu}$. The average data rate for the XR
data source is suggested to be in the range of $R_{\nu}=30-45\text{Mbps}$.
Furthermore, the baseline requirement for the packet success rate is set to
$99$% which translates to a packet error rate (PER) of $10^{-2}$. Note that,
this reliability constraint is with respect to the PDB.

###  II-C Transmission Frame Structure

As depicted in Figure 2, we assume a time division duplex structure where
$T_{f}$ denotes the duration of the transmission time interval (TTI). The
duration of the downlink transmission (from AP to XRDD) is denoted by
$T_{d}=\eta T_{f},0\leq\eta\leq 1$ and the $(1-\eta)T_{f}$ is dedicated for
uplink transmission (from sensors and XRDD to AP). We do not consider any re-
transmission schemes and thus, a video frame is either transmitted
successfully in a single TTI or dropped. With this assumption, the TTI
duration, $T_{f}$, represents the over-the-air latency and therefore PDB of
the packets.

As depicted in Figure 2, only a subset of all IBSs (referred to as active
IBSs) have packets to transmit at each TTI. This stems from the periodicity of
the traffic model and the fact that $T_{f}\leq T_{\nu}$. To determine if the
$n$-th IBS is active or not, we use a binomial random variable $b_{n}$. The
distribution of $b_{n}$ depends on the packet inter-arrival time and the frame
duration, i.e. $Pr\\{b_{n}=1\\}=\frac{T_{f}}{T_{\nu}}$. To prove this
probability, consider a sufficiently large time interval $T_{0}$. In this
interval, there would be $T_{0}/T_{f}$ TTIs and $T_{0}/T_{\nu}$ packets;
hence, out of $T_{0}/T_{f}$ TTIs, $T_{0}/T_{\nu}$ of them would have a packet
to transmit, i.e., $Pr\\{b_{n}=1\\}=\frac{T_{0}/T_{\nu}}{T_{0}/T_{f}}$.
Moreover, the set of the active IBSs for a frame is defined as
$\mathcal{A}=\\{n|b_{n}=1,\;1\leq n\leq N\\}$. We assume that this set is
known for the central controller.

##  III Sub-Band Allocation Algorithms

In this section, we first describe the resource allocation problem that will
be studied in this paper. Then, we will look into the sub-band allocation
algorithms that were analyzed.

###  III-A Problem description

The objective of this work is to study the reliability of the IBS systems, in
terms of PER, in supporting the video frame streams of XR applications in the
DL. As discussed earlier, we assume that the rendering process is done on the
AP of the IBS. The rendered video frames are transmitted to the XRDD for
viewing by the user. The SINR expression of $n$-th active IBS
($n\in\mathcal{A}$) over $k$-th sub-band is defined as

$\rho_{n}^{k}=\frac{a_{n}^{k}P_{n}|h_{nn}^{k}|^{2}}{\sum_{\begin{subarray}{c}i\in\mathcal{A}\\\
i\neq
n\end{subarray}}a_{i}^{k}P_{i}|h_{ni}^{k}|^{2}+\sigma_{w}^{2}},\;n\in\mathcal{A}$
(3)

where $P_{n}$ denotes the transmit power and
$\sigma_{w}^{2}=10^{(-204+NF+10\log_{10}B_{s})/10}$ is the additive white
Gaussian noise (AWGN) power. Furthermore, the binary indicator
$a_{i}^{k}\in\\{0,1\\}$ with $\sum_{k=1}^{K}a_{i}^{k}=1$ determines if the
$k$-th sub-band is allocated to $k$-th IBS or not. Using the Shannon capacity
formula, the throughput of $n$-th IBS over one TTI can be described as:

$\mathcal{R}_{n}=\eta
T_{f}B_{s}\sum_{k=1}^{K}a_{n}^{k}\log_{2}(1+\rho_{n}^{k}),\;n\in\mathcal{A}.$
(4)

Since no re-transmission is considered, a transmission would be successful if
the throughput is greater than or equal to the packet size, i.e.,

$\eta T_{f}B_{s}\sum_{k=1}^{K}a_{n}^{k}\log_{2}(1+\rho_{n}^{k})\geq S_{\nu}.$
(5)

This condition is used to evaluate the reliability performance of the sub-band
allocation algorithms, i.e., the percentage of packets that are successfully
transmitted. Using the constraint on sub-band allocation and denoting
$\rho_{th}=2^{\frac{S_{\nu}}{\eta T_{f}B_{s}}}$, the condition in (5)
translates to a condition on SINR for successful packet transmission, i.e.,

$\rho_{n}^{i}\geq\rho_{th},$ (6)

if the $i$-th sub-band is allocated to $n$-th active IBS.

###  III-B Reference Signal Transmission

We assume that each IBS is assigned orthogonal reference signals in each sub-
band in each TTI. Before starting the data transmission stage in each TTI, the
active IBSs transmit their assigned reference signal with a fixed power. Then,
the XRDDs measure the reference signal received powers (RSRP) on each sub-band
and send them to their corresponding APs. Depending on whether the sub-band
allocation algorithm is distributed or centralized, the APs either select a
sub-band or report the measured RSRPs to the centralized controller. We denote
the RSRP of the $n$-th XRDD from the $i$-th IBS on the $k$-th sub-band as
$q_{ni}^{k}$, and the RSRP from the AP of the same IBS as $p_{n}^{k}$. The
RSRP matrix of the $n$-th IBS is denoted by
$\textbf{Q}_{n}\in\mathcal{R}^{K\times N}$ where

$[\textbf{Q}_{n}]_{ki}=\begin{cases}q_{ni}^{k},&n\neq i,\;i\in\mathcal{A},\\\
p_{n}^{k},&n=i,\;i\in\mathcal{A},\\\ 0,&i\notin\mathcal{A}.\end{cases}$ (7)

In what follows, we will discuss how the centralized and distributed sub-band
allocation algorithms utilize the available information to allocate sub-bands
to each IBS.

###  III-C Centralized Sub-band Allocation

####  III-C1 Centralized Graph Coloring (CGC)

In this algorithm, a graph $\mathcal{G}$ is first defined by considering the
sets of vertices $\mathcal{V}$ and edges $\mathcal{E}$. The vertices represent
the IBSs and edges are added based on the available information. Then, using
the completed graph $\mathcal{G}(\mathcal{V},\mathcal{E})$, a vertex-coloring
algorithm is used to assign a color (equivalent to a sub-band) to each vertex.
The vertex-coloring assigns colors in a way so that no neighboring vertices
have the same color. Thus, a vertex must be connected to the IBSs with the
highest interference power so that they would be assigned different sub-bands.
One way of creating the edges is using the RSRPs received from the IBSs. For
this purpose, the average interference-to-signal-ratio (ISR) of the $n$-th IBS
from $l$-th IBS is calculated over the $K$ sub-bands:

$\bar{W}_{n}(l)=\frac{1}{K}\sum_{k=1}^{K}\frac{q_{nl}^{k}}{p_{n}^{k}}.$ (8)

After calculating the average ISRs, the $n$-th vertex will be connected to the
$K-1$ IBSs with the highest average ISRs. In addition, the averaging could be
done locally on the APs which will decrease the signaling overhead of the CGC
algorithm. Furthermore, the greedy vertex-coloring algorithms try to assign
the minimum number of colors to the vertices. To ensure that the graph,
$\mathcal{G}$, is colorable by the available number of colors (corresponding
to the number of sub-bands), we propose to add the edge $(i^{*},j^{*})$ to the
graph as follows,

$(i^{*},j^{*})=\underset{\underset{(i,j)\notin\mathcal{E}}{i,j\in\mathcal{A}}}{\text{arg}\,\text{max}\,\,}\bar{W}_{i}(j),$
(9)

if the number of assigned colors is less than $K$. Otherwise, we remove the
edge

$(i^{*},j^{*})=\underset{\underset{(i,j)\in\mathcal{E}}{i,j\in\mathcal{A}}}{\text{arg}\,\text{min}\,\,}\bar{W}_{i}(j).$
(10)

The process of adding or removing the edges is done sequentially until the
number of assigned colors equals the number of sub-bands.

###  III-D Sequential Iterative Sub-band Allocation (SISA)

The sequential iterative sub-band allocation (SISA) [9] focuses on minimizing
the sum ISR. The ISR level is calculated based on the RSRP described in the
Subsection III-A as follows:

$W_{n}(l,k)=\frac{q_{nl}^{k}}{p_{n}^{k}},\;1\leq l\leq N.$ (11)

This measurement matrix of size $N\times K$ is calculated at the AP and then
forwarded to the central controller. After receiving the ISRs from all IBSs,
the controller runs the sequential iterative sub-band allocation algorithm
which is described in Algorithm 1. This algorithm starts first with a random
allocation that we denote by the mapping function $\mathcal{C}^{0}(n),\forall
n\in\mathcal{A}$. Moreover,
$\mathcal{B}^{0}_{k}=\\{n|n\in\mathcal{A},\mathcal{C}^{0}(n)=k\\}$ denotes the
set of active users with the initial allocated sub-band of $k$. After
initialization, the algorithm goes through the active IBSs sequentially and
chooses the sub-band with the lowest mutual ISR according to the current
allocation which is described in line 5 of the algorithm. After selecting the
sub-band, the allocation is updated at line 6 before moving to the next IBS.
This process is carried out for $M$ iterations.

Algorithm 1 Sequential Iterative Sub-band Allocation

1:$W_{n}(l,k),n,l=1,\ldots,N,\;k=1,\ldots,K$

2:Initialize: $\mathcal{C}^{0}(n),\forall n\in\mathcal{A}$,
$\mathcal{B}^{0}_{k}=\\{n|n\in\mathcal{A},\mathcal{C}^{0}(n)=k\\}$

3:for $m=1,2,\ldots,M$ do

4: for $n=1,2,\ldots,|A|$ do

5: $t=|A|(m-1)+n$

6: $\mathcal{C}^{t}(n)=\underset{1\leq k\leq
K}{\text{arg}\,\text{min}}\,\,\sum_{\underset{l\neq
n}{l\in\mathcal{B}^{t-1}_{k}}}W_{n}(l,k)+W_{l}(n,k)$

7: $\mathcal{B}^{t}_{k}=\\{n|n\in\mathcal{A},\mathcal{C}^{t}(n)=k\\}$

8: end for

9:end for

10:$\mathcal{B}^{d}_{k}$, where $d=|A|\times M$

###  III-E Distributed Sub-band Allocation

####  III-E1 Greedy Selection

Greedy selection is the simplest distributed algorithm where IBSs after
receiving reference signals, independently select a $k^{*}$ as follows:

$k^{*}=\underset{k=1,...,K}{\text{arg}\,\text{min}\,\,}\sum_{l\in\mathcal{A}}\frac{q_{nl}^{k}}{p_{n}^{k}}.$
(12)

This algorithm would significantly reduce the signaling overhead. The overhead
can be further decreased by removing the reference signal transmission. Then,
the IBSs have to rely on the sensed aggregate interference level during the
data transmission stage.

####  III-E2 Sequential Greedy Selection

In this algorithm, there is a setup stage before data transmission during
which the active IBSs coordinate to select sub-bands. After the RS
transmission phase and based on a predefined order, the first IBS selects a
sub-band and broadcasts the selection to the other IBSs. Each subnetwork
stores a local copy of the sub-band selection sets $B_{k}$, denoting the set
of active IBSs that have selected the $k$-th sub-band, and updates them
according to the broadcast on sub-band selection. Then, the next IBS will
choose a sub-band based on the information it has received about the sub-band
selection and the ISR levels. The $n$-th IBS will choose the $k^{*}$ based on
the following criteria:

$k^{*}=\underset{k=1,...,K}{\text{arg}\,\text{min}\,\,}\sum_{l\in\mathcal{B}_{k}}\frac{q_{nl}^{k}}{p_{n}^{k}},$
(13)

This process is carried out $M$ times which is set beforehand. One important
note here is that since IBSs take turns selecting a sub-band and informing
others, the algorithm might not be practically feasible to be carried out in a
limited preamble time before data transmission, especially in dense scenarios.

Finally, we also evaluate the performance of random sub-band selection as the
worst-case baseline for the sub-band allocation algorithms. With this scheme,
the IBSs randomly select one of the sub-bands to transmit their data.

##  IV Performance Evaluation

###  IV-A Simulation Assumptions and Methodology

We have considered a $20\text{m}\times 20\text{m}$ environment where there are
$5$ attraction points as the parent points of the Thomas Cluster Process that
corresponds to $\lambda_{p}=12.5\times 10^{3}\,\text{points}/\text{km}^{2}$.
After deploying these parent points, the IBSs are dropped according to a
normal distribution with $\sigma_{n}=2\text{m}$. The total number of IBSs is
considered to be $\\{30,40,\ldots,100\\}$ which is equivalent to $\lambda$
from $7.5\times 10^{4}\text{IBS}/\text{km}^{2}$ to $25\times
10^{4}\text{IBS}/\text{km}^{2}$. The shape of the IBSs is assumed to be a
cylinder with a radius $R=0.25$m and a height of $H=1.9$m. To capture the
difference in the heights of people, the z-coordinate of the APs and XRDDs are
generated to be uniformly distributed between $0.9\text{m}-1.15\text{m}$ and
$1.65\text{m}-1.90\text{m}$, respectively.

Furthermore, a bandwidth of $B_{t}=1\text{GHz}$ with the carrier frequency of
$f_{c}=2.4\text{GHz}$ is considered which is divided into
$K=\\{5,10,15,\ldots,30\\}$ sub-bands. The TTI duration is assumed to be
$T_{f}=5\text{ms or }10\text{ms}$ to evaluate the effect of it on the
reliability. A summary of the simulation parameters is presented in Table I.
The large scale fading for intra-body channel is modeled as
$\beta_{nn}(dB)=8.6\log_{10}(d_{nn})+46.1+2\chi$ where the log-normal
shadowing factor is denoted by $\chi\sim\mathcal{N}(0,1)$ [16]. The 3-D
distance between the $n$-th XRDD and the corresponding AP is denoted by
$d_{nn}$ in meters. The large-scale fading for the inter-body channel is
assumed to follow the proposed model for D2D communications by 3GPP [17],

$\displaystyle\beta_{nl}(dB)=\begin{cases}38.8+16.9\log_{10}(d_{nl})+3\chi,&S_{nl}=LOS\\\
17.5+43.3\log_{10}(d_{nl})+4\chi,&S_{nl}=NLOS,\end{cases}$ (14)

for $n\neq l$. $d_{nl}$ denotes the 3-D distance of $n$-th XRDD and the $l$-th
AP in meters. The variable $S_{nl}\in\\{LOS,NLOS\\}$ determines the state of
the inter-body link. We assume that $S_{nl}$ follows the probability
distribution

$Pr\\{S_{nl}=LOS\\}=\exp(-\gamma d_{nl}^{2D}),$ (15)

where $\gamma=2\lambda R$ and $\lambda$ is the deployment density of IBSs in
the environment. Furthermore, $d_{nl}^{2D}$ is the 2-D distance of
transceivers in the xy-plane. This model is based on the consideration of the
human body (IBSs) with the shape of a cylinder as the blockages. The results
in the following section are calculated based on $5\times 10^{6}$ snapshots.

TABLE I: Values of Simulation Parameters

Parameter | Value
---|---
Deployment Area ($\text{m}^{2}$) | $20\text{m}\times 20\text{m}$
Number of IBSs | $\\{30,40,\ldots,100\\}$
Number of Clusters |  $5$, ($\lambda_{p}=12.5\times 10^{3}/\text{km}^{2}$)
Std. of IBSs Relative Location, $\sigma_{n}$ (m) | $2$
Minimum Distance of Clusters and IBSs (m) |  $d_{c}^{min}=4$, $d_{IBS}^{min}=0.5$
IBS Cylinder Dimensions (m) |  $R=0.25$, $H=1.9$
Total Bandwidth, $B_{t}$ (GHz) | $1$
Carrier Frequency, $f_{c}$ (GHz) | $2.4$
Number of Sub-bands, $K$ | $\\{5,10,15,\ldots,30\\}$
Transmission Frame Duration, $T_{f}$ (ms) | $5\text{ or }10$
Percentage of Downlink Transmission, $\eta$ | 0.8
Average Data Rate, $R_{\nu}$ (Mbps) | $45$
Packet Inter-arrival Time, $T_{\nu}$ (ms) | $16.67$
Transmit Power (dBm) | $0$
Noise Figure (dB) | 10

###  IV-B Simulation Results

Figure 3: CDF of SINR, for $N=60$, $K=10$, $T_{f}=5,10$ms. Figure 4: Packet
Error Rate (PER) versus the Number of IBSs ($N$), for $K=10$ and
$T_{f}=5\text{ms},10$ms. Figure 5: Packet Error Rate (PER) versus the Number
of Sub-bands ($K$), for $N=60$ and $T_{f}=5$ms.

Figure 3 depicts the cumulative distribution function (CDF) of SINR for the
sub-band allocation algorithms described in Section III for $T_{f}=5\text{ms
and }10\text{ms}$. As seen in this figure, the worst performing sub-band
allocation algorithm, Greedy, shows more than $10$dB gain compared to random
selection in the tail of the CDF. SISA has at least $10\text{dB}$ higher SINR
in the tail compared to the other algorithms. SG outperforms the SISA in the
high percentiles, but, it fails in the tail which is likely due to the impact
of the mutual ISR that is not considered. In addition, it is shown that lower
TTI duration yields higher SINR due to a lower number of active IBSs per TTI
which is evident from the SINR increase in $T_{f}=5$ms compared to
$T_{f}=10$ms.

Furthermore, Figure 4 shows the packet error rate (PER) versus the number of
IBSs. The PER is evaluated for two different $T_{f}$ values. The figure shows
that SISA and SG are the best-performing sub-band allocation algorithms due to
their superiority in mitigating interference. By considering the baseline PER
of $10^{-2}$ as suggested in [3], SISA and SG can support more than around
$N=70$ IBSs for $T_{f}$ is $5$ms. This is around 75% more than CGC. Second,
although higher $T_{f}$ decreases the threshold on SINR according to (6), it
also reduces the SINR by increasing the number of active IBSs per TTI, as seen
in Figure 3. The effect of a larger number of active IBSs overcomes the
advantage of having a lower SINR threshold which is seen in Figure 4 with
higher PER for $T_{f}=10$ms compared to $T_{f}=5$ms.

Figure 5 illustrates the impact of $K$, number of sub-bands, on PER for $N=60$
and $T_{f}=5\text{ms}$. It is shown that for all sub-band allocation
algorithms, PER decreases at first by increasing the $K$ which stems from the
higher degree of interference diversity, due to the higher number of sub-
bands. However, after a certain number of sub-bands, e.g., $K=15$ for SISA,
the PER begins to increase with the $K$. This is an indication that having
more narrower sub-bands is no longer beneficial. This is reasonable since a
higher number of sub-bands with fixed total bandwidth means smaller sub-band
bandwidth for data transmission and therefore, putting a more strict
constraint for successful packet transmission based on (6). This effect is
evident in a rapid increase of PER in Figure 5 for all sub-band allocation
algorithms.

###  IV-C Signaling Overhead of the sub-band allocation algorithms

Apart from the performance in mitigating interference, another aspect that
differentiates the sub-band allocation algorithms is the sensing and signaling
overhead that they incur. There are two sources for the signaling overhead of
the sub-band allocation algorithms: 1) transmission of the reference signals,
and 2) feedback of the measurements. Transmitting the reference signals has an
overhead in the order of $\mathcal{O}(NK)$ since the systems need $N$
orthogonal sequences that IBSs have to transmit over $K$ sub-band which is the
same for all the algorithms. CGC depends on the ISR levels reported by the
sub-networks to create the graph and assign the colors. Each IBS has to report
a vector of size $N$ to the centralized controller. Thus, the signaling
overhead of CGC is in order of $\mathcal{O}(N^{2})$. Since the number of sub-
bands, $K$ is typically much less than the number of IBSs, the overhead
incurred by CGC is dominated by $\mathcal{O}(N^{2})$.

For SISA, each IBS must report individual ISR levels on all the $K$ sub-bands
translating to signalling of a matrix of size $N\times K$. Overall, the
asymptotic signaling overhead of the SISA will be in the order of
$\mathcal{O}(N^{2})$. For the centralized algorithms, there is also an
overhead of reporting the allocation to the IBSs which is in order of $N$ and
is negligible compared to other parts of signaling overhead. Furthermore, the
greedy algorithm only requires the devices to report the sensed aggregate ISR
levels over sub-bands to their APs which results in the overhead in order of
$\mathcal{O}(N)$. The SG algorithm has a signaling overhead in the order of
$\mathcal{O}(NM)$ as each IBS must broadcast its sub-band selection $M$ times
to the whole network. In conclusion, SISA and CGC have higher overhead than SG
and greedy which is caused majorly by the report of measurements to the
central controller.

##  V Conclusions

We considered the problem of sub-band allocation for IBS supporting the XR
downlink video frames stream. We evaluated the performance of state-of-the-art
sub-band allocation methods and their capability to mitigate interference and
showed that these algorithms could help improve IBS network performance
compared to non-interference-aware algorithms. It is illustrated that the SISA
and SG algorithms have the best performance in terms of PER; however, SISA
incurs higher signaling overhead to the network which results in poor
scalability. In addition, the results show that a higher number of sub-bands
and larger TTI duration will not necessarily result in a better packet success
rate. Finally, we can conclude that although IBS can act as a platform for XR
applications, further improvements are needed in resource allocation
algorithms for supporting ultra-dense deployments.

## References

  * [1] W. Saad, M. Bennis, and M. Chen, “A Vision of 6G Wireless Systems: Applications, Trends, Technologies, and Open Research Problems,” IEEE Network, vol. 34, no. 3, pp. 134–142, 2020.
  * [2] M. Gapeyenko, V. Petrov, S. Paris, A. Marcano, and K. I. Pedersen, “Standardization of Extended Reality (XR) over 5G and 5G-Advanced 3GPP New Radio,” IEEE Network, vol. 37, no. 4, pp. 22–28, 2023.
  * [3] 3GPP, “Extended Reality (XR) in 5G (Release 18),” Tech. Rep. 26.928, 3rd Generation Partnership Project (3GPP), 2023.
  * [4] M. A. Uusitalo, P. Rugeland, M. R. Boldi, E. C. Strinati, P. Demestichas, M. Ericson, G. P. Fettweis, M. C. Filippou, A. Gati, M.-H. Hamon, et al., “6G Vision, Value, Use Cases and Technologies from European 6G Flagship Project Hexa-X,” IEEE Access, vol. 9, pp. 160004–160020, 2021.
  * [5] G. Berardinelli, P. Baracca, R. O. Adeogun, S. R. Khosravirad, F. Schaich, K. Upadhya, D. Li, T. Tao, H. Viswanathan, and P. Mogensen, “Extreme Communication in 6G: Vision and Challenges for ‘in-X’ Subnetworks,” IEEE Open Journal of the Communications Society, vol. 2, pp. 2516–2535, 2021.
  * [6] S. Movassaghi, M. Abolhasan, J. Lipman, D. Smith, and A. Jamalipour, “Wireless Body Area Networks: A Survey,” IEEE Communications Surveys and Tutorials, vol. 16, no. 3, pp. 1658–1686, 2014.
  * [7] R. Adeogun, G. Berardinelli, and P. E. Mogensen, “Enhanced Interference Management for 6G in-X Subnetworks,” IEEE Access, vol. 10, pp. 45784–45798, 2022.
  * [8] R. Adeogun, G. Berardinelli, and P. Mogensen, “Learning to Dynamically Allocate Radio Resources in Mobile 6G in-X Subnetworks,” in IEEE PIMRC, pp. 959–965, 2021.
  * [9] D. Li, S. R. Khosravirad, T. Tao, and P. Baracca, “Advanced Frequency Resource Allocation for Industrial Wireless Control in 6G Subnetworks,” in IEEE WCNC, pp. 1–6, 2023.
  * [10] R. Adeogun and G. Berardinelli, “Multi-Agent Dynamic Resource Allocation in 6G in-X Subnetworks with Limited Sensing Information,” Sensors, vol. 22, no. 13, 2022.
  * [11] R. Adeogun and G. Berardinelli, “Distributed Channel Allocation for Mobile 6G Subnetworks via Multi-Agent Deep Q-Learning,” in IEEE WCNC, pp. 1–6, 2023.
  * [12] X. Du, T. Wang, Q. Feng, C. Ye, T. Tao, L. Wang, Y. Shi, and M. Chen, “Multi-agent Reinforcement Learning for Dynamic Resource Management in 6G in-X Subnetworks,” IEEE Transactions on Wireless Communications, 2022.
  * [13] G. Berardinelli and R. Adeogun, “Hybrid Radio Resource Management for 6G Subnetwork Crowds,” IEEE Communications Magazine, vol. 61, no. 6, pp. 148–154, 2023.
  * [14] M. Haenggi, Stochastic Geometry for Wireless Networks.  Cambridge University Press, 2012.
  * [15] 3GPP, “Study on XR (Extended Reality) Evaluations for NR (Release 17),” Tech. Rep. 38.838, 3rd Generation Partnership Project (3GPP), 2021.
  * [16] K. Takizawa, T. Aoyagi, J. ichi Takada, N. Katayama, K. Yekeh, Y. Takehiko, and K. R. Kohno, “Channel Models for Wireless Body Area Networks,” in 30th Annual International Conf. of the IEEE Engineering in Medicine and Biology Society, pp. 1549–1552, 2008.
  * [17] 3GPP, “Study on LTE Device to Device Proximity Services,” Tech. Rep. 36.843, 3rd Generation Partnership Project (3GPP), 2014.