# Mapping Urban Population Growth from Sentinel-2 MSI and Census Data Using
Deep Learning: A Case Study in Kigali, Rwanda
††thanks: This work was supported in part by the Swedish National Space Agency
under Grant dnr 155/15, the Digital Futures under the grant for the EO-
AI4GlobalChange Project, and in part by ESA-China Dragon 5 Program under the
EO-AI4Urban Project.

Sebastian Hafner  KTH Royal Inst. of Technology
Stockholm, Sweden
<EMAIL_ADDRESS>Stefanos Georganos  KTH Royal Inst. of Technology
Stockholm, Sweden
<EMAIL_ADDRESS>Theodomir Mugiraneza  University of Rwanda
Kigali, Rwanda
<EMAIL_ADDRESS>Yifang Ban  KTH Royal Inst. of Technology
Stockholm, Sweden
<EMAIL_ADDRESS>
###### Abstract

To better understand current trends of urban population growth in Sub-Saharan
Africa, high-quality spatiotemporal population estimates are necessary. While
the joint use of remote sensing and deep learning has achieved promising
results for population distribution estimation, most of the current work
focuses on fine-scale spatial predictions derived from single date census,
thereby neglecting temporal analyses. In this work, we focus on evaluating how
deep learning change detection techniques can unravel temporal population
dynamics at short intervals. Since Post-Classification Comparison (PCC)
methods for change detection are known to propagate the error of the
individual maps, we propose an end-to-end population growth mapping method.
Specifically, a ResNet encoder, pretrained on a population mapping task with
Sentinel-2 MSI data, was incorporated into a Siamese network. The Siamese
network was trained at the census level to accurately predict population
change. The effectiveness of the proposed method is demonstrated in Kigali,
Rwanda, for the time period 2016–2020, using bi-temporal Sentinel-2 data.
Compared to PCC, the Siamese network greatly reduced errors in population
change predictions at the census level. These results show promise for future
remote sensing-based population growth mapping endeavors. Code is available on
GitHub111https://github.com/SebastianHafner/PopulationGrowthMapping_Kigali.git.

###### Index Terms:

Population mapping, Sub-Saharan Africa, Siamese network

##  I Introduction

The projections in the World Population Prospects 2022 report suggest that the
global population could reach 9.7 billion in 2050 [1]. At the forefront of the
anticipated population growth are countries of Sub-Saharan Africa. In light of
this, frequent updates of existing population data in that region are crucial,
particularly considering that knowledge of population distribution is a
necessary requisite for a wide range of applications. For example, population
distribution maps provide vital information for vaccination campaigns,
disaster response deployment, and urban mobility and transport planning.

In recent years, census-independent (i.e., bottom-up) population mapping using
deep learning and satellite imagery has shown promise in providing accurate
population estimates. For example, Doupe et al. [2] mapped population density
at 8 km spatial resolution in Tanzania and Kenya using a Convolutional Neural
Network (CNN) based on the VGG architecture and Landsat 7 imagery. Landsat 7
imagery and the VGG-net were also used by Robinson et al. [3] to predict
population counts in the United States at 1 km spatial resolution. Authors in
[4] proposed to fuse Landsat 8 optical data with Sentinel-1 radar data to
predict population density at 4.5 km spatial resolution for rural villages in
India and demonstrated that dual-branch fusion networks outperform uni-modal
networks. Sentinel-2 (S2) MultiSpectral Instrument (MSI) imagery was used by
Huang et al. [5] to map population distribution at 1 km spatial resolution for
the Atlanta, Georgia, and Dallas, Texas metropolitan areas in the United
States of America. Recently, Neal et al. [6] used WorldView-2 imagery for
estimating population in two districts of Mozambique using representation
learning. A ResNet was also used in [7] to map population in Sub-Saharan
African cities with multisource satellite imagery from Pleiades and S2.
Building footprints were further used to improve the geographical
transferability of models.

While deep learning-based population mapping from satellite imagery has gained
traction in recent years [6, 7, 5, 4, 3, 2], little attention has been paid to
population growth mapping with the exception of [8]. Using a ResNet and
Landsat 5 imagery, Zhuang et al. [8] performed population growth analysis in
China for the 1985-2010 period by mapping population distribution at 1 km
spatial resolution with a 5-year interval. However, analyzing population
growth by comparative analysis of independently produced population maps,
i.e., change detection by Post-Classification Comparison (PCC), is well-known
to suffer from the error propagation of the individual population maps. To
that end, we propose an end-to-end population growth mapping method to
overcome the error propagation of PCC in uni-temporal population maps. This
study is, up to the best of our knowledge, the first to map population growth
in an end-to-end fashion from satellite imagery.

##  II Study Area and Data

Kigali, the capital city and economic hub of Rwanda, was selected as the study
area. Kigali encompasses an area of approximately 730 km2. In 2012, Kigali had
a population of approximately 1.1 million and placed among the fastest-growing
cities in Africa [9]. In recent years, rapid urbanization resulted in the
conversion of major cropland areas into built-up areas in the urban fringe
zones of Kigali, which increased ecosystem service demands and negatively
affected the habitat for biodiversity service function [10].

S2 MSI imagery of Kigali for 2016 and 2020 was retrieved from Google Earth
Engine [11]. Specifically, cloud-free composites were generated by collecting
all S2 Level-1C (top-of-atmosphere) scenes acquired during the wet season of
the respective year. Thereafter, cloudy pixels (i.e., cloud probability $>$ 50
%) were masked for each scene, before the scenes were combined using median
compositing. The resulting cloud-free composites for 2016 and 2020 are
visualized in Figure 1(a) and Figure 1(b), respectively.

(a)

(b)

(c)

(d)

Figure 1: S2 MSI composites for (a) 2016 and (b) 2020, and (c) population
labels at grid level. (d) shows the data set splits.

Population census data at the level of designated census enumeration areas
were acquired for Kigali for the years 2016 and 2020 (161 administrative
polygons). These areas are corresponding to the smallest administrative
entities in Rwanda called villages. The data consist of number of population
(head counts) and were acquired from two institutions including Kigali city
One Stop center and the Local Administrative Entities Development Agency.
Using iterative merging, we aggregated the dataset into a smaller number of
units, to reflect a more realistic scenario regarding data availability, but
also to adapt to the needs of the experiment (i.e., 100 meter predictive
spatial resolution). Finally, the census units were randomly split into a
training, validation, and test set (60/20/20 split) (Figure 1(d)).

##  III Methodology

###  III-A Problem Setup

We consider two S2 MSI images that cover the same geographical area (Kigali)
but were acquired at two different times, $t_{1}$ and $t_{2}$. Furthermore, we
consider the census units constituting the City of Kigali, where each census
unit, $U$, contains an accurate count of the population, $Y$, for $t_{1}$ and
$t_{2}$. The goal is to train a network that accurately predicts the
population growth for a census unit $D$ ($=Y^{\rm t_{2}}-Y^{\rm t_{1}}$) from
the part of the S2 images $I^{\rm t_{1}}$ and $I^{\rm t_{2}}$ covering $U$.
However, each census unit has a unique non-rectangular shape and, therefore,
cannot be used directly as network input. A common way to deal with this is to
operate on a grid level by dividing the entire study area into patches that
constitute the census areas [7]. Consequently, census units are composed of a
varying number of patches (100 x 100 m). Therefore, the network input to
predict the population growth for a census unit is, in practice, the
collection of S2 patches, $x^{t1}$ and $x^{t2}$, constituting the census unit.

###  III-B Proposed Method

The proposed population growth mapping method consists of two stages: 1) an
encoder model is pretrained by mapping population at the grid level, and 2) a
Siamese network, incorporating the pretrained encoder, is trained at the
census level to map population growth.

#### Population Mapping at Grid Level

Our previous work demonstrated that an encoder based on the ResNet-18
architecture suffices to learn salient features from S2 MSI imagery for
population mapping [7]. The same architecture is employed in this work (Figure
2). Specifically, the first layer of the ResNet-18 encoder is replaced with a
3 x 3 conv layer with 4 input channels to accommodate the 10 m S2 bands (Band
2, 3, 4, and 8) as input, while the remaining conv blocks constituting the
encoder remain unchanged. The features extracted with the encoder are
converted to a population prediction, $p$, using a fully connected layer.
Finally, the ReLu activation function is used to constrain $p$ values to
positive numbers.

Hyper-parameters for training are tuned on the validation set using grid
search with 3 learning rates ($10^{-5}$, $10^{-4}$, and $10^{-3}$) and 2 batch
sizes (8, 16). AdamW is used as optimizer, and the training duration is set to
100 epochs with early stopping (patience 5) to prevent models from overfitting
to the training set. As in [7], flips (horizontal and vertical) and rotations
($k*90^{\circ}$, where $k\in\\{0,1,2,3\\}$) are applied to the training data
for data augmentation, and the Mean Square Error (MSE) loss (commonly known as
L2 loss) is used as loss function. L2 loss is defined as follows:
$L2=(y-p)^{2}$, where the true and predicted population value is denoted by
$y$ and $p$, respectively. An NVIDIA GeForce RTX 3090 graphics card is used
for training.

Figure 2: Diagram of the ResNet-18 model used for grid-level population
mapping.

#### Population Growth Mapping at Census Level

For population growth mapping, we incorporate the pretrained ResNet-18 encoder
into a Siamese network (Figure 3). Siamese networks consist of two encoders
with shared weights that are used to separately extract features from the
inputs, before deriving the change information from the combined features. Due
to their inherent suitability to detect differences, Siamese networks have
also become a popular architecture for change detection in bi-temporal pairs
of satellite images. In this work, the pretrained encoder is employed to
extract features on population count from both images separately. The pair of
bi-temporal features is then converted to a population growth prediction using
a fully connected layer. No activation function is applied to the output of
that layer to allow for negative growth predictions.

An important challenge of supervised population growth mapping is that bi-
temporal population counts are required for the derivation of growth labels.
While it is possible to accurately disaggregate a census to a grid, this
requires auxiliary data such as land cover maps or building footprints.
However, this data is often not available for both timestamps. Therefore, the
Siamese network is trained at the census level by adapting the weakly
supervised learning strategy proposed in [12]. Specifically, Metzger et al.
[12] trained a population mapping model using population count at the census
level as labels by comparing them to the aggregated model predictions (patch-
level) for corresponding census units. Likewise, we use the Siamese network to
predict population growth separately for all patches of a census unit, before
applying the loss to the sum of predicted growth, $D$, using $\Delta Y$ as
label. The training setup (i.e., hyper-parameter tuning, early stopping, and
data augmentations) is identical to that for population mapping. It should be
noted, however, that the pretrained encoder is frozen during training, meaning
that only the fully connected layer ($f_{\rm fc}$ in Figure 3) is trained.

Figure 3: Diagram of the proposed population growth mapping method consisting
of two pretrained ResNet-18 encoders, $f_{\rm en}$, with shared weights and a
fully connected layer, $f_{\rm fc}$. The network is trained at the census
level with frozen encoders.

###  III-C Accuracy Metrics

We make use of three commonly employed metrics in population studies [13],
namely the Root Mean Squared Error (RMSE), the Mean Absolute Error (MAE), and
the coefficient of determination (R2). RMSE and MAE are defined as follows:

$\text{RMSE}=\sqrt{\frac{\sum^{n}_{i=1}(y_{i}-p_{i})^{2}}{n}},\
\text{MAE}=\sqrt{\frac{\sum^{n}_{i=1}\left|y_{i}-p_{i}\right|}{n}},$ (1)

where $y$ and $p$ are true and predicted values, respectively, and $n$ is the
sample size. On the other hand, R2 is defined as 1 minus the fraction of the
residual sum of squares and the total variability of the data.

##  IV Results

Table I lists the quantitative population mapping results at the grid level
for 2020 and at the census level for 2016 and 2020. All three accuracy metrics
indicate that accurate population predictions were achieved at the grid level.
However, the aggregated results at the census level provide a stronger
validation since the census population counts are official data. While RMSE
and MAE values are not comparable between the grid and census level, the R2
values at the census level indicate good performance (0.70 +), although worse
than the performance achieved at the grid level (0.84). It is also apparent
that the obtained accuracy values for 2016 and 2020 are relatively similar.
Consequently, applying the model to new data from a different year had little
impact on model performance.

TABLE I: Quantitative population mapping results at the gird and census level for the test set. Level |  RMSE $\downarrow$ |  |  MAE $\downarrow$ |  |  $\text{R}^{2}$ $\uparrow$
---|---|---|---|---|---
2016 | 2020 |  | 2016 | 2020 |  | 2016 | 2020
Grid | - | 19 |  | - | 10 |  | - | 0.84
Census | 3,199 | 3,253 |  | 2,368 | 2,196 |  | 0.72 | 0.73

Figure 4 quantitatively compares the population growth predictions of (a) the
PCC with (b) the proposed end-to-end method. The former, PCC, performed
poorly, resulting in very high errors (RMSE = 1,471 and MAE = 1,082). In
contrast, the proposed method achieved satisfactory results with an RMSE of
202 and an MAE of 165. In terms of $\text{R}^{2}$ values, the results are more
similar, but better performance was also achieved by the proposed method (0.55
vs. 0.67). However, it is also apparent that the proposed method generally
underestimates population growth.

The qualitative population growth mapping predictions of the proposed method
are visualized in Figure 5(b), next to the ground truth in Figure 5(a).
Although the magnitude of growth was underestimated, the proposed method
picked up on the population growth that occurred on the outskirts of Kigali
(e.g., in the northeast and in the central south). However, the model failed
to detect population growth in the small census units of central Kigali for
which it predicted slightly negative growth values.

(a)

(b)

Figure 4: Population growth test results at the census level for (a) PCC and
(b) the proposed method.

(a)

(b)

Figure 5: Population growth maps for the test census units. (a) shows the
ground truth and (b) our predictions.

##  V Discussion and Limitations

We find the proposed method to be effective for population growth mapping from
S2 MSI imagery, especially compared to PCC. Our findings also emphasize that
salient features about population count can be learned from S2 imagery using a
ResNet model. These results are in line with [7].

Our work is also subject to several limitations. First of all, to train the
Siamese network, bi-temporal census data is required. However, census data,
let alone bi-temporal census data, is difficult to obtain in Sub-Saharan
Africa, or often not available at all [13]. Moreover, the S2 mission was
launched less than 8 years ago, while censuses are typically conducted every
10 years. Consequently, bi-temporal census data for time periods starting
after 2015 are largely unavailable. Another limitation of this work is that
population predictions are based on the presence of built-up areas, but the
land use of these areas may not be residential [3]. To overcome this, Neal et
al. [6] suggest including additional data modalities like, for example, night-
time light data. Our quantitative results in central Kigali (Figure 5(b)) also
suggest that densification of urban areas, and the subsequent increase in
population, may be challenging to accurately predict. Finally, further work is
needed to assess if the proposed method can accurately detect negative
population growth as a result of, for example, slum evictions.

##  VI Conclusion

In this paper, a population growth mapping method based on a Siamese network
is proposed and evaluated in Kigali, Rwanda for the time period 2016–2020.
Using S2 MSI data as input, the proposed method achieved satisfactory
population growth mapping results at the census level (RMSE = 202, MAE = 165,
$\text{R}^{2}$ = 0.67), and greatly outperformed PCC in terms of RMSE (-1,269)
and MAE (-917).

Our future work will extend the study area to other Sub-Saharan African
cities. Furthermore, we will investigate semi-supervised learning for Siamese
network training (e.g., [14]) to reduce the dependence on bi-temporal census
data.

## References

  * [1] United Nations Department of Economic and Social Affairs, Population Division,  “World population prospects 2022: Summary of results,”  Tech. Rep. UN DESA/POP/2022/TR/NO. 3, 2022.
  * [2] Doupe, P. et al,  “Equitable development through deep learning: The case of sub-national population density estimation,”  in Proceedings of the 7th Annual Symposium on Computing for Development, 2016, pp. 1–10.
  * [3] Robinson, C. et al,  “A deep learning approach for population estimation from satellite imagery,”  in Proceedings of the 1st ACM SIGSPATIAL Workshop on Geospatial Humanities, 2017, pp. 47–54.
  * [4] Hu, W. et al,  “Mapping missing population in rural india: A deep learning approach with satellite imagery,”  in Proceedings of the 2019 AAAI/ACM Conference on AI, Ethics, and Society, 2019, pp. 353–359.
  * [5] Huang, X. et al,  “Sensing population distribution from satellite imagery via deep learning: Model selection, neighboring effects, and systematic biases,”  IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 14, pp. 5137–5151, 2021.
  * [6] Neal, I. et al,  “Census-independent population estimation using representation learning,”  Scientific Reports, vol. 12, no. 1, pp. 1–12, 2022.
  * [7] Georganos, S. et al,  “A census from heaven: Unraveling the potential of deep learning and earth observation for intra-urban population mapping in data scarce environments,”  International Journal of Applied Earth Observation and Geoinformation, vol. 114, pp. 103013, 2022.
  * [8] Zhuang, H. et al,  “Mapping multi-temporal population distribution in china from 1985 to 2010 using landsat images via deep learning,”  Remote Sensing, vol. 13, no. 17, pp. 3533, 2021.
  * [9] National Institute of Statistics of Rwanda,  “The 2012 population and housing census results,”  Tech. Rep., 2012.
  * [10] Mugiraneza, T. et al,  “Monitoring urbanization and environmental impact in kigali, rwanda using sentinel-2 msi data and ecosystem service bundles,”  International Journal of Applied Earth Observation and Geoinformation, vol. 109, pp. 102775, 2022.
  * [11] Gorelick, N. et al,  “Google earth engine: Planetary-scale geospatial analysis for everyone,”  Remote sensing of Environment, vol. 202, pp. 18–27, 2017.
  * [12] Metzger, N. et al,  “Fine-grained population mapping from coarse census counts and open geodata,”  Scientific Reports, vol. 12, no. 1, pp. 20085, Nov 2022.
  * [13] Linard, C. et al,  “Population distribution, settlement patterns and accessibility across africa in 2010,”  PloS one, vol. 7, no. 2, pp. e31743, 2012.
  * [14] Hafner, S. et al,  “Urban change detection using a dual-task siamese network and semi-supervised learning,”  in IGARSS 2022-2022 IEEE International Geoscience and Remote Sensing Symposium. IEEE, 2022, pp. 1071–1074.


# Energy preserving boundary conditions in field theory

Manuel Asorey and Fernando Ezquerro  Centro de Astropartículas y Física de
Altas Energías, Departamento de Física Teórica. Universidad de Zaragoza,
E-50009 Zaragoza, Spain

###### Abstract

The dynamics of classical field theories is usually governed by field
equations, but when fields are constrained to bounded domains it is also
dependent on its boundary conditions. Usually boundary conditions are
constrained by the requirement of preserving the maximal symmetry of the
system. In the case of charged particles the symmetry is $U(1)$, but there are
many fields (e,g, electromagnetic fields) which are neutral and charge
conservation does not constraint its boundary conditions. In this paper we
explore the most general boundary conditions that preserve another symmetry
that all relativistic field theories do preserve: space-time translations. In
particular the families of boundary conditions of isolated systems which
preserve energy for scalar, electromagnetic and Yang-Mills field theories. We
point out the global properties of the space of all possible boundary
conditions of confined systems in two special domains. We also explore the
connection between energy preserving and charge preserving boundary
conditions.

##  1 Introduction

The new role of quantum boundary effects is boosting in condensed matter a new
era of quantum technologies. Indeed, the presence of plasmons and other
surface effects in metals and dielectrics [1], the appearance of edge currents
in the Hall effect [2]–[4], and the discovery of new edge effects in
topological insulators [5]–[10] and Weyl semiconductors [11] have very rich
potential implications.

Although boundary effects arise in quantum physics since the early days of the
theory, the role of boundary effects in quantum field theory have a much later
development.

Boundary effects also appear today as an essential ingredient in fundamental
physics. Since the discovery of the Casimir effect [12] they also arise behind
new quantum effects like Hawking radiation, Black Hole horizons effects,
topological defects, topology change [13]–[16] and holography in the AdS/CFT
correspondence.

The increasing relevant role of boundary effects is demanding a comprehensive
theory of boundary conditions. In spite of the fact that quite a lot of work
has been devoted to establish the foundations of the quantum theory, a
comprehensive theory of boundary conditions for quantum field theories is
still lacking. This gap was filled in the non-relativistic case by Asorey,
Ibort and Marmo [14] by using the unitarity principle for time evolution which
was further extended to quantum field theory [17, 18] . This first global
analysis was based on the preservation of unitarity for time evolution.
Another essential remark was made by Asorey, García-Alvarez and Muñoz-
Castañeda [19] who pointed out that the generalization for relativistic field
theories requires a change of the basic principles, from unitarity to the
preservation of the $U(1)$ symmetry which is responsible of electric charge
conservation. However, this principle does not apply to neutral or gauge
fields where the new approach does not provide any fundamental law to be
preserved.

In this paper we address the analysis of the theory of boundary conditions in
field theories confined in isolated domains based only on the requirement of
conservation of energy. This method applies to any bosonic or fermionic
relativistic field theory including neutral fields like gauge fields, with the
only exception of gravitation and topological field theories.

In section 2 we start our analysis with the neutral scalar fields in the case
of one boundary wall. In section 3 we study the case of charged scalar fields.
In section 4 we complete the discussion about charged scalar fields showing
the compatibility between the condition of charge conservation and energy
conservation. In section 5 we generalize this approach to Yang-Mills theories
developing the theory of boundary conditions which preserve energy and gauge
invariance. In section 6 we analyze the case of two different electromagnetic
active media separated by a planar surface. Finally, in section 7 we summarize
the applications and conclusions of this work.

##  2 Neutral scalar fields in half space

The dynamics of a real scalar field is governed by the Lagrangian density

$\mathcal{L}=\frac{1}{2}\partial_{\mu}\phi\partial^{\mu}\phi-\frac{1}{2}m^{2}\phi^{2}-V(\phi),$
(1)

where $V(\phi)$ is any arbitrary local potential function. The space-time
translation symmetry induces by Noether theorem four conservation laws

$\partial^{\mu}T_{\mu\nu}=0,$ (2)

where the energy-momentum tensor is given by

$T_{\mu\nu}=\frac{1}{2}\partial_{\mu}\phi\partial_{v}\phi-\frac{1}{2}\eta_{\mu\nu}\partial_{\alpha}\phi\partial^{\alpha}\phi+\frac{1}{2}m^{2}\phi^{2}\eta_{\mu\nu}+V(\phi)\eta_{\mu\nu},$
(3)

$\eta_{\mu\nu}$being the Minkowski metric. In particular when $\nu=0$ we get
the energy conservation law

$\partial^{\mu}T_{\mu 0}=\partial_{t}{\mathcal{E}}+\partial^{i}P_{i}=0,$

where

${\mathcal{E}}=\frac{1}{2}\partial_{t}\phi\partial_{t}\phi+\frac{1}{2}\partial_{i}\phi\partial_{i}\phi+\frac{1}{2}m^{2}\phi^{2}+V(\phi)$

is the energy density and

$P_{i}=\frac{1}{2}\partial_{0}\phi\,\partial_{i}\phi,$

the momentum density of the field (we assume that $c=1$ from now on). Thus,
for any bounded domain $\Omega$ with regular boundary $\partial\Omega$ the
variation of the intrinsic field energy is

$\frac{d}{dt}{E}_{\Omega}=\int_{\Omega}\partial_{t}{\mathcal{E}}\
d^{3}x=-\int_{\Omega}\partial^{i}P_{i}\
d^{3}x=-\int_{\partial\Omega}n^{i}P_{i}\ d^{2}x,$

where ${\mathbf{n}}=(n^{i})$ denotes the normal vector to the boundary surface
$\partial\Omega$.

Let us consider a simple case where $\Omega$ is just a half-space
$\Omega=\\{{\mathbf{x}}=(x^{1},x^{2},x^{3})|x^{3}\geq 0\\}$ whose boundary
$\partial\Omega=\\{{\mathbf{x}}=(x^{1},x^{2},0)\\}$ is the plane perpendicular
to the vector ${\mathbf{n}}=(0,0,-1)\in{\mathbb{R}}^{3}$. In that case the
conservation of energy implies that

$\frac{d}{dt}{E}_{\Omega}=\int_{\partial\Omega}T_{03}\
dx^{1}dx^{2}=\int\dot{\varphi}\,\varphi^{\prime}\ dx^{1}dx^{2}=0,$

where

$\dot{\varphi}=\partial_{t}\phi\Big{|}_{\partial\Omega}$

and

$\varphi^{\prime}=\partial_{3}\phi\Big{|}_{\partial\Omega}$

are the boundary values of the time derivative and the normal derivative of
the fields $\phi$ across the boundary $\partial\Omega$111In dimensions higher
than 2 there are some technical difficulties concerning the regularity of
boundary values [18] . In this paper we will restrict ourselves to cases of
regular boundary conditions . If we consider only homogeneous boundary
conditions which are invariant under translations along the boundary plane
$\partial\Omega$, the conservation of energy requires that

$\dot{\varphi}\,\varphi^{\prime}=0$ (4)

which has two solutions $\varphi^{\prime}=0$ and $\dot{\varphi}=0$. The first
solution corresponds to Neumann boundary condition which in string theory is
the usual boundary condition for open strings, whereas the second solution
includes the Dirichlet boundary conditions which corresponds to D-branes in
that theory [20] .

##  3 Charged scalar fields

Let us now consider the case of complex scalar fields with Lagrangian density

$\mathcal{L}=\frac{1}{2}\partial_{\mu}\phi^{\ast}\partial^{\mu}\phi-\frac{1}{2}m^{2}|\phi|^{2}-V(|\phi|),$
(5)

where $V(\phi)$ is any arbitrary local density potential function. Using the
same arguments as in the case of real scalar fields we get that energy
conservation requires the vanishing of

$\dot{\varphi^{*}}\varphi^{\prime}+\varphi^{*\prime}\dot{\varphi}=0.$ (6)

After the change of variables

$\psi_{1}=\begin{pmatrix}\varphi^{\prime}+\dot{\varphi}\\\
\varphi^{*\prime}+\dot{\varphi^{*}}\end{pmatrix}\hskip
28.45274pt\psi_{2}=\begin{pmatrix}\varphi^{\prime}-\dot{\varphi}\\\
\varphi^{*\prime}-\dot{\varphi^{*}}\end{pmatrix},$

the vanishing condition becomes

$\lvert\psi_{1}\rvert^{2}-\lvert\psi_{2}\rvert^{2}=4\left(\varphi^{*\prime}\dot{\varphi}+\dot{\varphi^{*}}\varphi^{\prime}\right).$

If we restrict ourselves to boundary conditions which are translation
invariant along the boundary the most general solution satisfies

$\left|\begin{pmatrix}\varphi^{\prime}+\dot{\varphi}\\\
\varphi^{*\prime}+\dot{\varphi^{*}}\end{pmatrix}\right|^{2}=\left|\begin{pmatrix}\varphi^{\prime}-\dot{\varphi}\\\
\varphi^{*\prime}-\dot{\varphi^{*}}\end{pmatrix}\right|^{2}$ (7)

and is given by

$\begin{pmatrix}\varphi^{\prime}+\dot{\varphi}\\\
\varphi^{*\prime}+\dot{\varphi^{*}}\end{pmatrix}=U\begin{pmatrix}\varphi^{\prime}-\dot{\varphi}\\\
\varphi^{*\prime}-\dot{\varphi^{*}}\end{pmatrix},$ (8)

where $U$ is an arbitrary 2x2 unitary matrix. Conjugation of equation (8)
leads to

$\begin{pmatrix}\varphi^{*\prime}+\dot{\varphi^{*}}\\\
\varphi^{\prime}+\dot{\varphi}\end{pmatrix}=U^{*}\begin{pmatrix}\varphi^{*\prime}-\dot{\varphi^{*}}\\\
\varphi^{\prime}-\dot{\varphi}\end{pmatrix};\quad\sigma_{1}\psi_{1}=U^{*}\sigma_{1}\psi_{2},\quad\psi_{1}=\sigma_{1}U^{*}\sigma_{1}\psi_{2}.$

This implies that the unitary matrix $U$ has to satisfy an extra condition

$U=\sigma_{1}U^{*}\sigma_{1}.$ (9)

The meaning of this restriction is that the unitary matrix $U$ has to belong
also to the $O(1,1)$ rotation group, because from (9) it follows that $U$
leaves the $(1,1)$ metric

$\sigma_{1}=\begin{pmatrix}0&1\\\ 1&0\end{pmatrix}$

invariant. i.e.

$U^{\perp}\sigma_{1}U=\sigma_{1}.$

Thus, the set of homogeneous local boundary conditions is given by the one-
dimensional subgroup of unitary matrices $G=O(1,1)\cap U(2)$222This is one of
three maximal compact subgroups of $U(2)$[21] which has two disjoint
components $G=G_{+}\cup G_{-}$ whose matrices only differ by the sign of their
determinant $\det U=\pm 1$. The component $G_{+}$ contains all the matrices of
the form

$U_{+}(a)=e^{\,ia\sigma_{3}}\qquad a\in[0,2\pi),$ (10)

that are continuously connected with the identity, whereas the other component
$G_{-}$ is given by all the matrices

$U_{-}(b)=\sigma_{1}e^{\,ib\sigma_{3}}\qquad b\in[0,2\pi),$ (11)

which cannot be continuously connected with the identity. The general solution
of the first type (10) is

$\varphi^{\prime}=-i\cot\frac{a}{2}\dot{\varphi},$

whereas

$\displaystyle\operatorname{Re}(\dot{\varphi}+\varphi^{\prime})+\operatorname{Re}(\dot{\varphi}-\varphi^{\prime})\cos
b=\operatorname{Im}(\dot{\varphi}-\varphi^{\prime})\sin b$ .
$\displaystyle\operatorname{Im}(\dot{\varphi}+\varphi^{\prime})-\operatorname{Im}(\dot{\varphi}-\varphi^{\prime})\cos
b=\operatorname{Re}(\dot{\varphi}-\varphi^{\prime})\sin b$

is the general solution of the second type (11).

Let us consider some particular cases of physical interest.

  * i)

$U_{D}=\mathbb{I}$: Static boundary conditions. This is a boundary condition
of the first type with $a=0$

$\dot{\varphi}=0.$

  * ii)

$U_{N}=-\mathbb{I}$: Neumann boundary conditions. This is a boundary condition
of the first type with $a=\pi$

$\varphi^{\prime}=0.$

  * iii)

$U_{\mathrm{c}}=\pm i\sigma_{3}$: Chiral boundary conditions. These are
boundary conditions of the first type with $a=\pm\frac{\pi}{2}$

$\varphi^{\prime}=\mp i\dot{\varphi}.$

  * iv)

$U_{\mathrm{t}}=\pm\sigma_{1}$: Twisted boundary conditions. These are
boundary conditions of the second type with $b=0$ or $b={\pi}$

$\displaystyle\operatorname{Im}\varphi^{\prime}$
$\displaystyle=0\qquad\operatorname{Re}\dot{\varphi}=0$
$\displaystyle\operatorname{Re}\varphi^{\prime}$
$\displaystyle=0\qquad\operatorname{Im}\dot{\varphi}=0.$

  * iv)

$U_{\mathrm{tc}}=\pm\sigma_{2}$: Twisted chiral boundary conditions. These are
boundary conditions of the second type with
$\beta=\frac{\pi}{2},\frac{3\pi}{2}$

$\operatorname{Re}\varphi^{\prime}=\mp\operatorname{Im}\varphi^{\prime}\qquad\operatorname{Im}\dot{\varphi}=\pm\operatorname{Re}\dot{\varphi}.\\\
$

###  3.1 Charged scalar fields two parallel plates

Let us consider a complex scalar field confined between two parallel plates
$\Omega=\\{(x^{1},x^{2},x^{3})\Big{|}-L\leq x^{3}\leq L\\}$. In this case
homogeneous local boundary conditions which preserve energy must satisfy that

$\varphi_{1}^{*\prime}\dot{\varphi}_{1}+\dot{\varphi_{1}^{*}}\varphi_{1}^{\prime}-\varphi_{2}^{*\prime}\dot{\varphi}_{2}-\dot{\varphi_{2}^{*}}\varphi_{2}^{\prime}=0,$
(12)

where $\varphi_{1}(x^{1},x^{2})=\phi(x^{1},x^{2},-L)$ and
$\varphi_{2}(x^{1},x^{2})=\phi(x^{1},x^{2},L)$ . Written in terms of the
following vectors

$H_{1}=\begin{pmatrix}\dot{\varphi}_{1}+\varphi_{1}^{\prime}\\\
\dot{\varphi}_{2}-\varphi_{2}^{\prime}\\\
\dot{\varphi_{1}^{*}}+\varphi_{1}^{*\prime}\\\
\dot{\varphi_{2}^{*}}-\varphi_{2}^{*\prime}\end{pmatrix}\hskip
14.22636ptH_{2}=\begin{pmatrix}\dot{\varphi}_{1}-\varphi_{1}^{\prime}\\\
\dot{\varphi}_{2}+\varphi_{2}^{\prime}\\\
\dot{\varphi_{1}^{*}}-\varphi_{1}^{*\prime}\\\
\dot{\varphi_{2}^{*}}+\varphi_{2}^{*\prime}\end{pmatrix}$

the restriction (12) reads

$\lvert H_{1}\rvert^{2}-\lvert
H_{2}\rvert^{2}=4\left(\dot{\varphi^{*}_{1}}\varphi^{\prime}_{1}+\varphi_{1}^{*\prime}\dot{\varphi}_{1}-\dot{\varphi^{*}_{2}}\varphi^{\prime}_{2}-\varphi_{2}^{*\prime}\dot{\varphi}_{2}\right)=0.$

This means that any solution has to be of the form

$H_{1}=UH_{2}$ (13)

with $U$ an unitary matrix of $U(4)$. There is another requirement that this
solution must satisfy. Indeed, if we conjugate (13) we get

$H_{1}^{*}=U^{*}H_{2}^{*}$

that implies

$H_{1}=\begin{pmatrix}0&{\mathbb{I}}_{2}\\\
{\mathbb{I}}_{2}&0\end{pmatrix}U^{*}\begin{pmatrix}0&{\mathbb{I}}_{2}\\\
{\mathbb{I}}_{2}&0\end{pmatrix}H_{2},$

which implies a further restriction on the unitary matrix

$U=\begin{pmatrix}0&{\mathbb{I}}_{2}\\\
{\mathbb{I}}_{2}&0\end{pmatrix}U^{*}\begin{pmatrix}0&{\mathbb{I}}_{2}\\\
{\mathbb{I}}_{2}&0\end{pmatrix}.$ (14)

The meaning of this restriction is that the unitary matrix $U$ has to belong
also to the $O(2,2)$ rotation group, because from (14) it follows that $U$
leaves the $(2,2)$ metric

$\begin{pmatrix}0&{\mathbb{I}}_{2}\\\ {\mathbb{I}}_{2}&0\end{pmatrix}$

invariant. i.e.

$U^{\perp}\begin{pmatrix}0&{\mathbb{I}}_{2}\\\
{\mathbb{I}}_{2}&0\end{pmatrix}U=\begin{pmatrix}0&{\mathbb{I}}_{2}\\\
{\mathbb{I}}_{2}&0\end{pmatrix}.$

Thus, the general solution of the homogeneous local boundary conditions is
given by the six-dimensional subgroup of unitary matrices $G=O(2,2)\cap U(4)$
which has two disjoint components $G=O(2,2)\cap U(4)=G_{+}\cup G_{-}$
distinguished by the sign of the determinant $\det U=\pm 1$, i.e. the
component $G_{+}$ contains all the matrices of the form

$U_{+}(\mathbf{a},\mathbf{b},\mathbf{c})={\mathrm{exp}}\,i\begin{pmatrix}a_{1}&b_{1}+ib_{2}&0&c_{1}+ic_{2}\\\
b_{1}-ib_{2}&a_{2}&-c_{1}-ic_{2}&0\\\ 0&-c_{1}+ic_{2}&-a_{1}&-b_{1}+ib_{2}\\\
c_{1}-ic_{2}&0&-b_{1}-ib_{2}&-a_{2}\end{pmatrix}$ (15)

with $\mathbf{a},\mathbf{b},\mathbf{c}\in\mathbb{R}^{6}$, that are
continuously connected with the identity. The other component is given by the
matrices of the form

$U_{-}(\mathbf{a},\mathbf{b},\mathbf{c})=\frac{1}{2}\begin{pmatrix}1&1&-1&1\\\
1&1&1&-1\\\ -1&1&1&1\\\
1&-1&1&1\end{pmatrix}U_{0}(\mathbf{a},\mathbf{b},\mathbf{c})$ (16)

that are disconnected from $G_{+}$. In summary, the homotopy group of $G$ is
$\pi_{0}(G)={\mathbb{Z}}_{2}$.

This group contains the solutions of the type considered in the previous case
for each of the plane boundaries. But the fact that there are two boundaries
gives rise to other remarkable boundary conditions like periodic boundary
conditions.

Particular cases of interest are

  * i)

$U_{N}={\mathbb{I}}_{4}$: Neumann boundary conditions. i.e. Neumann boundary
conditions for both walls $\varphi_{1}^{\prime}=\varphi_{2}^{\prime}=0$.

  * ii)

$U_{D}=-{\mathbb{I}}_{4}$: Static boundary conditions. i.e. static boundary
conditions for both walls $\dot{\varphi}_{1}=\dot{\varphi}_{2}=0$.

  * iii)

$U_{p}=\begin{pmatrix}\sigma_{1}&0\\\ 0&\sigma_{1}\end{pmatrix}:$ Periodic
boundary conditions connecting the two walls
$\dot{\varphi}_{1}=\dot{\varphi}_{2}$,
$\varphi_{1}^{\prime}=\varphi_{2}^{\prime}$, and

  * iv)

$U_{\mathrm{ap}}=\begin{pmatrix}-\sigma_{1}&0\\\ 0&-\sigma_{1}\end{pmatrix}$:
Antiperiodic boundary conditions connecting the two walls in a different
manner $\dot{\varphi}_{1}=-\dot{\varphi}_{2}$,
$\varphi_{1}^{\prime}=-\varphi_{2}^{\prime}$.

##  4 Compatibility between energy and charge conservations

In the case of complex fields there is another conserved quantity, the
electric charge. The charge conservation law provides another condition to be
preserved by boundary conditions [19]. In principle the families of boundary
conditions which preserve charge and energy are different. However, in the
case of charged scalar fields both families of boundary conditions are
compatible. This very relevant property is a consequence of the compatibility
of gauge transformations and space-time traslations. Indeed, the actions of
the $U(1)$ group of gauge transformations

$G(\alpha)\phi=e^{i\alpha}\phi$ (17)

and the group of translations $T_{4}$

$T_{4}\phi(x)=\phi(x-a)$ (18)

do commute. Moreover, a consequence of that property is that the Poisson
bracket of the charge density 333 $\Pi$ and $\Pi^{\ast}$ are the canonical
momenta
$\Pi=\frac{\partial\mathcal{L}}{\partial\dot{\phi}}=\dot{\phi^{*}}\hskip
14.22636pt\Pi^{*}=\frac{\partial\mathcal{L}}{\partial\dot{\phi^{*}}}=\dot{\phi}$
(19)

$\rho=\frac{i}{2}\left(\phi^{*}\dot{\phi}-\phi\dot{\phi^{*}}\right)=\frac{i}{2}\left(\phi^{*}\Pi^{*}-\Pi\phi\right)$

and the energy density

$\displaystyle{\mathcal{E}}$
$\displaystyle=\frac{1}{2}\left(\dot{\phi^{*}}\dot{\phi}+\nabla\phi^{*}\nabla\phi\right)+\frac{1}{2}m^{2}\phi^{2}+V(\phi)$
$\displaystyle=\frac{1}{2}\left(\Pi\Pi^{*}+\nabla\phi^{*}\nabla\phi\right)+\frac{1}{2}m^{2}\phi^{2}+V(\phi)$
(20)

vanishes, i.e.

$\displaystyle\\{\rho,{\mathcal{E}}\\}$ $\displaystyle=0.$

The conservation of electric charge is given by the continuity equation

$\partial_{t}\rho+\partial^{i}j_{i}=0,$

where

${\mathbf{j}}=\frac{i}{2}\left(\phi^{*}{\boldsymbol{\nabla}}\phi-({\boldsymbol{\nabla}}\phi^{*})\phi\right).$

In other words, charge conservation requires the vanishing of the electric
current flux through the boundary [19]

$-\int_{\Omega}\dot{\rho}\ d^{3}x=\int_{\Omega}\partial^{i}j_{i}\
d^{3}x=\int_{\partial\Omega}{\mathbf{j}}\,d{\boldsymbol{\sigma}}=\frac{i}{2}\int_{\partial\Omega}\left(\varphi^{*}{\boldsymbol{\nabla}}\varphi-({\boldsymbol{\nabla}}\varphi^{*})\varphi\right)d{\boldsymbol{\sigma}}$

which in the right half space case reduces to

$\int_{\Omega}\dot{\rho}\
d^{3}x=\frac{i}{2}\int\left(\varphi^{*}\varphi^{\prime}-\varphi^{*\prime}\varphi\right)dx^{1}dx^{2}.$

Thus, homogeneous local boundary conditions must satisfy

$\varphi^{*}\varphi^{\prime}-\varphi^{\ast}{}^{\prime}\varphi=0,$ (21)

and the most general boundary condition that satisfies this constraint (21) is
given by [14]

$\varphi+i\varphi^{\prime}=U_{c}(\varphi-i\varphi^{\prime}),$ (22)

in terms of an arbitrary unitary matrix $U_{c}$ of
$L^{2}(\mathbb{R}^{2})$444There is a technical subtlity associated to the fact
that in higher dimension the baundary values $\varphi$ of the fields $\phi$
are singular which can be solved with a slight modification of the theory [18]
. Now if $\phi$ is a monochromatic field

$\phi=e^{i\omega t}\chi(x^{1},x^{2}),$

we have that

$\dot{\phi}=i\omega\phi\hskip 14.22636pt\dot{\phi^{*}}=-i\omega\phi^{*}.$

Thus, for monochromatic fields the vanishing condition associated to the
conservation of energy (6) implies the conservation of charge (4)

$\dot{\varphi^{*}}\varphi^{\prime}+\varphi^{*\prime}\dot{\varphi}=-i\omega\left(\varphi^{*}\varphi^{\prime}-\varphi\varphi^{*\prime}\right)=0,$

and viceversa, the conservation of charge implies the conservation of energy.

However, the fact that $U$ is independent of $\omega$ does not guarantee that
boundary conditions which preserve the energy also preserve the electric
charge for multifrequency fields (22). Only boundary conditions of the type

$U=\begin{pmatrix}U_{c}&0\\\ 0&U_{c}^{\ast}\end{pmatrix},$

where $U_{c}$ has all eigenvalues $\pm 1$ define boundary conditions that
preserve energy and electric charge simultaneously. This type of boundary
conditions includes Dirichlet, Neumann or periodic boundary conditions [22] .

The application to gauge fields is straightforward as we shall see below.
However, the methods based on energy consevation fail when applied to
gravitation and other higher spin gauge fields fields [23].

##  5 Yang-Mills theory

Let us now consider the interesting and non-trivial case of Yang-Mills gauge
theories. The Lagrangian density is given by

$\mathcal{L}=-\frac{1}{16\pi}\text{tr }F_{\mu\nu}F^{\mu\nu}$ (23)

where
$F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}-[A_{\mu},A_{\nu}]$,
$A_{\mu}$ being the gauge fields $A_{\mu}=A_{\mu}^{a}T^{a}$ and $T^{a}$ the
generators of the Lie algebra. The symmetric energy-momentum tensor has the
form

$T_{\mu\nu}=\frac{1}{4\pi}\text{tr
}\left(\frac{1}{2}\eta_{\nu\beta}F_{\mu\alpha}F^{\alpha\beta}+\frac{1}{2}\eta_{\mu\beta}F_{\nu\alpha}F^{\alpha\beta}+\frac{1}{4}\eta_{\mu\nu}F_{\alpha\beta}F^{\alpha\beta}\right).$
(24)

Applying the translation symmetry as we did for the scalar case we have the
following conservation energy law

$\partial^{\mu}T_{\mu 0}=\partial_{t}\mathcal{E}+\partial^{i}P_{i}=0,$ (25)

for the energy density

$\mathcal{E}=\frac{1}{8\pi}\text{tr }(\textbf{E}^{2}+\textbf{B}^{2})$

and the non-abelian Poynting vector

$\textbf{P}=\frac{1}{4\pi}\text{tr }\left(\textbf{E}\times\textbf{B}\right),$

where $E_{i}=F_{0i}$ and $B_{i}=-\frac{1}{2}\epsilon_{ijk}F_{jk}$ are the
chromo-electric and chromo-magnetic fields for the non-abelian case. Thus, for
any bounded domain $\Omega$ with regular boundary $\partial\Omega$ we have

$\frac{d}{dt}E_{\Omega}=\int_{\Omega}\partial_{t}\mathcal{E}\
d^{3}x=-\int_{\Omega}\partial_{i}P^{i}\
d^{3}x=-\int_{\partial\Omega}n_{i}P^{i}\ d^{2}x,$

where $\textbf{n}=n^{i}$ is the normal vector to the boundary surface
$\partial\Omega$.

###  5.1 Gauge fields in half space

As in the scalar theory, we consider the half-space
$\Omega=\\{\textbf{x}=(x^{1},x^{2},x^{3})|x^{3}\geq 0\\}$ that has a boundary
$\partial\Omega=\\{\textbf{x}=(x^{1},x^{2},0)\\}$. The conservation of energy
is given in this case by

$\frac{d}{dt}E_{\Omega}=\int_{\partial\Omega}T_{30}\
dx^{1}dx^{2}=\int_{\partial\Omega}\text{tr
}(B_{2}E_{1}-B_{1}E_{2})dx^{1}dx^{2},$

Restricting ourselves to homogeneous local boundary conditions that are
invariant under translation along $\partial\Omega$, the requirement of energy
conservation leads to

$\text{tr }(B_{2}E_{1}-B_{1}E_{2})=0.$ (26)

In terms of the following auxiliary vectors

$H=\begin{pmatrix}E_{1}+iE_{2}+i(B_{1}+iB_{2})\\\
E_{1}-iE_{2}-i(B_{1}-iB_{2})\end{pmatrix}\hskip
28.45274ptG=\begin{pmatrix}E_{1}-iE_{2}+i(B_{1}-iB_{2})\\\
E_{1}+iE_{2}-i(B_{1}+iB_{2})\end{pmatrix},$

the condition (26) reads

$\text{tr }(\lvert H\rvert^{2}-\lvert G\rvert^{2})=-8\text{tr
}(E_{1}B_{2}-E_{2}B_{1}),$

where we use that the trace of the product is conmutative. Thus, for
homogeneous local boundary conditions the most general solution satisfying
that

$\lvert H\rvert^{2}=\lvert G\rvert^{2}$

is

$H=U_{b}\otimes U_{c}\ G,$ (27)

where $U_{b}$ is an unitary 2x2 matrix acting on the ${\mathbb{C}}^{2}$
chromoelectromagnetic vector fields $G$ and $U_{c}$ is an unitary matrix that
acts on the color component of the fields. Since $F_{\mu\nu}$ is covariant
under gauge transformations $\Phi$,
$F_{\mu\nu}\rightarrow\Phi^{-1}F_{\mu\nu}\Phi$, $E$ and $B$ are always gauge
covariant, but not gauge invariant and therefore the most general solution of
(27) is not gauge invariant

$\Phi^{-1}H\Phi=\Phi^{-1}U_{b}\otimes U_{c}\ G\Phi\neq U_{b}\otimes U_{c}\
\Phi^{-1}G\Phi.$ (28)

However, all physically consistent boundary conditions must be gauge
invariant. This implies that the only physically possible boundary conditions
are those where $U_{c}$ is an element of the center of the gauge group, i.e.
$U_{c}\in Z(G)$. In such a case the transformation (28) can be rewritten as

$\Phi^{-1}H\Phi=U_{b}\otimes U_{c}\ \Phi^{-1}G\Phi$

and the corresponding boundary conditions of gauge fields are gauge invariant.
In the case where the gauge group is $SU(N)$ we have that the center of the
group has the form

$Z(SU(N))=\\{e^{\frac{2\pi in}{N}}I_{N};n=0,1\ldots,N-1\\}.$

Using this expression we can write the boundary condition (27) as

$H=e^{\frac{2\pi in}{N}}U_{b}\otimes I_{N}\ G.$ (29)

This extra phase that the color adds can be reabsorbed into $U_{b}$ since it
does not affect the condition of being unitary. Thus, we can rewrite (29) as

$H=U\otimes I_{N}G,$

where $U=e^{\frac{2\pi in}{N}}U_{b}$ is an unitary 2x2 matrix. If we conjugate
this relation we have

$\sigma_{1}H=U^{*}\otimes I_{N}\sigma_{1}G,$

and we get that the unitary 2x2 matrix has to satisfy the extra condition

$U=\sigma_{1}U^{*}\sigma_{1}.$

Thus, we get the same type of matrices as for the charged scalar field in half
space with the two disjoint components (10) and (11). In the first case (10)
the local boundary conditions are given by

$\begin{pmatrix}E_{2}\\\
B_{2}\end{pmatrix}=\tan\frac{a}{2}\begin{pmatrix}E_{1}\\\ B_{1}\end{pmatrix},$
(30)

whereas in the second case (11) they are

$\begin{pmatrix}B_{1}\\\
B_{2}\end{pmatrix}=-\tan\frac{b}{2}\begin{pmatrix}E_{1}\\\
E_{2}\end{pmatrix}.$ (31)

Some interesting particular cases are:

  1. i)

$U_{+}(0)=\mathbb{I}$ $\Rightarrow$ $E_{2}=B_{2}=0$.

  2. ii)

$U_{+}(\pi)=-\mathbb{I}$ $\Rightarrow$ $E_{1}=B_{1}=0$.

  3. i)

$U_{+}(\pm\frac{\pi}{2})=\pm i\sigma_{3}$ $\Rightarrow$ $E_{2}=\pm E_{1}$ and
$B_{2}=\pm B_{1}$.

  4. ii)

$U_{-}(0)=\sigma_{1}$ $\Rightarrow$ $B_{1}=B_{2}=0$.

  5. iii)

$U_{-}(\pi)=-\sigma_{1}$ $\Rightarrow$ $E_{1}=E_{2}=0$.

  6. iv)

$U_{-}(\pm\frac{\pi}{2})=\pm\sigma_{2}$ $\Rightarrow$ $B_{1}=\mp E_{1}$ and
$B_{2}=\mp E_{2}$.

###  5.2 Boundary conditions for the Yang-Mills potentials

The above boundary conditions can be formulated in terms of chromo-
electromagnetic potentials.

In the case $U_{-}$ (31) we have the boundary conditions

$\displaystyle\partial_{2}A_{3}-A^{\prime}_{2}-[A_{2},A_{3}]=-\tan\frac{b}{2}(\dot{A}_{1}-\partial_{1}A_{0}-[A_{0},A_{1}]),$
(32)
$\displaystyle\partial_{1}A_{3}-A^{\prime}_{1}-[A_{1},A_{3}]=\tan\frac{b}{2}(\dot{A}_{2}-\partial_{2}A_{0}-[A_{0},A_{2}]),$
(33)

in terms of $A_{\mu}$, where

$\dot{A}=\partial_{t}A|_{\partial\Omega}$ (34)

and

$A^{\prime}=\partial_{3}A|_{\partial\Omega}.$ (35)

The above boundary conditions (32) and (33) can also be rewritten using the
covariant derivative $D_{\mu}A_{\nu}=\partial_{\mu}A_{\nu}-[A_{\mu},A_{\nu}]$,

$\displaystyle
D_{2}A_{3}-A_{2}^{\prime}=-\tan\frac{b}{2}\left(D_{0}A_{1}-D_{1}A_{0}-[A_{0},A_{1}]\right),$
$\displaystyle
D_{1}A_{3}-A_{1}^{\prime}=\tan\frac{b}{2}\left(D_{0}A_{2}-D_{2}A_{0}-[A_{0},A_{2}]\right).$

Choosing as gauge fixing condition

$A_{0}=0,\hskip
14.22636pt\partial_{1}A_{1}+\partial_{2}A_{2}+A^{\prime}_{3}=0,$

we can rewrite the conditions as

$\begin{pmatrix}A_{1}^{\prime}\\\ A_{2}^{\prime}\\\
A_{3}^{\prime}\end{pmatrix}=M\begin{pmatrix}A_{1}\\\ A_{2}\\\
A_{3}\end{pmatrix},$

where

$M=\begin{pmatrix}0&-\tan\frac{b}{2}D_{0}&D_{1}\\\
\tan\frac{b}{2}D_{0}&0&D_{2}\\\ -D_{1}&-D_{2}&0\end{pmatrix}$ (36)

is a selfadjoint differential operator that is gauge covariant and relates the
normal derivatives of the gauge fields with the very gauge fields in that
gauge. The expression (5.2) is reminiscent of the one obtained for scalar
fields where the matrix $M$ is related to the family of unitary matrices by
means of a Cayley transform [14].

In the case of the family of solutions (30) given by $U_{+}$ we have

$\displaystyle\dot{A}_{2}-\partial_{2}A_{0}-[A_{0},A_{2}]=\tan\frac{a}{2}(\dot{A}_{1}-\partial_{1}A_{0}-[A_{0},A_{1}]),$
$\displaystyle\partial_{2}A_{3}-A^{\prime}_{2}-[A_{2},A_{3}]=-\tan\frac{a}{2}(\partial_{1}A_{3}-A^{\prime}_{1}-[A_{1},A_{3}]),$

that can be rewritten using the covariant derivatives

$\displaystyle
D_{0}A_{2}-\partial_{2}A_{0}=\tan\frac{a}{2}(D_{0}A_{1}-\partial_{1}A_{0})$
$\displaystyle\partial_{2}A_{3}-D_{3}A_{2}=-\tan\frac{a}{2}(D_{1}A_{3}-A^{\prime}_{1}).$

Choosing a diferent gauge fixing condition

$A_{1}=0,\hskip
14.22636pt\partial_{t}A_{0}+\partial_{2}A_{2}+\partial_{3}A_{3}=0,$

these boundary conditions reduce to

$\partial_{2}\begin{pmatrix}A_{0}\\\ A_{2}\\\
A_{3}\end{pmatrix}=\begin{pmatrix}\tan\frac{a}{2}D_{1}&D_{0}&0\\\
-D_{0}&0&-D_{3}&\\\
0&D_{3}&-\tan\frac{a}{2}D_{1}\end{pmatrix}\begin{pmatrix}A_{0}\\\ A_{2}\\\
A_{3}\end{pmatrix}$

in terms of the selfadjoint differential operator that is gauge covariant.

##  6 Maxwell theory

The analysis of previous section includes the Maxwell theory since $U(1)$
gauge theories are simpler cases of Yang-Mills theories. In this case the
chromoelectromagentic fields become the standard electromagnetic fields, and
covariant derivatives are replaced by standard derivatives in the selfadjoint
differential operators that appear in the motion equations.

An interesting case of study is the analysis of the boundary conditions on the
interface between two different electromagnetic active media. Let us consider
a linear, non dispersive and isotropic case where $D=\epsilon E$ and $B=\mu
H$, with $\epsilon$ the dielectric permitivity and $\mu$ the magnetic
permeability of the material, and work on Gaussian units. In this type of
media the electromagnetic energy density is given by

$\mathcal{E}=\frac{1}{8\pi}(\textbf{E}\cdot\textbf{D}+\textbf{B}\cdot\textbf{H})$
(37)

and the Poynting vector by

$\mathbf{S}=\frac{c}{4\pi}\mathbf{E}\times\mathbf{H},$ (38)

where we recovered the dependency on the vacuum speed of light $c$. The energy
conservation law (25) becomes

$\partial_{t}\mathcal{E}+\nabla\textbf{S}=0.$ (39)

If we consider the case of two different electromagnetic active media in
$\Omega_{+}=\\{\textbf{x}=(x^{1},x^{2},x^{3})|x^{3}\geq 0\\}$ and
$\Omega_{-}=\\{\textbf{x}=(x^{1},x^{2},x^{3})|x^{3}\leq 0\\}$ which are
separated by the boundary $\partial\Omega=\\{\textbf{x}=(x^{1},x^{2},0)\\}$,
the condition for conservation of energy for both media is given by

$\frac{d}{dt}(E_{\Omega_{+}}+E_{\Omega_{-}})=\int_{\Omega_{+}}\partial_{t}\mathcal{E}\
d^{3}x\ +\int_{\Omega_{-}}\partial_{t}\mathcal{E}\
d^{3}x=\int_{\partial\Omega_{+}}S^{3}dx^{1}dx^{2}\
-\int_{\partial\Omega_{-}}S^{3}dx^{1}dx^{2},$

where $\partial\Omega_{+}$ is the boundary at the right side and
$\partial\Omega_{-}$ is the boundary at the left side. Inserting the explicit
expression the Poynting vector for both materials the condition has the form

$\frac{c}{4\pi\mu^{+}}\int_{\partial{\Omega_{+}}}\\!\\!\\!(E_{1}^{+}B_{2}^{+}-E_{2}^{+}B_{1}^{+})\
dx^{1}dx^{2}-\frac{c}{4\pi\mu^{-}}\int_{\partial\Omega_{-}}\\!\\!\\!(E_{1}^{-}B_{2}^{-}-E_{2}^{-}B_{1}^{-})\
dx^{1}dx^{2}=0$

where $\mu^{\pm}$ is the magnetic permeability in each media, and $E^{\pm}$
and $B^{\pm}$ are the electromagnetics field on each side. Considering only
homogeneous boundary conditions that are invariant under translations along
the boundary plane $\partial\Omega$, the energy conservation condition reduces
to

$\frac{1}{\mu^{+}}(E_{1}^{+}B_{2}^{+}-E_{2}^{+}B_{1}^{+})-\frac{1}{\mu^{-}}(E_{1}^{-}B_{2}^{-}-E_{2}^{-}B_{1}^{-})=0.$
(40)

We can define two auxiliary vectors

$\displaystyle
H=\begin{pmatrix}\frac{1}{\sqrt{\mu^{-}}}\left(E_{1}^{-}+iE_{2}^{-}+i(B_{1}^{-}+iB_{2}^{-})\right)\\\
\frac{1}{\sqrt{\mu^{+}}}\left(E_{1}^{+}-iE_{2}^{+}+i(B_{1}^{+}-iB_{2}^{+})\right)\\\
\frac{1}{\sqrt{\mu^{-}}}\left(E_{1}^{-}-iE_{2}^{-}-i(B_{1}^{-}-iB_{2}^{-})\right)\\\
\frac{1}{\sqrt{\mu^{+}}}\left(E_{1}^{+}+iE_{2}^{+}-i(B_{1}^{+}+iB_{2}^{+})\right)\end{pmatrix},$
$\displaystyle
G=\begin{pmatrix}\frac{1}{\sqrt{\mu^{-}}}\left(E_{1}^{-}-iE_{2}^{-}+i(B_{1}^{-}-iB_{2}^{-})\right)\\\
\frac{1}{\sqrt{\mu^{+}}}\left(E_{1}^{+}+iE_{2}^{+}+i(B_{1}^{+}+iB_{2}^{+})\right)\\\
\frac{1}{\sqrt{\mu^{-}}}\left(E_{1}^{-}+iE_{2}^{-}-i(B_{1}^{-}+iB_{2}^{-})\right)\\\
\frac{1}{\sqrt{\mu^{+}}}\left(E_{1}^{+}-iE_{2}^{+}-i(B_{1}^{+}-iB_{2}^{+})\right)\end{pmatrix},$

in terms of which the condition reads

$\lvert H\rvert^{2}-\lvert
G\rvert^{2}=8\left(\frac{1}{\mu^{+}}E_{1}^{+}B_{2}^{+}-\frac{1}{\mu^{+}}E_{2}^{+}B_{1}^{+}-\frac{1}{\mu^{-}}E_{1}^{-}B_{2}^{-}+\frac{1}{\mu^{-}}E_{2}^{-}B_{1}^{-}\right).$

Thus the general solution of the boundary condition is given by

$H=UG,$

where $U$ is a 4x4 unitary matrix. Gauge invariance is automatic since the
elements of $U(1)$ are simply a phase $e^{i\beta}$ which can be absorbed in
the unitary matrix. Conjugating this equation leads to

$\begin{pmatrix}0&\mathbb{I}\\\ \mathbb{I}&0\end{pmatrix}H=U^{*}\otimes
I_{N}\begin{pmatrix}0&\mathbb{I}\\\ \mathbb{I}&0\end{pmatrix}G,$

which gives an extra constraint for the unitary matrix

$U=\begin{pmatrix}0&\mathbb{I}\\\
\mathbb{I}&0\end{pmatrix}U^{*}\begin{pmatrix}0&\mathbb{I}\\\
\mathbb{I}&0\end{pmatrix}.$

Thus, we get the same 4x4 unitary matrix as for the scalar charged fields for
the two parallel plates, which are given by the two disjoint components (15)
and (16). Let see some examples of this boundary conditions:

  1. i)

$U_{+}\left(0,0,0,\frac{\pi}{2},0,-\frac{\pi}{2}\right)=\begin{pmatrix}0&I\\\
I&0\end{pmatrix}$. In this case the boundary magnetic fields vanish in both
sides $B_{1}^{+}=B^{+}_{2}=B_{1}^{-}=B_{2}^{-}=0$.

  2. ii)

$U_{+}\left(0,0,0,\frac{\pi}{2},0,\frac{\pi}{2}\right)=-\begin{pmatrix}0&I\\\
I&0\end{pmatrix}$. In this case the boundary electric fields vanish in both
sides $E_{1}^{+}=E_{2}^{+}=E_{1}^{-}=E_{2}^{-}=0$.

  3. iii)

$U_{+}\left(\frac{\pi}{2},\frac{\pi}{2},\mp\frac{\pi}{2},0,0,0\right)=\pm\begin{pmatrix}\sigma_{1}&0\\\
0&\sigma_{1}\end{pmatrix}$. In this case $E_{1}^{+}=\pm\alpha E_{1}^{-}$,
$E_{2}^{+}=\pm\alpha E_{2}^{-}$, $B_{1}^{+}=\pm\alpha B_{1}^{-}$ and
$B_{2}^{+}=\pm\alpha B_{1}^{-}$, which leads to periodic/antiperiodic boundary
conditions with a jump given by $\alpha=\sqrt{\frac{\mu^{+}}{\mu^{-}}}$.

These boundary conditions include those obtained in the case of one simple
wall, like in the first and second case, where simply the Poynting vector is
zero at each side of the boundary so there is no flux of energy through the
materials. The new phenomena due to the two different materials occurs in the
boundary conditions where the fluxes of energy from each side are not zero.
The constraint of energy conservation induces in this case a jump in the
values of the electromagnetic fields related to the different magnetic
permeabilties of the materials. This can be seen in the periodic and
antiperiodic boundary conditions showed above, where instead of just getting
the same values at the sides of the boundary there is a jump on them across
the boundary.

##  7 Conclusions

We have developed a consistent theory of boundary conditions for relativistic
field theories based on the conservation of fundamental quantities: charge and
energy. These are the two generic conservation laws of Nature. Boundary
conditions which preserve charge are quite well known [14]–[19] whereas those
that preserve the energy have not been explored so deeply. We have analyzed in
this paper the theory of boundary conditions from this perspective. We found
infinite families of energy preserving boundary conditions in scalar,
electromagnetic and non-abelian gauge theories. In some cases the boundary
conditions preserve both charge and energy. In the Maxwell and Yang-Mills
theory the boundary conditions are also gauge invariant. We also have shown
how this method can be used to describe boundaries between different material
media applying the conservation of energy.

Spinor fields have not been considered in this paper but the extension of the
analysis to this case is straightforward. However, the generalization for
gravitation and higher spin fields is not so simple. In the gravitational case
the difficulty resides on the fact that the vanishing condition required to
define the appropriate boundary conditions is automatically satisfied by any
solution of classical equation of motion. In this case a generalization of the
application of the above conservation laws is required.

Another avenue worth to explore is the extension of this analysis for
regularizations of field theories on the lattice. This is very interesting for
numerical approaches both for the classical dynamics of interacting field
theories and for the analysis of non-perturbative effects in the corresponding
quantum theories.

One of the interesting boundary effects associated to the boundary conditions
is the dependence of the non-perturbative Casimir energy in 2+1 gauge theories
on the boundary conditions [24][25]. The exponential decay of the Casimir
energy with the distance between two parallel plates is dependent on the
boundary conditions and points out to a new mass parameter much lower than the
glueball mass, opening the door to new interpretations of the confinement
mechanism,

## Acknowledgments

M. A. thanks to A.P. Balachandran, A. Ibort and G. Marmo for so many
discussions.

We are partially supported by Spanish Grant PGC2022-126078NB-C21 funded by
MCIN/AEI/ 10.13039/501100011033 and “ERDF A way of making EuropeGrant , the
DGA-FSE grant 2020-E21-17R of the Aragon Government and the European Union -
NextGenerationEU Recovery and Resilience Program on ’Astrofísica y Física de
Altas Energías’ CEFCA-CAPA-ITAINNOVA

## References

  * [1] J. Zhang, L. Zhang and W. Xu , J. Phys. D: Appl. Phys. 45 (2012) 113001
  * [2] R. B. Laughlin Phys. Rev. B 23 (1981) 5632.
  * [3] B. I. Halperin Phys. Rev. B 25 (4) (1982) 2185.
  * [4] Y. Hatsugai, Phys. Rev. Lett. 71 (1993) 3697.
  * [5] C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95 (2005) 146802.
  * [6] L. Fu and C. L. Kane, Phys. Rev. B 76 (2007) 045302.
  * [7] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82 (2010) 3045–3067.
  * [8] J. E. Moore, Nature 464 (2010) 194–198.
  * [9] X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys.83 (2011) 1057–1110.
  * [10] B. A. Bernevig and T. L. Hughes, Topological insulators and topological superconductors, Princeton U. Press (2013).
  * [11] M. Zahid Hasan, S.-Y. Xu, M. Neupane, in Topological Insulators: Fundamentals and Perspectives, Wiley and Sons, Weinheim, Germany (2015) 55–100.
  * [12] H. B. G. Casimir, Proc. K. Ned. Akad. Wet. 51 (1948) 793.
  * [13] A. P. Balachandran, G. Bimonte, G. Marmo and A. Simoni, Nucl. Phys. B446 (1995) 299.
  * [14] M. Asorey, A. Ibort and G. Marmo, Int. J. Mod. Phys. A 20 (2005) 1001
  * [15] A.D. Shapere, F. Wilczek, Z. Xiong, Models of topology change, arXiv:1210.3545 (2012)
  * [16] M. Asorey et al, Quantum Physics and Fluctuating Topologies: Survey, arXiv:1211.6882 [hep-th] (2012).
  * [17] M. Asorey, D. García. Álvarez and J. M. Muñoz-Castañeda, J. Phys. A 39 (21)(2006) 6127
  * [18] M. Asorey and J. M. Muñoz-Castañeda, Nucl. Phys. B 874 (2013) 852;
  * [19] M. Asorey, D. García-Alvarez, and J. M. Muñoz-Castañeda, Int. J. Geom. Meth. Mod. Phys. 12 (06) (2015) 1560004
  * [20] J. Polchinski, String Theory, Vol. 1 , Cambridge Univ. Press (2005).
  * [21] F. Antoneli, M. Forger and P. Gaviria, J. Lie Theory, 22 (4) (2012) 949-1024.
  * [22] M. Asorey, D. García-Alvarez, and J. M. Muñoz-Castañeda, J. Phys. A 40 (2007) 6767.
  * [23] J. Asorey, M. Asorey and D. García-Alvarez Phys. Rev. D 92, 103517 (2015)
  * [24] M.N. Chernodub, V.A. Goy, A.V. Molochkov, H.H Nguyen, Phys. Rev. Lett. 121 (2018) 191601.
  * [25] D. Karabali, V.P. Nair, Phys. Rev. D 98 (2018) 105009