# Super-resolving Dark Matter Halos using Generative Deep Learning

D. Schaurecker,11footnotetext: Corresponding author. Y. Li  J. Tinker  S. Ho
and A. Refregier

###### Abstract

Generative deep learning methods built upon Convolutional Neural Networks
(CNNs) provide great tools for predicting non-linear structure in cosmology.
In this work we predict high resolution dark matter halos from large scale,
low resolution dark matter only simulations. This is achieved by mapping lower
resolution to higher resolution density fields of simulations sharing the same
cosmology, initial conditions and box-sizes. To resolve structure down to a
factor of 8 increase in mass resolution, we use a variation of U-Net with a
conditional Generative Adversarial Network (GAN), generating output that
visually and statistically matches the high resolution target extremely well.
This suggests that our method can be used to create high resolution density
output over Gpc/h box-sizes from low resolution simulations with negligible
computational effort.

##  1 Introduction

Cosmological galaxy surveys, both current and planned, span increasingly large
volumes, eg. DESI [1], EUCLID [2], DES [3] and LSST [4]. In the previous
decade, surveys such as the Baryon Oscillation Spectroscopic Survey [5] have
used bright red galaxies as tracers of large scale structure. To make larger
maps, the current and next generation of surveys rely on utilizing fainter
galaxies. Spectroscopic surveys specifically are planning to observe galaxies
targeted by the presence of star formation within them. Such galaxies, some of
them Emission Line Galaxies (ELGs; [6]), are the cosmological workhorse of the
next decade. An important issue with creating bigger large-scale structure
with smaller galaxies is being able to simulate such maps computationally.
ELGs reside in halos at or below the mass scale of our own Milky Way [7, 8].
Thus, simulations of such galaxy samples require not just large volume but
also higher mass resolution than previously needed. A single simulation with
equal volume of the DESI ELG sample that properly resolves all halos in which
ELGs could form would require upwards of $9500^{3}$ particles, assuming a
volume of 122 Gpc and a resolution of halos at around $1000$ particles for
$10^{11}\ \textup{M}_{\Sun}$ halos. We often need a large suite of such
simulations to explore various parameters or build up statistics. This will
quickly become computationally expensive.

A promising solution to this problem is the use of deep learning networks.
Various works in recent years applied convolutional neural networks to
cosmological applications successfully, generally increasing small scale
structure information (eg. galaxy distributions in [9], mapping dark matter to
galaxies in [10, 9], extracting cosmological information while marginalizing
over baryonic effects in [11], or for general simulation output [12], [13] and
[14]). In the case of creating cosmological simulations spanning large
volumes, a simple first approach is to use DM-only density fields.

In this paper, we propose a novel method to predict the $z=0$ density field of
a high resolution simulation simply by giving it the $z=0$ density field of a
low resolution simulation of the same volume, run using same cosmologies and
initial conditions. High and low resolution here refers to the simulation
particle’s mass and the number of particles used inside the simulation volume.
Once a network is trained successfully it takes little computational effort to
produce new output from low resolution density fields. This would allow one to
run a computationally less expensive, lower resolution simulation spanning
Gpc/h in volume, and then using the resulting density field to predict the
structure at higher resolution using the CNN. The halos and subhalos in this
high-resolution output can be populated with galaxies in any number of ways,
as required by the specific goal of the mock galaxy catalog (see, [15] for a
review on the galaxy-halo connection). This paper provides a first proof of
work that uses density fields to super resolve dark matter halos. Previous
works either focus on other statistics such as void abundance [14], or use
particle displacement fields that assumes lattice pre-initial conditions (pre-
ICs, the Lagrangian particle positions) [16] and therefore cannot be applied
to many of the state-of-the-art simulations that use glass pre-ICs.

In particular, our presented method here can be applied to dark matter
simulations with glass initial conditions [17] (here after: IC) while previous
super-resolution techniques such as [18] cannot. Many state of the art
simulations use glass IC, such as Illustris project [19] while there is
currently no satisfactory way to interpolate those glass IC particles onto a
grid in order to compute the displacement field needed by previous super-
resolution methods such as [16]. Any kind of interpolation technique will
always create voxels that contain 0 particles, thus making displacement field
calculation impossible for that pixel. Rather than wasting computational
resources on new costly simulations for super-resolution tasks, we develop the
alternative Eulerian method to exploit the existing state-of-the-art
cosmological simulations that started with glass IC.

##  2 Method

###  2.1 Dataset

In this work, we used redshift=0 data from the Illustris-2-Dark (high-res) and
Illustris-3-Dark (low-res) dark matter only simulations [19] for training and
testing. Illustris-2-Dark has a box size of $75$ Mpc/h, and contains $910^{3}$
($m_{\text{dm}}=6\times 10^{7}\ \textup{M}_{\Sun}$, high-res) and
Illustris-3-Dark contains $455^{3}$ ($m_{\text{dm}}=4.8\times 10^{8}\
\textup{M}_{\Sun}$, low-res) dark matter particles. This results in both
simulation boxes having the same mean mass density, meaning if more particles
are simulated in the higher resolution simulation (Illustris-2-Dark in this
case) they have to be lighter. The cosmology for these simulations is
consistent with WMAP-9 measurements [20] and the glass tiled initial
conditions are described in detail in [21]. Similar work done by [16] used
displacement fields created from particle positions for training, and resolved
only higher mass halos (lowest physical halo mass resolved at $z=0$ is
$10^{11}\ \textup{M}_{\Sun}$, containing 120 tracer particles). However,
displacement fields are only naturally defined for lattice initial conditions,
and are difficult to produce from glass tiled initial conditions, making this
work ideal for simulations that do not have initial conditions that can be
turned into lattice like configuration.

Most state of the art hydro simulations (like Illustris) use glass ICs to
improve statistics at high redshift snapshots, thus the underlying DM-only
simulations are usually initialized using glass tiles, which are simply copied
next to each other to fill the simulation volume. There is no simple way to
interpolate those initial particles onto a grid in order to then use the
displacement field super-resolution method pioneered in [16]. Any kind of
interpolation technique will always have some residual pixels contain zero or
2 particles. Finding a general, automatic, way to fill those remaining pixels
with exactly one particle in a coherent way is very impractical if not
impossible, and has so far not been done. Additionally using displacement
fields also fixes the number of pixels a given simulation has to be divided
into. If a simulation contains for example $512^{3}$ particles, exactly one
pixel per spatial coordinate and particle needs to be created for training and
testing, complicating general applicability on any given simulation.

In our proposed method, we do however use grids to create density grids, but
we do not need to create a displacement field (which requires having at least
one particle per pixel in order to find the displacement field). Our deep
learning method uses convolution layers, which are typically applied imaging
data that has translational symmetry. Convolutional neural networks (CNNs)
generally work well with equidistantly gridded n-dimensional data. In our
case, a fixed spatial resolution is perfect for a convolutional neural
network. We prepare dark matter density maps, with a $2048^{3}$ grid, creating
a number count density field where each voxel of side-length $\approx 36.6$
kpc/h contains information about the number of particles that lie inside that
small region. Each of the data sets is then subsequently divided into eight
sub-cubes of equal size (37.5 Mpc/h on the side), from which six are used for
training, one for validation, and the last one for testing. This results in
the low resolution input (Illustris-3-Dark) and high resolution target
(Illustris-2-Dark) simulations being divided into $32768$ 3-dimensional
$64^{3}$ pixel cubes, of which an eighth is used for testing.

###  2.2 Data pre-processing

Neural networks benefit from certain input and target data distributions more
than others, in particular, they should not have a very long tail or an
extremely large dynamic range. This makes the following pre-processing steps
necessary for a successful and efficient training process.

####  2.2.1 Mesh painting

To achieve a data PDF shape which is continuous and roughly resembles a
Gaussian distribution, the Illustris simulations’ catalogs are painted onto a
mesh using a TSC (Triangular Shaped Cloud) window kernel, since most voxels
are empty, this helps to "smooth" out the integer number count values
especially in low density voxels which importantly helps smoothing out the
individual training cubes’ PDF curves and removes large discontinuities.

####  2.2.2 Normalization

The ranges of feature distributions vary a lot among different training cubes
and are also dependent on the training cube’s position inside the simulation.
To speed up and stabilize the training process, we normalize all training
cubes as we, for example, do not want to learn to predict certain features
better than others. The following simple normalization $N$ is used:

$N_{\text{log1p}}(x)=\text{log}\left(\epsilon+1+x\right),$ (2.1)

where the not strictly necessary factor $\epsilon$ is set to $\epsilon\equiv
10^{-7}$.

####  2.2.3 Cube division

The simulations painted on 3D meshes are too big to fit in a single GPU memory
all at once. We then choose to use sub-volumes with $112^{3}$ voxels as the
input low-res data to the networks, and $64^{3}$ as the target high-res data.
These include periodically padded regions of width of 24 voxels, and require
no further “paddings” once inside the networks. By doing so we fully preserve
the translational symmetry, which is crucial for generating outputs without
artifacts near their boundaries. Using the model described in 2.3, this
results in a generated output cube size of $64^{3}$ voxels.

###  2.3 Models and training

At the $36.6$ kpc/h pixel scale, dark matter in the low and high resolution
simulations are already spatially shifted from each other, due to the addition
of the high-frequency modes in the high-resolution simulation. The high
density regions can shift by a few pixels. This makes the simpler training
approach of supervised learning (e.g. using a mean squared error loss
function) impossible as the model does not have enough information to predict
this shift and will simply blur the output. Recent works (e.g. on displacement
fields not density fields in [16]) showed that an unsupervised approach to
deep learning is a possible pathway forward given the limitation described
here.

First proposed by Ian Goodfellow in 2014 [22] Generative Adversarial Networks
(GANs) were introduced as a new class of generative models. GANs are able to
produce output from input sampled from a random noise distribution and
generally do not rely on information from the input data to make predictions.
GANs consist of two different networks that are trained in parallel to each
other: the generator G and the discriminator D. For a given set of real
training samples, the generator learns to produce fake samples that look
similar to the real ones, while the discriminator tries to learn to
differentiate between real and fake samples. Formally (as in [22]): we can
define a mapping to data space $G(z,\theta_{g})$ where $z$ is taken from a
prior, and usually random, distribution $p_{z}(z)$, whereas $\theta_{g}$ are
G’s learnable parameters. Similarly, we define a mapping $D(x,\theta_{d})$
that outputs a single scalar representing the probability that $x$ was taken
from the data distribution ($p_{\text{data}}$) and not from $p_{g}$ (the
"fake" distribution produced by G). Training both D and G means maximizing the
probability of D correctly identifying fake and real data, while
simultaneously training G to minimize $\log\left(1-D(G(z)\right)$. This means
that D and G play a minmax game with value function:

$\displaystyle\min_{G}\max_{D}V(G,D)=$ $\displaystyle\mathbb{E}_{x\sim
p_{\text{data}}(x)}\left[\log D(x)\right]+$ (2.2)
$\displaystyle\mathbb{E}_{z\sim p_{z}(z)}\left[\log(1-D(G(z)))\right],$

where $\mathbb{E}$ denotes the expectation value, which in practice is
replaced by the average over a batch with stochastic gradient descent. As both
D and G improve in their training, the desired end result would be that G
produces fake samples that D is no longer to differentiate, thus having a 0.5
probability to guess correctly. This state is also called Nash-equilibrium
[23]. Of course, the training of actual models will stop before this happens,
but if the training was stable, G might be trained well enough to produce
valuable results before Nash-equilibrium is reached.

Nevertheless, in this work G needs to generate output not from a randomly
sampled input distribution, but from the low resolution training cubes and
from white noise fields in the noise layers. In order to help the
discriminator differentiate between real and fake data, the low-res input is
concatenated to the generated output and high-res target before passing the
data through D. The training process thus is no longer purely unsupervised,
and our method is a conditional GAN, an adapted version of the standard
approach first introduced by [24], which utilized a different architecture.

This results in the following loss functions for the discriminator and
generator in this work:

$\displaystyle\begin{aligned} \mathcal{L}_{D}=\
&-\mathbb{E}_{\hat{x},x}\left[\log D(\hat{x},x)\right]\\\
&-\mathbb{E}_{\hat{x},z}\left[\log(1-D(\hat{x},G(\hat{x},z)))\right]\\\
&-\gamma\ \mathbb{E}_{\hat{x},x}[\left\|\nabla
D(\hat{x},x)\right\|^{2}]\end{aligned}$ (2.3) $\displaystyle\begin{aligned}
\mathcal{L}_{G}=\
&-\mathbb{E}_{\hat{x},z}\left[\log(D(\hat{x},G(\hat{x},z)))\right],\end{aligned}$
(2.4)

where $x$ is sampled from the high-res distribution $p_{\text{high-res}}$,
$\hat{x}$ is sampled from the low-res distribution $p_{\text{low-res}}$ and
$z$ is sampled from a random white noise distribution $p_{z}$. The third term
in equation 2.3 is the $R_{1}$ regularization (introduced in [25]), which is
applied for real data only. The chosen penalty weight $\gamma=5$ is constant.

####  2.3.1 Models

##### Discriminator:

Training is not too sensitive on the exact discriminator architecture, as long
as the main idea of progressively extracting higher-level features by down-
sampling convolutions is followed. Figure 1 depicts a qualitative
approximation of the architecture as not all shapes and sizes are completely
accurate. It is important to understand that each convolution changes the
shape of the output feature map slightly, but since this would barely be
visible in the graph it is neglected here. Additionally, in this case the
boxes indicate the operator acting on the feature map of each stage. Meaning
that the actual data being passed through the network follows the path
indicated by the arrows and is reshaped by the operators along the way. For
the exact model architecture, please refer to the data availability section.

Figure 1: A plot showing the rather simple architecture of the discriminator
used for GAN training. Of course, every stage applies 3D operations on a 3D
input, but the visualization works better in 2D. The net’s stages are
indicated as: light orange: a 3D convolution; red: LeakyReLU activation;
green: Residual Block: consisting of two 3D convolutions each followed by a
LeakyReLU activation, another 3D convolution and a linear down-sampling layer.
At the end D outputs a $4^{3}$ voxel cube which is averaged over before being
passed to a sigmoid inside the BCELoss function. For further details please
refer to the data availability section.

##### Generator:

The model’s exact architecture is extremely important to a successful training
process. The model, itself not too complicated or large, is a shallower
adaptation of the often used U-Net (introduced in [26], adapted in eg. [27,
12] and [28]). As explained in [29], convolutional neural nets are built on
translational equivariance since their component stages only depend on
relative spatial positions as they operate on local input regions. Thus at all
blocks each voxel’s receptive field has the same size, independent on how
large or small the initial input is. This translational equivariance is very
important as it ensures that the net will produce good output on a variety of
different inputs. The most important change from a conventional U-Net is that
we replace the transposed upsampling convolution by a tri-linear interpolation
followed by a usual convolution (detailed implementation in figure 2 and the
data availability section, explained further in [30]). This greatly reduced
the checkerboard effect produced by the transposed convolutions, thus helping
the training process immensely.
In order to help the generator predict the highly non-linear structure of the
high-res target, uniformly sampled random noise is added before each
convolution. This provides the basis upon which the generator can build its
small scale prediction on. Similar to the hierarchical structure formation of
dark matter, the generator’s U-Net structure causes the layers producing small
scale structure to be conditioned on the large scale features that are passed
through G. This inherent similarity makes U-Nets an obvious choice for working
with cosmological simulation data.

Figure 2: The U-Net used as a generator in this work. Each stage’s output’s
number of channels is depicted by the small number below. The spatial
dimensions of the feature maps at each stage are approximately depicted by the
width and height of each stage. The net’s stages are indicated as (different
colors as for D in figure 1): light red: a 3D convolution; light orange: a
noise layer followed by another 3D convolution; dark orange / brown: LeakyReLU
activation; green: noise layer followed by a 3D (down-)convolution; blue:
noise layer followed by a linear up-sampling layer, a 3D convolution, a
LeakyReLU activation, another noise layer, another 3D convolution and lastly
another activation; cat: concatenation of two outputs from two layers along
the channel dimension. For more exact details please refer to the data
availability section.

####  2.3.2 Training

Each epoch of the training process was divided into training with the training
data set, followed by validation on the entire validation data set. To
increase the training set’s distribution variability, we augmented the dataset
by flipping and permuting axes of the training cubes. This helps the network
learn discrete rotations and parity symmetries which are inherent in
cosmological simulations. Such augmentation was not done on the validation set
for deterministic results at every epoch. In the first epoch of training, only
the generator was trained using the standard supervised approach using MSE
Loss. This helped the generator learn large scale shapes of the low-res input
and high-res target which was beneficial to the later epochs of GAN training.
With the AdamW optimizer [31] and small weight-decay ($\lambda=10^{-4}$), we
adopt different learning rates for the two networks: $1\times 10^{-5}$ for the
discriminator and $5\times 10^{-5}$ for the generator, with a ratio similar to
the one suggested in [32]. It is important to notice that these learning rates
are not universally applicable and will need to be adjusted for different
training applications to reach an equilibrium training state. Furthermore a
$R_{1}$ regularization gradient penalty scheme [25] was used every 16 batches
to penalize the discriminator upon deviating from equilibrium on real data
alone, thus increasing training stability.

##  3 Results

The sections below will summarize the quality of the network’s output after
training it for 600 epochs using the method described in section 2.3.2. Not
every epoch produces equally good results, although at later stages of
training the output quality is stable. This fluctuation is due to both
networks (D and G) hovering around a local loss minimum, only incrementally
improving. Most statistics are more or less unaffected by this, but for
example the halo two-point function for small mass halos is very sensitive to
the exact output. The best performing epoch of the last 10 is selected, in
this case: epoch 599. Best performing in this case refers to good accordance
with all target statistics. Epoch 599 was not the only one producing good
results across all statistics, but was subjectively selected for providing the
best results. It should also be mentioned that extending that run further
until epoch 1200 did not increase the quality of the output as the training
seemed to hover around a local minimum.

The testing dataset consists, as described in section 2.1, of one eighth of
the entire simulation region. Thus all testing cubes combine, stitched back
together, a volume of $37.5$ Mpc/h. All plots will show comparisons between
the Illustris-3 low-res input, the generator’s output and the high-res
Illustris-2 target, labeled as low-res, generated and high-res.

Some statistics (eg. two-point function, see below) require the data to be
stored as particle catalogs and not as number-count-field pixel values. In
order to go back to catalog format in the most simple and physical way, two
steps are taken:

First, all voxels are Poisson sampled as follows:

$x_{\text{poisson}}=\begin{cases}x&\text{if}\ r>x-\lfloor x\rfloor\\\
x+1&\text{else},\end{cases}$ (3.1)

where $r$ is uniformly drawn between $0$ and $1$. This results in the net
output’s voxels being converted to integers only. Secondly, if a voxel’s
Poisson sampled value is eg. 4, four particles are randomly placed inside the
voxel’s volume and depending on that voxel’s relative position inside the
simulation, a position-catalog entry is made for each particle.

This method obviously completely erases any information about particle
positions below the pixel-size scale, but this information is already lost
upon creation of the training data. Thus this kind of pixelization poses
fundamental limits and challenges to using number count or density fields for
this training method, as compared to e.g. displacement field training (which
however does not apply to simulations with glass initial conditions such as
Illustris here). This effect will also become more apparent in section 3.2,
while we look at statistical properties of the output.

###  3.1 Visual Comparison

Figure 3: left: low-res catalog, middle: generated output catalog, right:
high-res catalog. A projection on the z-axis of each catalog’s particles,
where particles which are part of FoF halos, are highlighted in orange. The
generated output plot is indistinguishable from the high-res plot by eye on
large scales. At smaller scales, the shapes and positions of small halos look
statistically consistent but vary a bit between the two as expected, whereas
they are completely missing in the low-res plot. The depicted testing-region
box is not periodic as it only spans one eighth in volume of the entire
corresponding Illustris simulation, which spans $(75\,\textrm{Mpc}/h)^{3}$.
This figure was created with the help of Yueying Ni.

The visual differences between the low-res input and generated output / high-
res target simulation become clear immediately in figure 3, especially when
looking at the smaller mass halos that are missing in the low-res simulation
data. As expected, the generated output match the high-res simulation on large
scales thanks to the conditional GAN training, while on small scales the
generated fine structures are different from that in the simulated ones but
look statistically consistent. For quantitative comparison, below we look at
various matter and halo statistics which provide a more solid and consistent
way to quantify model performances.

###  3.2 Statistics

There are numerous different statistical properties that can be looked at, but
most of them, eg. the power spectrum, don’t provide any new information, as
both the low-res input (Illustris-3) and high-res target (Illustris-2)
simulation are already at a quite high resolution and thus are non-
differentiable at this pixel-size scale. But, crucially, some statistics do
show clear differences, proving the performance of the trained network.

####  3.2.1 Power spectrum and 2PCF

Figure 4 depicts the particle 3D power spectrum and particle auto-correlation
function of the entire testing box, defined as in for example [33]. They are
compared to the power spectrum and two-point function calculated from the
"true" Illustris simulation particle catalogs over the same region. The power
spectrum $P(k)$ is calculated using a FFT over a cic-painted meshgrid whereas
the two-point function $\xi_{2}(x)$ is calculated by counting particle pairs
in binned distance-regions.

$\displaystyle\xi_{2}(x)$
$\displaystyle\equiv\left\langle\delta(\bm{x_{1}})\delta(\bm{x_{2}})\right\rangle,\quad
x=|\bm{x_{1}}-\bm{x_{2}}|$ (3.2) $\displaystyle\xi_{2}(x)$
$\displaystyle=\frac{1}{(2\pi)^{3}}\int d^{3}k\ P(k)e^{ikx},$ (3.3)

where the power $P(k)$ is simply defined as the Fourier transform of the two-
point correlation function and
$\delta(\bm{x_{i}})=\frac{\rho(\bm{x_{i}})}{\bar{\rho}}-1$ is the matter
overdensity. The effect of cic-compensation (as described in eq 18 of [34]),
which aims to reduce the statistical effect that painting particles onto a
mesh has on the power, is negligible and does not reduce the difference
between the pixelated and true power spectra.

Figure 4: left: The particle power spectrum comparison between the catalogs
obtained from Poisson sampling the training data and the "true" (unpixelized)
Illustris simulation catalogs. It is important to also compare to the "true"
Illustris statistics, as this will give an idea of the effects the
pixelization has on the low and high-res testing data. The divergence at small
scales due this pixelation becomes clearly visible. The low-res and and high-
res simulation’s power agree extremely well with each other up to the Nyquist
frequency, thus making the comparison futile as a way to prove the neural
net’s performance. right: The particle two-point correlation function up to
bin sizes of 2 Mpc/h, calculated using the Landy-Szalay estimator. As with the
power spectrum, the low and high resolution simulations coincide at scales
above the pixelsize of ${\sim}36.6$ kpc/h, only proving that the network can
at least mirror the input roughly in that regard. Small differences between
low-res input and generated output originate from the fact that the net does
not simply learn to copy the low-res input. The "true" Illustris two-point
correlations differ slightly from the testing data as an effect of the
pixelation.

For both, power and two-point function, the output agrees very well with the
low-res input / high-res target. The power differences at small physical
scales between testing data / generated output and the true simulation
catalogs, are caused by the pixelation and the ensuing loss of information
below pixel-size scales. The autocorrelation plot only covers a smaller range
of distances, thus doesn’t show this diversion. Since both the power and two-
point function of input and target are equal at these scales, it comes as no
surprise that the output will match these statistics reasonably well. It
should also be stated that even though the low-res input distribution matches
the high-res target’s in these two statistics, the generator does not produce
them exactly, as it does not simply learn to produce a one to one mapping
between input and target.

####  3.2.2 Halo mass function

All halos in figure 5 where calculated using nbodykit’s [35] FoF (Friends of
Friends) halo finder using a relative linking length of $0.2$. Again, the
pixelation effect becomes visible at smaller mass halos below and around
$10^{10}$ M☉, as the HMF found from pixelized data starts to deviate from the
halos found using the true simulation catalog. At smaller/larger pixel sizes
this effect will shift to smaller/larger halo masses.

Figure 5: The halo mass function obtained from FoF halos (linking length
$0.2$) inside the testing region. Illustris halos are taken from the "true"
high- and low-res $z=0$ simulation snapshots, which importantly show the
effect the randomization of particle positions inside each voxel has on FoF
halo finding. This effect becomes clearly visible for halos below $10^{10}$ M☉
as the output and target halo mass functions start to diverge from the "true"
Illustris-2 halos, whereas the "true" Illustris-3 halos still coincide with
the halos found from the low-res testing data as they are more massive and
thus unaffected by the pixelation.

####  3.2.3 Halo two-point function

The halo auto-correlation function is calculated by using the halo’s central
positions. Figure 6 and 7 shows the two point function of different halo mass
bins. Importantly for halos between $1\times 10^{10}$ M☉ and $4\times 10^{10}$
M☉, a clear difference between the low resolution and high resolution
simulation can be seen. The network manages to match the high resolution
simulation perfectly, proving it’s performance. It also matches the two-point
functions for higher mass halos, which is to be expected as those halos are
larger in size and easier to predict. The two-point function for halo masses
below $10^{10}$ M☉ becomes less meaningful, as the FoF halos found there will
be affected due to the pixelation. Nevertheless the two point function of the
(pixelized) generated output and high-res target still coincide with the two-
point function calculated from FoF halos found by using true Illustris-2
simulation particles.

Figure 6: The halo two-point function comparison between the testing data and
the "true" un-pixelated Illustris simulation FoF halos of the testing region,
calculated using the Landy-Szalay estimator. Only halos between $1\times
10^{10}$ M☉ and $6\times 10^{10}$ M☉ are plotted. The network correctly
predicts the high-res target’s (and also "true" Illustris-2) correlation of
halo positions, which importantly differs from the low-res input’s. This plot
is the most critical one, as it indicates that the method learns to predict
the high-resolution simulation and not simply up-scales the low-resolution
density fields. It also becomes clear from the comparison to the "true"
Illustris halo two-point correlation functions, that pixelation only has a
minor effect on the testing data in this regard. Figure 7: Halo two-point
functions of higher mass bins. Here low-res input and high-res target start to
coincide again so the net’s performance becomes less meaningful at these mass
scales.

###  3.3 Generalizing to another simulation

Ideally, the above model would be perfectly applicable to any low resolution
DM-only simulation’s density field of similar pixel size and find new halos.
The extend of this performance is tested in this section by using a set of the
Illustris-TNG DM-only simulations [36], which are the latest improvement on
the previous Illustris simulations. They use almost the same cosmology, but of
course are run using different initial conditions. By dividing both the
TNG-300-3-Dark (low-res) and TNG-300-2-Dark (high-res) simulations into
$5632^{3}$ density voxels, a similar pixelsize as before of ${\sim}36.6$ kpc/h
is fixed for the data. For memory efficiency, only an eighth of the entire
simulation volume of $205^{3}$ Mpc/h is used for testing again. Remarkably,
even though the particle masses are different compared to the Illustris
simulations (high-res TNG-300-2: $m_{\text{dm}}=5.6\times 10^{8}\
\textup{M}_{\Sun}$, low-res TNG-300-3: $m_{\text{dm}}=4.5\times 10^{9}\
\textup{M}_{\Sun}$), the network still manages to predict new halos
successfully which can again be seen in both the halo mass function and halo
two-point function in figure 8. This is an important indication that
convolutional neural networks are, once trained successfully, even applicable
outside their initial training and testing region. Of course once the pixel
size differs too much from the initial pixel size the network was trained on,
deteriorating results are to be expected.

Figure 8: left: The halo mass function of FoF halos of the TNG testing region,
similar to figure 5. Clearly the net performs weaker than on the Illustris
data, but still manages to predict new lower mass halos, compared to the low-
res input (TNG-300-3) simulation. Due to higher particle masses in the
simulations, only higher mass FoF halos, compared to the Illustris halo mass
function, are found. This also leads to the pixelation effect in figure 5 not
yet to take effect, as it only affects FoF halos below $10^{10}$ M☉ at this
pixelsize. right: Halo two-point functions of halos between $1\times 10^{11}$
M☉ and $4\times 10^{11}$ M☉. Here low-res input and high-res target differ
similarly to the Illustris simulations, but at higher halo masses, since the
simulation’s particle masses are different. The network makes the correct
prediction without ever training on the TNG data set, indicating that it is
applicable to a wide range of simulations. Again it becomes clear from the
comparison to the "true" TNG-300 halo two-point correlation functions, that
pixelation only has a minor effect on the testing data in this regard.

##  4 Conclusions

We apply GAN techniques to improve cosmological simulations to successfully
generate high resolution outputs from low resolution counterparts, mimicking
authentic high-res simulation results in all evaluated statistics (crucially
in figure 6).

This work pushes the resolution limit (from for example [18] and [37]), while
also being more generally applicable to general DM-only simulations than
displacement field training solutions, as it only relies on knowing the
particle positions at snapshot $z=0$. It provides an important first step to
building a full suite going from low-res dark matter simulations, to
ultimately producing ELG mock catalogs over Gpc in scale, by predicting new
halos at mass ranges relevant to ELGs.

In the case of simulations where initial conditions are set using a glass tile
configuration (like the Illustris simulations) it becomes almost impossible to
gain information about the dark matter displacement field, as the initial
particle positions cannot easily be placed on a grid together with their IDs.
This work allows to generate high-res density fields by only investing the CPU
time it takes to run the low-res simulation instead of actually running the
entire high-res simulation, as testing only took around 26 CPU hours for this
work. It is also naturally scaleable to much larger testing regions by simply
providing a larger low-res input.

A natural downside to using density fields for training is mainly the fact
that all information below pixel-size scales is lost naturally due to
pixelization. This poses natural limits to the lowest mass of physical halos
to be found when training on a given pixel-size. Nonetheless this work proves
that the method was successful down to scales where the low and high
resolution statistics start to deviate, and suggests that the same will also
hold up for training on smaller pixelsizes. This will also lead to a better
resolution of halos and subsequently also subhalos (identified using density
information), which will lead to the initial goal of creating large scale high
resolution ELG mocks from a low resolution simulation with almost no
additional computational effort.

##  5 Outlook

Future improvements of this method might make use of a smaller pixel-size, eg.
${\sim}18.3$ kpc/h, to improve smaller mass halo predictions and enable
subhalo predictions. In principle not much except for the learning rates will
have to be adjusted, although the amount of training data will increase
significantly, making the entire training process longer and more difficult to
tune. Additionally it will also be beneficial to train the network on multiple
redshifts, to make it more generally applicable to DM-only simulation output.
Recently similar work using displacement fields and super resolution methods
[16, 18] showed that training for higher redshift data is an easier task. This
also makes intuitively sense, as $z=0$ simulations went through the maximum
amount of non-linear structure development, making it harder to learn the
mapping accurately.

Another useful improvement to the training would be to add velocity fields to
the input, containing information about the particles’ bulk velocities inside
each pixel-volume, helping to infer halo field phase-space distributions. The
networks architecture supports training on multiple fields natively, so this
change could be easily executed.

It will also be interesting to provide the net with completely different
simulations, with different cosmological parameters and at similar mass
resolutions and pixel-scales during training, to gain a more universally
applicable network.

As we ultimately want to produce data comparable to observations, one will
also have to infer redshift-space distributions from our real-space data,
which will be subject of future work.

With the end goal in mind of creating a network that does not have to be
retrained and can easily enhance low resolution simulation density output,
this work provides an important foundation on the way to easily accessible
large scale ELG mock catalogs.

## Acknowledgements

We thank Yueying Ni for sharing her plotting scripts. This work was supported
in part through the NYU IT High Performance Computing resources, services, and
staff expertise. We would also like to thank AMD for the donation of critical
hardware and support resources from its HPC Fund, that also made this work
possible. The Flatiron Institute is supported by the Simons Foundation. JLT
acknowledges support of NSF grant number 2009291.

## Data Availability

All code and models used and created in this work are publicly available on
GitHub under the following url: https://github.com/dschaurecker/dl_halo.
Please refer to the ReadMe file in the dl_halo repository for a very in-depth
guide through the code and the process of training and testing. The training
code builds upon the map2map code repository, available on GitHub as well
(https://github.com/eelregit/map2map). It allows for general training of $n$
arbitrary input fields to $n$ arbitrary output fields using custom models,
normalizations and loss functions. Furthermore some utilities from nbodykit an
"Open-source, Massively Parallel Toolkit for Large-scale Structure" [35] were
used in pre-processing and statistical evaluation of this work.

## References

  * [1] DESI Collaboration, A. Aghamousa, J. Aguilar, S. Ahlen, S. Alam, L.E. Allen et al., _The desi experiment part i: Science,targeting, and survey design_ , 2016.
  * [2] R. Laureijs, J. Amiaux, S. Arduini, J.L. Auguères, J. Brinchmann, R. Cole et al., _Euclid definition study report_ , 2011.
  * [3] DES Collaboration, T.M.C. Abbott, M. Aguena, A. Alarcon, S. Allam, O. Alves et al., _Dark energy survey year 3 results: Cosmological constraints from galaxy clustering and weak lensing_ , 2021.
  * [4] Z. Ivezic, S.M. Kahn, J.A. Tyson, B. Abel, E. Acosta, R. Allsman et al., _Lsst: From science drivers to reference design and anticipated data products_ , _The Astrophysical Journal_ 873 (2019) 111.
  * [5] K.S. Dawson, D.J. Schlegel, C.P. Ahn, S.F. Anderson, . Aubourg, S. Bailey et al., _The baryon oscillation spectroscopic survey of sdss-iii_ , _The Astronomical Journal_ 145 (2012) 10.
  * [6] A. Raichoor, J. Comparat, T. Delubac, J.-P. Kneib, C. Yèche, K.S. Dawson et al., _The sdss-iv extended baryon oscillation spectroscopic survey: final emission line galaxy target selection_ , _Monthly Notices of the Royal Astronomical Society_ 471 (2017) 3955–3973.
  * [7] H. Guo, X. Yang, A. Raichoor, Z. Zheng, J. Comparat, V. Gonzalez-Perez et al., _Evolution of star-forming galaxies from z = 0.7 to 1.2 with eBOSS emission-line galaxies_ , _The Astrophysical Journal_ 871 (2019) 147.
  * [8] V. Gonzalez-Perez, J. Comparat, P. Norberg, C.M. Baugh, S. Contreras, C. Lacey et al., _The host dark matter haloes of [o ii] emitters at 0.5 < z < 1.5_, _Monthly Notices of the Royal Astronomical Society_ 474 (2017) 4024–4038.
  * [9] J.H.T. Yip, X. Zhang, Y. Wang, W. Zhang, Y. Sun, G. Contardo et al., _From dark matter to galaxies with convolutional neural networks_ , 2019.
  * [10] N. Kasmanoff, F. Villaescusa-Navarro, J. Tinker and S. Ho, _dm2gal: Mapping dark matter to galaxies with neural networks_ , 2020.
  * [11] F. Villaescusa-Navarro, B.D. Wandelt, D. Anglés-Alcázar, S. Genel, J.M.Z. Mantilla, S. Ho et al., _Neural networks as optimal estimators to marginalize over baryonic effects_ , 2020.
  * [12] S. He, Y. Li, Y. Feng, S. Ho, S. Ravanbakhsh, W. Chen et al., _Learning to predict the cosmological structure formation_ , _Proceedings of the National Academy of Sciences of the United States of America_ 116 (2019) 13825 [1811.06533].
  * [13] R.A. de Oliveira, Y. Li, F. Villaescusa-Navarro, S. Ho and D.N. Spergel, _Fast and accurate non-linear predictions of universes with deep learning_ , 2020.
  * [14] D. Kodi Ramanah, T. Charnock, F. Villaescusa-Navarro and B.D. Wandelt, _Super-resolution emulator of cosmological simulations using deep physical models_ , _Monthly Notices of the Royal Astronomical Society_ 495 (2020) 4227–4236.
  * [15] R.H. Wechsler and J.L. Tinker, _The connection between galaxies and their dark matter halos_ , _Annual Review of Astronomy and Astrophysics_ 56 (2018) 435–487.
  * [16] Y. Li, Y. Ni, R.A.C. Croft, T. Di Matteo, S. Bird and Y. Feng, _Ai-assisted superresolution cosmological simulations_ , _Proceedings of the National Academy of Sciences_ 118 (2021)  [https://www.pnas.org/content/118/19/e2022038118.full.pdf].
  * [17] S.D.M. White, _Formation and Evolution of Galaxies_ , in _Cosmology and Large Scale Structure_ , R. Schaeffer, J. Silk, M. Spiro and J. Zinn-Justin, eds., p. 349, Jan., 1996.
  * [18] Y. Ni, Y. Li, P. Lachance, R.A.C. Croft, T.D. Matteo, S. Bird et al., _Ai-assisted super-resolution cosmological simulations ii: Halo substructures, velocities and higher order statistics_ , 2021.
  * [19] M. Vogelsberger, S. Genel, V. Springel, P. Torrey, D. Sijacki, D. Xu et al., _Introducing the illustris project: Simulating the coevolution of dark and visible matter in the universe_ , _Monthly Notices of the Royal Astronomical Society_ 444 (2014) 1518 [1405.2921].
  * [20] C.L. Bennett, D. Larson, J.L. Weiland, N. Jarosik, G. Hinshaw, N. Odegard et al., _Nine-year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Final Maps and Results_ , _The Astrophysical Journal Supplement_ 208 (2013) 20 [1212.5225].
  * [21] M. Vogelsberger, S. Genel, V. Springel, P. Torrey, D. Sijacki, D. Xu et al., _Properties of galaxies reproduced by a hydrodynamic simulation_ , _Nature_ 509 (2014) 177 [1405.1418].
  * [22] I.J. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair et al., _Generative Adversarial Networks_ , _Proceedings - 2019 International Conference on Computer Vision Workshop, ICCVW 2019_ (2014) 3063 [1406.2661].
  * [23] J.F. Nash, _Equilibrium points in n-person games_ , _Proceedings of the National Academy of Sciences_ 36 (1950) 48 [https://www.pnas.org/content/36/1/48.full.pdf].
  * [24] M. Mirza and S. Osindero, _Conditional generative adversarial nets_ , 2014\.
  * [25] L. Mescheder, A. Geiger and S. Nowozin, _Which training methods for GANs do actually converge?_ , _35th International Conference on Machine Learning, ICML 2018_ 8 (2018) 5589 [1801.04406].
  * [26] O. Ronneberger, P. Fischer and T. Brox, _U-net: Convolutional networks for biomedical image segmentation_ , _Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics)_ 9351 (2015) 234 [1505.04597].
  * [27] E. Giusarma, M. Reyes Hurtado, F. Villaescusa-Navarro, S. He, S. Ho and C. Hahn, _Learning neutrino effects in Cosmology with Convolutional Neural Networks_ , _arXiv e-prints_ (2019) arXiv:1910.04255 [1910.04255].
  * [28] C. Chen, Y. Li, F. Villaescusa-Navarro, S. Ho and A. Pullen, _Learning the evolution of the universe in n-body simulations_ , 2020.
  * [29] E. Shelhamer, J. Long and T. Darrell, _Fully Convolutional Networks for Semantic Segmentation_ , _IEEE Transactions on Pattern Analysis and Machine Intelligence_ 39 (2017) 640 [1411.4038].
  * [30] A. Odena, V. Dumoulin and C. Olah, _Deconvolution and checkerboard artifacts_ , _Distill_ (2016) .
  * [31] I. Loshchilov and F. Hutter, _Decoupled weight decay regularization_ , 2019\.
  * [32] M. Heusel, H. Ramsauer, T. Unterthiner, B. Nessler and S. Hochreiter, _Gans trained by a two time-scale update rule converge to a local nash equilibrium_ , 2018.
  * [33] H. Mo, F. van den Bosch and S. White, _Probing the cosmic density field_ , in _Galaxy Formation and Evolution_ , pp. 262–318, Cambridge University Press (2010), DOI.
  * [34] Y.P. Jing, _Correcting for the alias effect when measuring the power spectrum using a fast fourier transform_ , _The Astrophysical Journal_ 620 (2005) 559–563.
  * [35] N. Hand, Y. Feng, F. Beutler, Y. Li, C. Modi, U. Seljak et al., _nbodykit: An open-source, massively parallel toolkit for large-scale structure_ , _The Astronomical Journal_ 156 (2018) 160.
  * [36] V. Springel, R. Pakmor, A. Pillepich, R. Weinberger, D. Nelson, L. Hernquist et al., _First results from the illustristng simulations: matter and galaxy clustering_ , _Monthly Notices of the Royal Astronomical Society_ 475 (2017) 676–698.
  * [37] D. Kodi Ramanah, T. Charnock and G. Lavaux, _Painting halos from cosmic density fields of dark matter with physically motivated neural networks_ , _Physical Review D_ 100 (2019) .


# Exceptional Points of Degeneracy with Indirect Bandgap Induced By Mixing
Forward and Backward Propagating Waves

Tarek Mealy, and Filippo Capolino Department of Electrical Engineering and
Computer Science, University of California, Irvine, CA 92697, USA
<EMAIL_ADDRESS>and<EMAIL_ADDRESS>
###### Abstract

We demonstrate that exceptional points of degeneracy (EPDs) are obtained in
two coupled waveguides without resorting to gain and loss. We show the general
concept that modes resulting from a proper coupling of forward and backward
waves exhibit EPDs of order two and that there the group velocity vanishes. We
verify our insight by using coupled mode theory and also by fullwave numerical
simulations of light in a dielectric slab coupled to a grating, when one
supports a forward wave whereas the other (the grating) supports a backward
wave. We also demonstrate how to realize photonic indirect bandgaps in guiding
systems supporting a backward and a forward wave, show its relations to the
occurrence of EPDs, and offer a design procedure.

###### Index Terms:

Exceptional point of degeneracy; Coupled mode theory; Parity time symmetry;
Degenerate band edge.

##  I Intruduction

An EPD is a point in the parameter space of a system at which the system’s
eigenvalues and eigenvectors coalesce [1, 2, 3, 4]. The term exceptional point
(EP) and the associated perturbation theory where discussed in the well known
Kato’s book in 1966 [4]. The phenomenon of degeneracy of both eigenvalues and
eigenvectors (polarization states), studied here, is a stronger degeneracy
compared to the traditional degeneracy of only two eigenvalues.

Non-Hermitian Hamiltonian can possess entirely real spectra when the system
obeys parity-time (PT) symmetry condition [5]. A system is said to be PT
symmetric if the PT operator commutes with the Hamiltonian [6, 7], where PT
operator applies a parity reflection and time reversal [5]. When the time
reversal operator is applied to physical systems, energy changes from damping
to growing and vice versa [8]. Based on this simple concept, two symmetrical
coupled waveguides with balanced gain and loss satisfy PT symmetry [9, 10,
11], where the individual application of each of the space or time reversal
would swap the gain and loss, therefore the simultaneous application of space
and time reversal operator to the system would end up with the same system.
The point separating the complex and real spectra regimes of PT-symmetric
Hamiltonians has been called exceptional point (EP) [4], also known as
transition point. Here, beside the mathematical aspects, we stress the role of
degeneracy, as implied also in [12], hence include the ’D’ in the EPD acronym.

In this paper, we present a class of two coupled waveguides where EPDs exist
without resorting to the presence of gain and loss. By using coupled mode
theory [13, 14], we show that two coupled waveguides, where one waveguide
supports forward propagation (i.e., where the phase and power propagate in the
same direction) whereas the other one supports backward propagation (i.e.,
where the phase and power propagate in opposite direction), experience a phase
transition as in the PT-symmetric case. We show the general conditions for
modes resulting from coupling two coupled waves to exhibit an EPD looking at
both the degenerate eigenvalues and eigenvectors. We show that the coupling of
two waves, that carry power in opposite directions, leads to an EPD and we
explain how this results in the vanishing of the group velocity of the
degenerate mode. We illustrate the concepts in a simple system made of two
coupled waveguides, i.e., a dielectric slab coupled to a grating, when one
supports a forward wave whereas the other (the grating) supports a backward
wave. Other general conditions that lead to exceptional degeneracies of two
modes in uniform waveguide were studied in [15] using a transmission line
approach. Finally, we relate the occurrence of EPDs to the presence of a
photonic indirect bandgap.

Figure 1: (a) Coupling between two electromagnetic waves whose complex
amplitudes are $A$ and $B$. Conditions that lead to EPDs are obtained by
introducing proper coupling between: (b) two waveguides with PT symmetry where
the two media have gain and loss supporting exponentially growing and
attenuating waves; (c) two waveguides with forward and backward propagating
waves, without resorting to PT symmetry (i.e., in this case the waveguides do
not have gain and loss). The waves with black arrows represent the directions
of propagation. The blue and red arrows represent the directions of power
flow.

##  II Second order EPD by mixing two waves

We consider two coupled electromagnetic waves as shown in Fig. 1(a). These two
waves are described by the complex time-domain notation

$\begin{array}[]{c}a(\mathbf{r},t)=A(z)f_{a}(\mathbf{\boldsymbol{\rho}})e^{i\omega
t},\\\ b(\mathbf{r},t)=B(z)f_{b}(\mathbf{\boldsymbol{\rho}})e^{i\omega
t},\end{array}$ (1)

where $A(z)$ and $B(z)$ are the complex amplitudes of waves along the $z$
direction, $\mathbf{r}=\boldsymbol{\rho}+z\hat{\mathbf{z}}$, and
$\boldsymbol{\rho}$ is the transverese coordinate.
$f_{a}(\mathbf{\boldsymbol{\rho}})$ and $f_{b}(\mathbf{\boldsymbol{\rho}})$
are the normalized modal field profile in the transverse direction for each
mode. When the two waves are uncoupled, i..e., when their waveguides are far
from from each other, the evolution of the amplitudes along the $z$ direction
is simply descried by $dA(z)/dz=-i\beta_{a}^{\prime}A(z)$ and
$dB(z)/dz=-i\beta_{b}^{\prime}B(z)$, where $\beta_{a}^{\prime}$ and
$\beta_{b}^{\prime}$ are the uncoupled propagation constants of each wave
(that is also an eigenmode of the structure since there is no coupling). The
solutions are $A(z)=A_{0}\exp(-i\beta_{a}^{\prime}z)$ and
$B(z)=B_{0}\exp(-i\beta_{a}^{\prime}z)$. The power carried by each wave in the
positive $z$ direction is given by $p_{a}(z)=\pm|A(z)|^{2}$ and
$p_{b}(z)=\pm|B(z)|^{2}$ where the sign depends on the type of the wave. The
sign is positive when the wave is forward, i.e., when the power is carried in
the same direction of wave propagation (i.e., when the phase and group
velocity have the same directions). The sign is negative when the waves
backward, i.e., when the phase propagates along the positive $z$ direction
whereas the power flows in the negative $z$ direction (i.e., when the phase
and group velocity have the opposite directions).

When coupling is introduced to those two waves, the system eigenmode is found
by solving the spatial-evolution equation that, based on coupled mode theory
[13, 14], is given by

$\dfrac{d}{dz}\left(\begin{array}[]{c}A(z)\\\
B(z)\end{array}\right)=-i\left(\begin{array}[]{cc}\beta_{a}&\kappa_{ab}\\\
\kappa_{ba}&\beta_{b}\end{array}\right)\left(\begin{array}[]{c}A(z)\\\
B(z)\end{array}\right),$ (2)

where $\beta_{a}$ and $\beta_{b}$ are “perturbed” propagation constants for
the coupled system, and $\kappa_{ab}$ and $\kappa_{ba}$ are the coupling
coefficients between the two modes. The relation between $\kappa_{ab}$ and
$\kappa_{ba}$ is determined by applying the power conservation principle. The
total power carried in the coupled structure is
$p_{t}(z)=|A(z)|^{2}\pm|B(z)|^{2}$ assuming wave $A$ is forward wave, and wave
$B$ to be either forward or backward, when taking the $+$ or the $-$ sign,
respectively. Thus, there are two possible scenarios: (i) “codirectional
coupling” when both waves carry power at the same direction and (ii)
“contradirectional coupling” when the two waves carry power in the opposite
direction [13].

When the system does not have gain and loss, conservation of energy states
that $dp_{t}(z)/dz=0$, and by using (2), one finds that the constraint
$\mathrm{Re}\left(AB^{*}\left(\kappa_{ba}\mp\kappa_{ab}^{*}\right)\right)=0$
should be satisfied. Therefore, we have $\kappa_{ab}=\kappa_{ba}^{*}$ in case
of codirectional coupling where the two waves are forward, and
$\kappa_{ab}=-\kappa_{ba}^{*}$ in case of contradirectional coupling where one
wave is a forward and the other one is backward [13].

The mixing of the two waves constitutes what is called the guiding system’s
eigenmode (some call it “supermode”) which is a weighted sum of the individual
guided waves. The eigenmode propagation constant is determined by solving the
characteristic equation of the coupled system in (2) assuming the wave
amplitudes to be in the form of $[A(z),B(z)]^{T}\propto e^{-ikz}$ which yields
$k^{2}-k(\beta_{a}+\beta_{b})+\left(\beta_{a}\beta_{b}-\kappa_{ab}\kappa_{ba}\right)=0.$
The characteristic equation has two solutions that are given by

$\begin{array}[]{c}k_{n}=\dfrac{\beta_{a}+\beta_{b}}{2}+(-1)^{n}\sqrt{\left(\dfrac{\beta_{a}-\beta_{b}}{2}\right)^{2}-\left(-1\right)^{p}\kappa^{2}},\end{array}$
(3)

where $\kappa=|\kappa_{ab}|$ and the indices $n=1,2$ denote the two modes of
the coupled system. Furthermore, $p=1$ and $p=2$ represent the case of
codirectional and contradirectional coupling, respectively. An EPD occurs when
two eigenmodes coalesce, i.e., $k_{1}=k_{2}=k_{e},$ with
$k_{e}=\left(\beta_{a}+\beta_{b}\right)/2$. This EPD occurs when
$\beta_{a}-\beta_{b}=2\sqrt{\left(-1\right)^{p}\kappa^{2}}$. At an EPD, the
eigenvectors must coalesce, and in this simple system their coalescence
follows from the coalescence of the eigenvalues. Indeed, the two eigenvectors
are
$\left[\begin{array}[]{cc}A_{n},&B_{n}\end{array}\right]^{T}=\left[\begin{array}[]{cc}1,&(k_{n}-\beta_{a})/\kappa_{ab}\end{array}\right]^{T}$and
it is easy to see that they coalesce when $k_{1}=k_{2}$.

The group velocity of the eigenmode with wavenumber $k_{n}$ is determined as
(assuming $k_{n}$ to be purley real)

$v_{g,n}=\frac{1}{d_{\omega}k_{n}}=\dfrac{k_{n}-k_{e}}{2k_{n}d_{\omega}\left(\beta_{a}+\beta_{b}\right)+2d_{\omega}\left(\beta_{a}\beta_{b}+\left(-1\right)^{p}\kappa^{2}\right)}.$
(4)

where $d_{\omega}\equiv d/d\omega$ denotes the derivative with respect to
angular frequency $\omega$. It is clear from the expression that $v_{g,1,2}=0$
when $k_{1}=k_{2}=k_{e},$i.e., exactly at the EPD. Next, we also show what
happens near the EPD.

Figure 2: The two wavenumbers of the guiding system versus the coupling
parameter $\kappa$, showing the existence of an EPD. Two cases are examined:
(a) coalescence of modes in PT-symmetrical waveguides, (b) coalescence of
modes obtained by coupling a forward wave (phase and group velocities have the
same direction) and a backward wave (phase and group velocities have opposite
directions). Both cases exhibit an EPD, represented by the bifurcation point.

###  II-A Codirectional coupling $(\kappa_{ab}=\kappa_{ba}^{*})$

For the case of codirectional coupling , $p=1,$ the EPD condition is
simplified to $\beta_{a}-\beta_{b}=\pm 2i\kappa$. The EPD condition puts a
constraint that the difference between the propagation constants of the
uncoupled waves has to be purely imaginary in order to exhibit an EPD. Thus,
we conclude that the EPD can never be obtained for any value of the coupling
parameter $\kappa$ in the case of a lossless/gainless system. If we resort to
a PT symmetry, as in Fig. 1(b), where the system has balanced gain and loss,
we have $\beta_{a}=\beta_{0}+i\alpha$ and $\beta_{b}=\beta_{0}-i\alpha$, and
an EPD is obtained when $\alpha=\kappa$ [9, 10, 11, 16, 17], and the
degenerate wavenumber is $k_{e}=\beta_{0}$.

For this case, the two propagation constants of the coupled system in the
vicinity of the EPD are $k_{1,2}=\beta_{0}\pm\sqrt{\kappa^{2}-\alpha^{2}}$ and
their derivatives are $d_{\omega}k_{1,2}=d_{\omega}\beta_{0}\pm\left(2\kappa
d_{\omega}\kappa-2\alpha d_{\omega}\alpha\right)/(k_{2}-k_{1})$. This is also
illustrated by determining the eigenvector of the degenerate eigenmode from
(2) (for codirectional coupling case where $p=1$) as

$\begin{array}[]{c}\left(\begin{array}[]{c}A_{e}(z)\\\
B_{e}(z)\end{array}\right)=\left(\begin{array}[]{c}1\\\
-ie^{-i\mathrm{arg}(\kappa_{ab})}\end{array}\right)e^{-ik_{e}z},\end{array}$
(5)

and one finds that the total power carried by the degenerate eigenmode is
$p_{t}(z)=|A_{e}(z)|^{2}-|B_{e}(z)|^{2}=0$ vanishes, in agreement with the
vanishing of the group velocity. In the vicinity of an EPD we have
$k_{1}\approx k_{2}$ and by neglecting the $d_{\omega}\beta_{0}$, usually
smaller than the other term, the group velocities of the two modes are
$v_{g,1,2}\approx\pm\left(k_{2}-k_{1}\right)/\left(2\kappa
d_{\omega}\kappa-2\alpha d_{\omega}\alpha\right)$ when $\kappa>\kappa_{e}$,
i.e., where the two eigenmodes are propagating with purley real wavenumbers.
Therefore we conclude that near an EPD in a PT-symmetric guiding system, the
two modes of the coupled system are phase-synchronized ($k_{1}\approx k_{2}$),
i.e., with almost identical phase velocity, but have opposite group velocity
($v_{g,1}\approx-v_{g,2}$), which eventually results on having a wave in the
guiding structure with vanishing group velocity when the system is exactly at
the EPD. As an example, for a system with PT symmetry where the uncoupled
waveguides have, respectively, a growing wave with $\beta_{a}=100+10i$ (1/m)
and an attenuating wave with $\beta_{b}=100-10i$ (1/m), an EPD is obtained
when the coupling parameter is $\kappa=10$ 1/m as shown in Fig. 2(a).

###  II-B Contradirectional coupling $(\kappa_{ab}=-\kappa_{ba}^{*})$

We consider coupling between forward wave with wavenumber $\beta_{a}>0$ and
$v_{g,a}>0$ ($v_{g,a}=d_{\omega}\beta_{a}$) and backward wave with wavenumber
$\beta_{b}>0$ and $v_{g,b}<0$ ($v_{g,b}=d_{\omega}\beta_{b}$). The two
propagating waves carry power in opposite directions and therefore they
exhibit contradirectional coupling, we use $p=2$ in Eq. (3). The EPD condition
($k_{1}=k_{2}$) for this case is simplified to $\beta_{a}-\beta_{b}=\pm
2\kappa$, which means that the difference between the propagation constants
should be purely real to have an EPD, which is possible for a
lossless/gainless system. Therefore, the EPD condition for this case is
satisfied through the proper design of the coupling parameters, i.e, when
$\kappa=\left|\beta_{a}-\beta_{b}\right|/2$ (we recall that $\kappa$ was
defined as purely real positive). This means that there are two possible EPD
conditions, $\beta_{a}-\beta_{b}=2\kappa$ and $\beta_{a}-\beta_{b}=-2\kappa$,
that may both occur when varying frequency. At those two frequencies one has
$\beta_{a}>\beta_{b}$ and $\beta_{b}>\beta_{a}$, respectively. A more detailed
discussion is provided later on when discussing the indirect bandgap.

In the vicinity of an EPD, the two propagation constants of the coupled system
are given by Eq. (3), and the derivatives of the two wavenumbers with respect
to the angular frequency are

$d_{\omega}k_{1,2}=\frac{1}{2}d_{\omega}\left(\beta_{a}+\beta_{b}\right)\pm\frac{1}{2}\dfrac{\left(\beta_{a}-\beta_{b}\right)d_{\omega}\left(\beta_{a}-\beta_{b}\right)-2\kappa
d_{\omega}\kappa}{(k_{2}-k_{1})}.$ (6)

When $\kappa<\kappa_{e}$, i.e., where the two eigenmodes are propagating with
purely real wavenumbers, in the vicinity of an EPD we have $k_{1}\approx
k_{2}$ and by neglecting the term $d_{\omega}\left(\beta_{a}+\beta_{b}\right)$
with respect to the second one, the group velocities of the two modes are

$v_{g,1,2}\approx\pm\frac{k_{2}-k_{1}}{\frac{1}{2}\left(\beta_{a}-\beta_{b}\right)\left(v_{g,a}^{-1}-v_{g,b}^{-1}\right)-2\kappa
d_{\omega}\kappa}.$ (7)

Therefore, we conclude that near an EPD, the coupled forward and backward
waves are syncronised in phase, i.e., $k_{1}\approx k_{2}$ but have nearly
opposite group velocity ($v_{g,1}\approx-v_{g,2}$), which eventually results
in having a wave in the guiding structure with vanishing group velocity when
the system is exactly at an EPD.

This is also illustrated by determining the eigenvector of the degenerate
eigenmode from (2) (for contradirectional coupling case where $p=2$) as

$\begin{array}[]{c}\left(\begin{array}[]{c}A_{e}(z)\\\
B_{e}(z)\end{array}\right)=\left(\begin{array}[]{c}1\\\
-e^{-i\mathrm{arg}(\kappa_{ab})}\end{array}\right)e^{-ik_{e}z},\end{array}$
(8)

and one finds that the total power carried by the degenerate eigenmode is
$P_{t}(z)=|A_{e}(z)|^{2}-|B_{e}(z)|^{2}=0$ vanishes, in agreement with the
vanishing of the group velocity. At the EPD, the system matrix is not
diagonalizable but rather similar to a $2\times 2$ Jordan matrix. The fields
in the two-waveguide system is represented using the degenerate and
generalized eigenvectors as

$\begin{array}[]{c}\left(\begin{array}[]{c}A(z)\\\
B(z)\end{array}\right)=\left(\begin{array}[]{c}1\\\
-e^{-i\mathrm{arg}(\kappa_{ab})}\end{array}\right)\left(u_{1}-jzu_{2}\right)e^{-ik_{e}z}\\\
\ \ \ \ \ \ \ \ \ \ \ \ \ \ +\left(\begin{array}[]{c}1\\\
ie^{-i\mathrm{arg}(\kappa_{ab})}/\kappa_{e}\end{array}\right)u_{2}e^{-ik_{e}z}.\end{array}$
(9)

where $u_{1}$ and $u_{2}$ are proper coefficients that depend on the system
excitation and boundary conditions.

The two waveguides with contradirectional coupling are schematically shown in
Fig. 1(c) and the dispersion diagram is in Fig. 2 where we see a forward wave
with $\beta_{a}=110$ (1/m) and backward wave with $\beta_{b}=90$ (1/m), and
the EPD is obtained at $k=k_{e}=100$ (1/m) , where $\kappa_{e}=10$ (1/m) as
shown in Fig. 2(b).

The contradirectional coupling case can be realized is in two possible
scenarios: (i) two modes exist in two separate waveguides where the first
waveguide supports a forward wave and the second waveguide supports a backward
wave and the coupling is introduced by bringing them near each other; (ii) two
modes exist in the same waveguide having periodicity where the one wave (e.g.,
the forward) has the fundamental Floquet harmonic equal to
$\beta_{1}=\beta_{0}$ and the other wave (e.g., the backward) has its 1st
harmonic Floquet harmonic equal to $\beta_{2}=-\beta_{0}+2\pi/d$, where $d$ is
the waveguide period, and the EPD is only possible at the band edge
$\beta_{0}=\pi/d$. An example belonging to the first scenario, where the EPD
is found in two coupled dielectric slab waveguides, is shown later on. The
second scenario instead exists in conventional periodic waveguides and it is
not further considered in this paper.

Figure 3: Example of modal EPD of order 2 between a forward wave (phase and
group velocities have the same direction) and backward wave (phase and group
velocities have opposite directions). (a) Two coupled Si layers where the top
one supports a forward wave (in red) while the bottom one is periodically
corrugated to support a backward wave (positive phase velocity and negative
group velocity, in blue). The red and blue arrows represent the direction of
power flow. (b) Dispersion relation showing the propagating eigenmodes when
the two waveguides are uncoupled (dashed) and when the are coupled (solid).
The dispersion of modes in the coupled waveguides show the existence of two
EPDs. The blue and red colors of the curves are related to the power flow
directions. Note also that because two EPDs are found, an indirect bandgap is
present between the upper and lower branches that can be designed ad-hoc. In
the shown case we have $\Delta k_{e}\equiv k_{e2}-k_{e1}<0$.

We present an example of a guiding system that supports two waves carrying
power in opposite directions and we show that it exhibits two EPDs. Consider
the guiding system made of a Si substrate (supporting the forward wave)
coupled to a Si grating waveguide as shown in Fig. 3, with dimensions
$w=p=h=70$ nm and $d=140$ nm. Silcon is modeled with a refractive index
$n_{\mathrm{Si}}=3.45$. We first show in Fig. 3 the dispersion of the two
wavenumbers $\beta_{a}^{\prime}$ and $\beta_{b}^{\prime}$ of the two uncoupled
waveguides ($s\to\infty$) as red dashed (uniform waveguide) and blue dashed
(grating waveguide). The figure shows that the structures support a forward
wave where the group velocity is positive
$v^{\prime}_{g,a}=1/(d_{\omega}\beta_{a}^{\prime})>0$ and a backward wave
where the group velocity is negative
$v^{\prime}_{g,b}=1/(d_{\omega}\beta_{b}^{\prime})<0$. In the same Fig. 3, we
show the dispersion of the two wavenumbers $k_{1}$ and $k_{2}$ of the coupled
guiding system, i.e., when the two waveguide are close to each other with a
gap of $s=70$ nm. The dispersion show the existence of two EPDs associated to
a wavenumber (momentum) displacement. The dispersion diagrams we show are for
modes with electric field polarized in the $y$ direction. The dispersion
diagrams have been found by using the finite element method-based eigenmode
solver implemented in CST Studio Suite, by numerically simulating only one
unit cell of the structure. The proposed condition allows to locate the EPDs
at band edges that are not necessarily at the center or at the edge of a
Brillioum zone, without using loss and gain. It also shows the capability to
engineer an indirect bandgap in these simple structures.

##  III Indirect bandgap in the contradirectional case

Figure 4: Schematic of a dispersion diagram showing the indirect bandgap that
results from two EPDs based on contradirectional coupling. The red-dashed line
represents the wavenumber of the forward wave $\beta_{a}^{\prime}$ whereas the
blue-dahsed line the one of the backward wave $\beta_{b}^{\prime}$, when the
waveguides are uncoupled. The coupling yields the two curves with two EPDs
that are labeled as EPD1 and EPD2. The indirect bandgap width is
$\Delta\omega_{IB}$. In the shown case $\Delta k_{e}\equiv k_{e2}-k_{e1}>0$.

In the contradirectional case where coupling occurs between a forward and
backward wave, as in Fig. 1(c), an indirect bandgap is possible and we show
here how it is formed. Since in this case one wave is forward and one is
backward, the uncoupled propagation constants $\beta_{a}^{\prime}$ and
$\beta_{b}^{\prime}$ have opposite slopes as schematically shown in Fig. 4
(see also dashed blue and red curves in Fig. 3), and an analogous trend is
expected for the parameters $\beta_{a}$ and $\beta_{b}$ of the coupled system.
By looking at the dispersion diagram shown in Fig. 4, the red-dashed curve is
the forward wave with wavenumber $\beta_{a}^{\prime}(\omega)$ and the blue-
dashed curve is the backward wave with wavenumber
$\beta_{b}^{\prime}(\omega)$. Assuming that the coupling it not so strong, one
may assume that $\beta_{a}(\omega)\approx\beta_{a}^{\prime}(\omega)$ and
$\beta_{b}(\omega)\approx\beta_{b}^{\prime}(\omega)$, at least in trend,
hence, in slope. The forward wave has positive slope,
$v_{g,a}=d\omega/d\beta_{a}>0$, whereas the backward wave has negative slope
$v_{g,b}=d\omega/d\beta_{b}<0$, and the dispersion curves for $\beta_{a}$ and
$\beta_{b}$ versus frequency should intersect at some frequency
($\beta_{a}=\beta_{b}$) because one wavenumber is increasing with frequency
whereas the other is decreasing with frequency; an example with a grating is
illustrated in Fig. 3 while a schemtic is in Fig. 4. Let us approximate the
dispersion curves locally, in the frequency range of interest, as straight
lines, i.e., $\beta_{a}(\omega)\approx a+v_{g,a}^{-1}\omega$ and
$\beta_{b}(\omega)\approx b+v_{g,b}^{-1}\omega$, where now $v_{g,a}$ and
$v_{g,b}$ are assumed to have the local fixed value. Furthermore, assuming,
for simplicity, that $\kappa$ is constant within the frequency range of
interest, one finds that EPDs occurs at two angular frequencies $\omega_{e1}$
and $\omega_{e2}$ such that
$\beta_{a}(\omega_{e1})-\beta_{b}(\omega_{e1})=-2\kappa$ and
$\beta_{a}(\omega_{e2})-\beta_{b}(\omega_{e2})=2\kappa$, where
$\omega_{e1}<\omega_{e2}$. Subtracting the previous two conditions leads to
the indirect bandgap determination

$\Delta\omega_{IB}\equiv\omega_{e2}-\omega_{e1}\approx\frac{4\kappa}{v_{g,a}^{-1}-v_{g,b}^{-1}}.$
(10)

Note that $v_{g,a}^{-1}-v_{g,b}^{-1}>0$, hence $\Delta\omega_{IB}>0$. The
bandgap width can be controlled by the slope of the two parameters $\beta_{a}$
and $\beta_{b}$; indeed, when $v_{g,a}^{-1}\approx-v_{g,b}^{-1}$, the
denominator of (10) is small and the bandgap is very wide, viceversa, the
bandgap is narrow when $v_{g,a}^{-1}$ is very different from $-v_{g,b}^{-1}$.
If we consider the dispersion of the coupling term $\kappa$, a more
complicated picture may arise that could be determined by the reasoning just
provided.

The degenerate wavenumbers at the two EPDs are
$k_{e,1}=\left(\beta_{a}(\omega_{e1})+\beta_{b}(\omega_{e1})\right)/2$ and
$k_{e,2}=\left(\beta_{a}(\omega_{e2})+\beta_{b}(\omega_{e2})\right)/2$. Using
the linear approximation formulas for the wavenumbers $\beta_{a}(\omega)$ and
$\beta_{a}(\omega)$, one finds that
$k_{e1}\approx[a+b+\left(v_{g,a}^{-1}+v_{g,b}^{-1}\right)\omega_{e1}]/2$ and
$k_{e2}\approx[a+b+\left(v_{g,a}^{-1}+v_{g,b}^{-1}\right)\omega_{e2}]/2$ . The
difference between the two degenerate wavenumbers is

$\Delta k_{e}\equiv
k_{e2}-k_{e1}\approx\frac{v_{g,a}^{-1}+v_{g,b}^{-1}}{2}\Delta\omega_{IB}.$
(11)

Therefore, it is necessary that $v_{g,a}^{-1}+v_{g,b}^{-1}\neq 0$ in order to
have indirect bandgap. When $|v_{g,a}^{-1}|>|v_{g,b}^{-1}|$ we get
$k_{e2}>k_{e1},$hence $\Delta k_{e}>0$, i.e., the EPD that occurs at the
smaller frequency $\omega_{e1}$ occurs also at the smaller degenerate
wavenumber $k_{e1}$ and this condition is depicted in Fig. 4. When
$|v_{g,a}^{-1}|<|v_{g,b}^{-1}|$, we get $k_{e1}>k_{e2},$ hence $\Delta
k_{e}<0$, i.e., the EPD that occurs at the smaller frequency $\omega_{e1}$
occurs at larger degenerate wavenumber $k_{e1}$ and this condition is depicted
in Fig. 3. Indeed, by looking at Fig. 3, one finds by naked eye that
$|v_{g,a}|>|v_{g,b}|$ , (absolute change with frequency of dashed red curve is
higher than the one of the dashed blue one), therefore
$|v_{g,a}^{-1}|<|v_{g,b}^{-1}|$ resulting in $k_{e1}>k_{e2}$ according to (11)
and the EPD at lower frequency, around 420 THz, occurs at higher degenerate
wavenumber of $k_{e1}=0.972\,\pi/d$ .

##  IV Conclusion

We have demonstrated that EPDs are not only obtained in PT-symmetric
waveguides but they are also obtained in two lossless and gainless waveguides
when they support forward and backward waves that are properly coupled. We
have shown a simple system that supports this condition made of a grating
(supporting a backward guided mode) coupled to a dielectric layer (supporting
the forward mode). We have elaborated that the scheme discussed here exhibits
a photonic indirect bandgap that can be controlled by changing the group
velocities of the forward and backward modes along with the coupling
coefficient. Two conditions may occur: the bandgap is associated either to a
positive or negative momentum difference between the two energy levels. The
finding in this paper can be useful to design systems with EPDs whose use is
of growing importance for enhancing light-matter interactions and nonlinear
photonic phenomena.

## Acknowledgment

This material is based on work supported by the Air Force Office of Scientific
Research award number FA9550-18-1-0355, and by the National Science Foundation
under award NSF ECCS-1711975. The authors are thankful to DS SIMULIA for
providing CST Studio Suite that was instrumental in this study. The authors
thank Dr. Hamidreza Kazemi, Dr. Mohamed Yehia Nada, and Dr. Ahmed Abdelshafy,
UC Irvine, for fruitful discussions on coupled mode theory.

## References

  * [1] M. I. Vishik and L. A. Lyusternik, “The solution of some perturbation problems for matrices and selfadjoint or non-selfadjoint differential equations i,” Russian Mathematical Surveys, vol. 15, pp. 1–73, jun 1960, doi: 10.1070/rm1960v015n03abeh004092.
  * [2] P. Lancaster, “On eigenvalues of matrices dependent on a parameter,” Numerische Mathematik, vol. 6, pp. 377–387, dec 1964, doi: 10.1007/bf01386087.
  * [3] A. P. Seyranian, “Sensitivity analysis of multiple eigenvalues,” Journal of Structural Mechanics, vol. 21, pp. 261–284, jan 1993, doi: 10.1080/08905459308905189.
  * [4] T. Kato, Perturbation Theory for Linear Operators.  Springer-Verlag New York Inc., New York, 1966, doi: 10.1007/978-3-662-12678-3.
  * [5] C. M. Bender and S. Boettcher, “Real spectra in non-Hermitian Hamiltonians having PT symmetry,” Physical Review Letters, vol. 80, no. 24, p. 5243, 1998.
  * [6] C. M. Bender, “Making sense of non-hermitian hamiltonians,” Reports on Progress in Physics, vol. 70, no. 6, p. 947, 2007.
  * [7] A. Mostafazadeh, “Exact pt-symmetry is equivalent to hermiticity,” Journal of Physics A: Mathematical and General, vol. 36, no. 25, p. 7081, 2003\.
  * [8] V. S. Asadchy, M. S. Mirmoosa, A. Díaz-Rubio, S. Fan, and S. A. Tretyakov, “Tutorial on electromagnetic nonreciprocity and its origins,” Proceedings of the IEEE, vol. 108, no. 10, pp. 1684–1727, 2020.
  * [9] A. Ruschhaupt, F. Delgado, and J. Muga, “Physical realization of-symmetric potential scattering in a planar slab waveguide,” Journal of Physics A: Mathematical and General, vol. 38, no. 9, p. L171, 2005.
  * [10] R. El-Ganainy, K. Makris, D. Christodoulides, and Z. H. Musslimani, “Theory of coupled optical PT-symmetric structures,” Optics Letters, vol. 32, no. 17, pp. 2632–2634, 2007.
  * [11] A. Guo, G. Salamo, D. Duchesne, R. Morandotti, M. Volatier-Ravat, V. Aimez, G. Siviloglou, and D. Christodoulides, “Observation of p t-symmetry breaking in complex optical potentials,” Physical Review Letters, vol. 103, no. 9, p. 093902, 2009.
  * [12] M. V. Berry, “Physics of nonhermitian degeneracies,” Czechoslovak Journal of Physics, vol. 54, no. 10, pp. 1039–1047, 2004.
  * [13] A. Yariv, “Coupled-mode theory for guided-wave optics,” IEEE Journal of Quantum Electronics, vol. 9, no. 9, pp. 919–933, 1973.
  * [14] A. Hardy and W. Streifer, “Coupled mode theory of parallel waveguides,” Journal of lightwave technology, vol. 3, no. 5, pp. 1135–1146, 1985.
  * [15] T. Mealy and F. Capolino, “General conditions to realize exceptional points of degeneracy in two uniform coupled transmission lines,” IEEE Transactions on Microwave Theory and Techniques, vol. 68, no. 8, pp. 3342–3354, 2020.
  * [16] Y. Li, X. Guo, L. Chen, C. Xu, J. Yang, X. Jiang, and M. Wang, “Coupled mode theory under the parity-time symmetry frame,” Journal of Lightwave Technology, vol. 31, no. 15, pp. 2477–2481, 2013.
  * [17] M. A. Othman and F. Capolino, “Theory of exceptional points of degeneracy in uniform coupled waveguides and balance of gain and loss,” IEEE Transactions on Antennas and Propagation, vol. 65, no. 10, pp. 5289–5302, 2017\.