# Numerical Study of Magnetic Island Coalescence Using Magnetohydrodynamics
With Adaptively Embedded Particle-In-Cell Model

Dion Li  Department of Nuclear Engineering and Radiological Sciences,
University of Michigan, Ann Arbor, Michigan 48109, USA  Yuxi Chen  Department
of Astrophysical Sciences, Princeton University, Princeton, New Jersey 08544,
USA  Princeton Plasma Physics Laboratory, Princeton University, Princeton, New
Jersey 08540, USA  Chuanfei Dong<EMAIL_ADDRESS>Princeton Plasma Physics
Laboratory, Princeton University, Princeton, New Jersey 08540, USA  Department
of Astrophysical Sciences, Princeton University, Princeton, New Jersey 08544,
USA  Department of Astronomy, Boston University, Boston, Massachusetts 02215,
USA  Liang Wang  Department of Astrophysical Sciences, Princeton University,
Princeton, New Jersey 08544, USA  Princeton Plasma Physics Laboratory,
Princeton University, Princeton, New Jersey 08540, USA  Gabor Toth  Department
of Climate and Space Sciences and Engineering, University of Michigan, Ann
Arbor, Michigan 48109, USA

###### Abstract

Collisionless magnetic reconnection typically requires kinetic treatments that
are, in general, computationally expensive compared to fluid-based models. In
this study, we use the magnetohydrodynamics with adaptively embedded particle-
in-cell (MHD-AEPIC) model to study the interaction of two magnetic flux ropes.
This innovative model embeds one or more adaptive PIC regions into a global
MHD simulation domain such that the kinetic treatment is only applied in
regions where kinetic physics is prominent. We compare the simulation results
among three cases: 1) MHD with adaptively embedded PIC regions, 2) MHD with
statically (or fixed) embedded PIC regions, and 3) a full PIC simulation. The
comparison yields good agreement when analyzing their reconnection rates and
magnetic island separations, as well as the ion pressure tensor elements and
ion agyrotropy. In order to reach a good agreement among the three cases,
large adaptive PIC regions are needed within the MHD domain, which indicates
that the magnetic island coalescence problem is highly kinetic in nature where
the coupling between the macro-scale MHD and micro-scale kinetic physics is
important.

††preprint: APS/123-QED

##  I  Introduction

Magnetic reconnection is a process that occurs on the order of Alfvénic
timescales in which the magnetic topology is rearranged and the magnetic
energy is converted into plasma kinetic or thermal energy Finn and Kaw (1977);
Bondeson (1983). Because of its broad applications in solar flares Forbes
(1991) and corona Chen  _et al._ (2021); Dong  _et al._ (2022), coronal mass
ejections (CMEs) Qiu  _et al._ (2007), Earth’s and planetary magnetospheres
Chen  _et al._ (2008); Le  _et al._ (2017); Phan  _et al._ (2018); Yang, Wang,
and Dong (2022), and also in laboratory plasmas Yamada  _et al._ (1994); Duan
_et al._ (2018); Raymond  _et al._ (2018), reconnection has been of immense
interest in recent years Ji  _et al._ (2022).

It has been well established that collisionless magnetic reconnection occurs
in the regime beyond ideal magnetohydrodynamics, where kinetic-scale physics
becomes prominent in the formation and evolution of the thin current sheets
near the reconnection sites Zweibel and Yamada (2009); Yamada, Kulsrud, and Ji
(2010). Thus, simulating collisionless magnetic reconnection typically
requires kinetic approaches such as the particle-in-cell (PIC) method. In
general, however, such kinetic simulations are computationally expensive and
are hence unable to efficiently solve large-scale problems involving
collisionless physics. In order to solve this issue with affordable
computational costs, two broad approaches have been proposed for large-scale
global simulations, i.e., the magnetohydrodynamics with embedded particle-in-
cell (MHD-EPIC) model Daldorff  _et al._ (2014); Tóth  _et al._ (2016); Chen
_et al._ (2017, 2019) and the multi-moment multi-fluid model Wang  _et al._
(2015, 2018); Dong  _et al._ (2019); Wang  _et al._ (2020); Jarmak  _et al._
(2020); Rulke, Wang, and Dong (2021). Recently, some progress has been made to
improve the fluid closure in the multi-moment multi-fluid model through
machine learning Qin  _et al._ (2022); Cheng  _et al._ (2023), and meanwhile,
MHD-EPIC has been improved to incorporate the feature of adaptively embedded
PIC regions (MHD-AEPIC) Shou  _et al._ (2021); Chen  _et al._ (2021); Wang,
Chen, and Tóth (2022), which offers flexibility for the PIC code to capture
the localized regions where kinetic physics is important.

It has been well demonstrated through a series of local studies Wang  _et al._
(2015); Ng  _et al._ (2015, 2017, 2019) that the multi-moment multi-fluid
model incorporating the higher-order moments is capable of reproducing some
critical aspects of the collisionless reconnection physics from fully kinetic
simulations. However, no such systematic local studies have been conducted
using the MHD-AEPIC model, which motivates this study to employ the MHD-AEPIC
model to investigate the magnetic island coalescence problem that is highly
kinetic in nature Stanier  _et al._ (2015); Ng  _et al._ (2015).

In an MHD-AEPIC simulation, only part of the simulation domain, where kinetic
effects are important, is simulated by the semi-implicit PIC code, Flexible
Exascale Kinetic Simulator (FLEKS) Chen  _et al._ (2021) and the rest of the
domain is handled by the MHD (or Hall MHD) model, Block-Adaptive-Tree
Solarwind Roe-type Upwind Scheme (BATS-R-US) Powell  _et al._ (1999); Tóth
_et al._ (2012). For the magnetic island coalescence problem, the embedded PIC
regions are applied to simulate the regions with strong current density. As we
will show later, the adaptive simulation comes in close agreement with the
simulation employing the fully kinetic PIC code, Object-oriented Simulation
Rapid Implementation System (OSIRIS) Fonseca  _et al._ (2002); Hemker (2015),
when analyzing their out-of-plane current densities, reconnection rates,
O-point separations, pressure tensor elements, and agyrotropy (a measure of
the deviation of the pressure tensor of a species from cylindrical symmetry
with respect to the direction of the local magnetic field Scudder and Daughton
(2008)).

In Sec. II, we discuss the model setup for the magnetic island coalescence
problem. In Sec. III, simulation results are presented and discussed by
comparing the full PIC simulation results with the outputs from the MHD-EPIC
and MHD-AEPIC models. Concluding remarks are given in Sec. IV.

##  II Model Setup and Methods

In this study, we set a Fadeev equilibrium Fadeev, Kvabtskhava, and Komarov
(1965) as the initial condition where the initial magnetic vector potential is
expressed as

$\displaystyle\textbf{A}=B_{0}\lambda\ln\big{[}\epsilon\cos(x/\lambda)+\cosh(y/\lambda)\big{]}\hat{\textbf{z}}$
(1)

where $B_{0}$ is the asymptotic magnetic field away from the x-axis,
$\epsilon=0.4$ is a measure of the island size, and the initial density
distribution is written as

$\displaystyle
n=\frac{n_{0}(1-\epsilon^{2})}{\big{[}\epsilon\cos(x/\lambda)+\cosh(y/\lambda)\big{]}^{2}}+n_{b}.$
(2)

where $n_{b}=0.2n_{0}$ is the background number density. We define the global
Alfvén time as $t_{A}=L_{x}/v_{A}$, where $L_{x}=4\pi\lambda$ is the length of
the simulation box, $\lambda=5d_{i0}$ is the half-width of the current sheet
(we also run additional simulations with $\lambda=10d_{i0}$; see Fig. 7),
$d_{i0}=(m_{i}/\mu_{0}n_{0}e^{2})^{1/2}$ is the ion inertial length, and
$v_{A}=B_{0}/(\mu_{0}n_{0}m_{i})^{1/2}$ is the Alfvén speed. The width of the
box is set as $L_{y}=2\pi\lambda$. In the full PIC simulation with OSIRIS, we
specify the ion and electron out-of-plane current densities as the following
relationship $J_{zi0}/J_{ze0}=T_{i0}/T_{e0}$, where $T_{i0}$ and $T_{e0}$ are
uniform. In the MHD-AEPIC simulations, the electron and ion bulk velocities
are initialized from the MHD current and momentum Daldorff  _et al._ (2014).
In the following simulations, we use the same configuration described in
Stanier  _et al._ (2015) and Ng  _et al._ (2015), setting $T_{i0}=T_{e0}$,
$m_{i}=25m_{e}$ (so $d_{i}=5d_{e}$), $t_{max}=2.5t_{A}$, and
$n_{0}k_{B}(T_{i0}+T_{e0})=B_{0}^{2}/2\mu_{0}$ (so $p=p_{mag}$). The
relationship between the electron plasma and cyclotron frequencies is given by
$\omega_{pe}=2\Omega_{ce}$ where
$\omega_{pe}=(n_{0}e^{2}/m_{e}\epsilon_{0})^{1/2}$ and
$\Omega_{ce}=|e|B_{0}/m_{e}$. The initial perturbation follows the same form
as in Daughton  _et al._ (2009a), with

$\displaystyle\delta B_{x}=\delta B_{0}\cos\bigg{(}\frac{2\pi
x}{L_{x}}\bigg{)}\sin\bigg{(}\frac{\pi y}{L_{y}}\bigg{)}$ (3a)
$\displaystyle\delta B_{y}=-\delta B_{0}\sin\bigg{(}\frac{2\pi
x}{L_{x}}\bigg{)}\cos\bigg{(}\frac{\pi y}{L_{y}}\bigg{)}$ (3b)

where we set the perturbation amplitude as $\delta B_{0}=0.1B_{0}$.

To simulate the merging of two magnetic islands, we first perform two
simulations using the MHD-AEPIC model, one with a large static, non-adaptive
(or fixed) PIC region and the other with adaptive PIC regions. We also run a
full PIC simulation using OSIRIS to validate and compare with the MHD-AEPIC
simulations. In all simulations, we have used periodic boundaries in the
horizontal ($x$) direction, and conducting electromagnetic boundaries and
reflecting particle boundaries in the vertical ($y$) direction.

The OSIRIS simulation domain consists of 2000$\times$1000 cells and each cell
has 256 particles per species. In the MHD-AEPIC simulations, the MHD grid
consists of 2048$\times$1024 cells. The effective PIC grid resolution and the
initial particle number per cell are the same as the OSIRIS simulation, and
the Gauss’s Law-satisfying Energy Conserving Semi-Implicit Method (GL-ECSIM)
Chen and Tóth (2019) is applied. In general, the MHD and PIC grid size does
not have to be the same Chen  _et al._ (2021) (see Appendix A). In the
magnetic island coalescence simulation, the two magnetic islands are
associated with strong current density, which suggests a notable separation
between the electron and ion velocities, and kinetic effects may play an
important role when two islands move toward each other. Outside the islands,
the current is generally weak, and neither the magnetic field nor plasma
quantities show a strong gradient, so an MHD model is sufficient to describe
this region. In the MHD-AEPIC run with adaptive PIC regions, the PIC regions
only cover the areas of $J\Delta x/B>0.01$, where the selection criterion
$J\Delta x/B$ is dimensionless. The threshold 0.01 is carefully chosen based
on numerical experiments such that the islands are well covered by the PIC
region. The selection criterion is calculated from the MHD results. We usually
run a pure MHD simulation first to estimate a proper threshold that would
determine the desired PIC regions, then we run short MHD-AEPIC simulations
with a coarse PIC grid to optimize the threshold. The small threshold for the
island coalescence problem is due to the fact that it is highly kinetic in
nature where the coupling between the macro-scale MHD and micro-scale kinetic
physics is important Stanier  _et al._ (2015); Ng  _et al._ (2019). In
general, the PIC region of an MHD-AEPIC simulation is selected based on the
nature of a problem. For comparison, we also perform a simulation using the
MHD-EPIC model with a relatively large static or fixed PIC region, which
covers the majority of the entire simulation domain.

In MHD-AEPIC and MHD-EPIC simulations, when Hall MHD is used (which is the
case for the current study), the Alfven and Whistler modes can travel through
the MHD-PIC boundaries smoothly Daldorff  _et al._ (2014). However, if there
are kinetic modes reaching the MHD-PIC boundaries, they may be artificially
reflected since they cannot propagate into the MHD regions. For such
applications, we should increase the PIC regions so that the kinetic modes
would have been damped near the PIC regions. The MHD-AEPIC model is initially
designed for global magnetosphere simulations with a relatively small kinetic
region. If the entire simulation domain is dominated by kinetic physics, a
full PIC or hybrid-PIC Dong  _et al._ (2021); Winske  _et al._ (2022) code is
a better choice.

We note that the boundary conditions, with the exception of periodic
boundaries, introduce numerical perturbations inevitably for all kinetic
codes, which is a natural consequence of numerical discretization. The
perturbations, however, would be reduced if the boundaries are far away from
the kinetic regions. We, therefore, usually choose PIC regions where the
solution is smooth at the MHD-PIC interface, such as the simulations presented
in this paper.

##  III Simulation Results

Figure 1:  Out-of-plane current density at $t=0.9t_{A}$ for the adaptive (top)
and fixed (middle) PIC region runs and the full-PIC run (bottom) with
$\lambda=5d_{i0}$, overlaid with the magnetic vector potential contours (black
dotted curves) and with the boundary between the MHD and PIC regions (green
solid curves). Figure 2:  Ratio of the active PIC regions over the entire
simulation domain for the MHD-AEPIC simulation. The regions that are covered
by the semi-implicit PIC code increase from $40\%$ to about $50\%$ during the
simulation.

Fig. 1 shows 2D plots of the out-of-plane current density at $t=0.9t_{A}$ from
the three different simulations, where the $x$ and $y$ coordinates have been
normalized over $d_{i0}$ and the current density over $|e|n_{0}c$. In each
panel of Fig. 1, the current sheet formation, which cannot be correctly
captured in a Hall MHD model for this problem Stanier  _et al._ (2015), are
indicated as regions of positive out-of-plane current densities near the
origin ($x,y=0$). Also illustrated in Fig. 1 are the boundaries (green curves)
between the MHD and PIC regions in both the adaptive-PIC-region and fixed-PIC-
region runs.

For the adaptive-PIC-region run, the variation of the active PIC cell number
inside the entire simulation domain over time is shown in Fig. 2. Initially,
about 40% of the simulation domain is run by PIC code, then the ratio
gradually increases to about 50%. In this simulation, the PIC regions are
defined as the areas with large current density values and they do not vary
significantly during the simulation. For other applications, such as modeling
the Earth magnetotail Wang, Chen, and Tóth (2022), the active PIC regions may
change dramatically. We note that the smallest granularity to turn on or off
PIC cells is a patch with two cells in each direction Chen  _et al._ (2021).
Since the MHD-AEPIC model adopts a semi-implicit PIC code, which requires an
iterative solver to update the electric field Chen and Tóth (2019), it
requires more computations per step. The computational efficiency also depends
on the implementation and problem properties. For the simulations presented in
Fig. 1, their computational cost (CPU hours) ratio is about $3:8:1$ among the
adaptive-PIC region case, the fixed-PIC region case, and the full-PIC case.
Given that FLEKS is a semi-implicit PIC code and the gird resolution between
the PIC and MHD regions can be different, one can reduce the grid resolution
of the PIC cells inside the MHD-AEPIC domain, which can lower the
computational cost of the MHD-AEPIC run and thus the MHD-AEPIC run becomes
computationally cheaper than the full-PIC case (see Appendix A). It is
noteworthy that in three-dimensional (3-D) cases, MHD-AEPIC can be
computationally even more efficient than a full PIC code.

Figure 3:  Reconnection rate as a function of time for the adaptive and fixed
PIC region runs and the full-PIC run with $\lambda=5d_{i0}$. Figure 4:
O-point separation of the coalescing magnetic islands as a function of time
for the adaptive and fixed PIC region runs and the full-PIC run with
$\lambda=5d_{i0}$.

An inspection of Fig. 1 reveals that the simulation results (such as the
current density and the magnetic flux) from three different runs present
nearly identical features and the system evolves at the same rate, indicating
the reconnection rates, $E_{R}$, from the three simulations are very similar.
In order to verify this idea, we explicitly compare the reconnection rates. We
normalize the reconnection rate to the maximum initial magnetic field between
the islands such that Stanier  _et al._ (2015)

$\displaystyle E_{R}=\frac{1}{v_{Am}B_{m}}\frac{\partial}{\partial
t}\big{(}A_{zX}-A_{zO}\big{)}$ (4)

where $v_{Am}=B_{m}/(\mu_{0}n_{0}m_{i})^{1/2}$, $A_{zX}$ is the out-of-plane
magnetic vector potential, evaluated at the $X$-point and $A_{zO}$ is
evaluated at the $O$-point. The time evolution of the reconnection rate for
the three simulations are plotted in Fig. 3, where a maximum reconnection rate
occurs around $t=0.9t_{A}$ for all the cases. After the reconnection rate
reaches a maximum, it decreases with time before reaching a constant value at
approximately $t=1.6t_{A}$. We note from Fig. 3 that the case with adaptive
PIC regions closely captures the behavior of $E_{R}$ for this problem, when
compared to the full-PIC simulation.

Since the O-point separation between the magnetic islands is largely dependent
on reconnection rates, we demonstrate the O-point separation as a function of
time in Fig. 4. In this plot, we have normalized the O-point separation by the
initial separation between the islands $L_{0}$. It is evident in Fig. 4 that
for all three cases the island separation decreases at an increasing rate
before $t=0.9t_{A}$, at which point the coalescence of the islands slows down,
reaching an approximately constant rate at $t=1.6t_{A}$. Both Figs. 3 and 4
show excellent agreement among the three simulation approaches, validating the
use of the MHD-AEPIC model.

At this point, we are interested in analyzing the effects of ion kinetics, as
its importance has been demonstrated previously Stanier  _et al._ (2015).
Stanier  _et al._ (2015) also pointed out that electron kinetics were not
crucial for the island coalescence problem. To do this, we write the
z-component of the normalized ion Ohm’s law as

$\displaystyle E_{z}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{d_{i}}{n}\bigg{[}\frac{\partial}{\partial
t}(nv_{iz})+\nabla\cdot(n\textbf{v}_{i}v_{iz})\bigg{]}+\frac{d_{i}}{n}\nabla\cdot\big{(}\bar{\bar{\textbf{P}}}_{i}\cdot\textbf{z}\big{)}$
(5)

where
$E_{z}^{\prime}=(\textbf{E}+\textbf{v}_{i}\times\textbf{B})\cdot\textbf{z}$ is
the z-component of the nonideal electric field, $\bar{\bar{\textbf{P}}}_{i}$
is the ion pressure tensor, and the resistivity is neglected for collisionless
reconnection. Near the X-point of the reconnection site, where the magnetic
field is small, pressure tensor effects become important. Specifically, Eq.(5)
demonstrates that the non-ideal electric field is heavily influenced by the
off-diagonal terms of the pressure tensor. Since it has been shown in Stanier
_et al._ (2015) that ion kinetic effects are of importance in this specific
reconnection setup, we focus our analysis on the ion pressure tensor. Fig. 5
compares the off-diagonal ion pressure tensor elements from the three cases
when the reconnection rate reaches its peak, once again showing spectacular
agreement among the cases. Fig. 5 also demonstrates that these two islands
should be covered by the PIC code, since the off-diagonal pressure tensor can
not be described by the isotropic MHD used here. In general, if the ion
pressure is anisotropic, but the off-diagonal terms are relatively small in
the local magnetic field coordinate system, MHD-AEPIC supports coupling with
an anisotropic MHD model Daldorff  _et al._ (2014), and the region that
requires PIC can be reduced. However, if the off-diagonal terms are prominent
in the local magnetic field coordinates, which is the case here (see Appendix
B), the PIC regions need to cover those areas, consistent with the previous
study using a similar approach Makwana, Keppens, and Lapenta (2018).

Figure 5:  Comparison of the off-diagonal ion pressure tensor elements at
$t=0.9t_{A}$ for the adaptive (top row) and fixed (middle row) PIC region runs
and the full-PIC run (bottom row) with $\lambda=5d_{i0}$. The left column
plots $P_{xy}$, the middle column shows $P_{xz}$, and $P_{yz}$ is illustrated
in the right column.

In addition to the off-diagonal ion pressure tensor terms, we also compare the
ion agyrotropies of the three simulations, employing the same formula for
agyrotropy as in Scudder and Daughton (2008). We note that agyrotropy is a
measure useful for characterizing the ion (or electron) diffusion region in
collisionless magnetic reconnection. We once again focus only on ion
agyrotropy as ion kinetics are required to capture the correct reconnection
rates and describe the global behavior of the system for the island
coalescence problem Stanier  _et al._ (2015). A comparison of the agyrotropy
from the three different simulations is shown in Fig. 6, indicating excellent
agreement among different cases, especially near the current sheet formation
where ion agyrotropy is of particular significance (indicated by the lighter
colors inside the ion diffusion region).

Figure 6:  Ion agyrotropy at $t=0.9t_{A}$ for the adaptive (top) and fixed
(middle) PIC region runs and the full-PIC run (bottom) with $\lambda=5d_{i0}$,
overlaid with the magnetic vector potential contours (in white).

In Fig. 7, we also present the reconnection rates for simulations with
$\lambda=10d_{i0}$ to demonstrate that reconnection in large systems becomes
slower. The islands are unable to coalesce on the first approach due to the
slower reconnection rates, and so bounce off each other Karimabadi  _et al._
(2011); Stanier  _et al._ (2015). All simulations reach nearly the same
maximum reconnection rate and show two peaks around $t=0.9t_{A}$ and
$t=1.2t_{A}$, respectively. Overall, all simulations show similar reconnection
processes.

##  IV Conclusion

It has been shown that simulating collisionless magnetic reconnection requires
kinetic approaches such as the PIC method. However, these simulation codes are
generally computationally expensive and hence make it difficult to efficiently
model large-scale global systems. In this study, we have presented an
alternative method by utilizing the MHD-AEPIC model that uses the PIC
treatment in regions where the current sheet formation is prominent and uses
the computationally cheap MHD model in the rest of the simulation domain. We
have shown that the case with adaptive PIC regions comes in close agreement
with the full PIC simulation when analyzing reconnection rates and O-point
separations as well as the ion pressure tensor elements and ion agyrotropy.
These results demonstrate that the MHD-AEPIC model may accurately and
efficiently simulate large-scale systems that involve collisionless
reconnection physics.

It should be noted that the magnetic island coalescence problem that we
studied here is highly kinetic in nature where the coupling between the MHD
and kinetic scales is important Daughton  _et al._ (2009b); Stanier  _et al._
(2015); Ng  _et al._ (2019), so relatively large PIC regions are needed within
the MHD domain. For large-scale simulations such as the case of the Earth’s
magnetosphere, Chen  _et al._ (2021) applied MHD-AEPIC to study the
magnetopause reconnection while treating the magnetotail reconnection using
Hall MHD, such that they drastically saved computational costs compared with
the full PIC magnetosphere simulations.

We also note that in magnetic confinement fusion, near the X-point where the
magnetic field is weak, gyrokinetic approximations begin to break down and it
becomes necessary to include full kinetic physics at these locations. The MHD-
AEPIC model may be useful for simulating FLARE Ji  _et al._ (2020) and NSTX-U
Ebrahimi and Raman (2016) at Princeton Plasma Physics Laboratory, MAST at the
UK Tanabe  _et al._ (2017) as well as TREX at Wisconsin Plasma Physics
Laboratory Olson  _et al._ (2021).

Figure 7:  Reconnection rate as a function of time for the adaptive and fixed
PIC region runs and the full-PIC run with $\lambda=10d_{i0}$. Figure 8:  Out-
of-plane current density at $t=0.9t_{A}$ for the adaptive PIC region run
(top), the adaptive PIC region run with half PIC grid resolution (middle), and
the full-PIC run (bottom) with $\lambda=5d_{i0}$, overlaid with the magnetic
vector potential contours (black dotted curves) and with the boundary between
the MHD and PIC regions (green solid curves). Figure 9:  Reconnection rate as
a function of time for the three cases shown in Fig.8. Figure 10:  Comparison
of the off-diagonal ion pressure tensor elements at $t=0.9t_{A}$ for the
adaptive (top row) and fixed (middle row) PIC region runs and the full-PIC run
(bottom row) in field-aligned coordinates with $\lambda=5d_{i0}$. The left
column plots $P_{12}^{\prime}$, the middle column shows $P_{13}^{\prime}$, and
$P_{23}^{\prime}$ is illustrated in the right column.

###### Acknowledgements.

This work was made possible by support from the Department of Energy for the
Summer Undergraduate Laboratory Internship (SULI) program. This work was
supported by the U.S. Department of Energy under contract number DE-
AC02-09CH11466 (through LDRD) and the DOE grant DE-SC0021205, the NASA grants
80NSSC21K1326 and 80NSSC22K0323, and the NSF grant AGS-2149787. Resources
supporting this work were provided by the NASA High-End Computing (HEC)
Program through the NASA Advanced Supercomputing (NAS) Division at Ames
Research Center. The authors acknowledge the OSIRIS Consortium, consisting of
UCLA and IST (Portugal) for the use of the OSIRIS 4.0 framework.

## Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

##  Appendix A MHD-AEPIC Simulations with Different PIC Grid Resolutions

Fig. 8 compares the out-of-plane current densities at $t=0.9t_{A}$ among the
adaptive-PIC-region case, the adaptive-PIC-region case (with half PIC grid
resolution), and the full-PIC simulation. The MHD grid resolution between the
two adaptive PIC cases is the same, which demonstrates the flexibility and
capability of the coupling algorithm between the MHD model and the semi-
implicit PIC code. As shown in Figs. 8 & 9, all the simulations show very
similar results. In general, a semi-implicit PIC code allows relaxed stability
constraints, therefore, it can use a coarser grid resolution than an explicit
PIC code.

##  Appendix B Off-Diagonal Pressure Tensor Elements in Field-Aligned
Coordinates

Fig. 10 depicts the off-diagonal ion pressure tensor elements in field-aligned
coordinates (or local magnetic field coordinates) for the three cases at
$t=0.9t_{A}$, showing that the off-diagonal terms are prominent in the local
magnetic field coordinates. The rotated pressure tensor elements
$P^{\prime}_{mn}$ in field-aligned coordinates are defined as follows,

$\displaystyle P^{\prime}_{mn}=(R^{T}PR)_{mn}=R_{im}P_{ij}R_{jn}$ (6)

where $P_{ij}$ are the non-rotated ion pressure tensor elements, the rotation
matrix is

and b is the unit vector along the local magnetic field direction.

## REFERENCES

## References

  * Finn and Kaw (1977) J. M. Finn and P. K. Kaw, “Coalescence instability of magnetic islands,” The Physics of Fluids 20, 72–78 (1977), DOI:10.1063/1.861709.
  * Bondeson (1983) A. Bondeson, “Linear analysis of the coalescence instability,” The Physics of Fluids 26, 1275–1278 (1983), DOI:10.1063/1.864287.
  * Forbes (1991) T. G. Forbes, “Magnetic reconnection in solar flares,” Geophysical $\&$ Astrophysical Fluid Dynamics 62, 15–36 (1991).
  * Chen  _et al._ (2021) Y. Chen, D. Przybylski, H. Peter, H. Tian, F. Auchère, and D. Berghmans, “Transient small-scale brightenings in the quiet solar corona: A model for campfires observed with Solar Orbiter,” Astronomy & Astrophysics 656, L7 (2021), arXiv:2104.10940 [astro-ph.SR] .
  * Dong  _et al._ (2022) C. Dong, L. Wang, Y.-M. Huang, L. Comisso, T. A. Sandstrom, and A. Bhattacharjee, “Reconnection-driven energy cascade in magnetohydrodynamic turbulence,” Science Advances 8, eabn7627 (2022).
  * Qiu  _et al._ (2007) J. Qiu, Q. Hu, T. A. Howard, and V. B. Yurchyshyn, “On the magnetic flux budget in low-corona magnetic reconnection and interplanetary coronal mass ejections,” The Astrophysical Journal 659, 758–772 (2007).
  * Chen  _et al._ (2008) L. J. Chen, A. Bhattacharjee, P. A. Puhl-Quinn, H. Yang, N. Bessho, S. Imada, S. Mühlbachler, P. W. Daly, B. Lefebvre, Y. Khotyaintsev, A. Vaivads, A. Fazakerley, and E. Georgescu, “Observation of energetic electrons within magnetic islands,” Nature Physics 4, 19–23 (2008).
  * Le  _et al._ (2017) A. Le, W. Daughton, L. J. Chen, and J. Egedal, “Enhanced electron mixing and heating in 3-D asymmetric reconnection at the Earth’s magnetopause,” Geophysical Research Letters 44, 2096–2104 (2017), arXiv:1703.10246 [physics.plasm-ph] .
  * Phan  _et al._ (2018) T. D. Phan, J. P. Eastwood, M. A. Shay, J. F. Drake, B. U. O. Sonnerup, M. Fujimoto, P. A. Cassak, M. Øieroset, J. L. Burch, R. B. Torbert, A. C. Rager, J. C. Dorelli, D. J. Gershman, C. Pollock, P. S. Pyakurel, C. C. Haggerty, Y. Khotyaintsev, B. Lavraud, Y. Saito, M. Oka, R. E. Ergun, A. Retino, O. L. Contel, M. R. Argall, B. L. Giles, T. E. Moore, F. D. Wilder, R. J. Strangeway, C. T. Russel, P. A. Lindqvist, and W. Magnes, “Electron magnetic reconnection without ion coupling in earth’s turbulent magnetosheath,” Nature 557, 202–206 (2018).
  * Yang, Wang, and Dong (2022) S.-D. Yang, L. Wang, and C. Dong, “Discovery of Double Hall Pattern Associated with Collisionless Magnetic Reconnection in Dusty Plasmas,” arXiv e-prints , arXiv:2206.08553 (2022), arXiv:2206.08553 [physics.plasm-ph] .
  * Yamada  _et al._ (1994) M. Yamada, F. M. Levinton, N. Pomphrey, R. Budny, J. Manickam, and Y. Nagayama, “Investigation of magnetic reconnection during a sawtooth crash in a high‐temperature tokamak plasma,” Physics of Plasmas 1, 3269–3276 (1994).
  * Duan  _et al._ (2018) S. C. Duan, M. X. Kan, B. Xiao, Q. Xu, L. Yang, G. H. Wang, X. B. Huang, and W. P. Xie, “Numerical modelling of inverse wire array z-pinch magnetic reconnection,” AIP Advances 8, 055018 (2018).
  * Raymond  _et al._ (2018) A. E. Raymond, C. F. Dong, A. McKelvey, C. Zulick, N. Alexander, A. Bhattacharjee, P. T. Campbell, H. Chen, V. Chvykov, E. Del Rio, P. Fitzsimmons, W. Fox, B. Hou, A. Maksimchuk, C. Mileham, J. Nees, P. M. Nilson, C. Stoeckl, A. G. R. Thomas, M. S. Wei, V. Yanovsky, K. Krushelnick, and L. Willingale, “Relativistic-electron-driven magnetic reconnection in the laboratory,” Physical Review E 98, 043207 (2018).
  * Ji  _et al._ (2022) H. Ji, W. Daughton, J. Jara-Almonte, A. Le, A. Stanier, and J. Yoo, “Magnetic reconnection in the era of exascale computing and multiscale experiments,” Nature Reviews Physics 4, 263–282 (2022).
  * Zweibel and Yamada (2009) E. G. Zweibel and M. Yamada, “Magnetic reconnection in astrophysical and laboratory plasmas,” Annual Review of Astronomy and Astrophysics 47, 291–332 (2009).
  * Yamada, Kulsrud, and Ji (2010) M. Yamada, R. Kulsrud, and H. Ji, “Magnetic reconnection,” Rev. Mod. Phys. 82, 603–664 (2010).
  * Daldorff  _et al._ (2014) L. K. S. Daldorff, G. Tóth, T. I. Gombosi, G. Lapenta, J. Amaya, S. Markidis, and J. U. Brackbill, “Two-way coupling of a global Hall magnetohydrodynamics model with a local implicit particle-in-cell model,” Journal of Computational Physics 268, 236–254 (2014).
  * Tóth  _et al._ (2016) G. Tóth, X. Jia, S. Markidis, I. B. Peng, Y. Chen, L. K. S. Daldorff, V. M. Tenishev, D. Borovikov, J. D. Haiducek, T. I. Gombosi, A. Glocer, and J. C. Dorelli, “Extended magnetohydrodynamics with embedded particle-in-cell simulation of Ganymede’s magnetosphere,” Journal of Geophysical Research (Space Physics) 121, 1273–1293 (2016).
  * Chen  _et al._ (2017) Y. Chen, G. Tóth, P. Cassak, X. Jia, T. I. Gombosi, J. A. Slavin, S. Markidis, I. B. Peng, V. K. Jordanova, and M. G. Henderson, “Global three-dimensional simulation of earth’s dayside reconnection using a two-way coupled magnetohydrodynamics with embedded particle-in-cell model: Initial results,” Journal of Geophysical Research: Space Physics 122, 10,318–10,335 (2017), DOI:10.1002/2017JA024186.
  * Chen  _et al._ (2019) Y. Chen, G. Tóth, X. Jia, J. A. Slavin, W. Sun, S. Markidis, T. I. Gombosi, and J. M. Raines, “Studying Dawn-Dusk Asymmetries of Mercury’s Magnetotail Using MHD-EPIC Simulations,” Journal of Geophysical Research (Space Physics) 124, 8954–8973 (2019), arXiv:1904.06753 [physics.space-ph] .
  * Wang  _et al._ (2015) L. Wang, A. H. Hakim, A. Bhattacharjee, and K. Germaschewski, “Comparison of multi-fluid moment models with particle-in-cell simulations of collisionless magnetic reconnection,” Physics of Plasmas 22, 012108 (2015), DOI:10.1063/1.4906063.
  * Wang  _et al._ (2018) L. Wang, K. Germaschewski, A. Hakim, C. Dong, J. Raeder, and A. Bhattacharjee, “Electron Physics in 3-D Two-Fluid 10-Moment Modeling of Ganymede’s Magnetosphere,” Journal of Geophysical Research (Space Physics) 123, 2815–2830 (2018).
  * Dong  _et al._ (2019) C. Dong, L. Wang, A. Hakim, A. Bhattacharjee, J. A. Slavin, G. A. DiBraccio, and K. Germaschewski, “Global Ten-Moment Multifluid Simulations of the Solar Wind Interaction with Mercury: From the Planetary Conducting Core to the Dynamic Magnetosphere,” Geophysical Research Letters 46, 11,584–11,596 (2019).
  * Wang  _et al._ (2020) L. Wang, A. H. Hakim, J. Ng, C. Dong, and K. Germaschewski, “Exact and locally implicit source term solvers for multifluid-Maxwell systems,” Journal of Computational Physics 415, 109510 (2020).
  * Jarmak  _et al._ (2020) S. Jarmak, E. Leonard, A. Akins, E. Dahl, D. R. Cremons, S. Cofield, A. Curtis, C. Dong, E. T. Dunham, B. Journaux, D. Murakami, W. Ng, M. Piquette, A. P. Girija, K. Rink, L. Schurmeier, N. Stein, N. Tallarida, M. Telus, L. Lowes, C. Budney, and K. L. Mitchell, “QUEST: A New Frontiers Uranus orbiter mission concept study,” Acta Astronautica 170, 6–26 (2020).
  * Rulke, Wang, and Dong (2021) E. Rulke, L. Wang, and C. Dong, “Three-dimensional, ten-moment multifluid simulation of the solar wind interaction with Asteroid 16 Psyche,” in  _AGU Fall Meeting Abstracts_ , Vol. 2021 (2021) pp. SM53C–08.
  * Qin  _et al._ (2022) Y. Qin, J. Ma, M. Jiang, C. Dong, H. Fu, L. Wang, W. Cheng, and Y. Jin, “Data-Driven Modeling of Landau Damping by Physics-Informed Neural Networks,” arXiv e-prints , arXiv:2211.01021 (2022), arXiv:2211.01021 [physics.plasm-ph] .
  * Cheng  _et al._ (2023) W. Cheng, H. Fu, L. Wang, C. Dong, Y. Jin, M. Jiang, J. Ma, Y. Qin, and K. Liu, “Data-driven, multi-moment fluid modeling of landau damping,” Computer Physics Communications 282, 108538 (2023).
  * Shou  _et al._ (2021) Y. Shou, V. Tenishev, Y. Chen, G. Tóth, and N. Ganushkina, “Magnetohydrodynamic with Adaptively Embedded Particle-in-Cell model: MHD-AEPIC,” Journal of Computational Physics 446, 110656 (2021).
  * Chen  _et al._ (2021) Y. Chen, G. Tóth, H. Zhou, and X. Wang, “Fleks: A flexible particle-in-cell code for multi-scale plasma simulations,” Earth and Space Science Open Archive (2021).
  * Wang, Chen, and Tóth (2022) X. Wang, Y. Chen, and G. Tóth, “Global magnetohydrodynamic magnetosphere simulation with an adaptively embedded particle-in-cell model,” Journal of Geophysical Research: Space Physics 127, e2021JA030091 (2022).
  * Ng  _et al._ (2015) J. Ng, Y. M. Huang, A. Hakim, A. Bhattacharjee, A. Stanier, W. Daughton, L. Wang, and K. Germaschewski, “The island coalescence problem: Scaling of reconnection in extended fluid models including higher-order moments,” Physics of Plasmas 22, 112104 (2015), DOI:10.1063/1.4935302.
  * Ng  _et al._ (2017) J. Ng, A. Hakim, A. Bhattacharjee, A. Stanier, and W. Daughton, “Simulations of anti-parallel reconnection using a nonlocal heat flux closure,” Physics of Plasmas 24, 082112 (2017), DOI:10.1063/1.4993195.
  * Ng  _et al._ (2019) J. Ng, A. Hakim, J. Juno, and A. Bhattacharjee, “Drift instabilities in thin current sheets using a two-fluid model with pressure tensor effects,” Journal of Geophysical Research: Space Physics 124, 3331–3346 (2019), DOI:10.1029/2018JA026313.
  * Stanier  _et al._ (2015) A. Stanier, W. Daughton, L. Chacón, H. Karimabadi, J. Ng, Y. M. Huang, A. Hakim, and A. Bhattacharjee, “Role of ion kinetic physics in the interaction of magnetic flux ropes,” Phys. Rev. Lett. 115, 175004 (2015), DOI:10.1103/PhysRevLett.115.175004.
  * Powell  _et al._ (1999) K. Powell, P. Roe, T. Linde, T. Gombosi, and D. L. De Zeeuw, “A solution-adaptive upwind scheme for ideal magnetohydrodynamics,” J. Comp. Phys. 154, 284–309 (1999).
  * Tóth  _et al._ (2012) G. Tóth, B. van der Holst, I. V. Sokolov, D. L. De Zeeuw, T. I. Gombosi, F. Fang, W. B. Manchester, X. Meng, D. Najib, K. G. Powell, Q. F. Stout, A. Glocer, Y.-J. Ma, and M. Opher, “Adaptive numerical algorithms in space weather modeling,” Journal of Computational Physics 231, 870–903 (2012).
  * Fonseca  _et al._ (2002) R. A. Fonseca _et al._ , “Osiris: A three-dimensional, fully relativistic particle in cell code for modeling plasma based accelerators,” Proceedings of The International Conference on Computer Science 2331, 342–351 (2002).
  * Hemker (2015) R. G. Hemker, “Particle-in-cell modeling of plasma-based accelerators in two and three dimensions,” arXiv:1503.00276. [Online]. Available: http://arxiv.org/abs/1503.00276 (2015).
  * Scudder and Daughton (2008) J. Scudder and W. Daughton, ““illuminating” electron diffusion regions of collisionless magnetic reconnection using electron agyrotropy,” Journal of Geophysical Research (Space Physics) 113, A06222 (2008), DOI:10.1029/2008JA013035.
  * Fadeev, Kvabtskhava, and Komarov (1965) V. M. Fadeev, I. F. Kvabtskhava, and N. N. Komarov, “An energy principle for hydromagnetic stability problems,” Nuclear Fusion 5, 202 (1965).
  * Daughton  _et al._ (2009a) W. Daughton, V. Roytershteyn, B. J. Albright, H. Karimabadi, L. Yin, and K. J. Bowers, “Influence of coulomb collisions on the structure of reconnection layers,” Physics of Plasmas 16, 072117 (2009a), DOI:10.1063/1.3191718.
  * Chen and Tóth (2019) Y. Chen and G. Tóth, “Gauss’s law satisfying energy-conserving semi-implicit particle-in-cell method,” Journal of Computational Physics 386, 632–652 (2019).
  * Dong  _et al._ (2021) C. Dong, A. Le, L. Wang, A. Stanier, B. Wetherton, W. Daughton, A. Bhattacharjee, J. Slavin, and G. DiBraccio, “Global Hybrid-VPIC Modeling of the Solar Wind Interaction with Mercury’s Dynamic Magnetosphere: Reconnection and Foreshock,” in  _AGU Fall Meeting Abstracts_ , Vol. 2021 (2021) pp. P25F–2211.
  * Winske  _et al._ (2022) D. Winske, H. Karimabadi, A. Le, N. Omidi, V. Roytershteyn, and A. Stanier, “Hybrid codes (massless electron fluid),” arXiv e-prints , arXiv:2204.01676 (2022), arXiv:2204.01676 [physics.plasm-ph] .
  * Makwana, Keppens, and Lapenta (2018) K. Makwana, R. Keppens, and G. Lapenta, “Two-way coupled MHD-PIC simulations of magnetic reconnection in magnetic island coalescence,” in  _Journal of Physics Conference Series_, Journal of Physics Conference Series, Vol. 1031 (2018) p. 012019.
  * Karimabadi  _et al._ (2011) H. Karimabadi, J. Dorelli, V. Roytershteyn, W. Daughton, and L. Chacón, “Flux pileup in collisionless magnetic reconnection: Bursty interaction of large flux ropes,” Physical Review Letters 107, 025002 (2011), DOI:10.1103/PhysRevLett.107.025002.
  * Daughton  _et al._ (2009b) W. Daughton, V. Roytershteyn, B. J. Albright, H. Karimabadi, L. Yin, and K. J. Bowers, “Transition from collisional to kinetic regimes in large-scale reconnection layers,” Physical Review Letters 103, 065004 (2009b), DOI:10.1103/PhysRevLett.103.065004.
  * Ji  _et al._ (2020) H. Ji, J. Yoo, J. Jara-Almonte, A. Goodman, Y. Ren, M. Yamada, A. Bhattacharjee, W. Fox, W. Daughton, A. Stanier, and The Flare Construction Team Team, “FLARE: a collaborative research facility to study magnetic reconnection and related phenomena,” in  _APS Division of Plasma Physics Meeting Abstracts_ , APS Meeting Abstracts, Vol. 2020 (2020) p. NO04.008.
  * Ebrahimi and Raman (2016) F. Ebrahimi and R. Raman, “Large-volume flux closure during plasmoid-mediated reconnection in coaxial helicity injection,” Nuclear Fusion 56, 044002 (2016).
  * Tanabe  _et al._ (2017) H. Tanabe, T. Yamada, T. Watanabe, K. Gi, M. Inomoto, R. Imazawa, M. Gryaznevich, C. Michael, B. Crowley, N. J. Conway, R. Scannell, J. Harrison, I. Fitzgerald, A. Meakins, N. Hawkes, K. G. McClements, T. O’Gorman, C. Z. Cheng, and Y. Ono, “Recent progress of magnetic reconnection research in the mast spherical tokamak,” Physics of Plasmas 24, 056108 (2017).
  * Olson  _et al._ (2021) J. Olson, J. Egedal, M. Clark, D. A. Endrizzi, S. Greess, A. Millet-Ayala, R. Myers, E. E. Peterson, J. Wallace, and C. B. Forest, “Regulation of the normalized rate of driven magnetic reconnection through shocked flux pileup,” Journal of Plasma Physics 87, 175870301 (2021).


# Hybrid Polar Encoding with Applications in Non-Coherent Channels

Mengfan Zheng
Department of Electronic and Computer Engineering
The Hong Kong University of Science and Technology
Email<EMAIL_ADDRESS>
###### Abstract

In coding theory, an error-correcting code can be encoded either
systematically or non-systematically. In a systematic encode, the input data
is embedded in the encoded output. Conversely, in a non-systematic code, the
output does not contain the input symbols. In this paper, we propose a hybrid
encoding scheme for polar codes, in which some data bits are systematically
encoded while the rest are non-systematically encoded. Based on the proposed
scheme, we design a joint channel estimation and data decoding scheme. We use
the systematic bits in the hybrid encoding scheme as pilots for channel
estimation. To mitigate the code rate loss caused by the pilots and to provide
additional error detecting capability, we propose a dynamic pilot design by
building connections between the systematic bits and non-systematic bits.
Simulation results show that the performance of the proposed scheme approaches
that of the traditional non-systematic polar coding scheme with perfect
channel state information (CSI) with the increase of SNR.

##  I Introduction

Polar codes [1] are the first family of capacity achieving channel codes and
have been adopted in the fifth generation wireless communication standards.
One of the key differences of polar codes from traditional channel codes is
their channel-specific design, indicating that the performance of polar codes
heavily relies on the knowledge of channel state information (CSI). Real-world
wireless communication channels are usually time-varying, and polar codes
designed for static channels may not be able to achieve the desired
performance.

The fading channel is a classic model for time-varying channels. In the fading
channel, the channel gain changes over time according to some distribution,
known as the channel distribution information (CDI). The block fading channel
models the case that the channel gain remains constant over a fixed time
interval $T_{c}$, called channel coherent time, and changes to another
independent value afterwards. Traditionally, the CSI is first estimated by
sending some known symbols, called pilots. Then data is encoded and
transmitted based on the estimated CSI. Such an approach is close to optimal
when the channel coherent time is very long. However, when the channel
coherent time is small, inserting pilots in the transmitted signal
significantly lowers the overall data rate. In this case, non-coherent
communication can be a better choice.

There have been some researches on polar coding for fading channels without
CSI in the literature. Polar coding for i.i.d. fading channels with only CDI
was considered in [2]. A multi-level polar coding scheme for non-coherent
block fading channels is proposed in [3], which is only suitable for block
fading channels with very small coherent time. Reference [4] proposed a
channel estimation scheme based on systematic polar codes by selecting pilot
symbols from the coded symbols. A list decoding scheme that can estimate the
channel coefficient and decode message simultaneously with the help of a few
pilots was proposed in [5], which requires a large list size to achieve good
performance. A pilot-free two-stage polar coding scheme that jointly estimates
the CSI and data based on the constraints imposed by the frozen bits was
proposed in [6], which can approach the performance of a coherent receiver
with perfect CSI.

All the aforementioned schemes are based on traditional non-systematic or
systematic polar codes. In this paper, we design a new polar coding scheme for
non-coherent communications over block fading channels based on a hybrid
encoding construction that combines non-systematic and systematic polar
encoding. In this construction, part of the message bits are encoded non-
systematically and the rest are encoded systematically. The idea of such a
hybrid polar encoding design has been proposed in [7] for the purpose of
reducing encoding and decoding resource consumption (called partial systematic
polar codes in that paper). In our proposed scheme, the systematically encoded
message bits serve as dynamic pilots which will be decoded first and used for
channel estimation. Then the rest message bits are decoded with the estimated
CSI. By exploiting the connection between the systematic and non-systematic
bits and with a properly designed two-phase decoder, we show by simulations
that the proposed scheme can approach the performance of the traditional non-
systematic polar coding scheme with perfect CSI as the SNR increases and can
even outperform it at high SNRs. We discuss the reason for this phenomenon in
details with our simulation results.

Notations: $[N]$ is the abbreviation for $\\{1,2,...,N\\}$. Vectors and
matrices are denoted by lowercase and uppercase boldface letters,
respectively. Vectors can also be denoted as
$x^{a:b}\triangleq[x_{a},x_{a+1},...,x_{b}]$ for $a\leq b$. For
$\mathcal{A}\subset[N]$ and $\mathcal{B}\subset[N]$, $x^{\mathcal{A}}$ denotes
the subvector $[x_{i}:i\in\mathcal{A}]$ of $x^{1:N}$,
$\mathbf{G}^{\mathcal{A}}$ is the submatrix of $\mathbf{G}$ consisting of rows
with indices in $\mathcal{A}$, and $\mathbf{G}^{\mathcal{A}\mathcal{B}}$ is
the submatrix of $\mathbf{G}$ consisting of elements $G_{ij}$, where
$i\in\mathcal{A}$ and $j\in\mathcal{B}$.

##  II Preliminaries

###  II-A System Model

We consider communications over a block fading channel with coherent time
$T_{c}$. Channel coding is performed across $B$ coherent blocks. The frame
size is then $N=BT_{c}$. At block $i$ ($i=1,...,B$), the channel is modelled
as

$\mathbf{y}^{(i)}=h_{i}\mathbf{x}^{(i)}+\mathbf{w}^{(i)},$ (1)

where $h_{i}\in\mathbb{R}$ is the channel gain at block $i$,
$\mathbf{x}^{(i)}=[x^{(i)}_{1},x^{(i)}_{2},...,x^{(i)}_{T_{c}}]$ is the
channel input,
$\mathbf{y}^{(i)}=[y^{(i)}_{1},y^{(i)}_{2},...,y^{(i)}_{T_{c}}]$ is the
channel output, and
$\mathbf{w}^{(i)}=[w^{(i)}_{1},w^{(i)}_{2},...,w^{(i)}_{T_{c}}]$ is the white
Gaussian noise, with $w^{(i)}_{j}\sim\mathcal{N}(0,\sigma^{2})$ for
$j\in[T_{c}]$. We assume BPSK modulation is used in this paper, i.e.,
$x^{(i)}_{j}\in\\{-1,1\\}$, although the proposed scheme can be easily
extended to higher order modulation scenarios. We study the non-coherent case
when neither the transmitter nor the receiver has the instantaneous CSI of the
channel. The CDI is assumed to be known by the receiver.

###  II-B Polar Codes

Polar codes are defined by the following polar transform [1]

$\mathbf{x}=\mathbf{u}\mathbf{G}_{N},$ (2)

where $\mathbf{u}=[u_{1},...,u_{N}]$ is the uncoded bit sequence,
$\mathbf{x}=[x_{1},...,x_{N}]$ is the encoded codeword and
$\mathbf{G}_{N}=\textbf{F}^{\otimes n}$. Here
$\textbf{F}=\begin{bmatrix}1&0\\\ 1&1\end{bmatrix}$ and $\otimes$ is the
Kronecker power. $\mathbf{G}_{N}$ can also be written as
$\mathbf{B}_{N}\textbf{F}^{\otimes n}$ with $\mathbf{B}_{N}$ being the bit-
reversal matrix. However, for systematic encoding it is more convenient to use
the generator matrix without $\mathbf{B}_{N}$.

####  II-B1 Non-Systematic Encoding

The non-systematic construction of polar codes partitions $\mathbf{u}$ into an
information set $\mathcal{I}$ and a frozen set $\mathcal{F}=\mathcal{I}^{c}$,
where $\mathcal{I}^{c}$ denotes the complementary set of $\mathcal{I}$.
Message bits are assigned to $\mathbf{u}^{\mathcal{I}}$ while
$\mathbf{u}^{\mathcal{F}}$ are assigned with some fixed value, such as 0. Then
the encoding process of (2) can be written as

$\mathbf{x}=u^{\mathcal{I}}\mathbf{G}_{N}^{\mathcal{I}}+u^{\mathcal{F}}\mathbf{G}_{N}^{\mathcal{F}}.$
(3)

####  II-B2 Systematic Encoding

In systematic encoding, message bits directly appear in $\mathbf{x}$. Suppose
$\mathbf{x}^{\mathcal{J}}$ ($\mathcal{J}\subset[N]$) are used to carry message
bits, where $|\mathcal{J}|=|\mathcal{I}|$. It is shown that if
$\mathbf{G}_{N}^{\mathcal{I}\mathcal{J}}$ is invertible [8], then

$u^{\mathcal{I}}=(x^{\mathcal{J}}-u^{\mathcal{I}^{c}}\mathbf{G}_{N}^{\mathcal{I}^{c}\mathcal{J}})(\mathbf{G}_{N}^{\mathcal{I}\mathcal{J}})^{-1}.$
(4)

Thus, to perform systematic polar encoding, one can first assign
$u^{\mathcal{I}^{c}}$ with the frozen value, compute $u^{\mathcal{I}}$
according to (4) and then obtain $\mathbf{x}$ according to (2). Note that the
encoding complexity of this approach is higher than that of non-systematic
encoding due to the matrix inversion operation. For low-complexity systematic
polar encoding algorithms, we refer the readers to [9, 10, 11].

####  II-B3 Decoding

Upon receiving $\mathbf{y}$, the receiver can use a successive cancellation
(SC) decoder to recover $\mathbf{u}$:

$\hat{u}_{i}=\begin{cases}u_{i},&\text{ if }i\in\mathcal{F}\\\
\arg\max_{u\in\\{0,1\\}}P_{U_{i}|\mathbf{Y},U^{1:{i-1}}}(u|\mathbf{y},\hat{u}^{1:{i-1}}),&\text{
if }i\in\mathcal{I}\end{cases}.$ (5)

For non-systematic encoding, the original information is retrieved from
$\hat{u}^{1:N}$. For systematic encoding, the original information is
retrieved from $\hat{x}^{1:N}=\hat{u}^{1:N}\mathbf{G}_{N}$.

To improve error performance, we can use the cyclic redundancy check
(CRC)-aided successive cancellation list (CA-SCL) decoding [12, 13]. The idea
of CA-SCL is to retain up to $L$ most probable paths during the SC decoding
process and use CRC to eliminate wrong paths.

###  II-C Hybrid Polar Encoding

In the hybrid polar encoding scheme, part of the message bits are mapped to
$\mathbf{u}$, while the rest are directly shown in $\mathbf{x}$. The hybrid
polar encoding goes as follows. Denote the non-systematic bit positions in
$\mathbf{u}$ by $\mathcal{I}_{ns}$ and the systematic bit positions in
$\mathbf{x}$ by $\mathcal{I}_{s}$. We assign frozen bits
$\mathbf{u}^{\mathcal{F}}$ to 0. Then we have

$\displaystyle x^{\mathcal{I}_{s}}$
$\displaystyle=u^{\mathcal{I}_{ns}}\mathbf{G}_{N}^{\mathcal{I}_{ns}\mathcal{I}_{s}}+u^{\mathcal{I}_{s}}\mathbf{G}_{N}^{\mathcal{I}_{s}\mathcal{I}_{s}},$
(6) $\displaystyle x^{\mathcal{I}_{ns}}$
$\displaystyle=u^{\mathcal{I}_{ns}}\mathbf{G}_{N}^{\mathcal{I}_{ns}\mathcal{I}_{ns}}+u^{\mathcal{I}_{s}}\mathbf{G}_{N}^{\mathcal{I}_{s}\mathcal{I}_{ns}},$
(7) $\displaystyle x^{\mathcal{F}}$
$\displaystyle=u^{\mathcal{I}_{ns}}\mathbf{G}_{N}^{\mathcal{I}_{ns}\mathcal{F}}+u^{\mathcal{I}_{s}}\mathbf{G}_{N}^{\mathcal{I}_{s}\mathcal{F}}.$
(8)

$x^{\mathcal{I}_{s}}$ and $u^{\mathcal{I}_{ns}}$ are used to carry
information. From (6) we have

$u^{\mathcal{I}_{s}}=(x^{\mathcal{I}_{s}}-u^{\mathcal{I}_{ns}}\mathbf{G}_{N}^{\mathcal{I}_{ns}\mathcal{I}_{s}})(\mathbf{G}_{N}^{\mathcal{I}_{s}\mathcal{I}_{s}})^{-1}.$
(9)

Then $x^{\mathcal{I}_{ns}}$ and $x^{\mathcal{F}}$ can be obtained according to
(9), (7) and (8).

##  III Joint Channel Estimation and Decoding with Hybrid Polar Codes

Figure 1: The proposed scheme.

Our proposed scheme is illustrated in Fig. 1. $\mathbf{s}$ is the message
vector to be sent and $\mathbf{v}$ is the output vector of the hybrid polar
encoder. The transmitted sequence in a frame is denoted by
$\mathbf{x}\triangleq[\mathbf{x}^{(1)},...,\mathbf{x}^{(B)}]$, which is an
interleaved version of $\mathbf{v}$. The received signal is denoted by
$\mathbf{y}\triangleq[\mathbf{y}^{(1)},...,\mathbf{y}^{(B)}]$. Let
$\mathbf{u}=\mathbf{v}\mathbf{G}_{N}$ and denote the non-systematic bit set,
systematic bit set and frozen bit set by $\mathcal{I}_{ns}$, $\mathcal{I}_{s}$
and $\mathcal{F}$, respectively. Suppose the code rate is $R$. Then the number
of message bits in each frame is $K=NR=BT_{c}R$. In our scheme, we use the
systematic bits to perform channel estimation. Thus, we also refer to them as
pilots in the sequel. The purpose of the interleaver is to distribute the
pilots uniformly among the coherent blocks. An example is shown in Fig. 2, in
which the pilots are mapped to the beginning of each coherent block.

Figure 2: Interleaver design.

A straightforward approach to apply hybrid polar coding for non-coherent
communications is to assign the systematic bits with some known value. The
decoder uses them to perform channel estimation and then recovers the message
embedded in the non-systematic bits with the estimated channel gains. The
drawback of this approach is obvious. Since the systematic bits does not carry
message, to maintain a target message transmission rate one has to increase
the rate of the polar code, which will deteriorate error performance.

###  III-A Dynamic Pilot Design

To utilise the systematic bits as pilots while not sacrificing the error
correcting capability, we propose the following dynamic pilot design. Let
$k_{s}=|\mathcal{I}_{s}|$ and denote
$\mathcal{I}=\mathcal{I}_{ns}\cup\mathcal{I}_{s}$. In the proposed scheme, we
require that the first $k_{s}$ indices in $\mathcal{I}$, denoted by
$\mathcal{I}^{(1:k_{s})}$, all belong to $\mathcal{I}_{ns}$. This is to ensure
that in an SCL decoder, the first $k_{s}$ non-frozen bits to be decoded (i.e.,
$u^{\mathcal{I}^{(1:k_{s})}}$) all belong to the non-systematic bits. At the
encoder side, a message vector $\mathbf{s}$ together with its CRC bits are
assigned to $u^{\mathcal{I}_{ns}}$. In the meantime, $v^{\mathcal{I}_{s}}$ is
chosen to be the same as the first $k_{s}$ bits in $u^{\mathcal{I}_{ns}}$,
denoted by $u^{\mathcal{I}_{ns}^{(1:k_{s})}}$, as illustrated in Fig. 3. In
this way, the values of pilots are not fixed, but dynamically change with the
message bits. The connection between the non-systematic and systematic bits
can be exploited by the decoder to perform channel estimation and also to
improve error performance.

At the receiver side, the decoder recovers the message in a two-phase manner.
In the first phase, it uses a CA-SCL decoder to find the most probable
$\hat{u}^{\mathcal{I}_{ns}^{(1:k_{s})}}$ with only the CDI of the channel.
Note that unlike conventional CA-SCL decoder, this decoding process terminates
as soon as the $k_{s}$-th non-frozen bit has been decoded, as we are only
interested in $u^{\mathcal{I}_{ns}^{(1:k_{s})}}$ in this phase. Since
$v^{\mathcal{I}_{s}}=u^{\mathcal{I}_{ns}^{(1:k_{s})}}$,
$\hat{u}^{\mathcal{I}_{ns}^{(1:k_{s})}}$ can be used for channel estimation.
This phase is referred to as the channel estimation phase in this paper. To
achieve lower decoding error rate in this phase,
$u^{\mathcal{I}_{ns}^{(1:k_{s})}}$ can consist of $k_{s}-m_{e}$
($m_{e}<k_{s}$) message bits and $m_{e}$ CRC bits of these message bits. In
the second phase, the decoder recovers the rest message bits with the
estimated channel gains and the knowledge of
$\hat{u}^{\mathcal{I}_{ns}^{(1:k_{s})}}$. We refer to this phase as the
message decoding phase.

Figure 3: Dynamic pilot design.

Besides playing the role of pilots for channel estimation, the systematic bits
in our proposed scheme can also serve as check bits in the message decoding
phase in addition to the CRC bits. For a decoding candidate
$\hat{\mathbf{u}}$, besides checking its validity using CRC, we can also
calculate $\hat{\mathbf{v}}=\hat{\mathbf{u}}\mathbf{G}_{N}$ and check whether
the corresponding $\hat{v}^{\mathcal{I}_{s}}$ coincides with the decoded
pilots $\hat{u}^{\mathcal{I}_{ns}^{(1:k_{s})}}$ in the channel estimation
phase. If not, this candidate is very likely to be wrong provided that
$\hat{u}^{\mathcal{I}_{ns}^{(1:k_{s})}}$ is correct.

###  III-B Encoding

Now we describe details of the encoding scheme. Suppose we use $m_{e}$ CRC
bits for the channel estimation phase and $m_{d}$ CRC bits for the message
decoding phase. The encoding procedures are summarized in Algorithm 1.

Algorithm 1 Encoding for Non-Coherent Communications

the message vector $\mathbf{s}$

the codeword to be transmitted $\mathbf{x}$

1\. The first $k_{s}-m_{e}$ message bits $s^{1:k_{s}-m_{e}}$ together with
their length-$m_{e}$ CRC bits $c_{e}^{1:m_{e}}$ are assigned to
$u^{\mathcal{I}_{ns}^{(1:k_{s})}}$ and $v^{\mathcal{I}_{s}}$.

2\. The rest message bits $s^{k_{s}-m_{e}+1:K}$ together with the
length-$m_{d}$ CRC bits $c_{s}^{1:m_{d}}$ of $\mathbf{s}$ are assigned to the
rest positions in $u^{\mathcal{I}_{ns}}$.

3\. The hybrid polar encoder’s output $\mathbf{v}$ is generated as described
in Sec. II-C.

4\. The final codeword $\mathbf{x}$ is generated with an interleaver $\pi$ by
$\mathbf{x}=\pi(\mathbf{v})$ as introduced previously.

###  III-C Decoding

In the channel estimation phase, the decoder uses a CA-SCL decoder with list
size $L_{e}$ to recover $u^{\mathcal{I}_{ns}^{(1:k_{s})}}$. Since the
correctness of the recovered pilots is vital for channel estimation and the
subsequent message decoding phase, we may choose a relatively large $L_{e}$.
As the number of pilots is much smaller than the length of the whole message,
the increase of complexity induced by choosing a large $L_{e}$ can be
controlled. The log-likelihood ratio (LLR) of $\mathbf{y}$ used in this
decoder is calculated assuming that the channel gain is constant and equals
$\mathbb{E}[h_{i}]$. Unlike conventional CA-SCL decoding, in this phase the
SCL decoding procedure terminates when the $k_{s}$-th non-frozen bit has been
decoded, and the most probable estimate
$\hat{u}^{\mathcal{I}_{ns}^{(1:k_{s})}}$ is selected with the help of CRC.
Then $\hat{u}^{\mathcal{I}_{ns}^{(1:k_{s})}}$ is used to estimate the channel
gains of each block, denoted by $\hat{h}_{i}$, with some standard channel
estimation method such as leat square (LS) or minimum mean square error
(MMSE).

In the message decoding phase, $u^{\mathcal{I}_{ns}^{(1:k_{s})}}$ is treated
as frozen bits with value $\hat{u}^{\mathcal{I}_{ns}^{(1:k_{s})}}$. The log-
likelihood ratios (LLRs) of $y_{j}$, $\pi^{-1}(j)\in\mathcal{I}_{s}$ (i.e.,
positions of pilots) are modified to infinity or minus infinity according to
the corresponding value of $\hat{u}_{i}$ ($i\in\mathcal{I}_{ns}^{(1:k_{s})}$).
The LLRs of the rest $y_{j}$ are calculated according to the estimated channel
gains. Then a conventional CA-SCL decoder with list size $L_{d}$ is used to
recover the rest message bits.

The overall joint channel estimation and decoding procedures are summarized in
Algorithm 2.

Algorithm 2 Joint Channel Estimation and Decoding

the channel output $\mathbf{y}$

the estimate of the transmitted message $\hat{\mathbf{s}}$

1\. Decode $\hat{u}^{\mathcal{I}_{ns}^{(1:k_{s})}}$ with $\mathbf{y}$ and the
CDI of the channel using a CA-SCL decoder of list size $L_{e}$.

2\. Estimate the channel gains $\hat{h}_{i}$, $i=1,...,B$ with
$\hat{u}^{\mathcal{I}_{ns}^{(1:k_{s})}}$.

3\. Update the LLRs of $y_{j}$ ($\pi^{-1}(j)\in\mathcal{I}_{s}$) by treating
them as known bits, update the LLRs of the rest $y_{j}$ according to the
estimated channel gains, set $u^{\mathcal{I}_{ns}^{(1:k_{s})}}$ as frozen bits
with value $\hat{u}^{\mathcal{I}_{ns}^{(1:k_{s})}}$, and then perform CA-SCL
decoding with list size $L_{d}$ to generate a list of candidates
$\mathcal{L}=\\{\hat{\mathbf{u}}_{i}\\}$, $i=1,...,L_{d}$.

4\. Find the most probable $\hat{\mathbf{u}}_{i}\in\mathcal{L}$, check whether
it can pass CRC. If yes, go to 5, otherwise eliminate it from $\mathcal{L}$
and repeat 4.

5\. Calculate $\hat{\mathbf{v}}_{i}=\hat{\mathbf{u}}_{i}\mathbf{G}_{N}$ and
see if $\hat{v}_{i}^{\mathcal{I}_{s}}=\hat{u}^{\mathcal{I}_{ns}^{(1:k_{s})}}$.
If not, eliminate $\hat{\mathbf{u}}_{i}$ from $\mathcal{L}$ and go back to 4,
otherwise use it to recover $\hat{\mathbf{s}}$ and terminate decoding.

##  IV Simulation Results

In this section, we compare the performance of the proposed scheme with
traditional polar codes. Since the frame error rate (FER) performance of
systematic and non-systematic polar codes are the same, we use the performance
of non-systematic polar codes as a benchmark. In our simulations, we assume
$h_{i}$ follows the Rayleigh distribution with PDF

$f(h_{i})=\frac{h_{i}}{\sigma_{h}^{2}}e^{-\frac{h_{i}^{2}}{2\sigma_{h}^{2}}},$
(10)

where $\sigma_{h}=1$. For traditional polar codes, we consider two cases,
i.e., the instantaneous CSI is known by the receiver and only the CDI is known
by the receiver. The decoding list size is 16 and the number of CRC bits is
16. For the proposed scheme, the decoding list size in the channel estimation
phase and that in the message decoding phase are both 16. 16 CRC bits are used
in total, with $m_{e}=12$ being allocated for the channel estimation phase and
$m_{d}=4$ for the message decoding phase. To see the decoding error rate of
the proposed scheme in the channel estimation phase, we also plot the error
rates when decoding $u^{\mathcal{I}_{ns}^{(1:k_{s})}}$ in the simulation
results.

Figure 4: Performance of the proposed scheme. $N=1024$, $R=1/2$, $T_{c}=64$,
$B=16$.

In the first example as shown in Fig. 4, the channel coherent time is
$T_{c}=64$ and the number of coherent blocks in a frame is $B=16$. Thus, the
frame length is $N=1024$. The number of systematic bits in the proposed scheme
is 48 (3 for each block). It can be seen that the performance of the proposed
scheme is similar to that of non-systematic polar codes with only CDI at low
SNRs, and gradually approaches and even outperforms that of non-systematic
polar codes with perfect CSI as the SNR increases.

It may look strange that the proposed scheme surpasses non-systematic polar
codes with perfect CSI at high SNRs. This can be explained as follows. As can
be seen from Fig. 4, the decoding error probability of pilots (denoted by
$P_{e}^{p}$ here) is much smaller than the overall FER. Suppose the pilots
have been decoded correctly and the estimated channel gains are almost perfect
(which is reasonable at high SNRs). Denote the error probability of the
message decoding phase by $P_{e}^{m}$. It can be expected that $P_{e}^{m}$ is
smaller than the FER of the non-systematic polar code with perfect CSI
(denoted by $P_{e}^{ns}$), since the two codes have the same rate but LLRs at
the pilot positions have been refined in the proposed scheme. The FER of the
proposed scheme can be upper bounded by

$\displaystyle P_{e}^{h}\leq P_{e}^{p}+(1-P_{e}^{p})P_{e}^{m}.$ (11)

When $P_{e}^{p}\ll P_{e}^{m}$, we have $P_{e}^{h}\approx
P_{e}^{m}<P_{e}^{ns}$. This explains why the proposed scheme is better than
non-systematic polar coding with perfect CSI at high SNRs.

We can also look at this problem from the viewpoint of concatenated codes. The
proposed scheme can be seen as a concatenation of traditional polar code with
a repetition code. Just like the concatenation with a CRC code can greatly
improve the performance of polar codes under SCL decoding, the concatenation
with a repetition code improves the performance of polar codes over non-
coherent block fading channels.

In the second example shown in Fig. 5, the channel coherent time is
$T_{c}=128$ and the number of coherent blocks in a frame is $B=8$. Thus, the
frame length is also $N=1024$. The number of systematic bits in the proposed
scheme is 48 (6 for each block). A similar result can be observed in this
example.

Figure 5: Performance of the proposed scheme. $N=1024$, $R=1/2$, $T_{c}=128$,
$B=8$.

##  V Conclusion

In this paper we have proposed a polar coding scheme for non-coherent
communications based on a hybrid encoding design that combines systematic and
non-systematic polar coding. Connections between the non-systematic and
systematic bits are exploited to perform joint channel estimation and
decoding. Simulations have shown that the proposed scheme can even achieve
better performance than traditional non-systematic polar codes with perfect
CSI over some block fading channels at high SNRs.

Although we have only considered BPSK modulation in this paper for the brevity
of our description, the idea can be readily extended to higher order
modulation and/or MIMO cases. We will investigate these scenarios in the
future journal version of this work.

## References

  * [1] E. Arıkan, “Channel polarization: A method for constructing capacity-achieving codes for symmetric binary-input memoryless channels,” _IEEE Transactions on Information Theory_ , vol. 55, no. 7, pp. 3051–3073, 2009.
  * [2] A. Bravo-Santos, “Polar codes for the rayleigh fading channel,” _IEEE Communications Letters_ , vol. 17, no. 12, pp. 2352–2355, 2013.
  * [3] M. Zheng, W. Chen, and C. Ling, “Polar coding for noncoherent block fading channels,” in _2018 10th International Conference on Wireless Communications and Signal Processing (WCSP)_ , 2018, pp. 1–5.
  * [4] L. Li, Z. Xu, and Y. Hu, “Channel estimation with systematic polar codes,” _IEEE Transactions on Vehicular Technology_ , vol. 67, no. 6, pp. 4880–4889, 2018.
  * [5] M. Xhemrishi, M. C. Coşkun, G. Liva, J. Östmanl, and G. Durisi, “List decoding of short codes for communication over unknown fading channels,” in _2019 53rd Asilomar Conference on Signals, Systems, and Computers_ , 2019, pp. 810–814.
  * [6] P. Yuan, M. C. Coşkun, and G. Kramer, “Polar-coded non-coherent communication,” _IEEE Communications Letters_ , vol. 25, no. 6, pp. 1786–1790, 2021.
  * [7] H. Jin and R. Liu, “Partial systematic polar coding,” in _International Conference on Communications and Networking in China_. Springer, 2018, pp. 137–146.
  * [8] E. Arıkan, “Systematic polar coding,” _IEEE Communications Letters_ , vol. 15, no. 8, pp. 860–862, 2011.
  * [9] G. Sarkis, I. Tal, P. Giard, A. Vardy, C. Thibeault, and W. J. Gross, “Flexible and low-complexity encoding and decoding of systematic polar codes,” _IEEE Transactions on Communications_ , vol. 64, no. 7, pp. 2732–2745, 2016.
  * [10] H. Vangala, Y. Hong, and E. Viterbo, “Efficient algorithms for systematic polar encoding,” _IEEE Communications Letters_ , vol. 20, no. 1, pp. 17–20, 2016.
  * [11] Chen _et al._ , “A low complexity encoding algorithm for systematic polar codes,” _IEEE Communications Letters_ , vol. 20, no. 7, pp. 1277–1280, 2016\.
  * [12] K. Niu and K. Chen, “CRC-aided decoding of polar codes,” _IEEE Communications Letters_ , vol. 16, no. 10, pp. 1668–1671, 2012.
  * [13] I. Tal and A. Vardy, “List decoding of polar codes,” _IEEE Transactions on Information Theory_ , vol. 61, no. 5, pp. 2213–2226, 2015.


# Formation of spiral dwarf galaxies: observational data and results of
numerical simulation 111Published: Khrapov S.S., Khoperskov A.V., Zaitseva
N.A., Zasov A.V., Titov A.V., Formation of spiral dwarf galaxies:
observational data and results of numerical simulation // St. Petersburg State
Polytechnical University Journal. Physics and Mathematics, 2023, v.16 (1.2),
pp.395–402. DOI: https://doi.org/10.18721/JPM.161.260
https://physmath.spbstu.ru/en/article/2023.64.60/

Sergey Khrapov1 Alexander Khoperskov1 Natalia Zaitseva3 Anatoly Zasov2,3
Alexander Titov1 1 Institute of Mathematics and Information Technology,
Volgograd State University, Volgograd, Russia
2 Lomonosov Moscow State University, Faculty of Physics, Moscow, Russia
3 Lomonosov Moscow State University, Sternberg Astronomical Institute, Moscow,
Russia<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Recent studies show the possibility of the formation of fairly regular and
global spiral patterns in dwarf galaxies (dS type). Our sample of observed
dwarf objects of this class also includes galaxies with a central stellar bar.
The analysis of the observational data provides a small rotation velocity and
a small disk component mass for dS galaxies, which is in poor agreement with
the spiral structure generation mechanism in isolated dwarfs due to the
development of disk gravitational instability. Numerical simulation of the
stellar–gaseous disks self-consistent dynamics imposes restrictions on the
stellar disk thickness and the maximum gas rotation velocity, at which the
gravitational mechanism of spiral formation can still be effective.

##  1 Introduction

Dwarf galaxies are small in size and mass compared to classical spiral
galaxies (S or SB types) and are usually considered structureless, irregular
objects (Irr). Some late-type dwarfs (Sd – Sm types) have a rotating stellar
disk without any regular and developed spiral structure. Such galaxies exhibit
flocculent type of spirals, which are discontinuous and consist of short
regions [1]. Gravitational instability in large massive disks is able to
provide a regular spiral pattern covering the entire disk, so Grand Design
spiral structure is common for Sa – Sc galaxy types [2, 3, 4, 5, 6].

Only a small part of dwarf galaxies shows a global, relatively regular spiral
pattern in their disk, and such objects belong to the fairly rare dS type [7,
8, 9]. The observations comparative analysis of normal S- and dS-galaxies
represents that such dwarfs are more than just smaller copies of large
objects, since their spectral characteristics are similar to Irr galaxies
[10]. The small size and mass of dwarf galaxies appear to be a theoretical
problem for the formation of extended spirals due to gravitational instability
[5, 6].

Here we consider the observed properties of the sample of dS galaxies compared
to dwarf galaxies without a regular spiral structure [7]. The numerical
simulation of the dynamics of dwarfs stellar disks with a rich gaseous
component makes it possible to determine the conditions under which
gravitational instability can generate sufficiently extended spiral patterns
in dS galaxies, which morphology is similar to Grand Design galaxies.

Figure 1:  Dwarf galaxies with a spiral pattern from our sample.

##  2 Sample of dS-galaxies and its properties

Our sample is limited to objects (usually type Sc – Irr) with the absolute
magnitude $M_{B}>-18^{m}$, the optical diameter $D_{25}<12$ kpc,
$m_{B}<15^{m}$, the inclination angle $i<75^{\circ}$ [7]. It is also important
to have images in different spectral ranges for deeper study (Figures 1, 2).
The logarithm of isolation index $\log(ii)$ characterizes the degree of
environmental influence [11] and we do not consider both Virgo, UMa, Fornax
clusters and peculiar galaxies with signs of strong interaction.

Figure 2:  Dwarf galaxies with a spiral pattern from our sample (continued).

Our sample of spiral dwarf galaxies includes 43 objects, which are compared
with the sample of dwarfs without spirals (119 objects of Sm and Irr types,
see detailed description in [7]). Figure 1 shows dwarf galaxies with bars and
two distinct arms. The images of SDSS 9 (the Sloan Digital Sky Survey), DECaLS
(the Dark Energy Camera Legacy Survey), DSS (the Digitized Sky Survey),
PanSTARRS (the Panoramic Survey Telescope and Rapid Response System), LEGA
(the DESI Legacy Imaging Surveys) [12] characterize the distributions of the
stellar components. The bottom rows in the fig. 1, 2 show the distributions of
gas and young stars according to GALEX data (the Galaxy Evolution Explorer).
Figure 2 demonstrates galaxies with a more complex spiral structure, where the
yellow lines highlight the positions of the spiral arms. Moreover, the spirals
in the stellar and gaseous components are in good agreement with each other.
Three-arm patterns indicate a rather massive dark halo.

Figure 3:  Positions of various dwarf galaxies: blue icons are spirals, green
squares are Sm-type, pink squares are irregular objects. a) In-plane
distribution of the “baryon mass – momentum” parameters compared to the
regression for isolated galaxies of the AMIGA sample [13] (black line),
colored lines show deviation 1$\sigma$,
$M_{barionic}=\Upsilon^{*}_{K}L_{K}+\eta M_{HI}$, where $\Upsilon^{*}_{K}=0.6$
[14] and $\eta=1.33$. b) Tully-Fisher relation for the K-band compared to
those obtained by Karachentsev et al. [15] for the sample of Local Volume
dwarf galaxies (black line, colored lines show 1$\sigma$ deviation).

Each galaxy in our two samples is characterized by the systemic velocity
$V_{sys}$, the diameter $D_{25}$, the maximum rotation velocity $V_{rot}$, the
HI mass $M_{HI}$, the estimate of the total gravitating mass $M_{dyn}$, the
luminosities $L_{K}$ according to the K-magnitudes of 2MASS catalog [16].
Statistical analysis gives the following conclusions (See also [7]).

— Dwarf galaxies with developed spiral structure are the most massive objects
in the sample.

— The distributions of dS galaxies and objects without spirals indicate the
absence of very significant differences for various pairs of parameters, for
example, $L_{K}-M_{dyn}$, $L_{K}-M_{HI}$, $D_{25}-M_{HI}$, $M_{dyn}-M_{HI}$,
$V_{rot}-M_{HI}$ and others (Figure 3).

— The HI mass in dS galaxies is, on average, about two times less than in
dwarf non-spiral galaxies, although the dynamical and photometric parameters
are close for both samples.

— The central stellar bar is found both in dS-type objects and in non-spiral
galaxies.

— Tidal influence is apparently not an essential factor of the formation in
considered galaxies.

— The proportion of baryonic matter in spiral dwarfs is, on average, lower
than in objects with irregular structure and in giant spiral galaxies, which
may indicate the influence of the dark halo on the formation of the spiral
patterns in dS-dwarfs.

— Remote dwarf galaxies follow the same Tully-Fisher relation as Local Volume
dwarfs, but their physical parameters are determined with a larger
uncertainties due to the low brightness of these objects, which significantly
increases the points scatter on the diagram.

Figure 4:  The dispersion dynamics of stellar vertical velocities in disk
models for different gravitational potential cutoff radii ($r_{c}$), the
number of particles in the disk ($N^{(*)}$) and the live halo ($N^{(h)}$). a)
$N^{(*)}=N^{(h)}=2^{18}$, $r=10^{-5}$ (red lines), $r=4\cdot 10^{-3}$ (green
lines), for different radii ($r_{d}$ is the radial exponential scale of
stellar disk); b) $N^{(*)}=N^{(h)}=2^{20}$, $r=10^{-5}$ (red lines), $r=4\cdot
10^{-3}$ (green lines), for different radii.

##  3 Numerical modeling of the galactic disk

The numerical models are based on the self-consistent dynamics of the N-body
gravitating system for the stellar disk and the gaseous component, which is
described by the hydrodynamic equations [3, 7]. We used direct method to
calculate self-gravity forces (each particle interacts with each other), which
is the most accurate modeling approach [17]. The GPUs application for parallel
code makes it possible to perform fairly fast numerical experiments with
$N=2^{20}-2^{23}$. The numerical model should ensure the collisionlessness of
the stellar component [7, 18], which is achieved by cutting off the Newtonian
potential at small radii $r_{c}$. Surface density of the stellar exponential
disk $\sigma(r)=\sigma_{0}\,\exp(-r/r_{d})$ characterized by the radial scale
$r_{d}$. We use dimensionless units of length $r=1\rightarrow 4$ kpc, velocity
$V=1\rightarrow 47$ km sec-1 and time $t=1\rightarrow 80$ Myr.

Figure 4 demonstrates the presence or absence of the stellar disk heating for
the corresponding values of $r_{c}$ and $N=N^{(*)}+N^{(h)}$ (where $N^{(*)}$
is the number of particles in the stellar disk, $N^{(h)}$ is the number of
particles that form dark live halo).

We see a noticeable linear increase in vertical velocity dispersion at a very
small cutoff radii due the absence of collisionlessness in such a model. The
value $r_{c}=0.004$ ensures the almost stationary behavior of the velocity
dispersions (green lines). The heating is stronger for a smaller values of the
number of particles $N^{(*)}$ and $N^{(h)}$ (comparison of Figure 4a and
Figure 4b).

Figure 5:  Examples of spiral patterns in numerical models of galaxies PGC
32390, PGC 29549, PGC 16617. Top row: stellar disk, bottom row: gaseous disk.

The figure 5 shows various spiral structures in dwarf galaxies models. We
obtain different stellar and gaseous disks morphology in numerical models by
varying the relative masses of stars, gas, and dark halo, as well as the
radial and vertical scales that determine the distributions of subsystem
parameters. The model structures in gaseous and stellar disks are close to the
observed patterns of the corresponding galaxies (See Figures 1, 2). Other
calculation examples are given in [7].

##  4 Conclusion

Figure 6:  Distribution of the Fourier harmonics maximum amplitude $A$ on the
plane of dimensionless parameters (where $M_{g}/M_{*}$ is the relative gas
mass, $z_{0}$ is the vertical disk scale, $r_{0}$ is the radial disk scale).

Photometric and kinematic observational data do not allow us to confidently
identify the factors that ensure the formation of a global spiral patterns in
dwarf galaxies, which are a rather rare phenomenon. There is only some gas
deficit in dS galaxies compared to dIrr objects.

We have studied the possibility of the global spiral structure formation in
numerical models of isolated dwarf galaxies due to the development of
gravitational instability in the stellar disk rich in gas. The presence of a
spiral pattern in dwarf models imposes some restrictions on the disks
thickness, the radial velocity dispersion profiles in stellar disk, the sound
speed in gas and the gas density. The results of numerical simulations show
that the maximum gas rotation velocity must be higher than 60 km sec-1 in
order to excite spiral waves with significant amplitude. Thicker stellar disk
requires more gas to form the spiral pattern.

### Acknowledgments

This research has made use of “Aladin sky atlas” developed at CDS, Strasbourg
Observatory, France. The work was supported by the Ministry of Science and
Higher Education of the Russian Federation (state assignment, project No.
0633-2020-0003, implementation of all numerical simulations) and by Russian
Foundation for Basic Research (project 20-02-00080 A, observational data
analysis).

## References

## References

  * [1] Mondal C, Subramaniam A, George K, Postma J E, Subramanian S and Barway S. 2021 _Astrophysical Journal_ 909 203
  * [2] Kormendy J 2011 _Secular evolution in disk galaxies_ ed J Falcón-Barroso and J H Knapen (Cambridge: Cambridge University Press)
  * [3] Khrapov S, Khoperskov A and Korchagin V 2021 _Galaxies_ 9 29
  * [4] Griv E, Jiang I G and Hou L G 2017 _Monthly Notices of the Royal Astronomical Society_ 468 3361–-67
  * [5] Dobbs C and Baba J 2014 _Publications of the Astronomical Society of Australia_ 31 e035
  * [6] Buta R 2013 Galaxy morphology _Secular Evolution of Galaxies_ ed J Falcón-Barroso and J H Knapen (Cambridge: Cambridge University Press) pp 155–258
  * [7] Zasov A V, Khoperskov A V, Zaitseva N A and Khrapov S S 2021 _Astronomy Reports_ 65 1215–32
  * [8] Magana-Serrano M A, Hidalgo-Gamez A M, Vega-Acevedo I and Castaneda H O 2020 _Revista Mexicana de Astronomia y Astrofisica_ 56 39–54
  * [9] Edmunds M G and Roy J R 1993 _Monthly Notices of the Royal Astronomical Society_ 261 L17-L19
  * [10] Hidalgo-Gamez A M 2004 _Revista Mexicana de Astronomía y Astrofísica_ 40 37–51
  * [11] Makarov D and Karachentsev I 2011 _Monthly Notices of the Royal Astronomical Society_ 412 2498–2520.
  * [12] Dey A et al. 2019 _Astronomical Journal_ 157 168
  * [13] Lisenfeld U et al. 2007 _Astronomy and Astrophysics_ 462 507–523
  * [14] McGaugh S S and Schombert J M 2014 _The Astronomical Journal_ 148 77
  * [15] Karachentsev I D, Kaisina E I, Kashibadze N and Olga G 2017 _The Astronomical Journal_ 153 6
  * [16] Huchra J P et al 2012 _Astrophys. J. Suppl. Ser._ 199 26
  * [17] Khrapov S S, Khoperskov S A and Khoperskov A V 2018 _Bulletin of the South Ural State University, Series: Mathematical Modelling, Programming and Computer Software_ 11 124–-136
  * [18] Smirnov A A and Sotnikova N Ya 2018 _Astronomical Society_ 481 4058–76