11institutetext: Institute of Astronomy, Faculty of Physics, Astronomy and
Informatics, Nicolaus Copernicus University,
ul. Grudziadzka 5, 87-100 Torun, Poland, 11email<EMAIL_ADDRESS>22institutetext: Jodrell Bank Centre for Astrophysics, School of Physics &
Astronomy, University of Manchester,
Oxford Road, Manchester M13 9PL, UK 33institutetext: School of Mathematical
and Physical Sciences, Macquarie University, Sydney, NSW 2109, Australia
44institutetext: Department of Geodesy, Faculty of Geoengineering, University
of Warmia and Mazury,
ul. Oczapowskiego 2, 10-719 Olsztyn, Poland

# Low velocity streams inside the planetary nebula H 2-18

A 3D photoionization and kinematical reconstruction††thanks: Based on
observations collected at the European Southern Observatory under ESO
programme 099.D-0386(A).

K. Gesicki 11 A. Zijlstra 2233 M. Hajduk 44 A. Iwanowska 11 K. Grzesiak 11 K.
Lisiecki 11 J. Lipinski 11

(Received ; accepted )

###### Abstract

Aims. Numerous planetary nebulae show complicated inner structures not
obviously explained. For one such object we undertake a detailed 3D
photoionization and kinematical model analysis for a better understanding of
the underlying shaping processes.

Methods. We obtained 2D ARGUS/IFU spectroscopy covering the whole nebula in
selected, representative emission lines. A 3D photoionization modelling was
used to compute images and line profiles. Comparison of the observations with
the models was used to fine-tune the model details. This predicts the
approximate nebular 3D structure and kinematics.

Results. We found that within a cylindrical outer nebula there is a hidden,
very dense, bar-like or cylindrical inner structure. Both features are co-
axial and are inclined to the sky by 40 deg. A wide asymmetric one-sided plume
attached to one end of the bar is proposed to be a flat structure. All nebular
components share the same kinematics, with an isotropic velocity field which
monotonically increases with distance from the star before reaching a plateau.
The relatively low velocities indicate that the observed shapes do not require
particularly energetic processes and there is no indication for the current
presence of a jet. The 3D model reproduces the observed line ratios and the
detailed structure of the object significantly better than previous models.

###### Key Words.:

planetary nebulae: general – planetary nebulae: individual: H 2-18 (PN G
006.3+04.4)

##  1 Introduction

The theory of mass loss at the late asymptotic giant branch evolutionary phase
is still in development and this is one of the most important missing factors
in the stellar evolution theory. The ejected matter often takes on interesting
axi- or a-symmetric shapes and the understanding of physical processes causing
these is far from clear. Theoretical approaches need some suggestions and
hints derived from observations. Such hints can be found from studies of
planetary nebulae (PN).

Figure 1: The HST image of the PN H 2-18 obtained in an H$\alpha$ filter. The
plot is in linear intensity scale. The angular dimensions are converted to
$10^{17}$ cm assuming that the PN is at the Galactic Bulge distance of 8 kpc
and centred at the star position. The image is rotated to easier compare with
computed structures.

A classical concept of a PN is a shell of gas ejected by a star at the end of
its life. Initially the shell is opaque and hidden from view. Only after the
shell expands sufficiently for the hot central core to shine through and
ionize the nearby material does the object enter the PN phase. This phase ends
when the ejecta expand to a large volume and very low density, or when the
remaining stellar core enters the White Dwarf cooling track and fades
significantly. Both timescales are comparable.

We study details of PNe by applying the so-called kinematical reconstruction
(sometimes called a ‘reverse engineering’). This derives the nebular structure
and kinematics from line profiles based on assumed velocity fields which are
verified through photoionization and emission line profiles computation
(Gesicki & Zijlstra, 2000). For this method to work and to reduce the inherent
ambiguities, both high quality images and high resolution spectroscopy are
needed in order to constrain the density and velocity spatial distributions.

Figure 2: The ARGUS/IFU monochromatic images with HST contours defined in Fig.
4. The high resolution HST image has been positioned to obtain the best
agreement with IFU images of obviously low resolution. Shown are examples of
H$\alpha$, [N ii] 658.3 nm, and He ii 468.6 nm lines.

Here we present new 2D spectroscopic observations and new 3D photoionization
modelling of a previously analyzed PN, allowing us to compare the new results
with the older ones, to verify the new approach and to indicate all its
advantages.

The planetary nebula H 2-18 (see Fig.1) was previously analyzed using the
kinematical reconstruction approach. Gesicki et al. (2014) applied the Torun
codes and derived a spherically symmetric model of this PN. Later, Gesicki et
al. (2016) used the same observational data but applied the pseudo-3D code py-
Cloudy. Both papers were based on the HST imaging and ESO/VLT spectroscopic
long-slit observations.

In 2017 we obtained ARGUS IFU 2D data of high spatial and spectral resolution
of this object. The observations covered selected spectral lines of high and
low excitation which probe different nebular regions and allowed for a more
precise kinematic and photoionization analysis. The collected data are
presented in the Sect. 2.

For modelling of this nebula we applied the publicly available code MOCASSIN
which computes photoionization of a 3D gaseous structure. To initiate the 3D
modelling and to compare the model outputs with the 2D spectra we wrote a
number of Python scripts. The modelling procedure and the obtained density
structure and velocity field are described in the Sect. 3.

The nebula H 2-18 appears to be complex and there are few published studies of
similar objects to which it can be compared. The generalization of our unique
results is also not easy. Some aspects of this are finally discussed in Sect.
4.

##  2 Observations

The ESO FLAMES/ARGUS instrument provided the high-resolution spectroscopy of
spatially resolved Galactic Bulge PN H 2-18 (PN G 006.3+04.4; 17:43:28.7,
$-$21:09:51.30 (J2000)). The nebula was observed in 2017, on Sep. 18
(wavelengths around 465 and 504 nm) and on Sep. 21 (627 and 651 nm). The wide
field of view is a mosaic of $22\times 14$ microlenses each of $0.52\times
0.52$ arcsec. The instrument setup was selected to cover at high spectral
resolution ($R\gtrsim 30\,000$) the regions around four lines: [O iii] 500.7
nm, H$\alpha$ (including [N ii] 658.3 nm), He ii 468.6 nm and [O i] 630.2 nm.
For this object at the time of observations, the line [O i] 630.2 nm happened
to be strongly affected by the overlying terrestrial atmospheric emission
which unfavourably overlapped in radial velocity, and so it was excluded from
the analysis.

The standard ESO reduction pipeline (EsoReflex) has been applied to the raw
data. Further data processing and plotting was performed using Python with
modules astropy, scipy and matplotlib. These concern the removal of bad pixels
by median filtering the data cube in the XY axes (sky plane) and smoothing the
spectral noise by the Savitzky-Golay filter along the wavelength axis. From
the original IFU field of $22\times 14$ pixels we extracted for the analysis a
square of $10\times 10$ pixels which covers the whole nebula and is shown in
the figures. The examples of the monochromatic images (i.e. integrated over
the full width of a single emission line) are shown in Fig. 2. They show
background emission increasing towards the right edges which was not flat-
fielded completely by the pipeline. We did not correct for this because it did
not affect the morphological and kinematical modelling. Absolute flux
calibration was not attempted, instead the line flux values needed for the
model fitting were taken from literature. In all figures we present the
observational data with the sky coordinates converted to $10^{17}$ cm
(assuming the PN distance of 8 kpc) and measured with respect to the position
of the central star to ease the comparison with models where the real physical
size is used. We do not present the [O iii] 500.7 nm monochromatic image
because it is very similar to H$\alpha$.

Previous HST observations in the H$\alpha$ filter show an overall axially
symmetric, elongated structure which in Fig. 1 is rotated so that it is
positioned vertically. The most intense, inner broad area is asymmetric with
strongest emission on the upper-left of the inner region. Outside of the main
cylinder near its middle, a low-intensity ears-like extensions can be seen
(oriented horizontally in the image) which could be a trace of a faint
equatorial ring.

Our new data confirm these features (Fig. 2 left panel) and show more than
that. The IFU data revealed an intriguing pair of small and bright blobs in
the [N ii] 658.3 nm image (Fig. 2 middle panel); their dominant emission is
axially symmetric with a weaker asymmetric component. This feature is
completely absent in He ii 468.6 nm (Fig. 2 right panel) and is barely visible
in H$\alpha$.

We see now that neither a spherical nor an elliptical model is sufficient to
describe the PN H 2-18 and we need to consider the inner component. Although
the elongation in the [N ii] image appears jet-like, because of the high
density we estimate and the low velocity, we shall rather call it bar-like or
cylindrical, leaving the case open.

##  3 3D model analysis

The photoionization modelling was performed with the publicly available
MOCASSIN code (Ercolano et al., 2003). Supplementary codes were written in
Python, these concern a construction of a 3D density distribution cube,
calculation of the emission line shapes derived from total emissivity with
assumed velocity field and emerging from a defined pixel area on the sky,
graphical presentation of the observed and modelled emissions, etc. The models
presented here were computed on a spatial grid of $121\times 121\times 121$
points which ensured sufficient spatial accuracy with reasonable time for
simulations and image processing.

We aimed at a model which, while being as simple as possible, fits the images
and simultaneously reproduces the line ratios and velocities.

###  3.1 Density structure

Figure 3: The density cross-section of the assumed model. The symmetry axis is
positioned vertically, all structures (except for the plume) are cylindrically
symmetrical. This graph shows the true spatial dimensions of the substructures
for a distance of 8 kpc. The plume is in the plane of the image and extends to
$\pm 2.4\times 10^{16}$ cm above and below it.

Figure 4: The observed image in the H$\alpha$ line (left panel) compared with
the one obtained from the assumed model (right panel). The axes are the same
as in Fig. 1. The contours (0.04, 0.09, 0.19, 0.4, 0.6 of maximum) were
selected to bring out the substructures, the same contours are shown in both
panels and in all other plots. For the model image the main nebular axis is
inclined to the sky plane at 40 deg and because of lack of background noise
the lowest contour is set at 0.01 level.

The components of the model discussed below are presented in Fig. 3 as a
cross-section of the density distribution taken along the presumed symmetry
axis. The values of the hydrogen number density are shown in the side bar,
while the size of the nebula in cm is based on assuming the Galactic Bulge
distance of 8 kpc.

The central region within a radius of $2.4\times 10^{16}$ cm was not
considered in the modelling and assumed empty. The ionizing source (central
star) was assumed a black-body with temperature $T_{\mathrm{bb}}$ of 62 kK
(previously it was estimated 52 kK).

The main nebular body shows nearly straight edges seen clearly in vertical
direction in the HST image in Fig. 1 (especially as the contours on the right
side in the left panel of Fig. 4). This cannot be the ionization boundary
because we do not see there the low-excitation [N ii] 658 nm emission. To
reproduce this feature we assumed a cylindrical outer structure similar to
that derived previously in Gesicki et al. (2016). The difference concerns the
edges which previously were spreading outwards and had higher density. Now we
found a reasonable fit assuming a plain cylinder fully filled with gas density
decreasing along the axis.

The nebular model structure and inclination was fine-tuned together with the
velocity field: both factors provide the emission lines profiles (Sect. 3.2)
which verify the model. The adopted, enhanced, constant density of the middle
part of the main cylinder is responsible for the high intensity of the
symmetric component of the H$\alpha$ emission profiles. The asymmetric parts
of emission profiles of H$\alpha$ and [N ii] (see Figs. 5, 6) originate in the
proposed inner bar and indicate that it is inclined to the plane of the sky at
an angle of 40 deg (previously it was estimated at 35 deg).

The two emission blobs seen in the low-excitation line [N ii] 658 nm (Fig. 2
middle panel) indicate that the ionization front is positioned just at that
location. To obtain an ionization front at a specified position there should
be enough material between this front and the ionizing source. The axial
symmetry of the observed feature guided our modelling towards a narrow and
dense bar-like structure. Sufficient for a rather good fit was the assumption
of a constant density with a value adopted to produce the ionization front at
both far ends. Within this setup the middle region of this bar reproduces the
He ii 468 nm unresolved central emission which requires intense ionization
flux. Were the [N ii] 658 nm lines formed in two separate distant blobs the
explanation of the unresolved He ii 468 nm emission would need a more
elaborated nebular structure.

To reproduce the asymmetric H$\alpha$ image we attached a plume to the bar-
like structure, defined as a quarter of a circle of fixed thickness and with
gas density decreasing away from the symmetry axis. This structure has the
same kinematics as the bar-like outflows, therefore it is natural to place it
near the axis and perpendicularly to the equatorial plane. The asymmetric
plume-like structure extends to $\pm 2.4\times 10^{16}$ cm above and below the
plane of the image shown in Fig. 3.

To complete the model we added an outer equatorial low density torus.

The total nebular mass of the proposed model is 0.09 M☉; the central dense bar
has a mass of 0.006 M☉ while the plume-like feature has a mass of 0.003 M☉.
The chemical composition adopted for the photoionization modelling was taken
from Chiappini et al. (2009) and supplemented by the latest data from Ventura
et al. (2017). The line ratios used for verification of the modelling were
taken from Górny et al. (2009).

In Table 1 we compare the new model parameters with our previous analyses from
2014 and 2016 and with selected, representative observational data adopted
from publications. We did not perform any automated search for a best-fit
model such as the genetic algorithm used previously with spherically symmetric
computations (Gesicki et al., 2006), as the model is now too complicated for
that. We started with the previously found basic parameters and searched for
an improvement (with trials and errors, not far from the starting values), all
the time verifying the agreement with the observed data. The nebular H$\beta$
flux was similarly well fitted with the older simpler codes, however details
of line ratios are significantly better reproduced with the new 3D model. This
in particular concerns the strong [O iii] 500 nm emission which is now excited
by the higher temperature of the central star. Previously, lower value of
$T_{\mathrm{bb}}$ was a compromise to fit high-excitation [O iii] 500 nm
simultaneously with low-excitation [N ii] 658 nm for spherical symmetry. This
conflict is now resolved by the spatial distribution of matter within the PN.

Table 1: Parameters of the models of PN H 2-18 obtained with different codes compared with observed values. The line intensities are on the scale H$\beta$=1 | m o d e l l e d | observedaa$a$dereddened line intensities from Górny et al. (2009)
---|---|---
| 2014bb$b$from Gesicki et al. (2014) | 2016cc$c$from Gesicki et al. (2016) | 2024 |
model version | sphere | ellipsoid | 3D |
ion. mass [M☉] | 0.27 | 0.051 | 0.096 |
kin. age [k yr] | 1.67 | 1.28 | 2.4 |
$T_{\mathrm{bb}}$ [kK] | 51 | 52 | 62 |
$L/L_{\sun}$ | 67 000 | 1 700 | 5 000 |
$\log F(\mathrm{H}\beta$) | $-11.5$ | $-11.6$ | $-11.5$ | $-11.5$
He ii 468 nm | 0.078 |  | 0.058 | 0.05
[N ii] 658 nm | 0.0026 | 0.07 | 0.097 | 0.078
[O i] 630 nm | 5$\times 10^{-8}$ | 8$\times 10^{-5}$ | 0.002 | 0.01
[O ii] 372 nm | 0.005 | 0.07 | 0.13 | 0.16
[O iii] 500 nm | 7.8 | 5.24 | 10.6 | 13.1
[Ne iii] 386 nm | 0.70 |  | 0.57 | 0.999
[S ii] 673 nm | 5$\times 10^{-6}$ |  | 0.013 | 0.018
[S iii 631 nm | 0.0002 |  | 0.009 | 0.01
[Cl iii] 551 nm | 0.007 |  | 0.081 | 0.0036
[Ar iv] 471 nm | 0.047 |  | 0.031 | 0.027

With the adopted parameters the resulting best fitting model produces the
image in the H$\alpha$ line in the plane of the sky presented in the right
panel of Fig. 4 compared with the HST image in the left panel. Both images are
rendered as an inverted colour map to better expose the faint structures and
with contours added to guide the eye. It is not perfectly fitted and it is not
unambiguous however it explains all the different observed characteristics and
is as simple as possible.

###  3.2 Velocity field

The ARGUS/IFU spectra allow much better for the velocity reconstruction than
previously. However the abundance of data required some ingenuity to present
clearly the 2D spectra and to compare observations with models. Fig. 5 shows
the presentation designed by us.

The left panel shows the observations, and the right panel the model results.
The panels cover a fragment of the sky plane, the corresponding HST image
(Fig. 1) is overlaid as thin gray contours; the coordinates are expressed in
units of $10^{17}$ cm.

Both panels are divided into $10\times 10$ small squares representing the IFU
pixels. Each square pixel shows the detailed emission profile of the H$\alpha$
line, with the horizontal axis expressed in km s-1 over the range $\pm$75 km
s-1 and the vertical axis extending from 0 to 1, normalized to the maximum
value of the whole image. To emphasize the value and direction of the velocity
we applied colour coding. As usual red means moving away from the observer
(red-shifted) and blue moving towards the observer. The higher the velocity
the more intense the colour; gray is for zero. We plotted the emission in the
form of bar-plots to show better the colours. The width of each bar is 10 km
s-1, the gray bar is centered on zero velocity.

In the observational panel on the left, we cut-off the noise at the 5% level.
In the model panel on the right we added a correction for the seeing. The
seeing was below the requested 0.9 arcsec for most of the observations, except
for the He ii 468 nm line. Inspecting the modelled emissions integrated over
the exact pixel size it became obvious that we should allow for additional
emission from some area around the given pixel which should mimic the actual
seeing. We found a satisfactory fit (comparable to observations) when
integration for the given pixel was extended by one pixel width (0.52 arcsec)
around.

In Fig. 6 the emission line [N ii] 658 nm is presented in the same setup as in
Fig. 5. Already at first glance we see that red colour dominates in the upper
half of the images, blue dominates in the lower half while the middle row is
approximately symmetric. It is obvious that we see elongated structure and if
we assume only expansion from the central source, then the upper part is
directed away from us while the lower part is expanding towards us. This also
clarifies that the polar axis corresponds to the vertical direction in the
image, something that is not obvious from the images only.

In Fig. 7 the emission line He ii 468 nm is presented in the same setup as in
Fig. 5. Although of lower quality this line is unsplit and limited to the
central pixels, in contrast to the [N ii] emission.

Figure 5: The detailed emission profiles of the H$\alpha$ line shown for each
of IFU pixels. The left panel shows the observed data, the right panel the
corresponding calculated model emissions. To emphasize the value and direction
of the velocity we applied colour coding where red means gas moving away from
the observer and blue moving towards the observer and the higher the velocity
the more intense the colour, gray is for zero. We plotted the emission in a
form of bar-plots, the width of each bar is 10 km s-1, the height (intensity)
is normalized to the maximum value over the whole image. The horizontal size
of each box extends at $\pm$75 km s-1, vertical from zero to unity.

Figure 6: The detailed emission profiles of the [N ii] 658 nm line shown for
each of IFU pixels. The presentation is the same as in Fig. 5

Figure 7: The detailed emission profiles of the He ii 468 nm line shown for
each of IFU pixels. The presentation is the same as in Fig. 5

The left panels of Figs. 5–7 show clearly that no extraordinary kinematics is
detected within H 2-18. The velocities are within normal ranges for PNe, and
there is no evidence for a fast-moving and/or collimated jet.

The data allows us to improve on the velocity model of Gesicki et al. (2016)
where the velocity was linearly increasing with radial distance. Now we have
spectral coverage over the whole nebula. Interestingly the velocities of the
different substructures identified previously, as seen in the H$\alpha$ line,
are similar to each other. This suggests an isotropic velocity field. Indeed,
attempts to apply different velocities to the different substructures failed
to improve the fit, in agreement with this. (There may be differences at
smaller scales: we are limited by the spatial resolution of the IFU data
cube.)

A radial gradient of the expansion is obviously present because the high
excitation He ii 468 nm line shows the lowest velocity while the low
excitation [N ii] 658 nm line has the highest. A very simple relation of
velocity linearly increasing with distance (homologous expansion) was not
satisfactory because it makes the [N ii] 658 nm line too broad in comparison
with observations. The region of the [N ii] line formation should rather have
a constant velocity. We adopted a velocity monotonically increasing with
radial distance, starting from 10 km s-1 at the inner edge and reaching an
asymptotic value of 60 km s-1 near the ends of the bar-like structure. In this
way we obtained a remarkable fit to the broad, split, and asymmetrical
H$\alpha$ emission and simultaneously to the narrow, fast and oppositely
directed N ii emissions and spatially unresolved, narrow, unsplit He ii 468 nm
line, as indicated by the high spectral resolution of the IFU data.

##  4 Discussion

###  4.1 New details and improved parameters

The collected IFU spectra provided new information about the Galactic Bulge
nebula H 2-18. Although of lower spatial resolution than HST they revealed a
feature proposed to be a bar-like structure embedded inside the main object
and co-axial with it. We call it a ‘bar’ rather than ‘jet’ because its
relatively low velocity is similar to that of the surrounding broader
cylinder. The adjoining single-sided plume of ejected gas inclined to the main
axis was earlier seen in the HST images. Its kinematics is also not different
from the rest of the nebula. The ears-like extensions seen in the HST image
are now simply explained as part of the outer equatorial torus.

The photoionization modelling in 3D allowed for improved estimation of stellar
and nebular parameters. The higher $T_{\mathrm{bb}}$ now better reproduces the
line ratios, the advantage of 3D density structure is obvious in this context.
The $T_{\mathrm{bb}}$ value indicates a not very old object, likely in the
middle of its way towards White Dwarf cooling track. The size of this PN
agrees with that: the kinematic age (derived from mass averaged velocity of 43
km s-1, lobes lengths of 0.15 pc, and the corrected formula of Gesicki et al.
(2014)) is 2400 yr.

###  4.2 Co-axial structures

Montoro-Molina et al. (2024) analysed the integral field spectroscopy of the
PN A 58 and derived a structure of wide bipolar outflows with an embedded
narrow co-axial collimated stream. Their kinematics (bipolar outflows at 280
km s-1) is different from ours but the morphology looks nearly the same
(certainly excluding the plume). A common envelope influence is suggested.
Their analyses were performed with the well known tool named
SHAPE111https://wsteffen75.wixsite.com/website which is focused on morphology
and does not consider the photoionization. A thin bar-like structure in the
polar direction is also seen in the JWST image of the planetary nebula NGC
3132 (De Marco et al., 2022), shown in Fig. 8. The equatorial ring is seen in
absorption against the bar, which indicates that the bar is located inside the
nebula. Both cases show that the co-axial inner and outer nebula in H 2-18 is
not a unique phenomenon. (A thin bar-like feature in JWST images of NGC 6720
(Wesson et al., 2024) is a different structure, as it is mainly seen as a gap
and is located in the equatorial direction: it is interpreted as the inner
edges of a wide bipolar flow.)

Figure 8: JWST image of NGC 3132 in the filter NIRcam F187N dominated by
Pa$\alpha$ emission, showing a bar-like, polar structure with some similarity
to that in H 2-18. Image from De Marco et al. (2022).

Comparison with published hydrodynamical calculations is hampered by the fact
that very few of them show the common envelope evolution outcome on a
planetary-nebula wide scale and are rendered as images (most of them are
limited to the solar radius or a.u. scale and they are usually presenting
intricately shaped density cross-sections). García-Segura et al. (2018)
present a common-envelope model that is extended into the PN phase. It assumes
a pre-existing equatorial density enhancement into which the common envelope
ejecta expand. It is restricted to axial symmetry with numerous
simplification, and is focused on explaining bipolar shapes while skipping
kinematics. Their Fig. 8 shows two examples of emission measures (a proxy for
images) for complicated elongated nebulae tilted by 40 deg which show
structures similar to those in H 2-18: the outer barrel shape, central
ellipsoidal enhancements, and symmetric extensions (ears). One of their models
shows blobs of gas piled-up at the symmetry axis which look similar to the
symmetric N ii emissions found by us. Interestingly two separate mass ejection
events are not required for the formation of both inner and outer lobes. A
follow-up study (García-Segura et al., 2022) concluded that lower mass stars
will produce more elongated PNe and the observed jets in PNe must be remnants
of early phases. Both conclusions might be applicable to our case.

###  4.3 Low velocity outflow

The spherical models of Gesicki et al. (2014) commonly revealed a velocity
increase towards the outer nebular edge. This is in agreement with spherically
symmetric hydrodynamical modelling of Schönberner et al. (2005). Here we
obtained a much more complex 3D structure, still with a similar velocity
field. We can speculate that the same acceleration mechanism known from
hydrodynamical models is acting here: the slowly spreading gas is accelerated
and heated by the increasing ionizing radiation of the gradually heating-up
central star.

One of the conclusions of García-Segura et al. (2018) was that the inner
nebulae can be ablated and photoevaporated from the excretion disk. One can
expect that such a process can result in a mild acceleration of the PN. A disk
wind, caused by the evaporation of the photoionized gas was our interpretation
of low velocity outflows in PN M 2-29 (Gesicki et al., 2010). Recently Icke
(2022) analysed a circumstellar ring irradiated by the central star and
obtained that star-driven evaporation would produce a cylindrically collimated
outflow. The simple model was intended to reproduce the very elongated PN Hen
3-401. Interestingly the cylindrical outflow appears to be relatively slow
(below 10 km s-1 near the irradiated disk) and presumably weakly accelerated
along the axis. In our Fig. 1 we do not see any obvious disk comparable to the
one seen in Hen 3-401, nevertheless some similarities exist. Our model of H
2-18 (shown in Fig. 3) has lengths comparable to Hen 3-401 and approximately
twice the width (Icke, 2022). The structure requires higher density in the
equatorial region which can be treated as an approximation to a true
equatorial disk or an equatorial density enhancement. Therefore, photo-
evaporation from this region with further photo-acceleration are a plausible
mechanism shaping H 2-18.

These models are not perfect matches for H 2-18: the lobes in Icke (2022) have
low mass, while for the models of García-Segura et al. (2018) most of the
common envelope mass does not reach escape velocities and remains bound. The
latter model also does not present the kinematics.

Circumbinary disks in PNe have low mass (e.g., De Marco et al., 2022) and may
have difficulty explaining the dense cylindrical structure in H 2-18. A
larger-sized equatorial density enhancement may however play a similar role.

###  4.4 Jets or not

Akashi & Soker (2021) propose that structures such as the ‘ears’ in H 2-18
form through polar jets, while García-Segura et al. (2018) form them by
equatorial gravitational focusing. In the case of H 2-18, the ‘ears’ are
presumably located in the equatorial plane so that the second mechanism is
more likely.

The cylindrical polar structure is more ambiguous. In García-Segura et al.
(2018), it forms through expansion of the common envelope ejecta into a pre-
existing equator-to-pole density gradient. Icke (2022) shows that the
cylindrical nebula of Hen 3-401 can originate from an evaporating disk (or
torus?), without requiring a jet. In contrast, Soker (2002) shows that
cylindrical structures in PNe can form through refocused jets. The paper uses
a jet velocity of 500 km s-1, as could come from a main sequence companion.
Such velocities are not seen in H 2-18, nor is a jet bow shock seen, but the
structure could have formed from a previous jet which is now extinct.

###  4.5 Binary central star

There is no direct evidence of a close binary star in H 2-18 centre,
nevertheless the axial symmetry of the dominant structures suggests the
possibility of a binary system shaping this PN. All the discussed
hydrodynamical models require binarity. However, none appear to depend on
common-envelope evolution. Wider binary interaction which avoids common
envelope can still produce equatorial density enhancement, which may suffice.
A binary (multiple?) system with a disk can be well hidden in the central
unresolved dense region.

###  4.6 Asymmetric plume

The one-sided plume-like structure can be compared with the one observed in
the PN M 2-29 (Gesicki et al., 2010), however in H 2-18 it expands a little
faster, it is seen at different angle, and is located inside the main nebula
although protruding a little out of it. Its flat structure with a nearly
constant density was proposed on the basis of the rather uniform intensity
seen in H$\alpha$ image; any plowing action of the ionization front is not
seen. The mass and brightness of the plume are comparable to those of the bar
which makes it difficult to disentangle both components. This simple model
does not reproduce all details so the true shape may be more complicated. The
models of García-Segura et al. (2018) show structures with a rough similarity
to the plume, although with point-symmetric multiplicity that is lacking in H
2-18.

The mass loss shaping of H 2-18 has in any case been a complex process. The
high density of the inner bar or cylinder, the plume and the lack of a
detected jet put significant constraints on its evolution.

###### Acknowledgements.

We thank Vincent Icke the referee for friendly and helpful comments. Help from
Roger Wesson during installing and running the MOCASSIN code is gratefully
acknowledged. We acknowledge financial support from the Nicolaus Copernicus
University through the University Centre of Excellence ‘Astrophysics and
astrochemistry’. Part of this work was supported by STFC through grant
ST/X001229/1. AAZ also acknowledges support from the Royal Society through
grant IES/R3/233287 and the University of Macquarie. This work made use of
Astropy:222http://www.astropy.org a community-developed core Python package
and an ecosystem of tools and resources for astronomy (Astropy Collaboration
et al., 2022).

## References

  * Akashi & Soker (2021) Akashi, M. & Soker, N. 2021, ApJ, 913, 91
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167
  * Chiappini et al. (2009) Chiappini, C., Górny, S. K., Stasińska, G., & Barbuy, B. 2009, A&A, 494, 591
  * De Marco et al. (2022) De Marco, O., Akashi, M., Akras, S., et al. 2022, Nature Astronomy, 6, 1421
  * Ercolano et al. (2003) Ercolano, B., Barlow, M. J., Storey, P. J., & Liu, X. W. 2003, MNRAS, 340, 1136
  * García-Segura et al. (2018) García-Segura, G., Ricker, P. M., & Taam, R. E. 2018, ApJ, 860, 19
  * García-Segura et al. (2022) García-Segura, G., Taam, R. E., & Ricker, P. M. 2022, MNRAS, 517, 3822
  * Gesicki & Zijlstra (2000) Gesicki, K. & Zijlstra, A. A. 2000, A&A, 358, 1058
  * Gesicki et al. (2006) Gesicki, K., Zijlstra, A. A., Acker, A., et al. 2006, A&A, 451, 925
  * Gesicki et al. (2014) Gesicki, K., Zijlstra, A. A., Hajduk, M., & Szyszka, C. 2014, A&A, 566, A48
  * Gesicki et al. (2016) Gesicki, K., Zijlstra, A. A., & Morisset, C. 2016, A&A, 585, A69
  * Gesicki et al. (2010) Gesicki, K., Zijlstra, A. A., Szyszka, C., et al. 2010, A&A, 514, A54
  * Górny et al. (2009) Górny, S. K., Chiappini, C., Stasińska, G., & Cuisinier, F. 2009, A&A, 500, 1089
  * Icke (2022) Icke, V. 2022, Galaxies, 10, 53
  * Montoro-Molina et al. (2024) Montoro-Molina, B., Tafoya, D., Guerrero, M. A., Toalá, J. A., & Santamaría, E. 2024, A&A, 684, A107
  * Schönberner et al. (2005) Schönberner, D., Jacob, R., Steffen, M., et al. 2005, A&A, 431, 963
  * Soker (2002) Soker, N. 2002, ApJ, 568, 726
  * Ventura et al. (2017) Ventura, P., Stanghellini, L., Dell’Agli, F., & García-Hernández, D. A. 2017, MNRAS, 471, 4648
  * Wesson et al. (2024) Wesson, R., Matsuura, M., Zijlstra, A. A., et al. 2024, MNRAS, 528, 3392


# Layered Decoding of Quantum LDPC Codes

Julien Du Crest111Univ. Grenoble Alpes, Grenoble INP, LIG, F-38000 Grenoble,
France (julien.du-crest@univ-grenoble-alpes.fr), Francisco Garcia-
Herrero222Department of Computer Architecture and Automatics (DACYA),
Complutense University of Madrid, Spain<EMAIL_ADDRESS>, Mehdi Mhalla333
Univ. Grenoble Alpes, CNRS, Grenoble INP, LIG, F-38000 Grenoble, France
(mehdi.mhalla@univ-grenoble-alpes.fr),
Valentin Savin 444 Univ. Grenoble Alpes, CEA-Léti, F-38054 Grenoble, France
<EMAIL_ADDRESS>and Javier Valls 555Instituto de Telecomunicaciones y
Aplicaciones Multimedia, Universitat Politecnica de Valencia, 46022 Valencia,
Spain<EMAIL_ADDRESS>

###### Abstract

We address the problem of doing message passing based decoding of quantum LDPC
codes under hardware latency limitations. We propose a novel way to do layered
decoding that suits quantum constraints and outperforms flooded scheduling,
the usual scheduling on parallel architectures. A generic construction is
given to construct layers of hypergraph product codes. In the process, we
introduce two new notions, $t$-covering layers which is a generalization of
the usual layer decomposition, and a modified scheduling called random order
scheduling. Numerical simulations show that the random ordering is of
independent interest as it helps relieve the high error floor typical of
message passing decoders on quantum codes for both layered and serial decoding
without the need of post-processing.

##  1 Introduction

A lot of work has been done in order to improve the decoding of quantum low-
density parity-check (qLDPC) codes using message-passing (MP) decoders. Most
of these works rely on the use of post-processing techniques [1, 2, 3], whose
feasibility is still to be demonstrated on actual hardware, due to the
stringent latency, power and scalability requirements of the quantum system. A
key attribute of MP decoding is the underlying scheduling, indicating the
order in which variable and check node messages are updated. This has been
subject to extensive research in the classical LDPC decoding literature, and
it has been shown that the MP scheduling may significantly impact the
convergence speed [4], the decoding performance (_e.g._ , in case of adaptive
scheduling strategies [5, 6, 7]), or the performance (_e.g._ , latency, area,
power-consumption) of the hardware design [8]. The vast majority of hardware
designs are based on partly-parallel architectures, implementing a layered
decoding scheduling, which can be considered as a de facto standard solution,
able to provide relevant complexity and performance advantages in most
applications [8].

For qLDPC codes, the MP decoding performance may depend even more on the
underlying scheduling, which can be most likely attributed to the code
degeneracy [2]. Moreover, some post-processing techniques may be highly
dependent on the MP decoding scheduling. For instance, the order statistics
decoding post-processing has been shown to provide very good performance when
a layered scheduling is used, but its performance may be drastically degraded
using a flooded (_i.e._ , fully parallel) scheduling [2].

To design an efficient partly parallel architecture implementing a layered
scheduling, one needs a layer decomposition of the parity check matrix. For
qLDPC codes this may be tricky, as they do not have an innate decomposition
into horizontal layers (as for instance in the case of classical quasi-cyclic
LDPC codes). To ensure a high degree of parallelism, it is also desirable to
have a decomposition into a minimal number of layers. In this paper, we first
give a generic construction of a minimal layer decomposition for hypergraph
product codes. Moreover, in an attempt to start bridging the gap between
hardware limitations and state of the art MP decoders, we propose two new
tools to implement layered decoding of qLDPC codes. The first is a
generalization of the notion of layer decomposition, consisting of a family of
$t$-covering layers, which can be seen as a layer decomposition of $t$
decoding iterations, and is aimed at increasing the parallelism degree of the
layered architecture. The second is a new scheduling called random order
scheduling, and is shown to significantly improve the decoding performance.
Our numerical simulations provide evidence that both could be used in the
future to meet hardware needs as they offer a good compromise of speed and
performance.

##  2 Preliminaries

###  2.1 Quantum Codes

Calderbank-Shor-Steane (CSS) codes are defined by two classical
$(m_{X}/m_{Z},n)$-parity check matrices $H_{X},H_{Z}$ satisfying
$H_{X}H_{Z}^{\perp}=0$. The dimension of the quantum code is
$n-\text{rank}(H_{X})-\text{rank}(H_{Z})$ and its minimum distance
$d=\\{\min|v|,v\in\ker H_{X}\backslash\operatorname{im}H_{Z}^{t}\cup\ker
H_{Z}\backslash\operatorname{im}H_{X}^{t}\\}$. One such class of CSS codes are
the hypergraph product codes (HPC), which given two classical codes $A$ and
$B$ give the quantum code $H_{X}=[A\otimes I,I\otimes B^{t}]$ and
$H_{Z}=[I\otimes B,A^{t}\otimes I]$ (see [9] for the construction and
parameters of the code).

In the following, we will only focus on decoding $Z$ errors using $H_{X}$. All
proofs are easily adapted to correcting $X$ errors.

###  2.2 MP decoding

MP decoders work by exchanging soft information between check and variable
nodes on the Tanner graph representation of the code, trying to converge to a
hard decision on the variable nodes that satisfy the syndrome. One crucial
factor is to decide in which order messages are exchanged and soft information
updated. There are three main decoding scheduling used classically: flooded,
in which messages are exchanged simultaneously and soft information updated in
parallel, serial, where the graph is updated sequentially going through all
the checks one by one, and layered, which lies in between, taking advantage of
checks that have a disjoint support to update them in parallel, essentially
doing a speed up of serial scheduling at no cost. For more details on
classical message passing, refer to [10].

###  2.3 Layered Scheduling

To avoid memory conflicts in a partly parallel architecture, implementing a
layered scheduling, the same memory slot should not be read/written to by two
different processing units at the same time. This motivates the following
definitions.

######  Definition 1.

A layer is a collection of check-nodes such that any two check-nodes have no
neighbouring variable-node in common.

######  Definition 2.

A layer decomposition $L_{0}\sqcup\dots\sqcup L_{k-1}$ is a partition of the
set of check-nodes into $k$ layers.

######  Definition 3.

A decomposition is said to be minimal if it is impossible to find a
decomposition in less layers.

A simple density argument is enough to state the following fundamental
inequality:

######  Lemma 1.

Any $k$ layers decomposition of a $(\\_,\delta)$-regular666A matrix is said to
be $(c,d)$-regular if every column is of weight $c$ and every row of weight
$d$. $(m,n)$-parity-check matrix satisfies $k\geq\frac{\delta m}{n}$.

######  Definition 4.

A decomposition is $\gamma$-balanced if

$\frac{|L_{i}|}{|L_{j}|}\leq\gamma,\quad\forall L_{i},L_{j}$

A decomposition will be said to be balanced if it is 1-balanced. Balanced
decompositions ensure an efficient use of hardware resources (check-node
processing units).

##  3 Hardware Requirements

In contrast to classical LDPC decoders, which prioritize optimizing
throughput, qLDPC ones must satisfy highly constrained values of latency to
avoid the backlog problem [11], which would lead to an exponential slowdown of
the quantum processor making the QEC implementation impractical.

Table 1: Latency approximation for the different architectures Parallel | Serial | Layered
---|---|---
$T_{min}^{(P)}\times 2\times it_{max}$ | $T_{min}^{(S)}\times(it_{max}/2)\times m$ | $T_{min}^{(L)}\times 2\times(it_{max}/2)\times k$

To illustrate the behavior of the decoder, the B1 code from [1] is taken as an
example (defined in Section 4). An MP decoder for the B1 code can achieve with
this architecture777All figures reported here come from our implementation of
an either fully parallel [12] or serial min-sum decoder architecture, with
exchanged messages quantized on 6 bits, on a Xilinx FPGA xcv095 board., a
clock period between 8 and 10ns which derives a latency between 480ns and
600ns at 30 iterations, that is close to the most constrained technology.
Taking this into account, the clock period usually can be reduced to 70% and
80% of the clock period obtained with the parallel version. The question is
that with this schedule and the derived architecture only one check node is
updated in a clock cycle, because of the sequential update of the messages.
Due to this, at least $m$ clock-cycles are required 888Assuming that due to
the reduced complexity of the units both the check node and the connected
variable nodes can be updated in parallel. to complete just one iteration of
the MP algorithm. Following the example of code B1, the clock period should be
between 5.6 and 7ns, but the total latency at 30 iterations would be between
74.26$\mu$s and 92.82$\mu$s. Even assuming a reduction of the number of
iterations to get similar performance to the flooded schedule, the range of
latencies would be out of the time budget of supercomputing qubits and
transmons. To meet the timing requirements, the clock period should be equal
to $\frac{5.6}{m/2}$=0.025ns, which is a maximum clock frequency of 40GHz.
This frequency cannot be achieved by any FPGA or ASIC, and on the other hand,
it would require a large power consumption that will cause another problem
with the power budget and the refrigeration system [13]. With the previous
examples, it is easy to conclude that the implementation of serial scheduling,
even if it has a better performance than the flooded one, it is not a
realistic solution when it comes to implementation. A trade-off solution
between both serial and flooded may be the layered schedule. If the number of
layers is small enough, the number of clock cycles per iteration will not grow
too much and the number of iterations will be usually reduced by two. Going
back to the B1 code example, assuming that in the worst case, the clock period
will be similar to the parallel architecture, and with a distribution of 3.5
layers999See Section 4 for the formal definition of fractional layer number.
the total latency for 30 iterations will be between 8x3.5x2x30=1.68$\mu$s and
10x3.5x2x30=2.10$\mu$s; and between 8x3.5x2x15=840ns and
10x3.5x2x15=1.05$\mu$s with 15 iterations, which is fairly close to the
constraints of superconducting qubits and meets the requirements of other
technologies.

As we will see in the following sections, the layered schedule will also
benefit from some non-negligible performance improvements, apart from the
reduction in the number of iterations, compared to the flooded schedule.

In table 1, we can find a summary of the total latency for different
architectures, where $T_{min}^{(P)}$, $T_{min}^{(S)}$ and $T_{min}^{(L)}$ are
the minimum clock periods achievable by the parallel, serial and layered
architectures respectively, and $it_{max}$ is the maximum number of iterations
configured in the parallel decoder. Note that we assume that the number of
iterations of serial and layered is usually half the number of iterations of
the parallel architecture [4].

##  4 Generic Constructions

###  4.1 Layered Construction for Hypergraph Product Codes

Consider an hypergraph product code defined by two matrices $A$ and $B$, such
that $H_{X}=[A\otimes I,I\otimes B^{t}]$ and $H_{Z}=[I\otimes B,A^{t}\otimes
I]$ [9]. There exist a layer construction from a layer decomposition of
$A,B,A^{t},B^{t}$.

######  Theorem 1.

Given minimal decompositions $A=A_{0}\sqcup\dots\sqcup
A_{k_{A}-1},B=B^{t}_{0}\sqcup\dots\sqcup B^{t}_{k_{B^{t}}-1}$, one can
construct a minimal decomposition of $H_{X}$ in $k=\max(k_{A},k_{B^{t}})$
layers.

A similar theorem can be stated for $A^{t},B$ and $H_{Z}$ with the same proof
techniques.

#### Construction

If $k_{A}\neq k_{B}$, without loss of generality, suppose $k_{A}<k_{B}$. The
first step is to add empty layers to A so that
$k^{\prime}_{A}=k^{\prime}_{B}=k$. That is let $A=A_{0}\sqcup\dots\sqcup
A_{k_{A}-1}\sqcup A_{k_{A}}\dots\sqcup A_{k}$ where
$A_{k_{A}}=\dots=A_{k}=\emptyset$. Let’s label each row of $H_{X}$ as

$a\varoast b:=[a\otimes e_{b},e_{a}\otimes b],\quad
a\in\text{rows}(A),b\in\text{rows}(B^{t})$

In the following we will denote by left the sub-matrix $[A\otimes I]$ and
right the sub-matrix $[I\otimes B^{t}]$. Create layers $L_{0}\dots L_{k-1}$
such that

$(a\varoast b)\in L_{i}\Leftrightarrow\exists j\quad a\in A_{j},b\in
B^{t}_{j+i\mod k}$ (1)

By definition, all checks belong to some layer, we now have to check that any
two checks in a given layer have disjoint variable nodes support. Suppose that
two checks $a\varoast b,a^{\prime}\varoast b^{\prime}$ belong to $L_{i}$. Case
A : $a=a^{\prime}$. They do not touch on the left thanks to the tensor product
with the identity. Furthermore, it means that $b\neq b^{\prime}$ but then both
belong to $B^{t}_{l}$ for some $l$, so they have disjoint support on the
right. Case B,C : $a\neq a^{\prime}$. They have disjoint support on the right
because of the tensor product with the identity. To show that they do not
intersect on the left, there are two cases : If $a$ and $a^{\prime}$ belong to
some $A_{l}$ (case B), then by definition they have disjoint support on the
left. If $a$ and $a^{\prime}$ belong respectively to $A_{l},A_{l^{\prime}}$
with $l\neq l^{\prime}$ (case C), then it means that $b\in B^{t}_{l+i\mod
k},b^{\prime}\in B^{t}_{l^{\prime}+i\mod k}$, two distinct classes. Hence even
though $a$ and $a^{\prime}$ might share variable nodes in $A$, they do not
intersect in the tensored version $A\otimes I$. Fig. 1 depicts a simple
example of the 3 cases.

1 |  |  | 1 |  |  | 1 | 1 |  |  |  |
---|---|---|---|---|---|---|---|---|---|---|---
| 1 |  |  | 1 |  | 1 |  |  |  |  |
|  | 1 |  |  | 1 |  | 1 |  |  |  |
\raisebox{-.9pt} {1}⃝ |  |  |  |  |  |  |  | \raisebox{-.9pt} {1}⃝ | \raisebox{-.9pt} {1}⃝ |  |
| 1 |  |  |  |  |  |  | 1 |  |  |
|  | 1 |  |  |  |  |  |  | 1 |  |
|  |  | \raisebox{-.9pt} {1}⃝ |  |  |  |  |  |  | \raisebox{-.9pt} {1}⃝ | \raisebox{-.9pt} {1}⃝
|  |  |  | 1 |  |  |  |  |  | 1 |
|  |  |  |  | 1 |  |  |  |  |  | 1

Case A : $a=a^{\prime}\implies b\neq b^{\prime},\exists l,\quad
b,b^{\prime}\in B^{t}_{l}$
---
Case \raisebox{-.9pt} {B}⃝ : $a\neq a^{\prime},\quad a,a^{\prime}\in A_{l}$
Case C : $a\neq a^{\prime},a\in A_{l},a^{\prime}\in A_{l^{\prime}}$

Figure 1: Small visualization example of proof cases where $A=B^{t}$

#### Minimality

The proof is by contradiction. Assume that there is a decomposition in less
than $k_{B^{t}}$ layers, then one could recover a decomposition for $B^{t}$ in
less than $k_{B^{t}}$ from a restriction to the $\\{a\varoast b,\quad\forall
b\\}$ positions for any given $a$. Any decomposition in less than $k_{A}$
layers would similarly give a decomposition for $A$ from the restriction of
$H_{X}$ to any $\\{a\varoast b,\quad\forall a\\}$ for a given $b$. Hence the
decomposition in $\max(k_{A},k_{B}^{t})$ is minimal for $H_{X}$.

Note that the construction is not unique, for example, Equation 1 can be
replaced by the following equation where $\sigma$ is any $k$-permutation,
although this is still not the most generic formula:

$(a\varoast b)\in L_{i}\Leftrightarrow\exists j,\quad a\in A_{\sigma(j)},b\in
B^{t}_{j+i\mod k}$ (2)

######  Theorem 2.

Given $k$-layerings for $A$ and $B^{t}$, respectively $\alpha$ and
$\beta$-balanced. Then $H_{X}$ is $\gamma$-balanced with :
(i) $\gamma<min(\alpha,\beta)$ if $\alpha,\beta>1$ (ii) $\gamma=1\qquad$
otherwise.

###### Proof.

(i) Let $a_{0}...a_{k-1}$ be the sizes of layers $A_{0}...A_{k-1}$, and
$b_{0}...b_{k-1}$ the sizes of the layers $B^{t}_{0}...B^{t}_{k-1}$. Then each
layer $L_{l}$ of $H_{X}$ has size $\sum a_{i}b_{i+l}$. We also suppose layers
of $A$ and $B^{t}$ are ordered from biggest to smallest hence $a_{0}=\alpha
a_{k-1}$ and $b_{0}=\beta b_{k-1}$.
The layer $L_{0}$ of size $\sum a_{i}b_{i}$ is the biggest layer, a classical
proof of that is by contradiction, using the fact that $\forall a\geq c,b\geq
d,\quad ab+cd\geq ad+bc$. The ratio between any other layer $L_{j}$ and
$L_{0}$ is smaller than $\beta$ since $\beta\sum a_{i}b_{i+j}>\sum
a_{i}b_{0}>\sum a_{i}b_{i}$ using the fact that $\forall b_{j},\quad\beta
b_{j}\geq b_{0}$ ( and similarly for $\alpha$ ).

(ii) Suppose $A$ is perfectly balanced. In that case, for any $b\in B^{t}$, it
will appear the same number of times in each layer $L_{i}$ since the
$a\varoast b,\forall a\in A$ will be equally balanced in the layers. Hence the
code will be balanced. The same argument holds if $B^{t}$ is perfectly
balanced.

∎

######  Corollary 2.1.

Given a $k_{A}$-layering or $A$, and a $k_{B}>k_{A}$ layering for $B^{t}$
$\beta$-balanced. Then $H_{X}$ is $\gamma$-balanced with :

$\gamma\leq\beta$

###### Proof.

Same as above, considering a $k_{B}$ layering for $A$ by adding empty layers.
This new layering is $\infty$-balanced. ∎

###  4.2 Random Ordering

We introduce a decoding technique called random ordering. This technique
consists of applying a random order on the layers’ application at each
decoding step. This is also generalized to serial decoding by considering that
each check belongs to its own layer (i.e. $k=m$). This seemingly anodyne step
helps to alleviate the error floor quite dramatically. In addition, further
simulations showed us that one does not even have to use a “good” pseudo-
random generator to generate the permutation, and this can be done with
virtually no cost using a simple congruent generator, a solution that is
hardware friendly.

###  4.3 $t$-Covering of Layers

For many codes, the theoretical bound on $k$ given by a density argument is
not tight. However, since for the quantum codes the number of layers is fixed
due to latency constraints, it is important to stay as close as possible to
the theoretical bound. We introduce a generalization of the layer
decomposition called a $t$-covering of ($k$) layers. We drop the requirement
that the layers should be disjoint, and only require that their union taken
with multiplicities should cover each check exactly $t$ times. In the
following, the parameters of a $t$-covers will be specified as $(t,k,\gamma)$,
giving the cover parameters and the balance of the layers. Note that when
using $t$-covers, the usual term of “iteration” becomes ambiguous because it
might be the case that the decoder stops while all the checks have not been
seen the same number of times. Since by pipelining the process, the syndrome
satisfaction could be checked after each layer application adding very low
latency, in the following we will often refer to the number of iterated layers
(but always specify it when we do so). To quickly compare a $t$-cover with
another or with a layer decomposition, it is useful to introduce the
fractional layer number as $\frac{k}{t}$, intuitively it captures the
“average” number of layers the decoder has to process to see each check once.
Finally, by concatenating $t$ times the matrix $H$, it is clear that the
density bound of lemma 1 applies to the fractional layer number. As a simple
application, for the code B1 given below, we found a (2,7,1)-cover,
$\frac{k}{t}=3.5$ . We could also find a (1,4,2) layer decomposition,
$\frac{k}{t}=4$, and the density bound gives us $\frac{k}{t}\geq 3$, since no
decomposition in $3$ layers is known for the B1 code, our $2$-cover sets a new
standard in decoding efficiency.

##  5 Applications on Particular Quantum Codes

###  5.1 C2 Code

The C2 code is a hypergraph product code generated from a single cyclic matrix
($A=B$) of generator polynomial $p(x)=1+x^{2}+x^{5}$ and length $l=31$. Since
this cyclic matrix (and its transpose) accepts a decomposition in $5$-layers,
using the technique from theorem 1, we can construct a 5 layer decomposition
for the C2 code. As said earlier about the balancing effect of the procedure,
the decomposition used for $A,B,A^{t},B^{t}$ is $(1,5,2)$-cover and it yields
a $(1,5,1.1)$-cover for C2. This shows the balancing effect of the procedure,
as we go from $\gamma=2$ to $\gamma=1.1$. Here are the layers used for the
quasicyclic matrix:

$A_{0}$ | 0 1 7 8 14 15 21 22
---|---
$A_{1}$ | 2 3 9 10 16 17 23 24
$A_{2}$ | 4 11 18 25 29
$A_{3}$ | 5 12 19 26 30
$A_{4}$ | 6 13 20 27 28

It should be noted that in order to improve the latency (at the cost of a more
complex construction), we were also able to create a $(224,961,1)$-cover of
C2, achieving a fractional layer number of 4.29 and giving close numerical
results.

###  5.2 B1 Code

(a) $B1[[882,24]]$ SP (b) $C2[[1922,50,16]]$ SP

(c) $B1[[882,24]]$ NMS (d) $C2[[1922,50,16]]$ NMS

Figure 2: Comparison of different decoders and scheduling on B1 and C2 codes
under Z-noise. In the simulations, we use a perturbated NMS, where each check
node message is multiplied by a normalization factor uniformly chosen at
random in {0.875,0.9275} at each iteration. This perturbation is important to
avoid an error-floor degradation.

The B1 code is a Generalized Hypergraph Product code (construction given in
[1], Appendix) : As such it shares similarities with the Hypergraph Product
Codes. Although we do not have a generic decomposition for this family of
codes, some ideas from the hypergraph product theorem apply when creating a
cover for the B1 code. The B1 code accepts a 2-cover in 7 layers, hence giving
a fractional layer number of 3.5. The layers are as follows :

$L_{i}=\\{i+7j\quad\forall j\\}\cup\\{3+i+7j\quad\forall j\\}$

This 2-cover comes from a decomposition of the quasi-cyclic matrix defined by
polynomial $p(x)=1+x+x^{6}$ in 7 layers $S_{0}\dots S_{6}$ such that
$S_{i}=\\{i+7j\quad\forall j\\}$ which is an obvious decomposition (albeit not
minimal) given the generating polynomial. Those layers have the property that
the union of any two layers $S_{i},S_{i+3\mod 7}$ is still a valid layer which
gives the basis for the layers of the code B1. In fact, those layers are
extended to the matrices of $H_{X},H_{Z}$ much in the fashion of what is done
with the hypergraph product but with a “twist” as the blocks are quasi-cyclic
shifts of identities instead of all identities so one has to be more careful
and cannot use the generic formula.

##  6 Numerical Results

Fig. 2 compares the different decoding techniques proposed on an HGP and
Hyperbicycle Codes under Z-noise. We consider Sum-Product (SP) and Normalized
Min-Sum (NMS) decoders, with serial, layered and flooded scheduling [10]. When
decoding classical codes, using serial scheduling yields a factor two
improvement in convergence speed over flooded scheduling. This is not the case
in our numerics on quantum codes, as the serial scheduling suffers from a high
error floor. This error floor can be virtually eliminated by using a random
ordering scheduling. For the flooded scheduling, the number of iterations used
is $i=128$, for serial $i=64$. For the layered scheduling of a
$(t,k,\\_)$-cover, the layer iteration number $i_{lay}=\lfloor
64\times\frac{k}{t}\rfloor$. Although we argued before that checking the
syndrome after each layer is essentially costless, to do a “fairer” comparison
with serial scheduling under random order scheduling we also tried checking
the syndrome only after $\lceil\frac{k}{t}\rceil$ layer iterations. We did not
include the numerics as the two curves match almost perfectly, making it a
non-issue.

On the B1 code, because it is a $t$-cover and not a layer decomposition, we
alter the random ordering scheduling a little bit to boost the performances by
requiring that the permutation is not chosen uniformly at random, and must
satisfy the additional constraint that two successive layers should not share
any check. These additional constraints help the decoder to converge faster as
processing the same check twice in a row in different layers would not change
its soft information.

##  7 Conclusion

We showed how to implement a layered scheduling for qLDPC codes to meet with
the hardware latency limitations. In our numerics, this decoder was more
efficient than what could be achieved using similar resources in flooded
scheduling which might make it the go to hardware option in the future. We
also show that the random order scheduling is a result interesting on its own,
as it can be applied to both serial and layered scheduling to alleviate the
high error floor of some codes without the need of a post-processing. It
should be noted that presently, the best decoders for those codes use some
kind of post-processing after message passing, something that was not studied
in this paper, as none of the known post-processings can meet the hardware
latency considerations. Knowing that our serial scheduling with random
ordering already achieves the performances of the Ordered Statistic Decoding
(OSD) post-processing on those codes101010See [2][Sec 4, Fig. 2], error
probability should be multiplied by 2/3 to compare the two since the error
model is depolarizing noise there.. Finding such hardware friendly post-
processing to use with our layered scheduling would be another step in the
direction we are aiming for.

## Acknowledgement

This work was supported by the QuantERA grant EQUIP (French
ANR-22-QUA2-0005-01, and Spain MCIN/AEI/10.13039/501100011033, grant
PCI2022-132922), and the Plan France 2030 (ANR-22-PETQ-0006) and by the
European Union “NextGenerationEU/PRTR”.

## References

  * [1] P. Panteleev and G. Kalachev, “Degenerate quantum LDPC codes with good finite length performance,” _Quantum_ , vol. 5, p. 585, 2021, arXiv:1904.02703.
  * [2] J. Du Crest, M. Mhalla, and V. Savin, “Stabilizer inactivation for message-passing decoding of quantum LDPC codes,” in _IEEE Information Theory Workshop (ITW)_. IEEE, 2022, pp. 488–493.
  * [3] N. Raveendran and B. Vasić, “Trapping sets of quantum LDPC codes,” _Quantum_ , vol. 5, p. 562, Oct. 2021.
  * [4] J. Zhang, Y. Wang, M. P. C. Fossorier, and J. S. Yedidia, “Iterative decoding with replicas,” _IEEE Transactions on Information Theory_ , vol. 53, no. 5, pp. 1644–1663, 2007.
  * [5] Y. Mao and A. H. Banihashemi, “Decoding low-density parity-check codes with probabilistic scheduling,” _IEEE Communications Letters_ , vol. 5, no. 10, pp. 414–416, 2001.
  * [6] V. Savin, “Iterative LDPC decoding using neighborhood reliabilities,” in _IEEE Int. Symp. on Inf. Theory (ISIT)_ , 2007, pp. 221–225.
  * [7] A. I. Vila Casado, M. Griot, and R. D. Wesel, “LDPC decoders with informed dynamic scheduling,” _IEEE Trans. on Communications_ , vol. 58, no. 12, pp. 3470–3479, 2010.
  * [8] E. Boutillon and G. Masera, “Hardware design and realization for iteratively decodable codes,” in _Channel coding: Theory, algorithms, and applications_ , D. Declercq, M. Fossorier, and E. Biglieri, Eds. Elsevier, 2014, pp. 583–642.
  * [9] J.-P. Tillich and G. Zémor, “Quantum ldpc codes with positive rate and minimum distance proportional to the square root of the blocklength,” _IEEE Transactions on Information Theory_ , vol. 60, no. 2, pp. 1193–1202, 2013, arXiv:0903.0566.
  * [10] V. Savin, “LDPC decoders,” in _Channel coding: Theory, algorithms, and applications_ , D. Declercq, M. Fossorier, and E. Biglieri, Eds. Elsevier, 2014, pp. 211–260.
  * [11] A. Holmes, M. R. Jokar, G. Pasandi, Y. Ding, M. Pedram, and F. T. Chong, “NISQ+: Boosting Quantum Computing Power by Approximating Quantum Error Correction,” in _Proceedings of the ACM/IEEE 47th Annual International Symposium on Computer Architecture_ , ser. ISCA ’20. IEEE Press, 2020, p. 556–569.
  * [12] J. Valls, F. Garcia-Herrero, N. Raveendran, and B. Vasić, “Syndrome-based min-sum vs OSD-0 decoders: FPGA implementation and analysis for quantum LDPC codes,” _IEEE Access_ , vol. 9, pp. 138 734–138 743, 2021.
  * [13] Y. Ueno, M. Kondo, M. Tanaka, Y. Suzuki, and Y. Tabuchi, “QECOOL: On-line quantum error correction with a superconducting decoder for surface code,” in _2021 58th ACM/IEEE Design Automation Conference (DAC)_. IEEE, dec 2021.