# Magnetic Fields Observed along the E-W Outflow of IRAS 16293-2422

Frankie J. Encalada Department of Astronomy, University of Illinois, 1002 West
Green St, Urbana, IL 61801, USA  Leslie W. Looney Department of Astronomy,
University of Illinois, 1002 West Green St, Urbana, IL 61801, USA  National
Radio Astronomy Observatory, 520 Edgemont Rd., Charlottesville, VA 22903 USA
Giles Novak Center for Interdisciplinary Exploration and Research in Astronomy
(CIERA), Northwestern University, 1800 Sherman Avenue, Evanston, IL 60201, USA
Department of Physics & Astronomy, Northwestern University, 2145 Sheridan
Road, Evanston, IL 60208, USA  Sarah Sadavoy Department for Physics,
Engineering Physics and Astrophysics, Queen’s University, Kingston, ON K7L
3N6, Canada  Erin G. Cox Center for Interdisciplinary Exploration and Research
in Astronomy (CIERA), Northwestern University, 1800 Sherman Avenue, Evanston,
IL 60201, USA  Fabio Pereira-Santos  Max Planck Institute for Astronomy,
Königstuhl 17, 69117 Heidelberg, Germany  Dennis Lee Center for
Interdisciplinary Exploration and Research in Astronomy (CIERA), Northwestern
University, 1800 Sherman Avenue, Evanston, IL 60201, USA  Department of
Physics & Astronomy, Northwestern University, 2145 Sheridan Road, Evanston, IL
60208, USA  Rachel Harrison Department of Astronomy, University of Illinois,
1002 West Green St, Urbana, IL 61801, USA  School of Physics and Astronomy,
Monash University, Clayton VIC 3800, Australia  Kate Pattle Department of
Physics and Astronomy, University College, Gower St., London WC1E 6BT, UK

(Accepted March 9, 2024)

###### Abstract

Magnetic fields likely play an important role in the formation of young
protostars. Multiscale and multiwavelength dust polarization observations can
reveal the inferred magnetic field from scales of the cloud to core to
protostar. We present continuum polarization observations of the young
protostellar triple system IRAS 16293-2422 at 89 $\mu$m using HAWC+ on SOFIA.
The inferred magnetic field is very uniform with an average field angle of
89${}^{\circ}\pm$23∘ (E of N), which is different from the $\sim$170∘ field
morphology seen at 850 $\mu$m at larger scales ($\gtrsim$ 2000 au) with JCMT
POL-2 and at 1.3 mm on smaller scales ($\lesssim$300 au) with ALMA. The HAWC+
magnetic field direction is aligned with the known E-W outflow. This alignment
difference suggests that the shorter wavelength HAWC+ data is tracing the
magnetic field associated with warmer dust likely from the outflow cavity,
whereas the longer wavelength data are tracing the bulk magnetic field from
cooler dust. Also, we show in this source the dust emission peak is strongly
affected by the observing wavelength. The dust continuum peaks closer to
source B (northern source) at shorter wavelengths and progressively moves
toward the southern A source with increasing wavelength (from 22 $\mu$m to 850
$\mu$m).

††journal: ApJ

##  1 Introduction

Magnetic fields are thought to be a critical component of the star formation
process on all scales, from the assembly of molecular cloud structure to
protostellar accretion (e.g., Pattle et al., 2023; Tsukamoto et al., 2023).
Nevertheless, the details of how magnetic fields affect star formation are not
well observationally constrained. The magnetic field is difficult to directly
observe in most regions. The most common approach is to use dust polarization
emission to infer the plane of the sky magnetic field orientation. This is
possible due to radiative alignment torques that tend to align the short axis
of elongated dust grains with the magnetic field (e.g., Lazarian, 2007). As a
result, thermal emission of dust grains will be polarized in a direction that
is perpendicular to the field. Many studies therefore rotate the polarization
angles by 90∘ to present the inferred plane-of-sky magnetic field morphology.
Although this technique works very well in most star forming regions and
spatial scales, when resolving the circumstellar disk other polarization
mechanisms, such as scattering, dominate (e.g., Kataoka et al., 2015; Yang et
al., 2016; Cox et al., 2018; Harris et al., 2018; Sadavoy et al., 2019). This
makes it very difficult to measure the magnetic field morphology on small
scales using dust polarization observations.

On the other hand, there have been a few surveys of the magnetic field in the
inner envelope of protostars comparing the average B-field axis with the
outflow axis (Hull et al., 2014; Galametz et al., 2018; Zhang et al., 2014).
Most recently, Huang et al. (2024) surveyed 61 protostars in Orion (with 56
detections) from the Class 0/I identified disk sources in the VANDAM survey
(Tobin et al., 2022). They found that $\sim$40% of the protostars exhibit
B-fields that are perpendicular to the source outflows on scales of 400-1000
au, with the remainder of the sources being consistent with random alignment.

The magnetic field may also play roles in the formation and/or evolution of
binaries and multiple systems. Main sequence stars exhibit a multiplicity rate
that increases with stellar mass: the mean frequency of stellar companions per
primary rises from 0.5 for solar-type main sequence stars to 2.1 for O-type
main sequence primaries (Moe & Di Stefano, 2017). In addition, the
multiplicity fraction is largest in the youngest populations of protostars
(Tobin et al., 2016; Encalada et al., 2021; Tobin et al., 2022), likely
decreasing as the systems age. To better understand binary formation, we must
observe the youngest protobinary systems during binary formation. For young
binary systems with separations $<$500 au, we expect the formation pathways
are likely a combination of disk fragmentation by gravitational instability
and turbulent fragmentation with migration, whereas in systems with
separations $>$1000 au the binary formation mechanism is likely dominated by
turbulent fragmentation (e.g., Adams et al., 1989; Lee et al., 2019; Padoan &
Nordlund, 2002; Moe & Kratter, 2018).

IRAS 16293-2422 is a well studied Class 0 protostellar triple system located
in the $\rho$-Ophiuchi star forming region inside the dark cloud L1689N
(Lynds, 1962) at a distance of 141 pc (Dzib et al., 2018). The IRAS 16293
system consists of two close binaries in the southeast (separated by 54 au),
sources A1 and A2, (Wootten, 1989; Maureira et al., 2020) and a single
protostar in the northwest, source B, (separated by 725 au) (Wootten, 1989;
Looney et al., 2000). Sources A and B are embedded in a large (6-8 $\times$
103 au) envelope (e.g., Schöier et al., 2002; Crimier et al., 2010) with a
bridge of material connecting them (e.g., Looney et al., 2000; Pineda et al.,
2012). Although there had been speculation that source B was a more evolved T
Tauri star (Stark et al., 2004), observations of an inverse P-Cygni infall
profile toward source B (Pineda et al., 2012) indicate that the differences in
the sources are more due to viewing angle, as source B is likely face-on (e.g.
Rodríguez et al., 2005) whereas the binary sources A1 and A2 are more edge-on
(Pineda et al., 2012).

With three protostars and their geometry, the system has both complicated
envelope and outflow structures. The envelope has complex chemistry (e.g.,
Jørgensen et al., 2016), and due to the multiplicity and the connecting
bridge, the envelope presents difficulties in modeling the source details
(Jacobsen et al., 2018). The system has large and impressive outflows that
even at their discovery were known to be multi-lobed and very complicated,
likely due to a multiple system (e.g., Fukui et al., 1986; Wootten & Loren,
1987; Walker et al., 1988; Mizuno et al., 1990). The multiple outflows, some
of which are observed only on large scales while others are observed only on
small scales, are well summarized in van der Wiel et al. (2019) and many
references within.

On the largest scales, there are two observed outflows: east (blue)-west (red)
(Fukui et al., 1986; Wootten & Loren, 1987; Walker et al., 1988; Mizuno et
al., 1990; Stark et al., 2004) and northeast (red)-southwest (blue) (Walker et
al., 1988; Mizuno et al., 1990; Stark et al., 2004). However at smaller
scales, the northeast-southwest outflow is not detected, implying that the
launching engine was quenched (van der Wiel et al., 2019). On the other hand,
the east-west outflow, which is driven by source A, is clearly detected on
smaller scales with interferometers (Girart et al., 2014; Yeh et al., 2008;
van der Wiel et al., 2019). An additional outflow detected only on the small
scale is aligned northwest (blue)-southeast (red) originating slightly north
of source A with a blue bowshock in the location of source B (Girart et al.,
2014; Kristensen et al., 2013). Although there are some observational hints of
an outflow from source B (e.g., Yeh et al., 2008), there are no observations
that show clear evidence of outflow emission from source B.

In this paper, we present far-infrared observations (89 $\mu$m) on the scale
of $\sim$1000 au toward IRAS 16293-2422. The observations used the High-
Resolution Airborne Wide-band Camera (HAWC+: Vaillancourt et al., 2007; Harper
et al., 2018) onboard the Stratospheric Observatory For Infrared Astronomy
(SOFIA). We use the 89 $\mu$m dust continuum observations to infer the
magnetic field in the envelope of IRAS 16293-2422, hereafter called IRAS
16293. The inferred magnetic field morphology is compared to 850 $\mu$m
observation from the POL-2 polarimeter on the James Clerk Maxwell Telescope
(JCMT) and to 1.3 mm observations from the Atacama Large
Millimeter/submillimeter Array (ALMA). In addition, we compare the peak flux
location variation with wavelength using archival observations from WISE,
Herschel, and SCUBA 2 on the JCMT.

The paper is organized as follows: in Section 2 we present the observations
and data reduction, Section 3 covers the results of the observations, Section
4 discusses the results, and Section 5 contains our conclusions.

##  2 Observations and Data Reduction

We observed IRAS 16293 on 2019 July 23rd, between 12:21:13 and 12:25:41 UTC,
using HAWC+ at Band C (89 $\mu$m). The observations were part of the Cycle 7
program PID 07_0147 (PI: Novak). The bandwidth was 16.9 $\mu$m, the angular
resolution (full width at half maximum; FWHM) was 7.8′′, and the field of view
for the total intensity and polarimetry was 4.2′ $\times$ 2.7′ and 2.1′
$\times$ 2.7′, respectively.

The HAWC+ instrument is a far-infrared detector with 5 broadband filters for
continuum, a rotating half-wave plate that modulates the incoming
polarization, and a grid that orthogonally separates the polarization.
Multiple cycles of standard nodding and chopping were used to build an image
within a few minutes (Hildebrand et al., 2000). The instrument typically nods
at 0.1 Hz and chops at 10 Hz. The four dithered positional observations are
used to construct a single imaging block, or dither set. Due to a lack of
bright visual guide stars in the region, there was a small offset in some of
the imaging blocks. To correct this, we shifted the fit map Gaussian peak
center to align with the archival map from Herschel’s Photodetector Array
Camera & Spectrometer (PACS) at 100 $\mu$m. No individual observation was
offset by more than half a beam.

The raw images are then processed through the HAWC+ data reduction pipeline,
which produces science ready continuum and polarization products. The details
of the process are summarized in Harper et al. (2018). Briefly, the pipeline
flat-fields the demodulated and usable chopped data, calculates the Stokes I,
Q, and U parameters from the combined nodded fluxes, and then corrects the
pointing and instrumental polarization. Calibrations are applied to the flux
via standard atmospheric opacity models using Neptune planetary observations.
Finally, the Stokes maps are combined via standard re-gridding and Gaussian
smoothing (Houde & Vaillancourt, 2007). We smoothed the final maps to increase
the signal to noise in the lower surface brightness regions away from the
central source, creating a final beam size of 11.7′′.

We compared our Stokes I, Q, and U maps across the different dither sets using
a $\chi^{2}$ analysis. A $\chi^{2}$ is calculated for each pixel, which is
then compare with the pipeline error result for that pixel. If the $\chi^{2}$
is larger, then we inflate the pipeline error (also see Chapman et al., 2013;
Novak, 2011). Following Cox et al. (2022), we perform the inflation by fitting
a parameterized $\chi^{2}$ with Stokes I intensity, allowing us to inflate
each pixel based on its continuum brightness. This was done to handle the
uncertainties more carefully in the the brightest regions, where intensity-
dependent errors dominate.

Finally, the percent polarization ($p$) and its error ($\sigma_{p}$) are
calculated for each pixel following the HAWC+ handbook (Clarke et al., 2020)
by

$p=100\sqrt{\left(\frac{Q}{I}\right)^{2}+\left(\frac{U}{I}\right)^{2}}$ (1)

and

$\sigma_{p}=\frac{100}{I}\sqrt{\begin{aligned}
\frac{1}{(Q^{2}+U^{2})}(Q\sigma_{Q}^{2}+U\sigma_{U}^{2}+2QU\sigma_{QU})+\\\
\left[\left(\frac{Q}{I}\right)^{2}+\left(\frac{U}{I}\right)^{2}\right]\sigma_{I}^{2}-2\frac{Q}{I}\sigma_{QI}-2\frac{U}{I}\sigma_{UI}\end{aligned}}$
(2)

where $\sigma_{Q}$ and $\sigma_{U}$ are the uncertainties in Stokes Q and U
and where $\sigma_{QU}$, $\sigma_{QI}$, and $\sigma_{UI}$ are the covariance
uncertainty terms. The polarization fraction is debiased ($p^{\prime}$),
following the HAWC+ handbook (Clarke et al., 2020) by
$p^{\prime}=\sqrt{p^{2}-\sigma_{p}^{2}}$. From all of these, the debiased
polarization intensity ($P^{\prime}=I\times p^{\prime}/100$) and its errors
are calculated. Although the uncertainties are calculated per pixel, the
median Stokes I RMS is $\sigma_{I}$ = 2.2 mJy/arsec2 and the median
polarization intensity RMS $\sigma_{p}$ = 0.22 mJy/arsec2.

For the final maps, we selected the debiased polarization vectors that met the
criteria of $P^{\prime}$/$\sigma_{p}$ $>$ 3, the debiased polarization
percentage was $<$ 50%, and the Stokes I flux values were at least $\times$10
the Stokes I RMS noise level. While this was done at every pixel in the image,
we present only enough vectors to be considered Nyquist sampled, for a total
of 45 polarization vectors.

We estimate the overall flux calibration uncertainty of the observations at
20%, but for the remainder of the paper, any flux uncertainty listed will only
consider statistical uncertainty.

###  2.1 Archival Data

To compare our polarization observations with other wavelengths, we used JCMT
POL-2 850 $\mu$m wavelength polarization observations at 14$\arcsec$
resolution from Pattle et al. (2021), and ALMA Band 6 (1.3 mm wavelength)
polarization observations from (Sadavoy et al., 2018), smoothed to 1$\arcsec$
resolution.

In addition, we also used Herschel 70, 100, and 160 $\mu$m PACS data111PACS
observing labels 1342205093, 1342205094, 1342227150, 1342227151 and 250, 350,
and 500 $\mu$m SPIRE data222SPIRE observing labels 1342205093, 1342205094
(Spectral and Photometric Imaging Receiver).

We included the 24 $\mu$m Spitzer MIPS data333MIPS observing label 4321536
(Multiband Imaging Photometer). All of these data were obtained from the
Herschel Science Archive444http://archive.esac.esa.int/hsa/whsa/. Lastly, the
22 $\mu$m WISE data555WISE observing coadd id: 2477m243_ac51 (Wide-field
Infrared Survey Explorer) from the NASA/IPAC Infrared Science
Archive666https://irsa.ipac.caltech.edu0/applications/wise/ were also
included.

##  3 Results

Figure 1: HAWC+ 89 $\mu$m continuum map overlaid with the inferred magnetic
field direction in black. The three protostellar peaks at 3 mm are shown in
white Maureira et al. (2020), source A1/A2 in the South and source B in the
North. The vectors were selected to show Nyquist sampling. The image has a
smoothed image beam size of 11.7′′. The polarization scale bar is in the top
right corner, the beam is in the bottom left corner, and a spatial scalebar of
10′′ (1400 au) is provided in the bottom right corner.

Figure 1 shows the HAWC+ Band C (89 $\mu$m) dust continuum observations of
IRAS 16293 with the inferred magnetic field polarization direction (i.e.
polarization rotated by 90 degrees) plotted over the 89 $\mu$m continuum map
as line segments, hereafter called vectors.

We fit the protostellar envelope continuum emission using a Gaussian. The
protostellar envelope is resolved with a deconvolved fit size of 10.3′′ $\pm$
0.4′′ by 8.7′′ $\pm$ 0.4′′ with a position angle (PA) of 108∘ $\pm$ 11∘. The
integrated fit flux is 1236 $\pm$ 25 Jy, and the peak flux is 1.25 $\pm$ 0.02
Jy/arcsec2.

The continuum Gaussian fit of the HAWC+ 89 $\mu$m observation is different
than Gaussian fits to the Herschel archival data. The PACS 100 $\mu$m source
fits are somewhat smaller (7.7′′ $\pm$ 0.4′′ by 4.8′′ $\pm$ 0.4′′ and a PA of
141∘ $\pm$5∘), whereas the PACS 70 $\mu$m source fits are a little larger
(11.9′′ $\pm$ 0.3′′ by 10.1′′ $\pm$ 0.3′′ and a PA of 113∘ $\pm$ 7∘). The
variation indicated in the Gaussian fits of the HAWC+ and 70 and 100 PACS
$\mu$m images likely arises from differences in the source morphologies (e.g.,
outflows, envelopes, and the bridge connecting source A and B) so close to the
overall emission peak near 100 $\mu$m. On the other hand, the PA from the
HAWC+ Band C fit is within 10∘ of the 850 $\mu$m dust emission (Pattle et al.,
2021), although the 850 $\mu$m core extension is measured on much larger
scales.

Figure 2: Our continuum map overlaid with our normalized, Nyquist sampled
inferred magnetic field vectors in black. In orange vectors, we show the 850
$\mu$m POL-2 normalized-vectors from Pattle et al. (2021) at the published
pixel scale of 12′′. The corresponding beams, 11.7$\arcsec$ and 14$\arcsec$
for our observation and POL-2, respectively, are in the bottom left corner.

We can compare the inferred magnetic field from other observations with
varying wavelength and spatial scales. Figure 2 shows the inferred magnetic
field from our 89 $\mu$m data from Figure 1 (black vectors) with the 850
$\mu$m data from POL-2 (orange vectors), which have 14$\arcsec$ resolution and
12$\arcsec$ pixels (Pattle et al., 2021), both wavelengths presented as
normalized vectors. The vectors align on the western side of the source (and
one vector in the southeast), but otherwise there is generally not much
agreement on the inferred magnetic field direction.

On the other hand, misalignment is also seen when comparing the core scale
(POL-2) with the cloud scale (Planck). As pointed out by Pattle et al. (2021),
there is a misalignment of the large scale Planck field with the POL-2 850
$\mu$m magnetic field near the IRAS 16293 protostars, which is also seen in
Figure 2 as a field shift near the protostars compared to the outer region. In
Figure 1, the inferred HAWC+ magnetic field is generally aligned EW. This is
offset by $\sim$65∘ from the overall $\sim$arcminute large-scale field
inferred from Planck observations of the Ophiuchus L1689 molecular cloud core,
which is 24∘ E of N (Pattle et al., 2021).

At higher spatial scale, Figure 3 compares the 89 $\mu$m HAWC+ data and 850
$\mu$m POL-2 data with high resolution polarization observations at 1.3 mm
with ALMA from Sadavoy et al. (2018), which are consistent with 880$\mu$m
polarization observations from Rao et al. (2009). The ALMA data have been
smoothed to 1$\arcsec$ resolution to better compare with our observations.

Note that the region mapped in Figure 3 is only slightly larger than the POL-2
and HAWC+ beam sizes, so it is difficult to do more than compare broadly.
However, in general, the vectors from the three observations are not well-
aligned. The HAWC+ and ALMA observations have some agreement to the northwest,
and the ALMA and POL-2 observations have some agreement in the bridge region
between source A and B. On the other hand, HAWC+ and ALMA observations are
anti-aligned in the South, while ALMA and POL-2 observations are anti-aligned
in the middle-eastern side.

Figure 3:  The HAWC+ 89 $\mu$m continuum image overlaid with the normalized
magnetic field vectors from HAWC+ in black and from POL-2 in orange. Both maps
show oversampled vectors at 4$\arcsec$ pixels to better match the ALMA scales.
The white vectors show the ALMA 1.3 mm polarization vectors, smoothed to
1$\arcsec$ resolution (white beam in the bottom left corner). The ALMA vectors
are not normalized, and a 5% scale bar is given in the top right corner. The
dotted contour shows the ALMA continuum at the 10$\sigma_{I}$ level, where
$\sigma_{I}$ = 1.72 mJy arcsec-1. The locations of the hierarchical triple
protostars (source A1/A2 and source B) at 3 mm (e.g., Maureira et al., 2020)
are labeled. The bridge is the region connecting the two sources

##  4 Discussion

Figure 4: Image of the IRAS 16293-2422 system with Gaussian-fit peak positions
observed by various telescopes/instruments. Background image is our HAWC+
Stokes I map. From right to left: the magenta marker is the Gaussian peak of
the WISE W4 filter at 22 $\mu$m; the white marker is Spitzer MIPS 24 $\mu$m
peak; light coral, pale green, and blue corresponding to Herschel PACS at 70,
160, and 100 $\mu$m; orange and chartreuse for SPIRE at 250 and 350 $\mu$m;
the pink marker corresponding to SCUBA-2 at 850 $\mu$m; and the cyan is SPIRE
at 500 $\mu$m. Errors represent pointing uncertainties collected from the
respective telescope’s manual, main publications, or relevant sources (WISE
(Wright et al., 2010), MIPS on Spitzer (MIPS Instrument And MIPS Instrument
Support Teams, 2011), PACS & SPIRE on Herschel (Sánchez-Portal et al., 2014),
SCUBA2 on JCMT (Pattle et al., 2015; Kirk et al., 2018)). Contours are 1.3 mm
continuum from ALMA (smoothed to 1′′) with levels at N$\sigma$ for N$\in$[10,
30, 100, 300] and $\sigma$ = 1.72 mJy arcsec-1 (Sadavoy et al., 2018).

###  4.1 Continuum Wavelength Dependence

When comparing the HAWC+ dust continuum emission Gaussian fits with the
Gaussian fits from Herschel, as well as the peak in the POL-2 observations, we
noted that the peak of the emission varied with wavelength– with the far
infrared wavelengths peaking between the protostellar source locations. As
expected, the flux peak dust continuum of the IRAS 16293-2422 system is
strongly dependent on observing wavelength, as can be seen in the SED from
Schöier et al. (2002) with a maximum around 100 $\mu$m, but it also seems that
the spatial location of the peak is very strongly dependent on observing
wavelength. By using the archival observations, we can compare this trend more
broadly.

We fit the different wavelength emission (e.g., WISE, Herschel, and SCUBA2
listed in Section 2) to Gaussians. Figure 4 shows the Gaussian fit peak to
wavelengths from 22 $\mu$m to 850 $\mu$m. There is a clear gradient from
northwest to southeast from short wavelengths to long wavelengths. In other
words, we detect a gradient in the cloud dust temperatures with wavelength.
This is likely due to differences in the morphology or detailed properties of
the source, possibly arising from the structure of the envelope or disks in
the binary system.

Another explanation of the shift with wavelength could be a difference in
evolutionary stage of the binary sources. As stated in §1, source B has been
argued to be more evolved than source A (Stark et al., 2004). As we probe
shorter and shorter wavelengths, we are most sensitive to the hot dust in the
less obscured source. However, the evolutionary stages of the protostars are
still uncertain with most recent work suggesting that the sources are likely
at about the same evolutionary stage or that source B is slightly less evolved
(van der Wiel et al., 2019), but the evolutionary ages are still a point of
contention.

One of the complications is that the two protostars have different
inclinations with respect to the line of sight. Source A, which is the tight
binary, is closer to edge-on (inclinations of 59∘ and 74∘; Maureira et al.,
2020) while source B is close to face-on (inclination $\sim$18∘; Zamponi et
al., 2021). The source A protostars are most likely Class 0 based on evidence
of them driving active outflows (e.g., Yeh et al., 2008), but the evolution of
source B is harder to verify due to its face-on inclination and lack of
outflow activity (cf., Yeh et al., 2008). The difference in inclination could
explain the observed peak shift seen at the shorter wavelengths, since the
shorter wavelengths are most sensitive to the hotter dust seen in the more
face-on disk, which is less obscured. In addition, source A is also surrounded
by a $\sim$100 au circumbinary disk with a PA of $\sim$50∘ (Maureira et al.,
2020), which could increase the obscuration for the binary. We suggest that
such continuum peak positional shifts with wavelength may be present in other
binary sources, requiring careful consideration of variation in peak position
at the shorter wavelengths.

Figure 5: Same as Figure 1 but overlaid with the blue and red 12CO 3-2 outflow
contour lines of Yeh et al. (2008) using contours at 3, 10, 30, and
50$\sigma$, where $\sigma$ is 1.4 Jy beam-1 for the blue and 2.3 Jy beam-1 for
the red.

###  4.2 Magnetic Fields

At the large scales in dense, elongated structures or filaments of molecular
clouds, the magnetic field’s alignment depends on column density: below a
critical column density, the magnetic field is preferentially oriented
parallel to the elongation of density structures, and above, the magnetic
field is preferentially perpendicular (Planck Collaboration et al., 2016).
However, the details of the alignment may vary across star formation regions
(Stephens et al., 2022). Using HAWC+ observations in L1688, it has been shown
that the transition occurs at a molecular hydrogen column density of
$\sim$1021.7 cm-2 (Lee et al., 2021). Magneto-hydrodynamic simulations have
suggested that the transition may coincide with the kinetic energy of the
gravitationally induced flows surpassing the magnetic energy (Chen et al.,
2016).

One of the goals of our observations was to compare the inferred magnetic
field in the far infrared dust continuum to other resolutions and scales that
impact the star formation process. The large scale Planck inferred magnetic
field in the L1689-N region is 24∘ E of N (Pattle et al., 2021), i.e. nearly
perpendicular to the large L1689/L1712 filament identified in Ladjelate et al.
(2020). These large-scale Planck fields are generally consistent with the 850
$\mu$m JCMT POL-2 inferred magnetic fields shown in Figure 2, although in
detail the POL-2 fields are less aligned with the overall Planck field at the
cloud core (e.g., see Figure 4 in Pattle et al., 2021). This is not surprising
as the overall Planck field seems to exhibit a small shift at the the core
that suggests unresolved morphology.

When one traces the magnetic field from the L1689-N core scale down to the
protostellar envelope scale with the three observations in Figure 3, we see
regions where the inferred magnetic fields are aligned but other regions where
the inferred magnetic fields are misaligned. The HAWC+ and ALMA observations
are aligned in the northwest, while the ALMA and POL-2 observations are well
aligned in the bridge region between source A and B. On the other hand, HAWC+
and ALMA observations are misaligned in the South, while ALMA and POL-2
observations are misaligned on the eastern side. However, as pointed out by
Pattle et al. (2021), the average field angles are consistent between ALMA and
POL-2: 166${}^{\circ}\pm$31∘ for POL-2 at the cloud core and
176${}^{\circ}\pm$54∘ and 130${}^{\circ}\pm$14∘ averaged around the protostars
only and the bridge between them only in the unsmoothed ALMA data in Sadavoy
et al. (2018), respectively. Whereas the HAWC+ observations have a very
different measured polarization angle: 89${}^{\circ}\pm$23∘, which is more in
line with the E-W outflow.

Indeed, one of the most striking features of the HAWC+ map is the uniform
field that is nearly completely EW in direction, which is very different
morphologically to the Planck, JCMT, and ALMA inferred magnetic field angles.
Although one could argue for large field morphological connections or trends
between the mapped fields in Figure 3, there is still not overall coherence
across wavelengths. The key to understanding these different wavelength
observations is that the observations are probing different optical depth and
dust temperatures (as seen clearly in Figure 3), making it difficult to piece
together the overall 3D magnetic field morphology without a better
understanding of the system details.

The IRAS 16293-2422 triple system has a complicated morphology with the
protostars, disks, a bridge feature, and multiple outflows observed on
different spatial scales (see Figure 1 in van der Wiel et al., 2019). To date,
the outflows have only been associated with source A, which is likely due to
the face-on inclination of source B or possibly a difference in evolutionary
stage. As discussed in §1, source A has two outflows on the 1000s of au scale:
one east-west (e.g., Yeh et al., 2008) and one northwest-southeast (e.g.,
Kristensen et al., 2013). Arguably, the east-west outflow dominates the system
on the cloud core scale (e.g., Mizuno et al., 1990; Stark et al., 2004), which
is also well detected in higher resolution interferometric observations (Yeh
et al., 2008). Figure 5 shows our inferred magnetic fields with the blue- and
red-shifted outflows as identified from CO (3-2) observations of Yeh et al.
(2008). Our magnetic field vectors are well aligned with this outflow
direction, with the exception of our vectors extending more in the south.
Nonetheless, based on this comparison, we posit that the HAWC+ observations
are tracing the warm dust from the E-W outflow cavity or cavity walls of IRAS
16293-2422 (Stark et al., 2004). This is in contrast to the average magnetic
field from the POL-2 and ALMA observations, which are approximately
perpendicular to the E-W outflow. On the other hand, this is consistent with
the Huang et al. (2024) survey, where $\sim$40% of the sources have average
magnetic fields perpendicular to the source outflow.

There are many morphological examples of the outflows of protostars seeming to
modify or shape the magnetic field of the core as traced by single dish
observations: Orion A filament (Pattle et al., 2017), NGC 2071IR in Orion B
(Lyo et al., 2021), and CB 54 (Pattle et al., 2022) or to modify or shape the
inner protostellar envelope as traced by interferometric observations: Serpens
SMM1 (Hull et al., 2017), B335 (Maury et al., 2018), Emb 8(N) (Le Gouellec et
al., 2019), and BHR 71 IRS2 (Hull et al., 2020).

Additionally, the short wavelength observational inferred magnetic field can
exhibit different morphologies compared to the longer wavelengths. In Orion,
for example, the HAWC+ shorter wavelengths (53 and 89 $\mu$m compared to 154
and 214 $\mu$m) is more aligned with the bipolar outflow structure in the
Becklin–Neugebauer/Kleinman–Low region, as traced by molecular tracers (Chuss
et al., 2019). The longer wavelengths are argued to be tracing the cooler dust
that is outside of the explosion influenced region.

In this case for IRAS 16293-2422, the 89 $\mu$m continuum emission is tracing
the warmer dust surrounding the outflow. The magnetic field morphology in this
region is dominated by the outflow and not the magnetic field in the cloud
core or the bridge, which is seen to dominate at longer wavelengths. This is
somewhat consistent with the continuum peak of IRAS 16293-2442 at $\sim$100
$\mu$m being near the center of the system so that the polarization is
dominated by the large opening angle cavity seen in Figure 5.

###  4.3 Possibility of Polarization by Dichroic Extinction

As the inferred magnetic direction in Figure 1 is different than what is
observed at 850 $\mu$m with POL-2, we should also examine other polarization
mechanisms that may account for the observations. The most likely other
polarization mechanism would be dichroic extinction, used in the optical or
infrared to infer magnetic fields from polarization observations of stars
extincted by aligned dust grains (e.g., Scarrott & Warren-Smith, 1989).
Although dichroic extinction depends on dust properties, optical depth, and
temperature gradients (also see, Hildebrand et al., 2000), we can consider the
effect in the typical case for young stellar objects when unpolarized, or
weakly polarized, emission is extincted by cooler foreground dust that has its
grains aligned by magnetic fields or other means. In that case, the dust
grain’s long axis will have more efficient extinction than the short axis, so
the light is now polarized along the dust’s short axis, which would result in
a 90∘ flip in the polarization, compared to the inferred magnetic field from
emission of aligned dust grains. The effect has been seen in far infrared dust
polarization observations, particularly in the Sagittarius B2 molecular cloud
(Novak et al., 1997; Dowell, 1997). The process will be most effective in
optically thick sources with a temperature gradient.

Although, we can not rule out dichroic extinction as the dominate mechanism
for polarization in IRAS 16293-2422 (see Zielinski et al., 2021, where they
argue that dichroic extinction polarization is not important in HAWC+
observations of B335), there are a few reasons to disfavor that explanation.
First, we do not see a strong signature of the 90∘ flip toward the outside of
the source, as is seen, for example, in observations of NGC 1333 IRAS4 (Ko et
al., 2020) and OMC-3/MMS 6 (Liu, 2021). In fact, with the exception of two
vectors, Figure 1 shows a very uniform inferred magnetic field. Of course such
a flip is not necessary, but if we flip all of our vectors by 90∘, then the
vectors on the west in Figure 2 that currently agree with the POL-2 and ALMA
vectors would disagree. Second, although a flip of 90∘ for the inferred
magnetic field would mean an average field of 179${}^{\circ}\pm$23∘, which is
closer to the bulk average magnetic fields of 166${}^{\circ}\pm$31∘ for POL-2
and 176${}^{\circ}\pm$54∘ for ALMA data, in detail the vectors in the inner
region of the source, see Figure 3 for example, would not better match with a
HAWC+ flip of 90∘. Third, the polarization fraction is lower toward the center
of our observations in Figure 1. Such polarization depression with increasing
density are commonly detected in magnetic field aligned dust grain
polarization observations (e.g., Chuss et al., 2019) and are not seen in
regions where the polarization is due to dichroic extinction (Novak et al.,
1997; Dowell, 1997; Ko et al., 2020). Finally, east-west magnetic fields have
also been detected along a few stellar sight lines in the L1689 cloud at
optical wavelengths (e.g., Vrba et al., 1976), suggesting that there are east-
west fields in some parts of the region. Overall, based on these four reasons,
the IRAS 16293-2422 polarization observations presented here are most likely
due to magnetic field alignment of dust grains.

##  5 Conclusions

We present the 89 $\mu$m continuum polarization emission toward the
protostellar system IRAS 16293-2422 using the HAWC+ polarimeter onboard SOFIA.
Our main conclusions are:

  1. 1.

We detect a uniform magnetic field in the inner region of IRAS 16293-2422 that
is aligned east-west. The average field angle in the HAWC+ observations is
89${}^{\circ}\pm$23∘, which is consistent with one of the known large scale
outflows. This is different from the average field angles at longer
wavelengths: Planck large-scale average field angle of 24∘, JCMT POL-2 average
field angle of 166${}^{\circ}\pm$31∘, and the ALMA average field angle near
the protostars of 176${}^{\circ}\pm$54∘.

  2. 2.

We posit that the magnetic field probed by the 89 $\mu$m continuum emission is
dominated by the outflow magnetic field, while the 850 $\mu$m dust emission is
dominated more by the core magnetic field.

  3. 3.

The continuum peak of the source varies significantly with wavelength, moving
from near Source B at the midIR to near Source A in the FIR. This is either a
consequence of evolutionary state of the two sources or more likely due to the
difference in inclination. The shorter wavelengths are more sensitive to the
warmer dust which is seen in Source B.

These observations suggest that magnetic fields in young protostars are likely
more complicated than simple pictures suggest with multiple field morphologies
dominating at various scales and in different structures. Any observational
constraints should include multiwavelength observations that sample multiple
scales and all the structures of star formation.

Based on observations made with the NASA/DLR Stratospheric Observatory for
Infrared Astronomy (SOFIA). SOFIA is jointly operated by the Universities
Space Research Association, Inc. (USRA), under NASA contract NNA17BF53C, and
the Deutsches SOFIA Institut (DSI) under DLR contract 50 OK 0901 to the
University of Stuttgart. L.W.L. acknowledges support from NSF AST-1910364 and
NSF AST-2307844. R.E.H acknowledges support from NSF AST-1910364. G.N. is
grateful for NASA support for this work, via award Nos. SOF 06-0116 and SOF
07-0147 issued by USRA to Northwestern University. E.G.C. acknowledges support
from the National Science Foundation through the NSF MPS Ascend Fellowship
Grant number 2213275. K.P. is a Royal Society University Research Fellow,
supported by grant no. URF\R1\211322. We thank Naomi Hirano for providing us
with the calibrated SMA CO observations in Figure 5. We thank John Tobin for
discussions on binarity and Zhi-Yun Li for discussions on magnetic field
alignment with column density.

## References

  * Adams et al. (1989) Adams, F. C., Ruden, S. P., & Shu, F. H. 1989, ApJ, 347, 959, doi: 10.1086/168187
  * Chapman et al. (2013) Chapman, N. L., Davidson, J. A., Goldsmith, P. F., et al. 2013, ApJ, 770, 151, doi: 10.1088/0004-637X/770/2/151
  * Chen et al. (2016) Chen, C.-Y., King, P. K., & Li, Z.-Y. 2016, ApJ, 829, 84, doi: 10.3847/0004-637X/829/2/84
  * Chuss et al. (2019) Chuss, D. T., Andersson, B. G., Bally, J., et al. 2019, ApJ, 872, 187, doi: 10.3847/1538-4357/aafd37
  * Clarke et al. (2020) Clarke, M., Berthoud, M., Kovács, A., Santos, F., & Novak, G. 2020, Guest Observer Handbook for HAWC+ Data Products, Rev. F.  https://www.sofia.usra.edu/sites/default/files/Instruments/HAWC_PLUS/Documents/hawc_data_handbook.pdf
  * Cox et al. (2018) Cox, E. G., Harris, R. J., Looney, L. W., et al. 2018, ApJ, 855, 92, doi: 10.3847/1538-4357/aaacd2
  * Cox et al. (2022) Cox, E. G., Novak, G., Sadavoy, S. I., et al. 2022, ApJ, 932, 34, doi: 10.3847/1538-4357/ac722a
  * Crimier et al. (2010) Crimier, N., Ceccarelli, C., Maret, S., et al. 2010, A&A, 519, A65, doi: 10.1051/0004-6361/200913112
  * Dowell (1997) Dowell, C. D. 1997, ApJ, 487, 237, doi: 10.1086/304577
  * Dzib et al. (2018) Dzib, S. A., Ortiz-León, G. N., Hernández-Gómez, A., et al. 2018, A&A, 614, A20, doi: 10.1051/0004-6361/201732093
  * Encalada et al. (2021) Encalada, F. J., Looney, L. W., Tobin, J. J., et al. 2021, ApJ, 913, 149, doi: 10.3847/1538-4357/abf4fd
  * Fukui et al. (1986) Fukui, Y., Sugitani, K., Takaba, H., et al. 1986, ApJ, 311, L85, doi: 10.1086/184803
  * Galametz et al. (2018) Galametz, M., Maury, A., Girart, J. M., et al. 2018, A&A, 616, A139, doi: 10.1051/0004-6361/201833004
  * Girart et al. (2014) Girart, J. M., Estalella, R., Palau, A., Torrelles, J. M., & Rao, R. 2014, ApJ, 780, L11, doi: 10.1088/2041-8205/780/1/L11
  * Harper et al. (2018) Harper, D. A., Runyan, M. C., Dowell, C. D., et al. 2018, Journal of Astronomical Instrumentation, 07, 1840008, doi: 10.1142/s2251171718400081
  * Harris et al. (2018) Harris, R. J., Cox, E. G., Looney, L. W., et al. 2018, ApJ, 861, 91, doi: 10.3847/1538-4357/aac6ec
  * Hildebrand et al. (2000) Hildebrand, R. H., Davidson, J. A., Dotson, J. L., et al. 2000, PASP, 112, 1215, doi: 10.1086/316613
  * Houde & Vaillancourt (2007) Houde, M., & Vaillancourt, J. E. 2007, PASP, 119, 871, doi: 10.1086/521109
  * Huang et al. (2024) Huang, B., Girart, J. M., Stephens, I. W., et al. 2024, ApJ, 963, L31, doi: 10.3847/2041-8213/ad27d4
  * Hull et al. (2020) Hull, C. L. H., Le Gouellec, V. J. M., Girart, J. M., Tobin, J. J., & Bourke, T. L. 2020, ApJ, 892, 152, doi: 10.3847/1538-4357/ab5809
  * Hull et al. (2014) Hull, C. L. H., Plambeck, R. L., Kwon, W., et al. 2014, ApJS, 213, 13, doi: 10.1088/0067-0049/213/1/13
  * Hull et al. (2017) Hull, C. L. H., Girart, J. M., Tychoniec, Ł., et al. 2017, ApJ, 847, 92, doi: 10.3847/1538-4357/aa7fe9
  * Jacobsen et al. (2018) Jacobsen, S. K., Jørgensen, J. K., van der Wiel, M. H. D., et al. 2018, A&A, 612, A72, doi: 10.1051/0004-6361/201731668
  * Jørgensen et al. (2016) Jørgensen, J. K., van der Wiel, M. H. D., Coutens, A., et al. 2016, A&A, 595, A117, doi: 10.1051/0004-6361/201628648
  * Kataoka et al. (2015) Kataoka, A., Muto, T., Momose, M., et al. 2015, ApJ, 809, 78, doi: 10.1088/0004-637X/809/1/78
  * Kirk et al. (2018) Kirk, H., Hatchell, J., Johnstone, D., et al. 2018, ApJS, 238, 8, doi: 10.3847/1538-4365/aada7f
  * Ko et al. (2020) Ko, C.-L., Liu, H. B., Lai, S.-P., et al. 2020, ApJ, 889, 172, doi: 10.3847/1538-4357/ab5e79
  * Kristensen et al. (2013) Kristensen, L. E., Klaassen, P. D., Mottram, J. C., Schmalzl, M., & Hogerheijde, M. R. 2013, A&A, 549, L6, doi: 10.1051/0004-6361/201220668
  * Ladjelate et al. (2020) Ladjelate, B., André, P., Könyves, V., et al. 2020, A&A, 638, A74, doi: 10.1051/0004-6361/201936442
  * Lazarian (2007) Lazarian, A. 2007, J. Quant. Spec. Radiat. Transf., 106, 225, doi: 10.1016/j.jqsrt.2007.01.038
  * Le Gouellec et al. (2019) Le Gouellec, V. J. M., Hull, C. L. H., Maury, A. J., et al. 2019, ApJ, 885, 106, doi: 10.3847/1538-4357/ab43c2
  * Lee et al. (2019) Lee, A. T., Offner, S. S. R., Kratter, K. M., Smullen, R. A., & Li, P. S. 2019, ApJ, 887, 232, doi: 10.3847/1538-4357/ab584b
  * Lee et al. (2021) Lee, D., Berthoud, M., Chen, C.-Y., et al. 2021, ApJ, 918, 39, doi: 10.3847/1538-4357/ac0cf2
  * Liu (2021) Liu, H. B. 2021, ApJ, 914, 25, doi: 10.3847/1538-4357/abf8b6
  * Looney et al. (2000) Looney, L. W., Mundy, L. G., & Welch, W. J. 2000, ApJ, 529, 477, doi: 10.1086/308239
  * Lynds (1962) Lynds, B. T. 1962, ApJS, 7, 1, doi: 10.1086/190072
  * Lyo et al. (2021) Lyo, A. R., Kim, J., Sadavoy, S., et al. 2021, ApJ, 918, 85, doi: 10.3847/1538-4357/ac0ce9
  * Maureira et al. (2020) Maureira, M. J., Pineda, J. E., Segura-Cox, D. M., et al. 2020, ApJ, 897, 59, doi: 10.3847/1538-4357/ab960b
  * Maury et al. (2018) Maury, A. J., Girart, J. M., Zhang, Q., et al. 2018, MNRAS, 477, 2760, doi: 10.1093/mnras/sty574
  * MIPS Instrument And MIPS Instrument Support Teams (2011) MIPS Instrument And MIPS Instrument Support Teams. 2011, 100, doi: 10.26131/IRSA488
  * Mizuno et al. (1990) Mizuno, A., Fukui, Y., Iwata, T., Nozawa, S., & Takano, T. 1990, ApJ, 356, 184, doi: 10.1086/168829
  * Moe & Di Stefano (2017) Moe, M., & Di Stefano, R. 2017, ApJS, 230, 15, doi: 10.3847/1538-4365/aa6fb6
  * Moe & Kratter (2018) Moe, M., & Kratter, K. M. 2018, ApJ, 854, 44, doi: 10.3847/1538-4357/aaa6d2
  * Novak (2011) Novak, G. 2011, in Astronomical Society of the Pacific Conference Series, Vol. 449, Astronomical Polarimetry 2008: Science from Small to Large Telescopes, ed. P. Bastien, N. Manset, D. P. Clemens, & N. St-Louis, 50
  * Novak et al. (1997) Novak, G., Dotson, J. L., Dowell, C. D., et al. 1997, ApJ, 487, 320, doi: 10.1086/304576
  * Padoan & Nordlund (2002) Padoan, P., & Nordlund, Å. 2002, ApJ, 576, 870, doi: 10.1086/341790
  * Pattle et al. (2023) Pattle, K., Fissel, L., Tahani, M., Liu, T., & Ntormousi, E. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 193, doi: 10.48550/arXiv.2203.11179
  * Pattle et al. (2015) Pattle, K., Ward-Thompson, D., Kirk, J. M., et al. 2015, MNRAS, 450, 1094, doi: 10.1093/mnras/stv376
  * Pattle et al. (2017) Pattle, K., Ward-Thompson, D., Berry, D., et al. 2017, ApJ, 846, 122, doi: 10.3847/1538-4357/aa80e5
  * Pattle et al. (2021) Pattle, K., Lai, S.-P., Di Francesco, J., et al. 2021, ApJ, 907, 88, doi: 10.3847/1538-4357/abcc6c
  * Pattle et al. (2022) Pattle, K., Lai, S.-P., Sadavoy, S., et al. 2022, MNRAS, 515, 1026, doi: 10.1093/mnras/stac1356
  * Pineda et al. (2012) Pineda, J. E., Maury, A. J., Fuller, G. A., et al. 2012, A&A, 544, L7, doi: 10.1051/0004-6361/201219589
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 586, A138, doi: 10.1051/0004-6361/201525896
  * Rao et al. (2009) Rao, R., Girart, J. M., Marrone, D. P., Lai, S.-P., & Schnee, S. 2009, ApJ, 707, 921, doi: 10.1088/0004-637X/707/2/921
  * Rodríguez et al. (2005) Rodríguez, L. F., Loinard, L., D’Alessio, P., Wilner, D. J., & Ho, P. T. P. 2005, ApJ, 621, L133, doi: 10.1086/429223
  * Sadavoy et al. (2018) Sadavoy, S. I., Myers, P. C., Stephens, I. W., et al. 2018, ApJ, 869, 115, doi: 10.3847/1538-4357/aaef81
  * Sadavoy et al. (2019) Sadavoy, S. I., Stephens, I. W., Myers, P. C., et al. 2019, ApJS, 245, 2, doi: 10.3847/1538-4365/ab4257
  * Sánchez-Portal et al. (2014) Sánchez-Portal, M., , Marston, A., et al. 2014, Experimental Astronomy, 37, 453, doi: 10.1007/s10686-014-9396-z
  * Scarrott & Warren-Smith (1989) Scarrott, S. M., & Warren-Smith, R. F. 1989, MNRAS, 237, 995, doi: 10.1093/mnras/237.4.995
  * Schöier et al. (2002) Schöier, F. L., Jørgensen, J. K., van Dishoeck, E. F., & Blake, G. A. 2002, A&A, 390, 1001, doi: 10.1051/0004-6361:20020756
  * Stark et al. (2004) Stark, R., Sandell, G., Beck, S. C., et al. 2004, ApJ, 608, 341, doi: 10.1086/392492
  * Stephens et al. (2022) Stephens, I. W., Myers, P. C., Zucker, C., et al. 2022, ApJ, 926, L6, doi: 10.3847/2041-8213/ac4d8f
  * Tobin et al. (2016) Tobin, J. J., Looney, L. W., Li, Z.-Y., et al. 2016, ApJ, 818, 73, doi: 10.3847/0004-637X/818/1/73
  * Tobin et al. (2022) Tobin, J. J., Offner, S. S. R., Kratter, K. M., et al. 2022, ApJ, 925, 39, doi: 10.3847/1538-4357/ac36d2
  * Tsukamoto et al. (2023) Tsukamoto, Y., Maury, A., Commercon, B., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Astronomical Society of the Pacific Conference Series, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 317
  * Vaillancourt et al. (2007) Vaillancourt, J. E., Chuss, D. T., Crutcher, R. M., et al. 2007, in Infrared Spaceborne Remote Sensing and Instrumentation XV, ed. M. Strojnik-Scholl (SPIE), doi: 10.1117/12.730922
  * van der Wiel et al. (2019) van der Wiel, M. H. D., Jacobsen, S. K., Jørgensen, J. K., et al. 2019, A&A, 626, A93, doi: 10.1051/0004-6361/201833695
  * Vrba et al. (1976) Vrba, F. J., Strom, S. E., & Strom, K. M. 1976, AJ, 81, 958, doi: 10.1086/111976
  * Walker et al. (1988) Walker, C. K., Lada, C. J., Young, E. T., & Margulis, M. 1988, ApJ, 332, 335, doi: 10.1086/166659
  * Wootten (1989) Wootten, A. 1989, ApJ, 337, 858, doi: 10.1086/167156
  * Wootten & Loren (1987) Wootten, A., & Loren, R. B. 1987, ApJ, 317, 220, doi: 10.1086/165270
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, The Astronomical Journal, 140, 1868, doi: 10.1088/0004-6256/140/6/1868
  * Yang et al. (2016) Yang, H., Li, Z.-Y., Looney, L., & Stephens, I. 2016, MNRAS, 456, 2794, doi: 10.1093/mnras/stv2633
  * Yeh et al. (2008) Yeh, S. C. C., Hirano, N., Bourke, T. L., et al. 2008, ApJ, 675, 454, doi: 10.1086/524648
  * Zamponi et al. (2021) Zamponi, J., Maureira, M. J., Zhao, B., et al. 2021, MNRAS, 508, 2583, doi: 10.1093/mnras/stab2657
  * Zhang et al. (2014) Zhang, Q., Qiu, K., Girart, J. M., et al. 2014, ApJ, 792, 116, doi: 10.1088/0004-637X/792/2/116
  * Zielinski et al. (2021) Zielinski, N., Wolf, S., & Brunngräber, R. 2021, A&A, 645, A125, doi: 10.1051/0004-6361/202039126


# CAPER: Coarsen, Align, Project, Refine
A General Multilevel Framework for Network Alignment

Jing Zhu  University of Michigan, Ann Arbor<EMAIL_ADDRESS>,  Danai Koutra
University of Michigan, Ann Arbor<EMAIL_ADDRESS>and  Mark Heimann
Lawrence Livermore Natl Laboratory<EMAIL_ADDRESS>
(2022)

###### Abstract.

Network alignment, or the task of finding corresponding nodes in different
networks, is an important problem formulation in many application domains. We
propose CAPER, a multilevel alignment framework that Coarsens the input
graphs, Aligns the coarsened graphs, Projects the alignment solution to finer
levels and Refines the alignment solution. We show that CAPER can improve upon
many different existing network alignment algorithms by _enforcing alignment
consistency across multiple graph resolutions: nodes matched at finer levels
should also be matched at coarser levels._ CAPER also accelerates the use of
slower network alignment methods, at the modest cost of linear-time coarsening
and refinement steps, by allowing them to be run on smaller coarsened versions
of the input graphs. Experiments show that CAPER can improve upon diverse
network alignment methods by an average of 33$\%$ in accuracy and/or an order
of magnitude faster in runtime.

††copyright: none††journalyear: 2022††copyright: usgovmixed††conference:
Proceedings of the 31st ACM International Conference on Information and
Knowledge Management; October 17–21, 2022; Atlanta, GA, USA††booktitle:
Proceedings of the 31st ACM International Conference on Information and
Knowledge Management (CIKM ’22), October 17–21, 2022, Atlanta, GA, USA††price:
15.00††doi: 10.1145/3511808.3557563††isbn: 978-1-4503-9236-5/22/10

##  1\. Introduction

Graphs or networks are foundational representations for relational structure
and their analysis is useful in innumerable scientific and industrial
applications. In many diverse tasks, such as recommendation across multiple
social networks, protein-protein interaction analysis, and database schema
matching (kazemi2016network, ), it is necessary to discover meaningful
correspondences between nodes in multiple networks. This general problem is
called network alignment.

Network alignment methods in general have two main limitations. First, they
may overfit to local structural similarity and fail to preserve higher-order
measures of matching consistency (conealign, ; refina, ). Second, especially
the most accurate methods tend to rely on solving challenging optimization
problems with high computational complexity, e.g. quadratic or cubic time in
the number of nodes in one of the input graphs (gwl, ; zhang2019kergm, ;
conealign, ).

We argue that multilevel network analysis is a powerful technique for
improving network alignment algorithms on both fronts. Accordingly, we design
the first general multilevel framework to pair with any network alignment
method, a four-step framework which we call CAPER: (1) Coarsening a graph into
multiple levels of varying coarseness, (2) Aligning at the coarsest level, and
(3) Projecting back to finer levels, and (4) Refining the solution at each
level. We can accelerate the use of slow network alignment algorithms by
running them on the smaller coarsened graphs, while refining the solutions at
multiple levels of structural resolution encourages greater consistency
between the local and global structure of matched nodes. Our contributions can
be summarized as follows:

  * •

General-Purpose Framework: We propose an intuitive multilevel framework
(CAPER) in which any network alignment method can be used.

  * •

Design Choices and Empirical Success: We propose and study specific design
choices and parameter settings that work well within CAPER. We provide code
and additional supplementary material at https://github.com/GemsLab/CAPER.

  * •

Study of Accuracy and Runtime Tradeoff: Through complexity analysis and
experiments, we show that CAPER is able to improve accuracy by 33$\%$ on
average across multiple datasets and/or is 10x faster runtime than baselines,
depending on the properties of the base methods employed.

Figure 1. An illustrative example of CAPER. The pink and blue nodes in the
leftmost figure have the same local structural similarity, so methods that
overfit the local structural similarity may misalign them. But with the help
of higher-order information (coarsened graphs in step 4), CAPER is able to
eventually correctly align the pink nodes as well as the blue nodes.

##  2\. Related Work

Graph Coarsening and Multilevel Methods. Graph coarsening (liu2018graph, ) is
the process of shrinking a large graph into a similar smaller one, such that
some properties or structures are preserved, e.g. spectral graph properties or
cliques. It has been used to accelerate many graph mining tasks, including
graph clustering (dhillon2007weighted, ), node embedding (mile, ; graphzoom, )
and graph neural networks (groupinn, ).

Network Alignment. We focus on unsupervised approaches requiring no known
matchings _a priori_. These can be categorized into two groups. (1) Classic
graph alignment approaches often formulate an optimization-based assignment
problem. FINAL (final, ) optimizes a topological consistency objective which
may be augmented with node and edge attribute information. MAGNA (magna, ),
applied to biological networks, uses genetic algorithms to evolve network
populations while maximizing proximity consistency criteria. More recently,
Zhang _et al._ (zhang2019kergm, ) leveraged kernel methods to solve the
quadratic assignment problem, but requires cubic computational complexity. (2)
Another line of work relies on embedding-based methods. REGAL (regal, )
matches structural node embeddings (jin2021toward, ; rossi2020proximity, )
that are directly comparable across networks. CONE-Align (conealign, ) uses
embeddings modeling proximity within each graph (rossi2020proximity, ) and
aligns the graphs’ embedding spaces with a subspace alignment procedure, while
GWL (gwl, ) solves a Gromov-Wasserstein optimization problem to jointly find
node embeddings and the graph matching. G-CREWE (gcrewe, ) uses graph
compression to accelerate the matching step of embedding-based alignment,
though the embedding step is performed on the entire input graphs.

Table 1. Comparing alignment meta-frameworks.

|  General |  Multiscale |  Improves |  Improves
---|---|---|---|---
|  |  |  accuracy |  runtime
MOANA (moana, ) | ✗ | ✓ | ✗ | ✓
Boosting (kyster2021boosting, ) | ✗ | ✗ | ✓ | ✗
RefiNA (refina, ) | ✓ | ✗ | ✓ | ✗
CAPER | ✓ | ✓ | ✓ | ✓

Network Alignment Meta-Frameworks. A few recent works have proposed meta-
frameworks to improve unsupervised network alignment algorithms. This includes
MOANA, the only other multilevel network alignment approach (moana, ). MOANA
uses multiresolution matrix factorization to accelerate FINAL (final, ) (it
produces negative-valued adjacency matrices that do not work with other
network alignment methods) at the cost of some accuracy. RefiNA (refina, )
makes the opposite tradeoff, enforcing greater local consistency to increase
accuracy of several base methods at the cost of adding additional runtime.
Another meta-framework (kyster2021boosting, ) studies how design choices of
recent embedding-based network alignment methods can be combined to increase
accuracy via boosting. Meanwhile, our approach inherits all these benefits, as
shown in Tab. 1.

##  3\. Preliminaries

Graphs. We consider two graphs $G_{1}$ and $G_{2}$ with nodesets
$\mathcal{V}_{1},\mathcal{V}_{2}$ and adjacency matrices
$\mathbf{A}_{1},\mathbf{A}_{2}$ containing edges between nodes. A graph
$G_{i}$ has a coarsened version $\tilde{G}_{i}$ with a smaller nodeset
$\tilde{\mathcal{V}}_{i}$ of $\tilde{n}<n$ nodes. Each node in the original
graph corresponds to a node in the coarsened graph, represented by an
assignment matrix $\mathbf{P}\in\\{0,1\\}^{n\times\tilde{n}}$. For clarity, we
drop the $\tilde{G}$ notation unless it is necessary to distinguish coarsened
and uncoarsened versions.

Alignment. An alignment between the nodes of two graphs can be represented by
a matrix $\mathbf{S}$, where $s_{ij}$ is the (real-valued or binary)
similarity between node $i$ in $G_{1}$ and node $j$ in $G_{2}$.

Problem Statement. Given two graphs $G_{1}$ and $G_{2}$ with meaningful node
alignments, but none known a priori, we seek to shrink them into coarsened
versions $\tilde{G}_{1}$ and $\tilde{G}_{2}$, and recover their alignment
$\mathbf{S}$ from $\tilde{\mathbf{S}}$ obtained by aligning their coarser
versions $\tilde{G}_{1}$ and $\tilde{G}_{2}$.

##  4\. Method

Next, we detail our CAPER framework, the first general-purpose multilevel
framework for unsupervised network alignment that can accommodate any base
network alignment approach. It consists of four steps that are carefully
designed in order to achieve higher accuracy and/or lower runtime compared to
its base alignment methods: Coarsen, Align, Project, Refine (CAPER). In Fig.
1, we provide an example of how CAPER can implicitly enforce higher-order
structural consistency that improves network alignment.

###  4.1. Graph Coarsening

Given an input graph $G_{i}$, we want to obtain a coarsened graph
$\tilde{G}_{i}$ using grouping-based coarsening methods. We leverage the
normalized heavy-edge matching (NHEM) heuristic (dhillon2007weighted, ) for
graph coarsening. This approach repeatedly combines pairs of adjacent nodes
into a supernode in decreasing order of degree-normalized edge weight (mile,
), which for edge ($u$,$v$) with weight $w_{uv}$ connecting nodes $u$ and $v$
with degrees $d_{u}$ and $d_{v}$ respectively is given by
$w_{uv}/\sqrt{d_{u}d_{v}}$, until no node is left uncombined or the uncombined
nodes do not have uncombined neighbors (isolated nodes). The resulting coarse
graph consists of these supernodes, which share an edge if any of the nodes in
one supernode shared an edge in the original graph with any of the nodes in
the other supernode.

Graph coarsening turns each input graph $G_{i}$ into a coarsened graph
$\tilde{G}_{i}$. We iteratively repeat this coarsening procedure up to
$\mathbf{L}$ times to produce a sequence of coarsened graphs
$\tilde{G}_{i}^{(0)},\ldots,\tilde{G}_{i}^{(L)}$, where the first level is the
input graph ($\tilde{G}_{i}^{(0)}=G_{i}$), and the coarsest (smallest) graph
is $\tilde{G}_{i}^{(L)}$. Assignments between nodes at consecutive levels
$\ell-1$ and $\ell$ are contained in a matrix $\mathbf{P}_{i}^{(\ell)}$ for
$\ell\in[1,\ldots,L]$.

###  4.2. Alignment of Coarsened Graphs

We can apply any unsupervised network alignment method to align the nodes of
the coarsest graphs $\tilde{G}_{1}^{(L)}$ and $\tilde{G}_{2}^{(L)}$ to produce
a matching $\mathbf{S}^{(\mathbf{L})}$. We observe that the coarsening
procedure sometimes generates slightly different numbers of nodes for the same
graph even if the input graphs have the same size, so the proposed formulation
must be able to handle graphs of different sizes. This can be done by adding
singleton nodes to the smaller graph (zhang2019kergm, ; conealign, ).

###  4.3. Projection

We project the alignment solution at the coarsest level $\mathbf{S}^{(\ell)}$
to a mapping between the nodes at the next finer level using the assignment
matrices:
$\mathbf{S}^{(\ell-1)}=\mathbf{P}_{1}^{(\ell)^{\top}}\mathbf{S}^{(\ell)}\mathbf{P}_{2}^{(\ell)}$.
Note that this solution is coarse, and all nodes in level $\ell-1$ mapped to
the same supernode in level $\ell$ will have the same match. Thus, we next
show how to use the finer graph structure to refine this coarse solution.

###  4.4. Soft Refinement

Recent work for refining network alignment (refina, ) operates on “hard”
initial solutions, where each node is mapped to at most one other node. Here,
we propose a new refinement operator that uses the “soft” initial alignments,
which better models the various strengths of several potential matches for
each node, as shown in Fig. 5. Given an initial soft (real-valued) alignment
$\mathbf{S}$, we iteratively apply the update rule
$\mathbf{S}=\text{NORMALIZE}(\mathbf{S}\circ\mathbf{A}_{1}\mathbf{S}\mathbf{A}_{2}+\epsilon)$,
where $\circ$ denotes Hadamard product, $\epsilon$ is a small positive minimum
matching score to any pair of nodes to prevent over-reliance on the initially
discovered matches (we set
$\epsilon=10^{-\lceil\log_{10}\max(n_{1},n_{2})\rceil}$) and NORMALIZE is a
single round of row-wise then column-wise normalization, as in (refina, ).

Table 2. Dataset statistics: These four datasets represent various phenomena
as shown in the description column.

Name | Nodes | Edges | Description
---|---|---|---
Arenas (koblenz, ) | 1,133 | 5,451 | communication network
Hamsterster (koblenz, ) | 2,426 | 16,613 | social network
Facebook (snapnets, ) | 4,039 | 88,234 | social network
Magna (magna, ) | 1,004 | 8,323 | protein-protein interaction

We iteratively apply this project-and-refine procedure between successive
levels until we arrive back at the input level, giving us the mapping between
nodes in the original graph.

###  4.5. Computational Complexity

We analyze the time complexity of CAPER as a function of the number of nodes
$n$ (to simplify notation we assume this is the same for both graphs), for
sparse graphs with $O(n)$ edges. Then the complexity of our CAPER framework is
$Lf_{\text{coarsen}}(n)+f_{\text{align}}(\frac{n}{2^{L}})+L\Big{(}f_{\text{project}}(n)+f_{\text{refine}}(n)$).
The coarsening time applied to each of $L$ levels, $f_{\text{coarsen}}(n)$, is
linear in the number of edges using heavy-edge matching (dhillon2007weighted,
), which is $O(n)$. Projection $f_{\text{project}}$ and refinement
$f_{\text{refine}}$ consist of matrix multiplications that, by maintaining a
sparse matching matrix, can also run in $O(n)$ time (refina, ).

Meanwhile, with NHEM shrinking the graph by approximately a factor of 2 at
each level (dhillon2007weighted, ), note that we are able to run the base
alignment step on a smaller graph, incurring a runtime of
$f_{\text{align}}(\frac{n}{2^{L}})$ as opposed to $f_{\text{align}}(n)$ by
applying the base alignment algorithm to the full input graphs. Thus, CAPER
can offer computational speedup particularly for slow base alignment methods,
where $f_{\text{align}}$ may be asymptotically large (such as $O(n^{3})$),and
the savings may outweigh the overhead of coarsening, projection, and
refinement.

##  5\. Experiments

We first describe our experimental setup and the datasets and baseline methods
used in our empirical analysis, and then show quantitative improvements from
CAPER and a closer ablation study.

Data. We use simulated and real alignment scenarios on graphs representing
various real-world phenomena (Tab. 2). Following prior works (regal, ; final,
; kyster2021boosting, ), we simulate a network alignment scenario with known
ground truth: a graph with adjacency matrix $\mathbf{A}$ is aligned to a noisy
permuted copy
$\mathbf{A}^{*}=\overline{\mathbf{S}}\mathbf{A}\overline{\mathbf{S}}^{\top}$
and $\overline{\mathbf{S}}$, for which we generate a random permutation matrix
$\overline{\mathbf{S}}$; we then randomly remove edges from $\mathbf{A}^{*}$
with probability $p\in[0.05,0.10,0.15,0.20,0.25]$. The MAGNA (magna, )
networks are protein-protein interaction (PPI) networks that are aligned to
versions of themselves with various percentages of low-confidence PPIs (edges)
added; thus, all edges in this graph represent real-world phenomena and we do
not need to synthesize an alignment scenario.

Baselines. We use (1) FINAL (final, ) and (2) REGAL (regal, ), which are
popular unsupervised network alignment methods that have usable public
codebases and represent different classes of techniques (optimization and node
embeddings), demonstrating the wide applicability of our framework. We also
use a more recent approach, (3) GWL (gwl, ), which combines optimization and
node embeddings, and achieves good accuracy but has slow runtime due to its
$O(n^{3})$ computational complexity. Moreover, we consider the post hoc
refinement method RefiNA applied to each of the network alignment methods: (4)
FINAL-RefiNA, (5) REGAL-RefiNA and (6) GWL-RefiNA. Additionally, we use (7)
MOANA (moana, ) as a baseline, the only other multilevel network alignment
method.

(a)

(a) Arenas

(b) Hamsterster

(c) Magna

(d) Facebook

Figure 2. Accuracy (solid lines) vs. different noise levels. CAPER outperforms
baselines, particularly as noise increases. Figure 3. Accuracy vs. runtime for
CAPER and RefiNA for 20% noise. CAPER yields better accuracy for FINAL and
REGAL by enforcing higher-order consistency. For GWL, CAPER runs up to 80x
faster because the alignment is run on smaller graphs.

For FINAL’s prior alignment information, we take the top
$k=\lfloor\log_{2}n\rfloor$ most similar nodes by degree for each node (regal,
; conealign, ). We set other parameters for REGAL (regal, ) and GWL (gwl, )
using the defaults recommended by the authors.

CAPER variants. We test variants of CAPER using each base alignment method:
CAPER(FINAL), CAPER(REGAL), and CAPER(GWL). We use 3 coarsening levels and 100
refinement iterations, as in (refina, ), to balance accuracy and computational
efficiency (we found that more refinement may increase performance if that is
desired and increased runtime is acceptable).

Evaluation. We measure alignment accuracy, or the proportion of correctly
aligned nodes, and runtime.

Table 3. Number of nodes in the coarsened graph after 2–4 levels of coarsening. Name |  2 |  3 | 4
---|---|---|---
Hamsterster | 1,288 | 702 | 418
Facebook | 2,078 | 1,078 | 572

###  5.1. Alignment Accuracy

Setup. In Fig. 2, we report the average accuracy and standard deviation (+
sign: standard deviation $>0.05$) over five trials for each setting, except
for Magna where we do not simulate alignments.

Results. While the existing multilevel alignment method, MOANA, has accuracy
below its single-level counterpart FINAL as expected, our multilevel
framework, CAPER, significantly outperforms different base alignment methods
as well as their single-level refined variants using RefiNA. Moreover, we can
see that CAPER is more robust to noise due to the multilevel consistency that
it encourages; this is especially notable for CAPER(REGAL) whose performance
is very stable even when the noise level increases.

###  5.2. Alignment Runtime

Setup & Evaluation. Due to GWL’s slow runtime, we only run it for one trial on
the largest Facebook dataset. Others are averaged over five trials in Fig. 3.

Results. For faster base methods such as FINAL and especially REGAL, our
improvements are mainly in accuracy (up to 50% higher accuracy); the
computational savings of performing the alignment on smaller graphs does not
outweigh the overhead of coarsening and refinement when the base alignment
method is fast. However, when the base alignment method is slow, as is the
case for GWL, our framework results in considerable computational savings
(5-80$\times$ faster), while largely preserving accuracy. Compared with MOANA,
CAPER achieves better accuracy.

###  5.3. Sensitivity Analysis

Number of levels. Figure 4 compares the performance of CAPER(REGAL) with
different numbers of coarsening levels on the Hamsterster and Facebook
datasets. The number of coarsening levels leads to a tradeoff between accuracy
and runtime: more coarsening leads to smaller coarsened graphs (Tab. 3) and
faster runtime at a cost of some accuracy. For our main experiments, we used 3
levels of coarsening for all datasets to balance this tradeoff, and could use
2 levels to achieve even higher accuracy.

(a) Hamsterster

(b) Facebook

Figure 4. Sensitivity to number of coarsening levels for CAPER(REGAL). In
general, 2 levels leads to highest accuracy. We use 3 levels for the best
accuracy/runtime tradeoff.

Hard vs. soft refinement. In Fig. 5, we see improvement from our refinement of
more expressive “soft” alignments (§ 4.4), most noticeably for the base method
REGAL. For FINAL, because its initial solution is less accurate on these
datasets, we used hard refinement when operating directly on its solution (at
the coarsest level) and soft refinement at subsequent levels. This also
explains the smaller gap in performance.

(a) Magna

(b) Facebook

Figure 5. Sensitivity to soft/hard refinement for CAPER(REGAL) and
CAPER(FINAL). Soft refinement works significantly better, especially for
accurate base methods.

##  6\. Conclusion

We describe the first general-purpose multilevel framework for unsupervised
network alignment. It works with various base network alignment algorithms,
making them more accurate and robust by incorporating multiscale graph
information, and accelerating the runtime by allowing them to operate on
smaller input graphs. However, not all coarsening methods work well. Some
recent spectral coarsening methods (graphzoom, ) will give clusters with zero
nodes and thus our multi-level alignment framework could fail. One possible
future direction is to characterize the effect of various coarsening methods
on multilevel network alignment.

## Acknowledgements

This work was performed under the auspices of the U.S. Department of Energy by
Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344,
Lawrence Livermore National Security, LLC. and was supported by the LLNL-LDRD
Program under Project No. 21-ERD-012. Jing Zhu was an intern at Lawrence
Livermore National Laboratory while working on this project. It was also
partially supported by the NSF under Grant No. IIS 1845491, and Amazon and
Facebook faculty awards.

## References

  * [1] Xiyuan Chen, Mark Heimann, Fatemeh Vahedian, and Danai Koutra.  Cone-align: Consistent network alignment with proximity-preserving node embedding.  In CIKM, 2020.
  * [2] Chenhui Deng, Zhiqiang Zhao, Yongyu Wang, Zhiru Zhang, and Zhuo Feng.  Graphzoom: A multi-level spectral approach for accurate and scalable graph embedding.  In ICLR, 2020.
  * [3] Inderjit S Dhillon, Yuqiang Guan, and Brian Kulis.  Weighted graph cuts without eigenvectors a multilevel approach.  IEEE transactions on pattern analysis and machine intelligence, 29(11):1944–1957, 2007.
  * [4] Mark Heimann, Xiyuan Chen, Fatemeh Vahedian, and Danai Koutra.  Refining network alignment to improve matched neighborhood consistency.  In SDM, 2021.
  * [5] Mark Heimann, Haoming Shen, Tara Safavi, and Danai Koutra.  Regal: Representation learning-based graph alignment.  In CIKM, 2018.
  * [6] Junchen Jin, Mark Heimann, Di Jin, and Danai Koutra.  Toward understanding and evaluating structural node embeddings.  ACM Transactions on Knowledge Discovery from Data (TKDD), 16(3):1–32, 2021.
  * [7] Ehsan Kazemi.  Network alignment: Theory, algorithms, and applications.  Technical report, EPFL, 2016.
  * [8] Jérôme Kunegis.  Konect: the koblenz network collection.  In WWW, 2013.
  * [9] Alexander Frederiksen Kyster, Simon Daugaard Nielsen, Judith Hermanns, Davide Mottin, and Panagiotis Karras.  Boosting graph alignment algorithms.  In Proceedings of the 30th ACM International Conference on Information & Knowledge Management, pages 3166–3170, 2021.
  * [10] Jure Leskovec and Andrej Krevl.  SNAP Datasets: Stanford large network dataset collection.  http://snap.stanford.edu/data, June 2014.
  * [11] Jiongqian Liang, Saket Gurukar, and Srinivasan Parthasarathy.  Mile: A multi-level framework for scalable graph embedding.  In ICWSM, 2021.
  * [12] Yike Liu, Tara Safavi, Abhilash Dighe, and Danai Koutra.  Graph summarization methods and applications: A survey.  ACM computing surveys (CSUR), 51(3):1–34, 2018.
  * [13] Kyle K Qin, Flora D Salim, Yongli Ren, Wei Shao, Mark Heimann, and Danai Koutra.  G-crewe: Graph compression with embedding for network alignment.  In Proceedings of the 29th ACM International Conference on Information & Knowledge Management, pages 1255–1264, 2020.
  * [14] Ryan A Rossi, Di Jin, Sungchul Kim, Nesreen K Ahmed, Danai Koutra, and John Boaz Lee.  On proximity and structural role-based embeddings in networks: Misconceptions, techniques, and applications.  ACM Transactions on Knowledge Discovery from Data (TKDD), 14(5):1–37, 2020.
  * [15] Vikram Saraph and Tijana Milenković.  Magna: maximizing accuracy in global network alignment.  Bioinformatics, 30(20):2931–2940, 2014.
  * [16] Hongteng Xu, Dixin Luo, Hongyuan Zha, and Lawrence Carin Duke.  Gromov-wasserstein learning for graph matching and node embedding.  In ICML, pages 6932–6941, 2019.
  * [17] Yujun Yan, Jiong Zhu, Marlena Duda, Eric Solarz, Chandra Sripada, and Danai Koutra.  GroupINN: Grouping-based interpretable neural network for classification of limited, noisy brain data.  In KDD, pages 772–782, 2019.
  * [18] Si Zhang and Hanghang Tong.  Final: Fast attributed network alignment.  In KDD, 2016.
  * [19] Si Zhang, Hanghang Tong, Ross Maciejewski, and Tina Eliassi-Rad.  Multilevel network alignment.  In The World Wide Web Conference, pages 2344–2354, 2019.
  * [20] Zhen Zhang, Yijian Xiang, Lingfei Wu, Bing Xue, and Arye Nehorai.  KerGM: Kernelized Graph Matching.  In NeurIPS19, pages 3330–3341, 2019.