# Droplet formation simulation using mixed finite elements111This article may
be downloaded for personal use only. Any other use requires prior permission
of the author and AIP Publishing. This article appeared in [1] and may be
found at https://doi.org/10.1063/5.0089752

Darsh K. Nathawani222Email address for correspondence<EMAIL_ADDRESS>Computational and Data-Enabled Science and Engineering, University at Buffalo,
Buffalo, NY 14260, USA  Matthew G. Knepley Department of Computer Science and
Engineering, University at Buffalo, Buffalo, NY 14260, USA

###### Abstract

Droplet formation happens in finite time due to the surface tension force. The
linear stability analysis is useful to estimate droplet size but fails to
approximate droplet shape. This is due to a highly non-linear flow description
near the point where the first pinch-off happens. A one-dimensional
axisymmetric mathematical model was first developed by Eggers and Dupont[2]
using asymptotic analysis. This asymptotic approach to the Navier-Stokes
equations leads to a universal scaling explaining the self-similar nature of
the solution. Numerical models for the one-dimensional model were developed
using the finite difference[2] and finite element method[3]. The focus of this
study is to provide a robust computational model for one-dimensional
axisymmetric droplet formation using the Portable, Extensible Toolkit for
Scientific Computation (PETSc). The code is verified using the Method of
Manufactured Solutions (MMS) and validated using previous experimental studies
done by Zhang and Basaran[4]. The present model is used for simulating pendant
drops of water, glycerol, and paraffin wax, with an aspiration of extending
the application to simulate more complex pinch-off phenomena.

Droplet formation simulation using mixed finite elements111This article may be
downloaded for personal use only. Any other use requires prior permission of
the author and AIP Publishing. This article appeared in [1] and may be found
at https://doi.org/10.1063/5.0089752

##  1 Introduction

Singularity in free surface flows is a crucial problem for time-accurate
simulations. A mathematical treatment, built on the self-similarity of the
pinch-off region, has made numerical simulations tractable. Our study is
concerned with the formation of pendant droplets. Rayleigh was the first to
demonstrate that droplet formation occurs in finite time due to the force of
surface tension acting against inertia [5]. The pinch-off dynamics of a
pendant drop has received the most attention from both mathematicians and
scientists. Considering the pendant drop as a fluid column, the surface
tension force and the gravitational force are initially balanced. The drop
becomes heavier as more fluid is added from the top. Eventually gravity
extends the drop, increasing the surface energy. In order to minimize this
energy, the radius of the fluid column shrinks, forming a neck region, what we
see as the action of surface tension. At some finite time, the radius becomes
zero at some location and the drop separates from the fluid column. This
location in the fluid column is called the pinch-off point or singularity.

Fluid motion in the immediate vicinity of the singularity is driven by very
high velocity gradients generated by the surface tension, inertial, and
viscous forces. In fact, the solution in this region is self-similar, meaning
that it does not depend on the initial or boundary conditions, but has a
universal character. After the first pinch-off, a long neck recoils back with
high velocity. This induces surface perturbations in the column and can lead
to further breakup into smaller satellite droplets, a phenomenon also observed
in liquid bridges and decaying jets. A now-classic treatment of the governing
dynamics of singularities as well as analysis of the self-similarity for these
cases is given by Eggers [6, 7]. Linear stability analysis can accurately
approximate droplet size, but fails to approximate the shape of a droplet [8].
Moreover, even the higher order analysis is not able to explain the shape of
the drop near singularity [8, 9]. One-dimensional analysis has proved useful
for circular liquid jets [10, 11] and pendant drops [12]. On the experimental
side, studies have examined the pinch-off dynamics of pendant drops [13, 14],
as well as characterized droplet dynamics in terms of non-dimensional
parameters [4].

A full one-dimensional mathematical model was constructed by Eggers and Dupont
[2] using the asymptotic expansion of the Navier-Stokes equations in
cylindrical coordinates. They used a finite difference scheme to discretize
the equations, and simulated both a pendant drop and a decaying jet. They
verified that one-dimensional treatment can accurately simulate the pinch-off
dynamics for a fluid in a quiescent background. However, their computational
approach could only simulate up to the first pinch-off. Other computational
models using one and two-dimensional analysis were explored by Ambravaneswaran
et al. using finite elements to discretize the problem [3]. They investigated
the effect of volume flow rate on droplet and neck shape. Moreover, they were
able to simulate satellite drops. However, they were not able to validate the
satellite drop simulations, as they had the primary droplet. Their two-
dimensional simulations support the conclusion that for axisymmetric droplet
formation, one-dimensional computational models are much faster and reliably
accurate. Ambravaneswaran et al. also propose a hybrid 1D-2D computational
model, and matching between 1D and 2D domains is also explored [15]. Other
numerical approaches like the Volume-of-Fluid (VOF) method were tested, but
either failed to accurately capture features of the flow, such as micro-
threads, or were limited to only a certain range of fluid parameters, such as
viscosity [16, 17]. Additionally, the VOF method requires at least a 2D
domain, which increases the computational cost.

Formation of droplets can be seen in many scientific and industrial
applications, a few examples being ink-jet printing [18], spray cooling [19],
and droplet entrainment in annular flow [20]. Moreover, droplet formation has
a pivotal role in fuel entrainment and burning in hybrid rockets [21], which
is the motivation to pursue this study. We present a computational model that
is both verified with MMS and validated with previous experimental work. We
extend the results with paraffin wax simulations to build a base for the
future work on droplets in a shear force environment.

##  2 Mathematical and Computational Model

In this section, we consider the Navier-Stokes momentum equation in
cylindrical coordinates for an axisymmetric fluid column, as treated by Eggers
and Dupont [2]. The fluid is considered incompressible with density $\rho$ and
kinematic viscosity $\nu$. Assuming no swirling motion, we consider the flow
only in radial and axial directions. A schematic of a pendant drop in
cylindrical coordinates is shown in Fig. 1. The surface of the droplet, which
is defined as variable $h(z,t)$, is moving with the velocity. Therefore, the
model equation for $h$ is given by

$\displaystyle\frac{\partial h}{\partial t}+u_{z}\frac{\partial h}{\partial
z}=u_{r}|_{r=h}$

Approaching the pinch-off point, the radial contraction is much faster than
the axial expansion. Therefore, considering the radius $r$ as an asymptotic
parameter, the axial velocity ($u_{z}$) and pressure ($p$) are expanded
asymptotically in even order terms to satisfy symmetry. Then using the
continuity equation, the radial velocity is derived.

$\displaystyle u_{z}$ $\displaystyle=u_{0}+u_{2}r^{2}+\ldots,$ $\displaystyle
u_{r}$ $\displaystyle=-\frac{\partial u_{0}}{\partial
z}\frac{r}{2}-\frac{\partial u_{2}}{\partial z}\frac{r^{3}}{4}-\ldots,$
$\displaystyle p$ $\displaystyle=p_{0}+p_{2}r^{2}+\ldots$

To introduce the surface tension force ($\gamma$) in the governing equations,
the following force balance is considered.

$\displaystyle\mathbf{\hat{n}}\sigma\mathbf{\hat{n}}$
$\displaystyle=-\gamma\left(\nabla\cdot\mathbf{\hat{n}}\right)$
$\displaystyle\mathbf{\hat{n}}\sigma\mathbf{\hat{t}}$ $\displaystyle=0$

Here, $\sigma$ is a stress tensor, $\mathbf{\hat{n}}$ is a unit outward
normal, $\mathbf{\hat{t}}$ is a unit tangent, and
$\nabla\cdot\mathbf{\hat{n}}$ is a mean curvature. The above force balance
explains that the normal stress is balanced by the surface tension and the
tangential stress is zero.

Using the force balance and the leading order terms in $r$ from the expansion,
we simplify the momentum equation. The advecting surface equation is already
in leading order. Dropping the subscripts, the governing equations for a one-
dimensional axisymmetric fluid column are given by

$\displaystyle\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial
z}+\frac{\gamma}{\rho}\frac{\partial(\nabla\cdot\mathbf{\hat{n}})}{\partial
z}-\frac{3\nu}{h^{2}}\frac{\partial}{\partial z}\left(h^{2}\frac{\partial
u}{\partial z}\right)-g$ $\displaystyle=0$ (1) $\displaystyle\frac{\partial
h}{\partial t}+u\frac{\partial h}{\partial z}+\frac{h}{2}\frac{\partial
u}{\partial z}$ $\displaystyle=0$ (2)

where, the mean curvature term $\nabla\cdot\mathbf{\hat{n}}$ is given by

$\displaystyle\nabla\cdot\mathbf{\hat{n}}=\left[\frac{1}{h\left(1+\frac{\partial
h}{\partial z}^{2}\right)^{1/2}}-\frac{\frac{\partial^{2}h}{\partial
z^{2}}}{\left(1+\frac{\partial h}{\partial z}^{2}\right)^{3/2}}\right]$ (3)

Here, the curvature term is not approximated to the leading order because it
has been shown that using the full curvature term better captures the
singularities [7].

Figure 1: Schematic of a pendant droplet.

Equations (1) and (2) govern a pendant drop under gravitational forcing. These
are solved for $u$ and $h$ using a finite element discretization. However, the
highest order derivative is of the third order, which is problematic for our
$C^{0}$ continuous element scheme. The approximation for this term will be
discontinuous across element interfaces. We could handle this using a
discontinuous Galerkin (DG) scheme, but instead, we choose a mixed-element
formulation, inspired by Ambravaneswaran et al. [3], in which we explicitly
discretize the axial derivative of the radius $h$ (or slope), $s=\partial
h/\partial z$, so that

$\displaystyle\nabla\cdot\mathbf{\hat{n}}=\left[\frac{1}{h\left(1+s^{2}\right)^{1/2}}-\frac{\frac{\partial
s}{\partial z}}{\left(1+s^{2}\right)^{3/2}}\right]$ (4)

The mixed finite element formulation is given by

$\displaystyle\int_{\Omega}q\left[\frac{\partial u}{\partial
t}+u\frac{\partial u}{\partial
z}+\frac{\gamma}{\rho}\frac{\partial(\nabla\cdot\mathbf{\hat{n}})}{\partial
z}-\frac{3\nu}{h^{2}}\frac{\partial}{\partial z}\left(h^{2}\frac{\partial
u}{\partial z}\right)-g\right]d\Omega$ $\displaystyle=0$ (5)
$\displaystyle\int_{\Omega}v\left[\frac{\partial h}{\partial
t}+u\frac{\partial h}{\partial z}+\frac{1}{2}h\frac{\partial u}{\partial
z}\right]d\Omega$ $\displaystyle=0$ (6)
$\displaystyle\int_{\Omega}w\left[s-\frac{\partial h}{\partial
z}\right]d\Omega$ $\displaystyle=0$ (7)

Here, $q$, $v$, and $w$ are test functions, and the mean curvature is defined
by Eq. (4). The third and fourth terms in Eq. (5) are simplified by performing
integration by parts, which makes the weak form with the highest order
derivative to first order. The equations after integration by parts take the
following form.

$\displaystyle\int_{\Omega}q\left[\frac{\partial u}{\partial
t}+u\frac{\partial u}{\partial z}-\frac{6\nu}{h}\frac{\partial h}{\partial
z}\frac{\partial u}{\partial
z}+\frac{\gamma}{\rho}\left\\{-\frac{s\frac{\partial s}{\partial
z}}{h\left(1+s^{2}\right)^{3/2}}-\frac{s}{h^{2}\left(1+s^{2}\right)^{1/2}}\right\\}-g\right]d\Omega$
$\displaystyle+\int_{\Omega}\nabla q\left[3\nu\frac{\partial u}{\partial
z}+\frac{\gamma}{\rho}\frac{\frac{\partial s}{\partial
z}}{\left(1+s^{2}\right)^{3/2}}\right]d\Omega-\int_{\Gamma}\nabla
q\left[3\nu\frac{\partial u}{\partial
z}+\frac{\gamma}{\rho}\frac{\frac{\partial s}{\partial
z}}{\left(1+s^{2}\right)^{3/2}}\right]d\Omega$ $\displaystyle=0$ (8)
$\displaystyle\int_{\Omega}v\left[\frac{\partial h}{\partial
t}+u\frac{\partial h}{\partial z}+\frac{1}{2}h\frac{\partial u}{\partial
z}\right]d\Omega$ $\displaystyle=0$ (9)
$\displaystyle\int_{\Omega}w\left[s-\frac{\partial h}{\partial
z}\right]d\Omega$ $\displaystyle=0$ (10)

Initially, the velocity is zero and the curvature profile is a hemisphere as
it minimizes surface energy. The inlet radius $h_{0}$ is fixed depending on
the nozzle radius and the inflow velocity $u_{0}$ is constant. The radius at
the tip of the droplet, at length $L(t)$, is zero. The set of Eq. (8)-(10) are
then solved using a continuous Galerkin formulation subject to the following
constraints.
Initial conditions:

$\displaystyle h$ $\displaystyle=\sqrt{h_{0}^{2}-z^{2}}$ $\displaystyle s$
$\displaystyle=-\frac{z}{\sqrt{h_{0}^{2}-z^{2}}}\ \text{for}\ (0\leq
z<L_{0}),\qquad s|_{L_{0}}=-C$ $\displaystyle u$ $\displaystyle=0$

where C is a large negative number. In our implementation, we use $-10$.
However, the code was tested with larger values and the results were
unchanged.
Boundary conditions:

$z=0$ | $h=h_{0}$ | $u=u_{0}$
---|---|---
$z=L(t)$ | $h=0$ | $u=\frac{dL}{dt}$

The length of the drop L(t) can be calculated as a part of the solution as
explained by Ambravaneswaran et al. [3] by calculating the volume of the drop,
which can then be used to calculate the velocity at the tip. However, this
results in a dense row in the Jacobian, so we instead produce $L(t)$ by self-
consistent iteration.

Initially, we are given $u(t)$ and $h(t)$, including the velocity at the end
of the droplet. We use that velocity to predict $L(t+dt)$ that is
$L(t)+u_{tip}*dt$, giving us our boundary condition, and we extend our mesh to
this length. We then solve our existing system (Eq. (8)-(10)) for $u(t+dt)$
and $h(t+dt)$. This allows us to calculate the droplet volume by integrating
$h(t+dt)$ along the length. This must match the volume from the last time-step
augmented by the amount of liquid flowing, which is $4\pi h^{2}_{0}u_{0}dt$.
The difference between the calculated volume and the theoretical volume is
used for self-consistency in adjusting the length $L(t+dt)$. We use bisection
to arrive at a consistent length $L$ for this new time step. This adaptation
loop is done when the conservation of volume is satisfied to a given tolerance
(we use 0.1%).

The one-dimensional mesh, representing the domain $0\leq z\leq L(t)$, moves as
the length $L(t)$ changes. We first update the position of the last vertex and
then move the remaining vertices to even out the cell lengths. The
interpolation of the discrete field representation between these two meshes
can be achieved using the Galerkin projection [22]. Galerkin projection is
optimal for the $L_{2}$ norm, which measures energy. Alternatively, we could
replace this with a volume constraint during interpolation. As shown in Fig.
2, the re-meshing is done using the calculated $L(t)$ between each time step.
The Galerkin projection is then used to interpolate the solution to the new
mesh. For instance, if the solution on the old mesh is $u^{old}$ and on the
new mesh is $u^{new}$, then the interpolation is done as follows:

$\displaystyle u^{new}$ $\displaystyle=\eta^{new}_{i}\left(u^{old}\right)$
$\displaystyle u^{new}$
$\displaystyle=\eta^{new}_{i}\left(\sum_{k}u^{old}_{k}\phi^{old}_{k}\right)$
$\displaystyle u^{new}$
$\displaystyle=\sum_{k}u^{old}_{k}\eta^{new}_{i}\left(\phi^{old}_{k}\right)$
where, $\displaystyle\eta^{new}_{i}\phi^{old}_{k}$
$\displaystyle=\sum_{q}w_{q}\phi^{old}_{k}(x_{q})$

Here, $\phi$ represents the basis, $\eta$ represents the dual basis, $x_{q}$
are the quadrature points and $w_{q}$ are the weights on the quadrature
points.

Figure 2: Scaling of the mesh according to $L(t)$.

  * •

PreStep:

    * –

Calculate $V_{t+1}=V_{t}+V_{dt}$ ($V_{0}$ is known)

    * –

Calculate $l_{t+1}=l_{t}+(u_{tip}*dt)$

    * –

Scale the mesh and Project the solution

    * –

Label the cells to refine

  * •

TSSolve: Solve the equations

  * •

TSAdapt:

    * –

Calculate $V_{new}$

    * –

If($V_{new}=V_{t+1}$) $\Rightarrow$ Accept the solution

    * –

Else

      * -

Reject the solution

      * -

Adapt the length using $V_{new}$ and $V_{t+1}$

      * -

Scale the mesh and Project the solution

      * -

GoTo: TSSolve

  * •

PostStep: Refine the labeled cells

Algorithm 1 Self consistency using TSAdapt

Calculating the length, scaling the mesh, projecting the solution and the
length adaptation (when volume lost is more than a specified threshold) can be
merged into a self-consistent loop as shown in Algorithm 1. The neck requires
sufficiently refined mesh to capture the singularity. Hence, we start with a
coarse mesh and refine it as we approach the singularity. The elements are
labeled for refinement based on the radius and velocity gradients. Before the
next time step, the labeled elements are then refined if necessary. This
adaptive mesh refinement is also included in the algorithm [23, 24]. We use
the Portable, Extensible Toolkit for Scientific Computation (PETSc) [25, 26]
to set up and solve the system using time-stepper (TS) object [27]. The self-
consistent algorithm is set up using TSAdapt functionality of TS. We use a
direct solver using LU factorization.

##  3 Results and discussion

In this section, we discuss the verification and validation of the numerical
model presented in the previous section. Then, we explore the pinch-off
dynamics in paraffin wax.

###  3.1 Verification and Validation

Table 1: Physical properties of the materials used for the simulations. Material |  $\rho$ (kg/m3) |  $\nu$ (m2/s) |  $\gamma$ ***Surface tension value at liquid-air interface. (N/m)
---|---|---|---
Water | $998$ | $1.0\times 10^{-6}$ | $0.0728$
$85\%$ Glycerol | $1223$ | $9.2\times 10^{-5}$ | $0.066$
Paraffin wax †††Wax properties are at 100∘ C. | $760$ | $4.2\times 10^{-6}$ | $0.028$

Before proceeding with the computational model, it is vital to perform a
verification test. Verification is a mathematical exercise that can be used to
examine the error evaluation done by the implemented code. One elegant method
for code verification is the Method of Manufactured Solution (MMS) [28]. This
is a very straightforward method, where we simply pick a non-trivial solution
and add the source term into the equation generated by applying the operator
on the solution. This way we know the exact solution of the modified equation
(original equation with the source term). Then the error evaluation for this
modified equation must be zero. The computational model we use is verified
using the MMS. The MMS helps to eliminate coding errors and is also useful to
test the discretization for problems with unknown exact solutions. Figure 3
shows log-log plot of $L_{2}$ norm of the error, $||u_{fe}-u_{mms}||_{L_{2}}$,
where $u_{fe}$ is finite element solution and $u_{mms}$ is the manufactured
solution. The velocity ($u$) and radius ($h$) are discretized using third-
order polynomials, whereas slope ($s$) is discretized using second-order
polynomials. The error reduces by order four for $u$ and $h$, and order three
for $s$, verifying the correct implementation of the numerical model. The
error evaluation is done on a moving mesh with the scaling factor of $1.0001$.
The MMS solution is evaluated every time step on a scaled mesh and compared to
the solution.

Figure 3: Discretization errors in velocity ($u$), radius ($h$) and slope
($s$) for different mesh sizes using manufactured solution. Figure 4:
Evolution of non-dimensional length in the time approaching the pinch-off for
(a) water and (b) $85\%$ glycerol. The inlet radius $h_{0}=0.0016$ m and
inflow rate $=$ 1 mL/min.

For the validation of our computational model, we use the experimental results
by Zhang and Basaran [4]. One crucial parameter to validate is the evolution
of the droplet length with time. Because the calculation of the length is
implicitly involved in the numerical model as explained in Algorithm 1. Figure
4 illustrates the comparison of the numerical simulation result with the
experimental data. The length ($L$) is non-dimensionalized by the inlet radius
($h_{0}$). The time axis shows the time distance from the pinch-off. The
comparison presented is for water and ($85\%$) glycerol solution. The
viscosity, density, and surface tension are given in Table 1. Initially, the
droplet evolves slowly because the surface tension force is stronger than the
effect of gravity. As we add more fluid, the drop becomes heavier and
eventually, the gravitational force surpasses the surface tension. The surface
tension starts to decrease the radius from the middle section, trying to
minimize the surface energy. The length evolves much faster after the necking
begins, suggesting an increase in the advection due to increasing velocity
gradients. The primary droplet is separated when the radius is zero. The
visualization of the droplet evolution after necking is also shown for water
and ($85\%$) glycerol. Each droplet profile is attached to the point in the
plot that corresponds to the time away from the pinch-off. The glycerol
solution shows more elongation approaching the pinch-off time, which explains
the effect of viscosity. The strong viscous effect allows glycerol droplet to
have a long neck. In case of the water droplet, the surface tension force is
much more dominant compared to the viscous forces approaching the singularity,
resulting a shorter neck.

Figure 5: Evolution of non-dimensional radius of the narrowest location in
time approaching the pinch-off for (a) water and (b) $85\%$ glycerol. The
inlet radius $h_{0}=0.0016$ m and inflow rate $=$ 1 mL/min.

Another important parameter for validation is the evolution of a minimum
radius in time away from the pinch-off. Figure 5 illustrates numerical results
for the evolution of this parameter, compared with the experimental profiles
for water and ($85\%$) glycerol solution. The profiles for both water and
glycerol show similar evolution until a point where the dynamics start to
become self-similar. From this point, the glycerol profile shows the influence
of high viscosity by more elongation. The minimum radius decreases slower,
which is evident by the long tail at the end of the profile. For the water
droplet, the radial shrinkage is much faster due to the small viscosity
compared to the glycerol solution. The numerical profile for glycerol agrees
with the experiments very well. For the water, the profile shows a small
amount of delay between the numerical and experimental results. This is due to
the initially added artificial viscosity to increase diffusive behavior, which
is inspired by the computational model by Eggers and Dupont [2]. We increase
the viscosity value initially for the low-viscosity materials since their
highly convective nature introduces surface fluctuations. We reduce this added
diffusion to zero well before the necking begins. Hence, the numerical profile
for the water droplet starts to agree with the experiments right where the
necking begins. This added diffusion shows no impact on the length evolution
or the pinch-off location at all.

###  3.2 Pinch-off dynamics of paraffin wax

Paraffin wax is explored as one of the potential candidate fuels for hybrid
propellant rockets [29]. In the combustion chamber of a hybrid rocket, a solid
paraffin wax form a liquid layer on its burning surface. This liquid layer,
under a high shear forcing, shows hydrodynamic instabilities that lead to the
formation of droplets. These droplets are then entrained in the flow. Here, we
explore the pinch-off dynamics of a pendant paraffin wax droplet using our
computational model.

Figure 6: The curvature profile and the velocity derivative profile when the
pinch-off happens. The inlet radius $h_{0}=0.0016$ m and inflow rate $=$ 1
mL/min.

The curvature is a chief attribute in understanding the pinch-off dynamics in
droplet formation. Figure 6 shows a paraffin wax droplet profile at pinch-off
alongside a curvature profile. Approaching the singularity, the curvature
starts to increase. But a finite time curvature blow-up happens at the pinch-
off location, where the radius is zero. As explained in the introduction
section, the fluid in close vicinity of this singular point is driven by very
high-velocity gradients. The motion is independent of the initial and boundary
conditions. This feature is mathematically described by a self-similar
solution. The velocity derivative profile on the right in Fig. 6 shows high
gradients close to the singularity. Moreover, the velocity gradients change
signs at the singular point location. This suggests that the fluid above and
below the pinch-off point moves in opposite directions very quickly. This
recoil provokes surface instabilities leading to the satellite drop formation.

However, it is also evident that the fluid motion approaching the pinch-off is
highly convective in nature. The finite element model can approximate the
solution but when the truncated terms are getting larger, the solution becomes
unstable. Hence, the numerical scheme is augmented with a stabilization
technique when the fluid viscosity is low. There are many options to consider
for stabilization [30]. We used the Streamline Upwinding (SU) and Streamline
Upwinding Petrov Galerkin (SUPG) and found that this 1D problem can become
stable with just the SU method. The SUPG method is better suited for problems
with cross-convection in 2D or 3D. Also, the SUPG method regularizes the
strong form residual that contains second-order derivatives, which can be
problematic for $C^{0}$ continuous elements. However, the SU scheme just adds
an artificial diffusion into the system and is over-diffusive in nature [31].
Therefore, we also decrease the artificial diffusivity as we refine the mesh
adaptively to avoid adding too much diffusion. This stabilization was only
enabled for low-viscosity fluids, like paraffin wax, water, etc. High-
viscosity fluids like glycerol can be handled without any stabilization.

##  4 Conclusion

A one-dimensional numerical model is reliable for simulating droplet
formation. The present model is validated using experimental pendant drops of
water and glycerol, which were excellently matched with simulation results.
The SU stabilization method was implemented since the low-viscosity fluids,
like water and paraffin wax are highly convective. The velocity and curvature
profile was shown for paraffin wax at the pinch-off time. The velocity was
approaching infinity at the pinch-off location with the opposite signs
suggesting that the pinch-off is followed by the recoil and then satellite
drop formation.

The current computational model can be extended to capture the motion after
pinch-off and satellite drops as well. Understanding the volume of satellite
droplets is useful in accurately predicting the regression rate of the fuel in
hybrid rockets. Also, the present model can be improved for droplet formation
in a shear environment. The droplet formation in a turbulent environment may
not be axisymmetric since the turbulence sets the initial droplet profile
[32]. For non-symmetric interfaces, where no swirling motion assumption fails,
two- or three-dimensional approaches are better, which is also recommended by
Ambraneswaran [3]. However, the droplet shape may not be assumed to be
axisymmetric, the solution in the singularity region is still self-similar and
follows the analogous dynamics.

##  5 Acknowledgement

Funded by the United States Department of Energy’s (DoE) National Nuclear
Security Administration (NNSA) under the Predictive Science Academic Alliance
Program III (PSAAP III) at the University at Buffalo, under contract number
DE-NA0003961.

This work was partially supported by the Department of Energy Office of
Science Award DE-AC02-0000011838.

## References

  * [1] Darsh K. Nathawani and Matthew G. Knepley  “Droplet formation simulation using mixed finite elements”  In _Physics of Fluids_ 34.6, 2022, pp. 064105  DOI: 10.1063/5.0089752
  * [2] Jens Eggers and Todd F Dupont  “Drop formation in a one-dimensional approximation of the Navier–Stokes equation”  In _Journal of fluid mechanics_ 262 Cambridge University Press, 1994, pp. 205–221
  * [3] Bala Ambravaneswaran, Edward D Wilkes and Osman A Basaran  “Drop formation from a capillary tube: Comparison of one-dimensional and two-dimensional analyses and occurrence of satellite drops”  In _Physics of Fluids_ 14.8 American Institute of Physics, 2002, pp. 2606–2621
  * [4] Xiaoguang Zhang and Osman A Basaran  “An experimental study of dynamics of drop formation”  In _Physics of fluids_ 7.6 American Institute of Physics, 1995, pp. 1184–1203
  * [5] Lord Rayleigh  “On the instability of jets”  In _Proceedings of the London mathematical society_ 1.1 Wiley Online Library, 1878, pp. 4–13
  * [6] Jens Eggers  “Nonlinear dynamics and breakup of free-surface flows”  In _Reviews of modern physics_ 69.3 APS, 1997, pp. 865  URL: https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.69.865
  * [7] Jens Eggers  “Theory of drop formation”  In _Physics of Fluids_ 7.5 American Institute of Physics, 1995, pp. 941–953
  * [8] KC Chaudhary and LG Redekopp  “The nonlinear capillary instability of a liquid jet. Part 1. Theory”  In _Journal of Fluid Mechanics_ 96.2 Cambridge University Press, 1980, pp. 257–274
  * [9] KC Chaudhary and T Maxworthy  “The nonlinear capillary instability of a liquid jet. Part 2. Experiments on jet behaviour before droplet formation”  In _Journal of fluid mechanics_ 96.2 Cambridge University Press, 1980, pp. 275–286
  * [10] DB Bogy  “Drop formation in a circular liquid jet”  In _Annual Review of Fluid Mechanics_ 11.1 Annual Reviews 4139 El Camino Way, PO Box 10139, Palo Alto, CA 94303-0139, USA, 1979, pp. 207–228
  * [11] HC Lee  “Drop formation in a liquid jet”  In _IBM Journal of Research and Development_ 18.4 IBM, 1974, pp. 364–369
  * [12] LE Cram  “A numerical model of droplet formation”  In _Computational Techniques and Applications: CTAC-83_ Elsevier, 1984, pp. 182–188
  * [13] D Howell Peregrine, G Shoker and A Symon  “The bifurcation of liquid bridges”  In _Journal of Fluid Mechanics_ 212 Cambridge University Press, 1990, pp. 25–39
  * [14] XD Shi, Michael P Brenner and Sidney R Nagel  “A cascade of structure in a drop falling from a faucet”  In _Science_ 265.5169 American Association for the Advancement of Science, 1994, pp. 219–222
  * [15] SF Kistler and Laurence Edward Scriven  “Coating flow theory by finite element and asymptotic analysis of the Navier-Stokes system”  In _International Journal for Numerical Methods in Fluids_ 4.3 Wiley Online Library, 1984, pp. 207–229
  * [16] John R Richards, Antony N Beris and Abraham M Lenhoff  “Drop formation in liquid–liquid systems before and after jetting”  In _Physics of fluids_ 7.11 American Institute of Physics, 1995, pp. 2617–2630
  * [17] Xiaoguang Zhang  “Dynamics of growth and breakup of viscous pendant drops into air”  In _Journal of colloid and interface science_ 212.1 Elsevier, 1999, pp. 107–122
  * [18] Brian Derby  “Inkjet printing of functional and structural materials: fluid property requirements, feature stability, and resolution”  In _Annual Review of Materials Research_ 40 Annual Reviews, 2010, pp. 395–414
  * [19] Jungho Kim  “Spray cooling heat transfer: The state of the art”  In _International Journal of Heat and Fluid Flow_ 28.4 Elsevier, 2007, pp. 753–767
  * [20] C Berna, A Escrivá, JL Muñoz-Cobo and LE Herranz  “Review of droplet entrainment in annular flow: Characterization of the entrained droplets”  In _Progress in Nuclear Energy_ 79 Elsevier, 2015, pp. 64–86
  * [21] MA Karabeyoglu and Brian J Cantwell  “Combustion of liquefying hybrid propellants: Part 2, stability of liquid films”  In _Journal of Propulsion and Power_ 18.3, 2002, pp. 621–630
  * [22] PE Farrell and JR Maddison  “Conservative interpolation between volume meshes by local Galerkin projection”  In _Computer Methods in Applied Mechanics and Engineering_ 200.1-4 Elsevier, 2011, pp. 89–100
  * [23] Michael Lange, Lawrence Mitchell, Matthew G. Knepley and Gerard J. Gorman  “Efficient mesh management in Firedrake using PETSc-DMPlex”  In _SIAM Journal on Scientific Computing_ 38.5, 2016, pp. S143–S155  DOI: 10.1137/15M1026092
  * [24] Joseph G. Wallwork, Matthew G. Knepley, Nicolas Barral and Matthew D. Piggott  “Parallel Metric-Based Mesh Adaptation in PETSc using ParMmg”  In _SIAM International Meshing Roundtable Workshop 2022_ , 2021, pp. 1–5  URL: http://imr.sandia.gov/papers/imr31/2026_imr31RJ_Wallwork.pdf
  * [25] Satish Balay et al.  “PETSc/TAO Users Manual”, 2021
  * [26] Satish Balay et al.  “PETSc Web page”, https://petsc.org/, 2021  URL: https://petsc.org/
  * [27] Shrirang Abhyankar et al.  “PETSc/TS: A Modern Scalable DAE/ODE Solver Library”, 2014
  * [28] Patrick J Roache  “Code verification by the method of manufactured solutions”  In _J. Fluids Eng._ 124.1, 2002, pp. 4–10
  * [29] MA Karabeyoglu, D Altman and Brian J Cantwell  “Combustion of liquefying hybrid propellants: Part 1, general theory”  In _Journal of Propulsion and Power_ 18.3, 2002, pp. 610–620
  * [30] Jean Donea and Antonio Huerta  “Finite element methods for flow problems”  John Wiley & Sons, 2003
  * [31] Alexander N Brooks and Thomas JR Hughes  “Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations”  In _Computer methods in applied mechanics and engineering_ 32.1-3 Elsevier, 1982, pp. 199–259
  * [32] Daniel J Ruth, Wouter Mostert, Stéphane Perrard and Luc Deike  “Bubble pinch-off in turbulence”  In _Proceedings of the National Academy of Sciences_ 116.51 National Acad Sciences, 2019, pp. 25412–25417


A Polar Subcode Approach to Belief Propagation List Decoding

Marvin Geiselhart, Ahmed Elkelesh, Jannis Clausius and Stephan ten Brink

Institute of Telecommunications, Pfaffenwaldring 47, University of  Stuttgart, 70569 Stuttgart, Germany

This work is supported by the German Federal Ministry of Education and Research (BMBF) within the project Open6GHub (grant no. 16KISK019).

MLmaximum likelihood
BPbelief propagation
BPLbelief propagation list
LDPClow-density parity-check
BERbit error rate
BPSKbinary phase shift keying
AWGNadditive white Gaussian noise
LLRlog-likelihood ratio
MAPmaximum a posteriori
BLERblock error rate
SCLsuccessive cancellation list
SCsuccessive cancellation
BI-DMCBinary Input Discrete Memoryless Channel
CRCcyclic redundancy check
CA-SCLCRC-aided SCL
CA-BPLCRC-aided belief propagation list
BECBinary Erasure Channel
BSCBinary Symmetric Channel
3GPP3rd Generation Partnership Project
5G5th generation mobile communications
eMBBenhanced Mobile Broadband
CNcheck node
VNvariable node
GenAlgGenetic Algorithm
CSIChannel State Information
OSDordered statistic decoding
MWPC-BPminimum-weight parity-check BP
FFGForney-style factor graph
MBBPmultiple-bases belief propagation
URLLCultra-reliable low-latency communications
DMCdiscrete memoryless channel
SGDstochastic gradient descent
PR-CA-BPpermuted relaxed CRC-aided BP
R-CA-BPrelaxed CRC-aided BP
R-CA-BPLrelaxed CRC-aided BPL
PR-CA-BPLpermuted relaxed CRC-aided BPL
P-CA-BPLpermuted CRC-aided BPL
P-CA-BPpermuted CRC-aided BP
LTAlower triangular affine
GAgeneral affine
AEDautomorphism ensemble decoding
SPAsum-product algorithm
PACpolarization adjusted convolutional

Permutation decoding gained recent interest as it can exploit the symmetries of a code in a parallel fashion. Moreover, it has been shown that by viewing permuted polar codes as polar subcodes, the set of usable permutations in permutation decoding can be increased. We extend this idea to pre-transformed polar codes, such as CRC-aided polar codes, which previously could not be decoded using permutations due to their lack of automorphisms. Using BP-based subdecoders, we showcase a performance close to CA-SCL decoding. The proposed algorithm outperforms the previously best performing iterative CA-BPL decoder both in error-rate performance and decoding latency.

§ INTRODUCTION

Polar codes were first introduced by Arıkan as a class of codes that achieve the capacity of binary memoryless channels [1] using the SC decoding algorithm. In the short block length regime, the performance of stand-alone polar codes (under optimal ML decoding) is however not favorable. It could however be shown that a high-rate outer code, such as a CRC, significantly improves the error-correcting performance of polar codes under SCL decoding [2]. As a consequence, CRC-concatenated polar codes have been selected for the control channel in the 5G standard [3].
Despite its good error-correcting performance, CA-SCL decoding is inherently sequential and hard-output. Consequently, the decoding latency is comparably high and the algorithm is not well-suited for iterative receivers (e.g., iterative detection and decoding loops).

An alternative to SC decoding is the iterative BP algorithm [4, 5], which is parallel by nature and, thus, better suited for high-throughput/low-latency applications. Unfortunately, its error-rate performance is generally worse than that of CA-SCL decoding.

Over the last couple of years, numerous improvements to BP decoding of polar codes have been proposed. Notable is the decoding over multiple permuted factor graphs [6], resulting in variations of the BP algorithm (BPL [7]) and SC decoding [8]. A version of BP decoding with an outer LDPC code has been proposed [9], where the Tanner graph of the outer code has been attached to the left side of the polar factor graph. We note that, unlike in SCL decoding, it turned out to be very difficult to effectively use the error-correcting capability of an outer CRC code in iterative decoding, as its Tanner graph contains many short cycles. Therefore, all known iterative decoding schemes of CRC-aided polar codes require additional steps to reach reasonable performance.
In [10], the edges are weighted to mitigate the effect of short cycles in the CRC Tanner graph. In [11], it has been proposed to reduce the number of edges in the factor graph and use an ensemble of stage-permuted polar BP decoders.
Lastly, the method of “polar relaxation” has been proposed to reduce the number of cycles in the combined factor graph [12].

At the same time, there have been advances in automorphism-based decoding, so-called AED. In AED, multiple permuted versions of the received sequence are decoded in parallel and the most promising candidate is chosen as the final codeword estimate. The used permutations come from the automorphism group of the code. AED has been successfully applied to codes with rich symmetries, such as cyclic codes [13, 14], quasi-cyclic codes [15], RM codes [16] and stand-alone polar codes (i.e., polar codes without the CRC-aid) [17].

Recently, in [18], it has been shown that for polar codes, additional permutations (outside the code's automorphism group) can be used in ensemble decoding. The key idea is to observe that such permutations result in a polar subcode [19].

In this work, we extend the polar subcode approach to iterative decoding. The main contributions of this paper are:

* We generalize the work of [18] to start with arbitrary polar subcodes, i.e., also CRC-aided polar codes.

* We propose to apply the polar subcode permutation scheme to BPL decoding of polar-CRC codes and report both a better error-rate performance and a decreased latency compared to stage-permutation based CA-BPL decoding [11].

* We show that a combination of this ensemble decoding scheme with polar relaxation further enhances the performance without any additional overhead.

§ PRELIMINARIES

Regular factor graph

Relaxed factor graph

Factor graph representations of a (8,5) polar code.

§.§ Polar Codes

Polar codes are a class of binary block codes with blocklength $N=2^n$ constructed from the $ n $-fold application of the polar transformation [1], resulting in the polar transformation matrix $
\mathbf{G}_N = \left[\begin{smallmatrix}
1 & 0 \\
1 & 1
\end{smallmatrix}\right]^{\otimes n}$.
Depending on the desired rate $R$, $k$ of the $N$ inputs (so-called synthetic channels) to the polar transformation carry the information, while the remaining $ N-k $ inputs are frozen bits and are set to zero. The set of information bit indices is denoted by $\mathcal{A}$, while the set of frozen indices is $\mathcal{A}^c$. The selection of $\mathcal{A}$ is subject to the code design and out of scope of this paper.

A graphical representation of polar codes is its factor graph, as depicted in Fig. <ref> for the (8,5) polar code with information set $\mathcal{A}=\{3,4,5,6,7\}$. Frozen bits are indicated by red dots.

§.§ Relaxation

Polar code relaxation refers to the removal of vertical connections from the factor graph that do not contribute to polarization [20]. In particular, a connection is removed if it connects either two frozen bits or two information bits. Fig. <ref> shows the relaxed factor graph of a (8,5) polar code. The removed connections are indicated by the dashed boxes.

In the case of stand-alone polar codes, relaxation does not change the set of codewords (but potentially the encoding). This means, relaxation can also be applied on only the receiver side, for example to reduce the decoding complexity or latency. Moreover, lower bit error-rates are reported in [20]. Mathematically, we denote the relaxed polar transformation matrix as $\mathbf{G}_N^{(R)}$, which is dependent on the set of frozen bits.

§.§ Polar Subcodes

While polar codes are proven to achieve the channel capacity as the block length goes to infinity ($N \to \infty$), they suffer from poor distance properties in the finite length regime.
A common method of enhancing the performance of polar codes is to concatenate with an outer code, which removes the low-weight codewords from the polar code [2]. The resulting code is a subcode of the polar code [19].
A polar subcode is defined by the dynamic freezing constraint matrix $\mathbf{V} \in \mathbb{F}_2^{(N-k)\times N}$. Rows in $\mathbf{V}$ with a single `1' correspond to (hard) frozen bits, and rows with multiple `1s' are so-called dynamic frozen bits. The parity-check matrix of the concatenated code is given by

\begin{equation}
\mathbf{H} = \mathbf{V}\cdot \mathbf{G}_N^T.
\end{equation}

Conversely, there exists a pre-transform matrix $\mathbf{W} \in \mathbb{F}_2^{k\times N}$, such that the whole polar subcode is generated by

\begin{equation}
\mathbf{G} = \mathbf{W}\cdot \mathbf{G}_N.
\end{equation}

Note that since $\mathbf{G}_N$ is invertible, any binary linear block code of length $N$ can be viewed as a polar subcode [19]. In the following we denote the uncoded message by $\mathbf{u}$, the pre-transformed vector by $\mathbf{v}=\mathbf{u}\mathbf{W}$ and the codeword by $\mathbf{c}=\mathbf{v}\mathbf{G}_N$.

§.§ Affine Permutations

In permutation-based polar decoding, codeword bits $c_i$ (and their channel LLR $L_{\mathrm{ch},i}$) are re-arranged according to affine permutations. An affine permutation is a permutation on $ \{ 0, \dots, N-1\} $, mapping the (least significant bit first) binary representations $ \mathbf{z} $ according to

\begin{equation}
\mathbf{z}' = \mathbf{A}\mathbf{z} + \mathbf{b} \mod 2,
\end{equation}

with $\mathbf{A}\in \mathbb{F}_2^{n\times n}$ and $\mathbf{b}\in \mathbb{F}_2^{n}$. There are three important groups of affine permutations in the context of polar codes:

* If $\mathbf{A}$ is only required to be invertible, we have the GA group of order $n$, denoted by $\operatorname{GA}(n)$. It is well known that $\operatorname{GA}(n)$ is the automorphism group of RM codes, i.e., the code remains the same when symbols are permuted according to any of these permutations [21].

* If $\mathbf{A}$ is invertible and lower-triangular, we have the LTA group of order $n$, denoted by $\operatorname{LTA}(n)$. $\operatorname{LTA}(n)$ is a subset of the automorphism group of polar codes [22].

* If $\mathbf{A}$ is a permutation matrix and $\mathbf{b}=\mathbf{0}$, the resulting permutation group corresponds to stage-shuffling of the polar factor graph [23]. We denote this group by $\Pi(n)$.

§.§ Belief Propagation Decoding of Polar Codes

Iterative BP decoding of polar codes is conventionally performed over the polar factor graph. In each iteration, $\mathbf{R}$-messages propagate from the left side to the right side and $\mathbf{L}$-messages vice versa through the stages of the graph. In each stage, a vertical connection in the (decoding) factor graph is called a processing element and its SISO MAP rule is applied in update. The $\mathbf{L}$-messages are initialized with the channel LLR and the frozen positions of the $\mathbf{R}$-messages with saturated LLR. Bit estimates for both $\hat{\mathbf{v}}$ and $\hat{\mathbf{c}}$ an be calculated by summing the corresponding $\mathbf{L}$- and $\mathbf{R}$-messages. A stopping condition based on
$\mathbf{\hat{v}}\mathbf{G}_N=\mathbf{\hat{c}}$ is employed to terminate the algorithm before reaching $N_{\text{it,max}}$ iterations. Due to the SISO nature of the algorithm, it can be connected to the Tanner graph of an outer code “turbo-like decoder” (see Fig. <ref>). For details on the BP algorithm, we refer the reader to [24], [11].

§.§ Belief Propagation List Decoding

To enhance the performance of BP decoding, an ensemble decoder, coined BPL, has been proposed in [7]. The idea is to perform BP on an ensemble of $L$ different stage-permuted factor graphs in parallel. Each sub-decoder contributes a codeword candidate to a list, from which the most likely valid candidate is picked as the decoder output. An extension towards CRC-aided polar codes is proposed in [11].

§ DECODING GRAPHS OF POLAR SUBCODES

CA-BP factor graph

Permuted (P)-CA-BP factor graph

Permuted and relaxed (PR)-CA-BP factor graph

Different factor graph representations of a (8,5) polar code concatenated with an outer CRC-3 code with $g(x)=x^3+x+1$. Note that the outer code factor graph is different in all three cases.

In this section we derive methods for obtaining different factor graphs of polar subcodes for ensemble decoding based on affine permutations and relaxation.

§.§ Post-Transformations

Given an invertible binary $N\times N$  matrix $\mathbf{T}$, we can write the generator matrix of any polar subcode as

\begin{equation}
\mathbf{W} \mathbf{G}_N \mathbf{T} \mathbf{T}^{-1} = \mathbf{W} \mathbf{G}_N \mathbf{T} \mathbf{G}_N^{-1} \mathbf{G}_N
\mathbf{T}^{-1} = \mathbf{W}_T \mathbf{G}_N \mathbf{T}^{-1}.
\end{equation}

In other words, the encoding is equal to that of a different polar subcode followed by the post-transformation $\mathbf{T}^{-1}$. The pre-transform matrix and dynamic freezing constraint matrix of the new polar subcode are

\begin{align}
\mathbf{W}_T &= \mathbf{W} \mathbf{G}_N \mathbf{T}\mathbf{G}_N \quad \text{and}\\
\mathbf{V}_T &= \mathbf{V} \left(\mathbf{G}_N \mathbf{T}^{-1}\mathbf{G}_N\right)^\mathrm{T},
\end{align}

respectively, where we used $\mathbf{G}_N^{-1}=\mathbf{G}_N$. A particular interesting case is when $\mathbf{T}$ is a permutation matrix, as this can be easily implemented at the first decoding stage as a permutation of the received sequence.

Example 1:
Fig. <ref> shows the factor graph of a (8,5) polar code concatenated with a 3-bit CRC (i.e., CRC-3) described by the generator polynomial $g(x)=x^3+x+1$. It can be interpreted as a polar subcode with dynamic freezing constraint matrix

\begin{equation}
\mathbf{V} = \left[\begin{array}{cccccccc}
\end{array}\right].
\end{equation}

Observe that the first three rows correspond to the frozen bits, while the remaining three rows contain the parity-check matrix of the CRC code.
Now, consider the LTA permutation $\pi=(1,6,3,4,5,2,7,0)$ corresponding to the pair

\begin{equation}
\mathbf{A} = \left[\begin{array}{ccc}
\end{array}\right], \; \mathbf{b} = \left[\begin{array}{c}
\end{array}\right].
\end{equation}

Thus, its post-transformation matrix is

\begin{equation}
\mathbf{T}^{-1} = \left[\begin{array}{cccccccc}
\end{array}\right].
\end{equation}

We then find the transformed dynamic freezing constraint matrix

\begin{equation}
\mathbf{V}_T = \mathbf{V} \left(\mathbf{G}_N \mathbf{T}^{-1}\mathbf{G}_N\right)^\mathrm{T} = \left[\begin{array}{cccccccc}
\end{array}\right].\label{eq:vt}
\end{equation}

Note that the first three columns remain frozen, which can be easily seen after a few row operations. Moreover, additional row-operations can be applied to the dynamic frozen bits, reducing the number of edges in the factor graph, as depicted in Fig. <ref>. The ensemble variant of this decoder is denoted by P-CA-BPL-$L$ decoding.

§.§ Polar Subcode Relaxation

If the same outer code is used together with a relaxed polar code, the resulting concatenated code is generally different than without relaxation. Hence, to comply with a non-relaxed encoder, the outer code must be altered to account for relaxation. The idea is to include the relaxed (i.e., removed) partial polar transformations in the outer code. Those removed polar transformations can be collected in the matrix $\mathbf{R}$, fulfilling

\begin{equation}
\mathbf{R} \cdot \mathbf{G}_N^{(R)} = \mathbf{G}_N.
\end{equation}

For a given dynamic freezing constraint matrix $\mathbf{V}$, the equivalent matrix for the relaxed polar code equates to

\begin{equation}
\mathbf{V}_R = \mathbf{V} \cdot \mathbf{R}^{\mathrm{T}}. \label{eq:relax}
\end{equation}

Similarly, $\mathbf{R}$ can also be applied to a permuted polar code, i.e., its  matrix $ \mathbf{V}_T$.
It should be noted that this variation of relaxation in general (i.e., for outer codes with dense parity-check matrix) retain the decoding complexity reduction as experienced in stand-alone polar codes, as the outer code usually remains high-density when multiplying by $\mathbf{R}$.

Example 1 (cont'd):
As the inner polar code is the same as in Fig. <ref>, relaxation removes the same processing elements. We find the relaxation matrix

\begin{equation}
\mathbf{R}= \left[\begin{array}{cccccccc}
\end{array}\right].
\end{equation}

Using Eq. (<ref>) and a few row operations, we find the permuted relaxed dynamic freezing constraint matrix for (<ref>) to be

\begin{equation}
\mathbf{V}_{T,R} = \left[\begin{array}{cccccccc}
\end{array}\right].\label{eq:vtr}
\end{equation}

Note that bit 5 turned from a dynamic frozen bit into a hard frozen bit. Fig. <ref> shows the permuted relaxed factor graph.

A decoder using this factor graph is called R-CA-BP, and its ensemble variant R-CA-BPL-$L$. If combined with an affine permutation, we refer to it as PR-CA-BPL-$L$ decoding.

§ RESULTS

For all shown results, we selected the information/frozen-bit positions according to the 5G standard [3] and the outer code is a CRC-8 described by the generator polynomial $g(x) = x^8 + x^5 + x^4 + x^3 + 1$. Similar to [11], we apply a heuristic density reduction to the dynamic freezing constraint matrices to be better suited for the SPA. For all shown ensemble decoders, we do not optimize the set of permutations that are used, i.e., we select random permutations from the respective groups. However, optimization methods (e.g., genetic algorithm-based as proposed in [11]) are applicable and expected to provide additional gains. In all simulations, we assume that the codeword is BPSK modulated and transmitted over an AWGN channel.

§.§ Permutation Subgroups

Error-rate performance of a single PR-CA-BP decoder using 100 iterations for the (128,64)-Polar-CRC8 code, random permutations.

While in [18] it is shown that permutations from all subsets of the GA group does not degrade the distance of the transformed code, the performance of the iterative decoding scheme is still highly dependent on the type of permutation. This is due to the changing number and position of the hard frozen bits (i.e., the effective polar code design) and the transformed dynamic frozen bits. To illustrate this, we plot the average performance of a single PR-CA-BP decoder using different permutations in Fig. <ref>. For every simulated frame, a new random permutation is selected from each permutation subgroup. As we can see, the LTA permuted decoder has the best performance, similar to the default (i.e., no) permutation. Stage shuffle permutations $\Pi$ deteriorate the performance on average by 0.2 dB at a BLER of $10^{-2}$, and GA permutations lose on average another 0.3 dB. This effect can be explained by the fact that LTA permutations, which are naturally automorphisms of the inner polar code, do not change the set of frozen/information bits. Hence, we expect all LTA permuted decoders to perform similarly to the unpermuted decoder, while permutations which are not automorphisms may deteriorate the effective polar code design. Based on this result, it is reasonable to limit the permutations used in the ensemble/list decoder to LTA permutations.[The limitation to LTA is for the sake of simplicity; the same argument holds for the full automorphism group of polar codes [17, 25].]

§.§ Ensemble Decoding Error-Rate Performance

Error-rate performance of BPL and SCL decoding for the (128,64)-Polar-CRC8 code. All iterative decoders use $N_\mathrm{it,max}=200$ iterations.

Error-rate performance of BPL and SCL decoding for the (256,128)-Polar-CRC8 code. All iterative decoders use $N_\mathrm{it,max}=200$ iterations.

Fig. <ref> compares the BLER performance of different iterative decoding algorithms with state-of-the-art SCL decoding ($L=8$) for the (128,64) CRC-aided polar code. All ensemble decoders use $L=8$ constituent decoders. At a BLER of $10^{-3}$, the proposed PR-CA-BPL scheme outperforms a single R-CA-BP decoder by more than 0.5 dB. Compared to CA-BPL (with stage shuffle permutations), the richer LTA permutations result in a gain of 0.4 dB. Roughly half of this gain is attributed to the use of relaxation. Interestingly, relaxation has a greater effect on LTA permutations than on $\Pi$ permutations, as their un-relaxed counter parts have equal performance. Lastly, PR-CA-BPL closes the gap to SCL decoding to 0.1 dB at a BLER of $10^{-4}$.
In Fig. <ref>, we show the same comparison for the longer, (256,128) CRC-aided polar code.
We see similar results as for $N=128$ with a gain of ensemble decoding over R-CA-BP decoding of roughly 0.4 dB.

§.§ Latency

Decoding latency comparison of iterative and CA-SCL decoding. All iterative decoders are limited to $N_\mathrm{it,max}=200$ iterations and use early stopping. As the blocklength $N$ increases, the latency gains due to the usage of iterative decoders become more pronounced.

Next, we compare the decoding latency of the different iterative decoding schemes and CA-SCL decoding. We assume a fully parallel implementation of each decoding algorithm, i.e., all operations that do not depend on each other can be executed at the same time and a single stage update takes one clock cycle. For iterative decoders, we assume a stopping criterion based on $\mathbf{\hat{v}}\mathbf{G}_N=\mathbf{\hat{c}}$ and the CRC check is implemented.
The average latency of a CA-BP decoder is the average number of iterations $N_\mathrm{it,avg}$ times the latency of a single iteration, consisting of $n$ left-to-right updates, $n$ right-to-left updates and two cycles for the outer SPA update. Hence, we have the average latency

\begin{equation}
\bar{\tau}_{\mathrm{CA-BP}} = \operatorname{E}\left[(2n+2)\cdot N_\mathrm{it}\right] = (2n+2)\cdot N_\mathrm{it, avg}.
\end{equation}

For the ensemble variants (i.e., x-CA-BPL), we have to wait until all $L$ decoders have converged. Hence, we need to consider the average maximum number of iterations. Additionally, we assume an extra cycle is required to select the final candidate. The latency is therefore

\begin{equation}
\bar{\tau}_{\mathrm{CA-BPL}} = \operatorname{E}\left[(2n+2)\cdot \max_{i=1,\dots,L}\{N_{\mathrm{it},i}\}+1\right].\label{eq:lat_bpl}
\end{equation}

Lastly, for CA-SCL decoding, there are $2N-1$ stage updates required, as well as $k$ additional cycles for metric sorting, and a final cycle for candidate selection [26]. Therefore, we have

\begin{equation}
\bar{\tau}_{\mathrm{CA-SCL}} = 2N+k = (2+R)N.\label{eq:lat_scl}
\end{equation}

Fig. <ref> shows the latency of the compared decoders with respect to the SNR. As expected, the single R-CA-BP decoder has the lowest latency. Moreover, we can see that the proposed LTA permuted decoders generally have a significantly lower latency than stage-shuffle permutations used in CA-BPL. In contrast, relaxation only has a very minor impact on the convergence speed. The proposed PR-CA-BPL algorithm outperforms CA-SCL in terms of latency in the medium-to-high SNR regime by a big margin. Additionally, as can be seen from Eq. (<ref>), Eq. (<ref>) and Fig. <ref>, the latency benefits of BPL-based decoders over SCL is larger for longer blocklengths $N$.

§ CONCLUSION

In this work, we extend the CA-BPL decoding algorithm to use LTA permutations rather than stage-shuffles. Combined with factor graph relaxation, the proposed PR-CA-BPL algorithm outperforms CA-BPL both in error-rate performance and complexity. Overall, an error-rate performance similar to CA-SCL can be reached, at a significantly lower decoding latency in the operating regime. Possible extensions of this work include the optimization of the used permutations (either for latency or error-rate performance) and the application of the proposed method to other pre-transformed polar codes, such as Polar-LDPC [9] codes or PAC codes, which can be described in the same polar subcode framework.

[1]
E. Arıkan, “Channel Polarization: A Method for Constructing
Capacity-Achieving Codes for Symmetric Binary-Input Memoryless Channels,”
IEEE Trans. Inf. Theory, vol. 55, no. 7, pp. 3051–3073, Jul. 2009.

[2]
I. Tal and A. Vardy, “List Decoding of Polar Codes,” IEEE Trans. Inf.
Theory, vol. 61, no. 5, pp. 2213–2226, May 2015.

[3]
“Technical Specification Group Radio Access Network,” 3GPP, 2018, TS
38.212 V.15.1.1. [Online]. Available:

[4]
E. Arıkan, “A Performance Comparison of Polar Codes and Reed-Muller
Codes,” IEEE Commun. Lett., vol. 12, no. 6, pp. 447–449, Jun. 2008.

[5]
——, “Polar Codes: A Pipelined Implementation,” Proc. 4th ISBC,
pp. 11–14, 2010.

[6]
N. Hussami, S. B. Korada, and R. Urbanke, “Performance of Polar Codes for
Channel and Source Coding,” in IEEE Inter. Symp. Inf. Theory (ISIT),
June 2009, pp. 1488–1492.

[7]
A. Elkelesh, M. Ebada, S. Cammerer, and S. ten Brink, “Belief Propagation
List Decoding of Polar Codes,” IEEE Commun. Lett., vol. 22, no. 8,
pp. 1536–1539, Aug. 2018.

[8]
S. A. Hashemi, N. Doan, M. Mondelli, and W. J. Gross, “Decoding
Reed-Muller and Polar Codes by Successive Factor Graph Permutations,” in
IEEE 10th Inter. Symp. on Turbo Codes Iterative Inf. Process. (ISTC),
Dec. 2018.

[9]
J. Guo, M. Qin, A. G. i Fàbregas, and P. H. Siegel, “Enhanced Belief
Propagation Decoding of Polar Codes through Concatenation,” in IEEE
Inter. Symp. Inf. Theory (ISIT), June 2014, pp. 2987–2991.

[10]
N. Doan, S. A. Hashemi, E. N. Mambou, T. Tonnellier, and W. J. Gross, “Neural
Belief Propagation Decoding of CRC-Polar Concatenated Codes,” in IEEE
Inter. Conf. on Commun. (ICC), May 2019.

[11]
M. Geiselhart, A. Elkelesh, M. Ebada, S. Cammerer, and S. ten Brink,
“CRC-Aided Belief Propagation List Decoding of Polar Codes,” in
IEEE Inter. Symp. Inf. Theory (ISIT), Jun. 2020.

[12]
A. S. Tan, O. Sahin, and S. Hong, “CRC Aided Short-Cycle Free BP Decoding for
Polar Codes,” in 2021 International Symposium on Networks, Computers
and Communications (ISNCC), 2021, pp. 1–6.

[13]
T. Hehn, J. B. Huber, O. Milenkovic, and S. Laendner, “Multiple-bases
belief-propagation decoding of high-density cyclic codes,” IEEE
Transactions on Communications, vol. 58, no. 1, pp. 1–8, 2010.

[14]
C. Chen, B. Bai, X. Yang, L. Li, and Y. Yang, “Enhancing Iterative Decoding
of Cyclic LDPC Codes Using Their Automorphism Groups,” IEEE Trans.
Commun., vol. 61, no. 6, pp. 2128–2137, 2013.

[15]
M. Geiselhart, M. Ebada, A. Elkelesh, J. Clausius, and S. ten Brink,
“Automorphism Ensemble Decoding of Quasi-Cyclic LDPC Codes by Breaking
Graph Symmetries,” IEEE Communications Letters, 2022.

[16]
M. Geiselhart, A. Elkelesh, M. Ebada, S. Cammerer, and S. ten Brink,
“Automorphism Ensemble Decoding of Reed–Muller Codes,” IEEE Trans.
Commun., 2021.

[17]
——, “On the automorphism group of polar codes,” in 2021 IEEE
International Symposium on Information Theory (ISIT), 2021, pp. 1230–1235.

[18]
M. Kamenev, “Improved permutation-based successive cancellation decoding of
polar codes,” in 2021 IEEE 22nd International Conference of Young
Professionals in Electron Devices and Materials (EDM), 2021, pp. 110–115.

[19]
P. Trifonov and V. Miloslavskaya, “Polar subcodes,” IEEE Journal on
Selected Areas in Communications, vol. 34, no. 2, pp. 254–266, 2016.

[20]
M. El-Khamy, H. Mahdavifar, G. Feygin, J. Lee, and I. Kang, “Relaxed polar
codes,” IEEE Transactions on Information Theory, vol. 63, no. 4, pp.
1986–2000, 2017.

[21]
F. J. MacWilliams and N. J. A. Sloane, The Theory of Error-Correcting
Codes, ser. North-Holland Mathematical Library.1em plus 0.5em minus
0.4emNorth-Holland Pub. Co., 1977, no. 16.

[22]
M. Bardet, V. Dragoi, A. Otmani, and J. Tillich, “Algebraic
Properties of Polar Codes From a New Polynomial Formalism,” in IEEE
Inter. Symp. Inf. Theory (ISIT), 2016, pp. 230–234.

[23]
N. Doan, S. A. Hashemi, M. Mondelli, and W. J. Gross, “On the
Decoding of Polar Codes on Permuted Factor Graphs,” in IEEE Global
Commun. Conf. (GLOBECOM), Dec. 2018.

[24]
A. Elkelesh, S. Cammerer, M. Ebada, and S. ten Brink, “Mitigating Clipping
Effects on Error Floors under Belief Propagation Decoding of Polar Codes,”
in Inter. Symp. Wireless Commun. Syst., Aug. 2017.

[25]
Y. Li, H. Zhang, R. Li, J. Wang, W. Tong, G. Yan, and Z. Ma, “The complete
affine automorphism group of polar codes,” in 2021 IEEE Global
Communications Conference (GLOBECOM), 2021, pp. 01–06.

[26]
A. K. Balatsoukas Stimming, “Hardware implementation aspects of polar
decoders and ultra high-speed LDPC decoders,” Ph.D. dissertation, EPFL,