# The Discovery of the Zeeman Effect in 38 GHz Class II Methanol Masers

E. Momjian  National Radio Astronomy Observatory
P. O. Box O, Socorro, NM, USA  A. P. Sarma  Department of Physics &
Astrophysics, DePaul University
2219 N. Kenmore Ave, Byrne 211, Chicago IL 60614, USA
Received 2023 August 8; revised 2023 September 6; accepted 2023 September 26,
for publication in ApJ

###### Abstract

Magnetic fields likely play an important role in star formation, but the
number of directly measured magnetic field strengths remains scarce. We
observed the 38.3 and 38.5 GHz Class II methanol (CH3OH) maser lines toward
the high mass star forming region NGC 6334 F for the Zeeman effect. The
observed spectral profiles have two prominent velocity features which can be
further decomposed through Gaussian component fitting. In several of these
fitted Gaussian components we find significant Zeeman detections, with
$zB_{\rm los}$ in the range from 8 to 46 Hz. If the Zeeman splitting factor
$z$ for the 38 GHz transitions is of the order of $\sim$1 Hz mG-1, similar to
that for several other CH3OH maser lines, then magnetic fields in the regions
traced by these masers would be in the range of 8-46 mG. Such magnetic field
values in high mass star forming regions agree with those detected in the
better-known 6.7 GHz Class II CH3OH maser line. Since Class II CH3OH masers
are radiatively pumped close to the protostar and likely occur in the
accretion disk or the interface between the disk and outflow regions, such
fields likely have significant impact on the dynamics of these disks.

Interstellar magnetic fields (845), Star forming regions (1565), Star
formation (1569), Astrophysical masers (103), Interstellar molecules (849),
Interstellar medium (847), High resolution spectroscopy (2096),
Spectropolarimetry (1973)

††facilities: VLA††software: AIPS (Greisen, 1990)

##  1 Introduction

A complete understanding of the role of magnetic fields in star formation
remains elusive, even though major strides have been made in understanding how
magnetic fields impact the star formation process (Pattle et al. 2023;
Tsukamoto et al. 2023). Masers offer the opportunity to observe star forming
regions at high angular resolution, because they are bright and compact
sources (Richards et al., 2020). In particular, Class II methanol (CH3OH)
masers are known to form close to protostars, pumped by infrared radiation
from the protostar itself. Ellingsen et al. (2018) observed 38.3 and 38.5 GHz
Class II CH3OH maser transitions at high angular resolution ($2\hbox
to0.0pt{\arcsec\hss}{.}\,3\times 1\hbox to0.0pt{\arcsec\hss}{.}\,4$) with the
Australia Telescope Compact Array (ATCA). One of the objects in their sample,
NGC 6334 F, shows strong maser lines ($>100$ Jy) at these frequencies.
Consequently, this prompted us to target NGC 6334 F for the Zeeman effect,
which is the most direct method for measuring magnetic fields in star forming
regions.

There are two classes of CH3OH masers: Class I and Class II (Menten 1991).
Class I CH3OH masers in star forming regions are known to be collisionally
pumped in outflows (see, e.g., Leurini et al. 2016). Class II CH3OH masers,
meanwhile, are pumped by infrared radiation and are therefore located close to
protostars (Cragg et al. 2005). They are exclusively associated with high mass
star forming regions (Ellingsen 2006). The most observed and well known Class
II transitions are at 6.7 GHz and 12.2 GHz (see, e.g., Nguyen et al. 2022, and
references therein), but about 20 different Class II CH3OH maser transitions
have been observed to date (Breen et al. 2019, and references therein). Maser
emission in the Class II CH3OH lines at 38.3 and 38.5 GHz was discovered by
Haschick et al. (1989) through single-dish observations. Ellingsen et al.
(2018) carried out radio interferometric observations of the 38.3 and 38.5 GHz
methanol masers toward a sample of 11 high mass star forming regions that host
strong 6.7 GHz CH3OH masers and detected 38.3 GHz transitions toward seven
sources in their sample, and 38.5 GHz transitions toward six sources. They
found that these transitions arise from the same location as the strong 6.7
GHz Class II CH3OH maser transitions, although they are less co-spatial with
the 6.7 GHz maser spots compared to the 12.2 GHz masers.

NGC 6334 is a well-known molecular cloud complex and star-forming region
located at a distance of 1.3 kpc (Chibueze et al. 2014). Strong infrared and
mm sources, hypercompact (HC) and ultracompact (UC) H II regions, and powerful
outflows reveal evidence of ongoing high mass star formation activity in this
region (Rodríguez et al. 2007; André et al. 2016; Brogan et al. 2016; Hunter
et al. 2021). At 4.9 GHz, Rodriguez et al. (1982) found six discrete continuum
sources in NGC 6334. These six sources lie along a ridge of emission parallel
to the Galactic Plane and were named A-F in order of increasing Right
Ascension (R.A.), as shown in Figure 1. Source F, the target of the
observations reported in this paper, is highly obscured at optical
wavelengths, but prominent in cold dust traced by 850 $\mu$m submm continuum
observations (Matthews et al. 2008). The dust clump associated with source F
was determined by Matthews et al. (2008) to have a mass of $2000~{}M_{\odot}$,
indicating that it is a large reservoir of material for star formation. Also
coincident with source F is the far-infrared (FIR) source I (Emerson et al.
1973; Loughran et al. 1986). In the CO $J=3-2$ transition, Zhang et al. (2014)
observed a blueshifted outflow toward the southwest of the FIR source I, and a
redshifted outflow toward the northeast. Hunter et al. (2006) discovered four
1.3 mm continuum sources toward NGC 6334 I; these were later renamed as
MM1-MM4 by Brogan et al. (2016). Of these, the ultracompact (UC) H II region
MM3 overlaps with the peak of NGC 6334 F. An outburst in the mm emission from
MM1, accompanied by simultaneous flaring in multiple maser species, was
ascribed to episodic accretion in this source (Hunter et al. 2021, and
references therein).

Figure 1: Contour image of the 1.67 GHz continuum from NGC 6334 taken from
Sarma et al. (2000), showing the sources A-F (Rodriguez et al. 1982). The
observations reported in this work are toward source F.

In this paper, we report the detection of the Zeeman effect in the 38.3 GHz
and 38.5 GHz Class II CH3OH masers toward NGC 6334 F. The details of the 2017
observations and data reduction are given in § 2. The results are presented in
§ 3, and discussed in § 4. In § 5, we state our conclusions.

##  2 Observations and Data Reduction

The observations of the Class II CH3OH maser emission lines $6_{2}\rightarrow
5_{3}$ A- at 38.3 GHz and $6_{2}\rightarrow 5_{3}$ A+ at 38.5 GHz toward NGC
6334 F were carried out with the Karl G. Jansky Very Large Array (VLA)111The
National Radio Astronomy Observatory (NRAO) is a facility of the National
Science Foundation, operated under cooperative agreement by Associated
Universities, Inc. on 2021 March 23 and 2021 April 7. Each of these observing
sessions was 2 hr long. The VLA was in its most compact (D) configuration with
a maximum baseline length of 1 km.

Table 1: Parameters for VLA Observations  Parameter | 38.3 GHz Value | 38.5 GHz Value
---|---|---
Date | 2021 Mar 23 & Apr 7
Configuration | D
R.A. of field center (J2000) | 17${}^{\text{h}}$ 20${}^{\text{m}}$ 53$\fs$370
Dec. of field center (J2000) |  $-$35° 47′ 02″. 0
Total bandwidth (MHz) | 4
No. of channels | 1024
Channel spacing (km s-1) | 0.0305
Approx. time on source (hr) … | 2.83
Rest frequency (GHz) | 38.293270 | 38.452629
FWHM of synthesized beam | $5\,\hbox to0.0pt{\arcsec\hss}.\,33\times 1\,\hbox to0.0pt{\arcsec\hss}.\,53$ | $5\,\hbox to0.0pt{\arcsec\hss}.\,60\times 1\,\hbox to0.0pt{\arcsec\hss}.\,50$
| P.A. = 6.99° | P.A. = 6.63°
Line rms noise (mJy beam-1) aaThe line rms noise was measured from the Stokes $I$ image cube using maser line free channels. | 10.5 | 12.0

The Wideband Interferometric Digital ARchitecture (WIDAR) correlator was
configured to deliver 4 MHz sub-bands with dual polarization products (RR, LL)
and 1024 spectral channels per sub-band. The resulting channel spacing was
3.90625 kHz, which corresponds to $\sim 0.0305$ km s-1 at the observed
frequencies. In addition to NGC 6334 F, the source J1331+3030 (3C 286) was
observed as the absolute flux density scale calibrator. The uncertainty in the
flux density calibration at the observed frequencies, accounting for various
observational parameters (e.g., weather, reference pointing, and elevation
effects), is expected to be up to 10%. We also observed the source
J1720$-$3552 as the complex gain calibrator for part of each observing session
to derive the absolute position of the masers in the target source. The phase
referencing cycle time was six minutes; 5 minutes on the target and 1 minute
on the complex gain calibrator. All the data reduction steps, including
calibration, imaging, and deconvolution, were carried out independently for
each observing session using the Astronomical Image Processing System (AIPS;
Greisen 1990). After Doppler correcting the spectral-line data of each
transition and observing session independently, the spectral channel with the
brightest maser emission in each was split off, and self-calibrated first in
phase, then in both phase and amplitude, and imaged in iterative cycles. The
final self-calibration solutions were then applied to the corresponding full
spectral-line $uv$ data sets of NGC 6334 F from each session. We note that
AIPS calculates the Stokes parameter $I$ as the average of the right circular
polarization (RCP) and left circular polarization (LCP), so that $I$ = (RCP +
LCP)/2, whereas Stokes V is calculated by AIPS as half the difference between
RCP and LCP, so that $V$ = (RCP $-$ LCP)/2; henceforth, all values of $I$ and
$V$ are based on this implementation in AIPS. Also, we note that RCP is
defined here in the standard radio convention, in which it is the clockwise
rotation of the electric vector when viewed along the direction of wave
propagation. Table 1 summarizes the parameters of the VLA observations and the
corresponding synthesized beamwidths at full width half maximum (FWHM) and
other parameters for each observing session.

##  3 Results

We detected Class II CH3OH masers at 38.3 and 38.5 GHz toward the position
located at R.A. (J2000) = 17${}^{\text{h}}$ 20${}^{\text{m}}$ 53$\fs$37, Decl.
(J2000) = $-$35° 47′ 01″. 40 in NGC 6334 F. The spectral profile of the Class
II methanol maser line detected at 38.3 GHz toward this position is shown in
Figure 2. It has two prominent spectral features in velocity. We have labeled
them as I and II in the figure, and separated them by a vertical dashed line.
Visually, spectral feature I looks like it could be fit by a single Gaussian
component, but it did require a second shallower and broader component to fit
the line wing on the side nearer to $-11$ km/s. We have labeled these as
component I and I-s respectively; despite the risk of confusion, we have
chosen this nomenclature to emphasize that the spectral feature I in Figure 2
is almost a single-component Gaussian. The peak intensity of components I and
I-s, together with the velocity at line center and the Full Width at Half
Maximum (FWHM) linewidth, are listed in Table 2 and shown in Figure 3.
Henceforth, we have tried to be clear from the context whether we are
referring to spectral feature I shown in Figure 2, or component I that is
listed in Table 2 and shown in Figure 3. The peak intensity of component I
listed in Table 2 is 64.49 Jy beam-1 and its center velocity is $-11.24$ km
s-1. It is narrow, with a linewidth of 0.18 km s-1. Component I-s is
shallower, with about 1/4 the peak intensity of component I. It is also
broader, with a FWHM linewidth of 0.327 km s-1. Note that even though
component I-s extends into the velocity space of component II, it was not used
in the fit for component II. Meanwhile, at 38.5 GHz we observe a spectral
profile that resembles the 38.3 GHz profile. As listed in Table 2, the
intensity of component I at 38.5 GHz is similar to that at 38.3 GHz; likewise
for component I-s. Within the errors, the center velocities and FWHM
linewidths at 38.3 and 38.5 GHz are the same for components I and I-s.

Figure 2: Observed spectral profile (black histogram-like line) toward NGC 6334 F at 38.3 GHz showing two prominent spectral features in velocity, which we have labeled as I and II. The dashed vertical line separates these two features. The blue curve to the left of the dashed vertical line is the resultant profile obtained by summing components I and I-s listed in Table 2 and shown in Figure 3, and the blue curve to the right of the dashed vertical line is the resultant profile obtained by summing components IIa, IIb, and IIc. The red dashed curve represents the residuals from the fit.  Table 2: Maser Line Parameters from Gaussian Fits  Frequency |  |  |  | Intensity |  | Center Velocity |  | Velocity Linewidthaa The velocity linewidth was measured at FWHM.
---|---|---|---|---|---|---|---|---
(GHz) |  | Component |  | (Jy beam-1) |  | (km s-1) |  | (km s-1)
38.3 |  | I |  | $64.49\pm 5.41$ |  | $-11.24\pm 0.02$ |  | $0.182\pm 0.020$
38.5 |  | I |  | $67.34\pm 3.89$ |  | $-11.25\pm 0.02$ |  | $0.181\pm 0.019$
38.3 |  | I-s |  | $16.62\pm 1.37$ |  | $-10.99\pm 0.04$ |  | $0.327\pm 0.205$
38.5 |  | I-s |  | $18.98\pm 3.63$ |  | $-10.97\pm 0.09$ |  | $0.356\pm 0.219$
38.3 |  | IIa |  |  $100.51\pm 8.54$1 |  | $-10.59\pm 0.02$ |  | $0.223\pm 0.021$
38.5 |  | IIa |  | $96.81\pm 5.14$ |  | $-10.60\pm 0.02$ |  | $0.228\pm 0.020$
38.3 |  | IIb |  | $78.33\pm 5.98$ |  | $-10.75\pm 0.02$ |  | $0.205\pm 0.020$
38.5 |  | IIb |  | $73.21\pm 3.55$ |  | $-10.76\pm 0.02$ |  | $0.202\pm 0.019$
38.3 |  | IIc |  | $35.67\pm 5.22$ |  | $-10.84\pm 0.04$ |  | $0.387\pm 0.057$
38.5 |  | IIc |  | $37.18\pm 2.86$ |  | $-10.83\pm 0.02$ |  | $0.403\pm 0.037$

Three Gaussian components were fitted to spectral feature II at 38.3 GHz; we
have designated them as IIa, IIb, and IIc. The intensity, center velocity, and
FWHM linewidth of these three components are listed in Table 2 and shown in
Figure 3. The peak intensities of components IIa and IIb, 100.51 and 78.33 Jy
beam-1 respectively, are larger than that of component I listed in Table 2,
and their center velocities are $-10.59$ and $-10.75$ km s-1. Like component
I, both component IIa and IIb are narrow, with FWHM linewidths of 0.22 km s-1
and 0.21 km s-1 respectively. Meanwhile, component IIc has lower intensity
(35.67 Jy beam-1), broader FWHM linewidth (0.39 km s-1), and a center velocity
of $-10.84$ km s-1.

Figure 3: Same spectral profile (black histogram-like line) as in Figure 2,
but now showing the individual Gaussian components fitted to the observed
profile at 38.3 GHz. The solid red and dashed blue curves show components I
and I-s listed in Table 2. The magenta, solid blue, and green curves show
components IIa, IIb, and IIc respectively. Note that for visual clarity, the
entire dashed blue and solid green curves are shown, even though the dashed
blue curve is used only in the fit for spectral feature I shown in Figure 2,
and the solid green curve is used only in the fit for feature II.

The resultant Gaussian profile obtained by summing component I and I-s, and
separately summing IIa, IIb, and IIc is shown in Figure 2, together with the
residuals from the fit. Again, the situation was similar at 38.5 GHz, where
components IIa and IIb had higher intensity than component I (see Table 2).
Within the errors, the center velocities and linewidths of IIa and IIb at 38.5
GHz were the same as those at 38.3 GHz. Component IIc also had lower intensity
and broader linewidth than IIa and IIb at 38.5 GHz. In the approach that was
followed, as described above, the two velocity features in both the 38.3 and
38.5 GHz maser transitions were treated separately and independently when
fitting Gaussian components to their profiles. This was necessitated by the
simultaneous optimization of the Gaussian fits to the Stokes $I$ spectra and
their scaled derivatives to the Stokes $V$ spectra. This process also took
into account the minimization of the residuals while using as few Gaussian
components as possible.

The observed spectral profiles shown in Figure 2 and 3 are the Stokes $I$
profiles at 38.3 GHz. The Stokes $I$ and $V$ profiles toward NGC 6334 F at
38.3 GHz are shown together in Figure 4. Following our usual procedure (see,
e.g., Momjian & Sarma 2017), we fit the Stokes $V$ profile to the derivative
of the Stokes $I$ profile and a scaled replica of the $I$ profile using the
equation (Troland & Heiles 1982; Sault et al. 1990):

$V=aI+\frac{b}{2}\,\frac{dI}{d\nu}$ (1)

The fit parameter $a$ is included to account for small calibration errors in
RCP versus LCP and we obtained $a\lesssim 10^{-3}$. The fit parameter
$b=zB_{\rm los}$, where $z$ is the Zeeman splitting factor and $B_{\rm los}$
is the line-of-sight magnetic field strength. We fitted independently for
spectral features I and II shown in Figure 2; the dashed vertical line in
Figure 4 separates the channels that we included in the fits for each of these
two features. The fits were carried out using the AIPS task ZEMAN (Greisen
2017), which allows for multiple Gaussian components in Stokes $I$ to be
fitted simultaneously to Stokes $V$ for different values of $b$ as was done,
for example, for the three components IIa-IIc. Figure 5 shows the individual
Gaussian components in the upper panel, and the derivatives of each component
scaled by the fitted value of $zB_{\rm los}$ in the lower panel. Unlike in
Figure 3, components I and I-s are shown only to the left of the dashed
vertical line, and components IIa-IIc are shown to the right of the dashed
vertical line; that is, confined to the velocity space in which they were
summed to obtain the blue profiles superposed on the black histogram-like
profiles in the upper and lower panels of Figure 4. Figure 6 is the equivalent
of Figure 4 and Figure 7 is the equivalent of Figure 5, but at 38.5 GHz.

The results of the fitting of Stokes $V$ using equation (1) are given in Table
3. Although Lankhaar et al. (2018) published values of the Zeeman splitting
factor for several CH3OH maser lines, values of $z$ for 38.3 and 38.5 GHz
lines are not available, so we will leave our results in terms of $zB_{\rm
los}$ in Hz. We obtained three significant detections and two upper limits.
Generally, we consider a detection significant if the signal-to-noise ratio in
$zB_{\rm los}$ is 3-$\sigma$ or higher. Here, we impose a stronger criterion;
we will consider the detection significant only if the signal-to-noise ratio
in $zB_{\rm los}$ is 3-$\sigma$ or higher at both 38.3 GHz and 38.5 GHz.
Component I listed in Table 2 shows a significant detection of $-45.10\pm
3.40$ Hz at 38.3 GHz, and $-46.08\pm 3.63$ Hz at 38.5 GHz. By convention, a
negative value for $B_{\rm los}$ implies that the line-of-sight magnetic field
is pointing toward the observer. Component IIa and IIb also show significant
detections, and $B_{\rm los}$ traced by component IIa points toward the
observer, whereas $B_{\rm los}$ traced by component IIb points away from the
observer. Component IIa shows a significant detection of $-7.48\pm 1.81$ Hz at
38.3 GHz, and $-12.88\pm 2.55$ Hz at 38.5 GHz. Finally, component IIb shows a
significant detection of $16.26\pm 2.59$ Hz at 38.3 GHz, and $20.07\pm 3.96$
Hz at 38.5 GHz. Although a value of $z$ is not available for these two
transitions, values published by Lankhaar et al. (2018) for many prominent
CH3OH transitions are near 1 Hz mG-1, and $B_{\rm los}$ values of 8-46 mG
appear reasonable in such regions.

Figure 4: Stokes $I$ (upper panel, black histogram-like line) and Stokes $V$ (lower panel, black histogram-like line) profiles toward NGC 6334 F at 38.3 GHz. The blue curve in the upper panel is the same as that shown in, and described in the caption to, Figure 2. The blue curve superposed on the Stokes $V$ profile in the lower panel is the sum of the solid red and dashed blue curves shown in the lower panel of Figure 5 to the left of the dashed vertical line, and the sum of the magenta, solid blue, and green curves in the lower panel of Figure 5 to the right of the dashed vertical line.  Table 3: Zeeman Effect Measurements  Frequency |  |  |  | z$B_{\rm los}$
---|---|---|---|---
(GHz) |  | Component |  | (Hz)
Significant Detections
38.3 |  | I |  | $-45.10\pm 3.40$
38.5 |  | I |  | $-46.08\pm 3.63$
38.3 |  | IIa |  |  $-7.48\pm 1.81$
38.5 |  | IIa |  | $-12.88\pm 2.55$
38.3 |  | IIb |  | $\phantom{-}16.26\pm 2.59$
38.5 |  | IIb |  | $\phantom{-}20.07\pm 3.96$
Upper Limits
38.3 |  | I-s |  | $-15\pm 49$
38.5 |  | I-s |  | $-27\pm 41$
38.3 |  | IIc |  | $\phantom{-}68\pm 10$
38.5 |  | IIc |  | $\phantom{-}34\pm 14$
Figure 5: Stokes $I$ (upper panel, black histogram-like line) and Stokes $V$
(lower panel, black histogram-like line) profiles toward NGC 6334 F at 38.3
GHz, as in Figure 4. The solid red and dashed blue curves in the upper panel
are Gaussian components I and I-s, respectively. The magenta, solid blue, and
green curves are components IIa, IIb, and IIc, respectively. The fit
parameters for all five components are listed in Table 2. The curves
superposed on Stokes $V$ in the lower panel are the derivatives of the
corresponding colored curves in the upper panel, scaled by the fitted value of
$zB_{\rm los}$ for each curve listed in Table 3.

##  4 Discussion

We observed the 38.3 and 38.5 GHz Class II CH3OH maser transitions toward NGC
6334 F with the VLA for the Zeeman effect. At both 38.3 and 38.5 GHz, the
observed spectral profiles contain two prominent velocity features. We fitted
one with a narrow Gaussian component and a shallow broad component, which we
have labeled as components I and I-s respectively (Table 2 and Figure 3). The
other was fitted with two narrow Gaussian components labeled IIa and IIb, and
a lower intensity broad component, which we have labeled as IIc. Our observed
maser profiles are consistent with the Australia Telescope Compact Array
(ATCA) observations of Ellingsen et al. (2018). However, it is worth noting
that their velocity resolution was 0.3 km s-1, a factor of $\sim$10 coarser
compared to ours. They observed our two prominent velocity features (which we
have labeled as I and II in Figure 2) to be blended in velocity space, but
with a significant asymmetry on the left side of the profile (near $-11$ km
s-1), consistent with the presence of our spectral feature I. Of interest is
that the single dish observations of Haschick et al. (1989) with 0.065 km s-1
velocity resolution show feature II but only weak emission at the velocity
corresponding to feature I. This could be due to variability, but there is no
way to directly compare the two observations taken with such different
instruments so many years apart. The base of the spectral profile observed by
Haschick et al. (1989) at both 38.3 and 38.5 GHz indicates the presence of a
broad emission component, consistent with the shallower and broader components
I-s and IIc that we have fitted (Figure 3). In both the single dish
observations of Haschick et al. (1989) and the ATCA observations of Ellingsen
et al. (2018), the 38.5 GHz spectral profile has higher intensity than the
38.3 GHz profile, consistent with our VLA observations (Table 2).

We have obtained significant detections of the Zeeman effect in components I,
IIa, and IIb (Table 3), with $zB_{\rm los}$ in the range 8-46 Hz. Obtaining
the line-of-sight magnetic field values from these is possible only if we know
the value of the Zeeman splitting factor $z$ for the 38.3 and 38.5 GHz
transitions. Lankhaar et al. (2018) published values of $z$ for a wide range
of CH3OH maser transitions, but not for the 38 GHz transitions. If, however,
the values of $z$ for these transitions are near 1 Hz mG-1, as it is for many
prominent transitions of CH3OH, then we would get $B_{\rm los}$ values of the
same order as $zB_{\rm los}$ in these regions. Values in the range 8-46 mG
appear reasonable for high mass star forming regions. For example, from the
statistical analysis of a flux-limited sample of 6.7 GHz Class II CH3OH
masers, Surcis et al. (2022) found line-of-sight magnetic fields in the range
9-40 mG. It is worth noting that the CH3OH maser transitions result from a
series of hyperfine transitions, each of which has a different $z$. Therefore,
the magnetic field values would be different depending on which hyperfine
transition, or combination of hyperfine transitions, is responsible for the
maser transition. Also, it is worth noting that the reversal in sign of
$B_{\rm los}$ from component IIa to IIb is usually interpreted by astronomers
as a field reversal in the regions traced by these masers. However, Lankhaar
et al. (2018) have proposed a different scenario in which the change in sign
is caused by the population inversion of two different hyperfine transitions.
Values in the range 8-46 mG are also in general agreement with field strengths
of 4.4-6.4 mG measured in 1.6 and 6.0 GHz OH maser lines in NGC 6334 F by
Chanapote et al. (2019).

Figure 6: As in Figure 4, but at 38.5 GHz.

Following standard procedure in reporting the Zeeman effect, we discuss now
whether the detected signal could be caused by instrumental effects, or by
processes other than the Zeeman effect. It is well known that a velocity
gradient across an extended source could mimic a Stokes $V$ signal caused by
the Zeeman effect. However, masers are compact sources confined to a narrow
velocity range. Therefore, it is unlikely that the observed Stokes $V$ signal
reported in this paper is due to such a velocity gradient. Moreover, the close
agreement of the detected $zB_{\rm los}$ value at 38.3 and 38.5 GHz, i.e., two
transitions observed in independent spectral windows, gives us confidence that
this is not a spurious detection. Processes other than the Zeeman effect that
could contribute to structure in the Stokes $V$ profile include, for masers
with strong linear polarization, changes in the orientation of the magnetic
field along the line of sight that could cause rotation of the linear
polarization vectors to produce circular polarization (Wiebe & Watson 1998).
In order to optimize our observations for detection of the Zeeman effect, we
did not propose for full polarization observations. We note, however, that for
a wide array of Class II CH3OH maser transitions, Breen et al. (2019) found
linearly polarized emission in the range 1.5-7.5%. If the 38 GHz Class II
CH3OH maser transitions reported in this paper have similar levels of linear
polarization, then it is unlikely that the observed Stokes $V$ could be caused
by the rotation of linear polarization vectors due to changes in the magnetic
field orientation along the line of sight. Another possible non-Zeeman origin
comes from maser radiation scattering off foreground molecules that can
enhance antisymmetric spectral profiles in Stokes $V$ (Houde 2014).

Figure 7: As in Figure 5, but at 38.5 GHz.

In such cases, if the Stokes $V$ were considered to be caused solely by the
Zeeman effect, then the resulting magnetic field strengths would be too large,
as demonstrated for SiO masers by Houde (2014). Unless $z$ for 38 GHz
transitions is determined to be extremely low so that the values of $zB_{\rm
los}$ we have detected translate to very high magnetic fields, it appears
unlikely that this effect is of consequence for our observations. Finally, if
the maser stimulated emission rate $R$ is larger than the frequency shift due
to the Zeeman effect $g\Omega$, which is the product of the Landé g-value for
the upper state of the transition and the gyrofrequency, $\Omega$, of the
electron, a rotation of the axis of symmetry for the molecular quantum states
could also cause circular polarization (Vlemmings et al. 2011). In our
observations toward NGC 6334 F, we have $g\Omega\approx 10$ s-1. The
stimulated emission rate $R$, taken from Vlemmings et al. (2011), is

$R\simeq\frac{AkT_{b}\Delta\Omega}{4\pi h\nu}$ (2)

where $k$ and $h$ are the Boltzmann and Planck constants, respectively, $\nu$
is the frequency of the maser transition, $A$ is the Einstein coefficient,
equal to $4.726\times 10^{-8}$ s-1 for the $6_{2}\rightarrow 5_{3}$ A-
transition at 38.3 GHz, and $4.829\times 10^{-8}$ s-1 for the
$6_{2}\rightarrow 5_{3}$ A+ transition at 38.5 GHz, $T_{b}$ is the maser
brightness temperature, and $\Delta\Omega/4\pi\approx 10^{-3}$ is a
conservative estimate for the maser beaming angle (Nesterenok 2016). The lower
limit on $T_{b}$ from our observations is $10^{6}$ K. Using these values,
equation (2) gives $R\approx 10^{-5}$ s-1, which in turn implies that $R\ll
g\Omega$. Therefore it is unlikely that a rotation of the axis of symmetry is
causing the splitting that results in the observed Stokes $V$ profile.
Moreover, such an effect would cause an intensity-dependent polarization, but
component IIa with higher intensity than IIb has a lower $zB_{\rm los}$ than
that detected in IIb.

Class II CH3OH masers are exclusively found in high mass star forming regions
and are located close to the protostar, unlike Class I CH3OH masers which are
located in outflows. However, there is as yet no consensus on the structures
that they trace in these protostellar environments. One scenario puts them in
the accretion disk itself, whereas an alternative places them in the
intermediate region between the outflow and the disk (Sanna et al. 2010; Goddi
et al. 2011; Sugiyama et al. 2014; Bartkiewicz et al. 2020). We can use our
calculated value for $B_{\rm los}$ to compare the magnetic energy density with
other relevant quantities in NGC 6334 F. The magnetic energy density is given
by $B^{2}/8\pi$, where $B^{2}=3B_{\rm los}^{2}$, as determined by Crutcher
(1999) on statistical grounds; we note that this is strictly only valid for an
ensemble of measurements. For $B_{\rm los}$=8 mG reported in Section 3, which
is the lowest magnetic field value we measure, the magnetic energy density is
equal to $7.6\times 10^{-6}$ erg cm-3. The kinetic energy density is a
relevant quantity with which to compare this value; it is given by
$(3/2)\,mn\sigma^{2}$, where $m=2.8~{}m_{p}$, with $m_{p}$ being the proton
mass, the numerical factor of 2.8 also accounting for 10% helium. The particle
density is $n=10^{6}$ cm-3, from modeling of several Class II CH3OH maser
lines in NGC 6334 F, including the 38 GHz transitions (Cragg et al. 2001). The
velocity dispersion is given by $\sigma=\Delta v/(8\ln 2)^{1/2}$. Rather than
using the narrower $\Delta v$ from maser observations, which may not be
indicative of the thermal motions in the clump of gas where the maser
transition arises, we use the broader value of 3 km s-1, from the velocity
gradient in HC3N and NH3 emission lines toward NGC 6334 F (Jackson et al.
1988; Bachiller & Cernicharo 1990). This gives a kinetic energy density equal
to $1.1\times 10^{-7}$ ergs cm-3. This means that the magnetic energy density
in the region traced by the masers in NGC 6334 F is at least comparable to, if
not larger than, the kinetic energy density. If the 38 GHz CH3OH masers occur
in a rotating accretion disk, another relevant quantity of interest is the
rotational energy density, given by $(1/2)\,I\omega^{2}/V$, where $\omega$ is
the angular velocity and $V$ is the volume of the disk. The moment of inertia
$I$ is of the order $MR^{2}$; $M$ is the mass and $R$ is the radius of the
disk. Using a rotational velocity of $v=4$ km s-1 for methanol maser disks
from Norris et al. (1998), where $v=R\omega$, the expression for the
rotational energy density then becomes $(1/2)\,(M/V)\,v^{2}$, where
$M/V=2.8\,m_{p}\,n$, with $m_{p}$ being the proton mass, and $n$ the particle
density with a value of $10^{6}$ cm-3, as discussed above. The derived
rotational energy density is then about $3.7\times 10^{-7}$ ergs cm-3. Once
again, the magnetic energy density is at least comparable to, if not greater
than, the rotational energy density. Therefore, if the 38 GHz Class II CH3OH
masers are located in the accretion disk around the protostar, then the
magnetic field likely plays a significant role in shaping the dynamics in that
accretion disk.

##  5 Conclusion

We observed the 38.3 and 38.5 GHz Class II CH3OH maser transitions toward the
high mass star forming region NGC 6334 F for the Zeeman effect. Both
transitions have similar spectral profiles, each with two prominent spectral
features. We fitted one of these two features with a single narrow Gaussian
and a shallow broad component, with FWHM linewidths of 0.18 km s-1 and 0.33 km
s-1 respectively at 38.3 GHz. These components were labeled I and I-s
respectively. We fitted the other prominent spectral feature with two narrow
Gaussian components labeled IIa and IIb, and a broad, lower intensity
component IIc; at 38.5 GHz, component IIc with a FWHM linewidth of 0.403 km
s-1 is almost twice as wide as component IIb, which has a FWHM linewidth of
0.202 km s-1. The center velocities of these components are quite close,
ranging from $-10.59$ km s-1 to $-11.25$ km s-1.

We have obtained significant detections of the Zeeman effect in components I,
IIa, and IIb. Values of $zB_{\rm los}$ in these masers range from 8 to 46 Hz,
and there is a change in sign from component IIa to IIb. The Zeeman splitting
factor $z$ for these 38 GHz transitions is not known, but if it is of the
order $\sim$1 Hz mG-1 as it is for a range of CH3OH maser transitions, then
the magnetic fields in the regions traced by these masers would be in the
range 8-46 mG. There is a reversal in the sign of $zB_{\rm los}$ from
component IIa to IIb, usually interpreted as a reversal in $B_{\rm los}$ from
one component to another; alternatively, it may be a result of different
hyperfine transitions being responsible for the maser transitions. Magnetic
fields in the range 8-46 mG agree well with fields detected in the better-
known 6.7 GHz CH3OH transition. Class II CH3OH masers are known to form close
to the protostar in accretion disks or in the interface region between the
accretion disk and outflow. Using the lowest magnetic field value of $B_{\rm
los}$= 8 mG from Table 3, we find that the magnetic energy density is at least
comparable to, if not greater than, the kinetic energy density and the
rotational energy density in the disk. Such fields may exert significant
influence on the dynamics of the accretion disk and likely play an important
role in the star formation process.

We thank an anonymous referee for insightful comments that helped to improve
the paper. APS acknowledges a Summer Research Grant from the University
Research Council (URC) at DePaul University.

## References

  * André et al. (2016) André, P., Revéret, V., Könyves, V., et al. 2016, A&A, 592, A54
  * Bachiller & Cernicharo (1990) Bachiller, R. & Cernicharo, J. 1990, A&A, 239, 276
  * Bartkiewicz et al. (2020) Bartkiewicz, A., Sanna, A., Szymczak, M., et al. 2020, A&A, 637, A15
  * Breen et al. (2019) Breen, S. L., Sobolev, A. M., Kaczmarek, J. F., et al. 2019, ApJ, 876, L25
  * Brogan et al. (2016) Brogan, C. L., Hunter, T. R., Cyganowski, C. J., et al. 2016, ApJ, 832, 187
  * Chanapote et al. (2019) Chanapote, T., Asanok, K., Dodson, R., et al. 2019, MNRAS, 482, 1670
  * Chibueze et al. (2014) Chibueze, J. O., Omodaka, T., Handa, T., et al. 2014, ApJ, 784, 114
  * Cragg et al. (2001) Cragg, D. M., Sobolev, A. M., Ellingsen, S. P., et al. 2001, MNRAS, 323, 939
  * Cragg et al. (2005) Cragg, D. M., Sobolev, A. M., & Godfrey, P. D. 2005, MNRAS, 360, 533
  * Crutcher (1999) Crutcher, R. M. 1999, ApJ, 520, 706
  * Ellingsen (2006) Ellingsen, S. P. 2006, ApJ, 638, 241
  * Ellingsen et al. (2018) Ellingsen, S. P., Voronkov, M. A., Breen, S. L., et al. 2018, MNRAS, 480, 4851
  * Emerson et al. (1973) Emerson, J. P., Jennings, R. E., & Moorwood, A. F. M. 1973, ApJ, 184, 401
  * Goddi et al. (2011) Goddi, C., Moscadelli, L., & Sanna, A. 2011, A&A, 535, L8
  * Greisen (1990) Greisen, E. W. 1990, Acquisition, Processing and Archiving of Astronomical Images, 125
  * Greisen (2017) Greisen, E. W. 2017, AIPS Memo, 118 http://www.aips.nrao.edu/aipsmemo.html
  * Haschick et al. (1989) Haschick, A. D., Baan, W. A., & Menten, K. M. 1989, ApJ, 346, 330
  * Houde (2014) Houde, M. 2014, ApJ, 795, 27
  * Hunter et al. (2006) Hunter, T. R., Brogan, C. L., Megeath, S. T., et al. 2006, ApJ, 649, 888
  * Hunter et al. (2021) Hunter, T. R., Brogan, C. L., De Buizer, J. M., et al. 2021, ApJ, 912, L17
  * Jackson et al. (1988) Jackson, J. M., Ho, P. T. P., & Haschick, A. D. 1988, ApJ, 333, L73
  * Lankhaar et al. (2018) Lankhaar, B., Vlemmings, W., Surcis, G., et al. 2018, Nature Astronomy, 2, 145
  * Leurini et al. (2016) Leurini, S., Menten, K. M., & Walmsley, C. M. 2016, A&A, 592, A31
  * Loughran et al. (1986) Loughran, L., McBreen, B., Fazio, G. G., et al. 1986, ApJ, 303, 629
  * Matthews et al. (2008) Matthews, H. E., McCutcheon, W. H., Kirk, H., et al. 2008, AJ, 136, 2083
  * Menten (1991) Menten, K. M. 1991, ApJ, 380, L75
  * Momjian & Sarma (2017) Momjian, E. & Sarma, A. P. 2017, ApJ, 834, 168
  * Nesterenok (2016) Nesterenok, A. V. 2016, MNRAS, 455, 3978
  * Nguyen et al. (2022) Nguyen, H., Rugel, M. R., Murugeshan, C., et al. 2022, A&A, 666, A59
  * Norris et al. (1998) Norris, R. P., Byleveld, S. E., Diamond, P. J., et al. 1998, ApJ, 508, 275
  * Pattle et al. (2023) Pattle, K., Fissel, L., Tahani, M., et al. 2023, Astronomical Society of the Pacific Conference Series, 534, 193
  * Richards et al. (2020) Richards, A. M. S., Sobolev, A., Baudry, A., et al. 2020, Advances in Space Research, 65, 780
  * Rodriguez et al. (1982) Rodriguez, L. F., Canto, J., & Moran, J. M. 1982, ApJ, 255, 103
  * Rodríguez et al. (2007) Rodríguez, L. F., Zapata, L. A., & Ho, P. T. P. 2007, ApJ, 654, L143
  * Sanna et al. (2010) Sanna, A., Moscadelli, L., Cesaroni, R., et al. 2010, A&A, 517, A78
  * Sarma et al. (2000) Sarma, A. P., Troland, T. H., Roberts, D. A., et al. 2000, ApJ, 533, 271
  * Sault et al. (1990) Sault, R. J., Killeen, N. E. B., Zmuidzinas, J., et al. 1990, ApJS, 74, 437
  * Sugiyama et al. (2014) Sugiyama, K., Fujisawa, K., Doi, A., et al. 2014, A&A, 562, A82
  * Surcis et al. (2022) Surcis, G., Vlemmings, W. H. T., van Langevelde, H. J., et al. 2022, A&A, 658, A78
  * Troland & Heiles (1982) Troland, T. H. & Heiles, C. 1982, ApJ, 252, 179
  * Tsukamoto et al. (2023) Tsukamoto, Y., Maury, A., Commercon, B., et al. 2023, Astronomical Society of the Pacific Conference Series, 534, 317
  * Vlemmings et al. (2011) Vlemmings, W. H. T., Torres, R. M., & Dodson, R. 2011, A&A, 529, A95
  * Wiebe & Watson (1998) Wiebe, D. S. & Watson, W. D. 1998, ApJ, 503, L71
  * Zhang et al. (2014) Zhang, Q., Qiu, K., Girart, J. M., et al. 2014, ApJ, 792, 116


# NUS-IDS at FinCausal 2021: Dependency Tree in Graph Neural Network for
Better Cause-Effect Span Detection

Fiona Anting Tan, See-Kiong Ng
Institute of Data Science
National University of Singapore, Singapore
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Automatic identification of cause-effect spans in financial documents is
important for causality modelling and understanding reasons that lead to
financial events. To exploit the observation that words are more connected to
other words with the same cause-effect type in a dependency tree, we construct
useful graph embeddings by incorporating dependency relation features through
a graph neural network. Our model builds on a baseline BERT token classifier
with Viterbi decoding, and outperforms this baseline in cross-validation and
during the competition. In the official run of FinCausal 2021, we obtained
Precision, Recall, and F1 scores of $95.56\%$, $95.56\%$ and $95.57\%$ that
all ranked 1st place, and an Exact Match score of $86.05\%$ which ranked 3rd
place.

##  1 Introduction

We worked on the shared task of FinCausal 2021 (Mariko et al., 2021) that aims
to identify cause and effect spans in financial news. This task builds on the
previous shared Task 2 (Mariko et al., 2020) by introducing additional
annotated data.

##### Contributions:

We propose a solution to include dependency relations in a sentence to improve
identification of cause-effect spans. We do so by representing dependency
relations in a graph neural network. Our model is an extension of a baseline
BERT token classifier with Viterbi decoding (Kao et al., 2020), and
outperforms this baseline in cross-validation and test settings. This
improvement also holds for two BERT pretrained language models that were
experimented with. During the competition, we ranked 1st for Precision,
Recall, and F1 scores and 3rd for Exact Match score.

##### Organization:

Section 2 outlines our approach for this task. Section 3 introduces the task
dataset and evaluation datasets. Our results are presented and discussed in
Section 4 while Section 5 concludes with some future directions.

##  2 Our Approach

In this section, we outline our approach 111Our source code is available at
https://github.com/tanfiona/CauseEffectDetection.. Additional architectural
and experimental details are provided in the Appendix.

###  2.1 Framing the Task

Given an example document, which could be one or multiple sentence(s) long,
the task is to identify the cause and effect substrings. We converted the span
detection task into a token classification task, similar to many state-of-the-
art methodologies for span detection (Pavlopoulos et al., 2021) and Named
Entity Recognition (Lample et al., 2016; Tan et al., 2020) tasks. Figure 1
demonstrates an example sentence that has its Cause (C) span highlighted in
green, and Effect (E) span highlighted in orange, while all other spans are
highlighted in grey. The sentence was tokenized and subsequently aligned
against the target token labels. We included the BIO format (Begin, Inside,
Outside) (Ramshaw and Marcus, 1995) in our labels to better identify the start
of spans. Thus, we have five labels: B-C, I-C, B-E, I-E and O.

Figure 1: Illustrative training example (ID: 0477.00020) with Cause (C) and
Effect (E) spans highlighted in green and orange respectively. We include
Begin (B), Inside (I) and Outside (O) prefixes to create 5 labels in our token
classification task.

###  2.2 Baseline

Kao et al. demonstrated that their BIO tagging scheme with a Viterbi decoder
(Viterbi, 1967) that utilised BERT-encoded document representations is useful
for this cause-effect span detection task. Their model topped the competition
last year across all metrics. Therefore, we adapted their pipeline and
proposed distinct additions highlighted later in Section 2.3 for improved
performance. In the immediate subsections, we motivate the benefits in
retaining the key components of Kao et al.’s model, and highlight any
differences in our approach.

####  2.2.1 BERT Embeddings

We employed the Bidirectional Encoder Representations from Transformers (BERT)
(Devlin et al., 2019) for its tokenizer and encoder model, fine-tuned on our
task. We used pretrained language models from Huggingface (Wolf et al., 2020).
Apart from bert-base-cased, we also used bert-large-cased for improved
performance.

####  2.2.2 Viterbi Decoder

The Viterbi decoding algorithm is only applied during evaluation. Since the
true cause-effect spans are consecutive sequences, but the token
classifications from the neural network could produce non-consecutive
sequences, the Viterbi algorithm serves as a forward error correction
technique and is an important element for the success of this pipeline.

####  2.2.3 Parts-of-Speech

Kao et al. (2020) showed that Parts-of-Speech (POS) did not improve their
model performance. We reconfirm this finding later in Section 4.2. However, we
found POS features to be useful inclusions for our proposed model.

Figure 2: Dependency-tree represented as directed graphs, with Cause, Effect,
and Other spans highlighted in green, orange and grey respectively.

###  2.3 Dependency Tree

Our key contribution is the inclusion of dependency tree relations into the
neural network for improved cause-effect token classification. Dependencies in
text can be mapped into a directed graph representation, where nodes represent
words and edges represent the dependency relation. Figure 2 shows some example
sentences, highlighted by their Cause and Effect labels. We notice that there
is a tendency for words of the same cause-effect label to be more connected in
these graphs and thus wish to incorporate these information into the model as
features.

Figure 3 reflects our neural network model, with the addition of our GNN
module that produces graph representations, which are then concatenated with
the BERT and POS embeddings and fed into a linear layer for token
classification. The following subsections describe the GNN module further.

Figure 3: Neural network pipeline. Appendix A.1 outlines this further. Model |  | Precision | Recall | F1 | ExactMatch
---|---|---|---|---|---
A. bert-base-cased
1 | Baseline |  | 95.62 | 95.90 | 95.76 | 88.18
2 | \+ Node features, BiLSTM |  | 95.42 | 95.95 | 95.68 | 88.07
3 | Baseline w/ POS |  | 95.39^ | 96.02 | 95.70 | 88.24
4 | \+ Node features |  | 95.66 | 95.81 | 95.73 | 88.20
5 | \+ BiLSTM |  | 95.46 | 95.99 | 95.72 | 88.06
6 | \+ Node features, BiLSTM (Proposed) |  | 95.73 | 96.05^ | 95.89 | 88.35
B. bert-large-cased
7 | Baseline |  | 95.64 | 96.01 | 95.82 | 87.65
8 | \+ Node features, BiLSTM |  | 93.75* | 95.56 | 94.61* | 86.44
9 | Baseline w/ POS |  | 94.46 | 96.08 | 95.22 | 87.28
10 | \+ Node features |  | 95.61 | 95.98 | 95.79 | 88.24
11 | \+ BiLSTM |  | 95.60 | 95.78^ | 95.69 | 87.94
12 | \+ Node features, BiLSTM (Proposed) |  | 95.72 | 96.11 | 95.91 | 88.35
Table 1: Average evaluation results over cross-validation sets from 5 random seeds, each with 3 folds. _Notes._ Scores are reported in percentages (%). Best score per Panel per column is bolded. Baseline models are our replications of the models introduced by Kao et al. (2020). For each Panel A and B, Paired T-test of the models was conducted against Row 1 and 7 respectively, with statistical significance indicated by: ***$<0.05$, **$<0.10$, *$<0.15$, ^$<0.20$. Model |  | Precision |  | Recall |  | F1 |  | ExactMatch
---|---|---|---|---|---|---|---|---
Baseline |  | 93.47 |  | 93.42 |  | 93.65 |  | 80.25
Proposed (bert-base-cased) |  | 94.24 |  | 94.21 |  | 94.37 |  | 83.23
Proposed (bert-large-cased) |  | 95.56 |  | 95.56 |  | 95.57 |  | 86.05
Best Score |  | 95.56 |  | 95.56 |  | 95.57 |  | 87.77
Our Ranking |  | 1st |  | 1st |  | 1st |  | 3rd
Table 2: Results over test sets submitted to Codalab (As of 01 September
2021). _Notes._ Scores are reported in percentages (%). Best score per column
is bolded. The models of the first three rows corresponds to Rows 1, 6, and 12
in Table 1 for CV respectively.

####  2.3.1 Document to Graph

Each example was represented as a directed graph, where nodes are token
features (the concatenation of BERT and POS embeddings), while edges are
directed connections pointing head to tail tokens 222If head or tail words are
split into multiple tokens, each head (sub) piece is connected to each tail
(sub) piece. based on dependency tree parsing 333Stanza (Qi et al., 2020) was
used for dependency parsing..

####  2.3.2 Graph Neural Network

Our graph neural network (GNN) comprised of two graph convolutional layers for
message passing across dependency relations that are two steps apart.
Specifically, we used SAGEConv operator (Hamilton et al., 2017) for its
ability to include node features and generate embeddings by aggregating a
node’s neighbouring information. To capture the long-term contextual
dependencies in both forward and backwards order of the original sentence, we
added a bi-directional long short-term memory (BiLSTM) layer (Hochreiter and
Schmidhuber, 1997) onto the graph embeddings.

##  3 Data

###  3.1 Task Dataset

Our main dataset is the FinCausal 2021 dataset (Mariko et al., 2021)
comprising of $2393$ train and $638$ competition test examples.

###  3.2 Evaluation

We evaluated our proposed models in cross-validation (CV) against the Viterbi
BERT model (Kao et al., 2020) that achieved first place in the previous run of
this shared task. We refer to this model as the Baseline in subsequent
sections. To check if our proposed models has statistically significant
improvements from the Baseline, we adapted Dietterich (1998)’s approach using
a 3-fold CV setup for 5 iterations, each iteration initialised with a random
seed 444The 5 random seeds used were $916,703,443,229,585$. This gives us
$5*3=15$ sets of evaluation results to perform Paired T-Tests for statistical
significance.

To obtain test predictions for submission on Codalab
555https://competitions.codalab.org/competitions/33102, we used a new
$seed=123$ to train on the full train data and applied the model onto the
unseen competition test data for online submission.

##  4 Results and Analysis

Table 1 reflects the model performances in CV, while Table 2 reflects the
results during the competition. In both cases, we demonstrate that our model
(Proposed) surpasses the Baseline.

For the CV setting, Proposed (Row 6) obtained $95.73\%$ Precision (P),
$96.05\%$ Recall (R), $95.89\%$ F1 and $88.35\%$ exact match (EM) scores.
These exceed the Baseline (Row 1) by $0.11\%$, $0.15\%$, $0.13\%$ and $0.18\%$
respectively 666We did not obtain P-values ($<5\%$) of statistical
significance when comparing Proposed against Baseline.. For test setting,
large performance increments were observed: Proposed achieved P/R/F1/EM scores
of $94.24\%$, $94.21\%$, $94.37\%$ and $83.23\%$, which are improvements from
Baseline by a magnitude of $0.77\%$, $0.78\%$, $0.72\%$ and $2.98\%$
respectively 777We were unable to run repeated iterations to conduct Paired
T-Tests for significance testing in competition test sets as we do not have
access to the true labels..

###  4.1 Size of Pretrained Models

The two sections of Table 1 shows that the inclusion of dependency-based graph
embeddings improved performance against Baseline irregardless of the
pretrained BERT model choice.

Between the two investigated BERT models, bert-large-cased outperforms bert-
base-cased by a small amount in CV and by a significant amount in testing
across metrics. With bert-large-cased, the Proposed model achieved P/R/F1/EM
scores of $95.56\%$, $95.56\%$, $95.57\%$ and $86.05\%$ during the
competition, which are improvements from Baseline by a magnitude of $2.09\%$,
$2.14\%$, $1.92\%$ and $5.80\%$ respectively 888We did not upload a Baseline
model using bert-large-cased, which would have allowed for a fairer
comparison, due to time and upload constraints..

###  4.2 Features and Layers

Table 1 also includes results from CV experiments where we removed components
of our model. In this subsection, we discuss the importance of each component
in the context of bert-base-cased, but note that similar findings persisted in
the bert-large-cased.

##### POS:

Inclusion of the POS into the Baseline led to mixed outcomes across the four
metrics against the Baseline (Row 3 vs Row 1). However, adding POS features in
Proposed improved performance in all metrics (Row 6 vs Row 2).

##### Node features:

We reran a model that takes in nodes with no features (i.e. all nodes are
represented by “$1$”). The results from this model corresponds to Row 5. No
obvious improvements from Baseline (Row 1) was found, however, the performance
was consistently worse off than Proposed for all metrics (Row 6), suggesting
that informative node features are important to include in the GNN.

##### BiLSTM:

Comparing our proposed model with (Row 6) and without (Row 4) the BiLSTM layer
suggested the layer was an important addition. Our hypothesis is that the
BiLSTM helps to align the graph embeddings into a sequential manner
corresponding to the original token order. A simple example is the punctuation
full-stop “.” in Figure 2. The target label of the full-stop in these cases
coincides with the immediate label of the word before it. However, our
dependency tree attributed the full-stop as a tail of another word far away
from it in the sentence (E.g. “said” $\rightarrow$ “.”).

##  5 Conclusions and Future Work

We have demonstrated the benefits of including dependency tree features as
graph embeddings in a neural network model for better cause-effect span
detection. A key caveat of our approach, which requires further research, is
that dependency parsing occurs within sentences, resulting in disconnected
graphs for examples with multiple sentences 999Preliminary experiments to tie
coreferential entities together to link dependency across sentences did not
produce fruitful results.. Another future work is to apply our cause-effect
detection model trained on financial texts onto other domains (E.g. academic
journals) to study its generalizability.

## References

  * Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2019.  BERT: Pre-training of deep bidirectional transformers for language understanding.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , pages 4171–4186, Minneapolis, Minnesota. Association for Computational Linguistics.
  * Dietterich (1998) Thomas G. Dietterich. 1998.  Approximate statistical tests for comparing supervised classification learning algorithms.  _Neural Comput._ , 10(7):1895–1923.
  * Hamilton et al. (2017) William L. Hamilton, Zhitao Ying, and Jure Leskovec. 2017.  Inductive representation learning on large graphs.  In _Advances in Neural Information Processing Systems 30: Annual Conference on Neural Information Processing Systems 2017, December 4-9, 2017, Long Beach, CA, USA_ , pages 1024–1034.
  * Hochreiter and Schmidhuber (1997) Sepp Hochreiter and Jürgen Schmidhuber. 1997.  Long Short-Term Memory.  _Neural Computation_ , 9(8):1735–1780.
  * Kao et al. (2020) Pei-Wei Kao, Chung-Chi Chen, Hen-Hsen Huang, and Hsin-Hsi Chen. 2020.  NTUNLPL at FinCausal 2020, task 2:improving causality detection using Viterbi decoder.  In _Proceedings of the 1st Joint Workshop on Financial Narrative Processing and MultiLing Financial Summarisation_ , pages 69–73, Barcelona, Spain (Online). COLING.
  * Lample et al. (2016) Guillaume Lample, Miguel Ballesteros, Sandeep Subramanian, Kazuya Kawakami, and Chris Dyer. 2016.  Neural architectures for named entity recognition.  In _Proceedings of the 2016 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies_ , pages 260–270, San Diego, California. Association for Computational Linguistics.
  * Mariko et al. (2020) Dominique Mariko, Hanna Abi-Akl, Estelle Labidurie, Stephane Durfort, Hugues De Mazancourt, and Mahmoud El-Haj. 2020.  The financial document causality detection shared task (FinCausal 2020).  In _Proceedings of the 1st Joint Workshop on Financial Narrative Processing and MultiLing Financial Summarisation_ , pages 23–32, Barcelona, Spain (Online). COLING.
  * Mariko et al. (2021) Dominique Mariko, Hanna Abi Akl, Estelle Labidurie, Hugues de Mazancourt, and Mahmoud El-Haj. 2021.  The Financial Document Causality Detection Shared Task (FinCausal 2021).  In _The Third Financial Narrative Processing Workshop (FNP 2021)_ , Lancaster, UK.
  * Pavlopoulos et al. (2021) John Pavlopoulos, Jeffrey Sorensen, Léo Laugier, and Ion Androutsopoulos. 2021\.  SemEval-2021 task 5: Toxic spans detection.  In _Proceedings of the 15th International Workshop on Semantic Evaluation (SemEval-2021)_ , pages 59–69, Online. Association for Computational Linguistics.
  * Qi et al. (2020) Peng Qi, Yuhao Zhang, Yuhui Zhang, Jason Bolton, and Christopher D. Manning. 2020\.  Stanza: A python natural language processing toolkit for many human languages.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics: System Demonstrations_ , pages 101–108, Online. Association for Computational Linguistics.
  * Ramshaw and Marcus (1995) Lance Ramshaw and Mitch Marcus. 1995.  Text chunking using transformation-based learning.  In _Third Workshop on Very Large Corpora_.
  * Tan et al. (2020) Chuanqi Tan, Wei Qiu, Mosha Chen, Rui Wang, and Fei Huang. 2020.  Boundary enhanced neural span classification for nested named entity recognition.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 34, pages 9016–9023.
  * Viterbi (1967) A. Viterbi. 1967.  Error bounds for convolutional codes and an asymptotically optimum decoding algorithm.  _IEEE Transactions on Information Theory_ , 13(2):260–269.
  * Wolf et al. (2020) Thomas Wolf, Lysandre Debut, Victor Sanh, Julien Chaumond, Clement Delangue, Anthony Moi, Pierric Cistac, Tim Rault, Remi Louf, Morgan Funtowicz, Joe Davison, Sam Shleifer, Patrick von Platen, Clara Ma, Yacine Jernite, Julien Plu, Canwen Xu, Teven Le Scao, Sylvain Gugger, Mariama Drame, Quentin Lhoest, and Alexander Rush. 2020.  Transformers: State-of-the-art natural language processing.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing: System Demonstrations_ , pages 38–45, Online. Association for Computational Linguistics.

##  Appendix A Appendix

###  A.1 Architecture of Proposed Model

An example with $n$ tokens can be represented by the vector of tokens
$\underline{w}=(w_{0},w_{1},...,w_{n})^{\mathsf{T}}$, where $w$ refers to a
BERT input token. A single token in this vector is represented by $w_{i}$,
where $i$ denotes the location index of the token within the example. The
tokenized example also has a matrix of POS features
$V=(\underline{v_{0}},\underline{v_{1}},...,\underline{v_{n}})^{\mathsf{T}}$,
where each $\underline{v_{i}}\in{\rm I\\!R}^{d_{POS}}$ refers to a one-hot
encoding across all POS tags, and so $d_{POS}$ reflects the number of possible
POS tags.

To obtain BERT representations ($r$), we run the tokenized sequence vector
through the BERT encoder ($T$) to obtain $r$ for each token $i$. Thus, we have
that,

$\displaystyle\underline{{r}_{i}}=T_{w}(\underline{w})_{i},$
$\displaystyle\underline{{r}_{i}}\in{\rm I\\!R}^{d_{BERT}}$ (1)

where $d_{BERT}$ refers to the output dimension for BERT encoder.

A dropout layer is represented by $\delta\sim Bernoulli(\rho)$, where $\rho$
refers to the dropout probability. We apply the dropout layer onto the BERT
representations as follows,

$\displaystyle\tilde{\underline{{r}_{i}}}=\delta*\underline{{r}_{i}},$
$\displaystyle{r}_{i}\in{\rm I\\!R}^{d_{BERT}}$ (2)

We combine the BERT and POS representations of each token together by a simple
concatenation along the feature dimension.

$\displaystyle\underline{{r_{2}}_{i}}=[\tilde{\underline{{r}_{i}}},\underline{v_{i}}],$
$\displaystyle\underline{{r_{2}}_{i}}\in{\rm I\\!R}^{d_{BERT}+d_{POS}}$ (3)

Our GNN model generates graph embedding based on these concatenated features,
which are then concatenated together to arrive at our final embeddings.

$\displaystyle\underline{{r_{3}}_{i}}=GNN(\underline{{r_{2}}_{i}}),$
$\displaystyle\underline{{r_{3}}_{i}}\in{\rm I\\!R}^{d_{GNN}}$ (4)
$\displaystyle\underline{{r_{4}}_{i}}=[\underline{{r_{2}}_{i}},\underline{{r_{3}}_{i}}],$
$\displaystyle\underline{{r_{4}}_{i}}\in{\rm I\\!R}^{d}$ (5)

$d_{GNN}$ refers to the output dimension for the last layer in our GNN module
introduced in Section 2.3. That is, if BiLSTM is opted, then $d_{GNN}$ refers
to the size of the output dimension of the BiLSTM layer. If not, it refers to
the output dimension for the second SAGEConv layer. Our final representations
have a feature dimension size of $d=d_{BERT}+d_{POS}+d_{GNN}$.

Next, we run our combined embeddings through a linear layer to obtain
predicted probabilities per token. $c$ refers to the number of classes to be
predicted.

$\displaystyle\underline{{o}_{i}}=\underline{{r_{4}}_{i}}*W+\underline{b_{i}},$
$\displaystyle W\in{\rm I\\!R}^{d\times c},$
$\displaystyle\underline{{o}_{i}},\underline{{b}_{i}}\in{\rm I\\!R}^{c}$ (6)

Cross entropy loss was used during training to optimize model weights. In
evaluation, the logits ran through a Viterbi decoder for adjustment. Finally,
all logits ran through an argmax function to obtain the predicted class that
had the highest probability.

In our implementation, we set $d_{BERT}=768$, $d_{POS}=51$, $d_{GNN}=512$ and
$c=5$.

###  A.2 Replication Checklist

  * •

Hyperparameters: Our pretrained BERT models were initialized with the default
configuration from Huggingface (Wolf et al., 2020). To train our model, we
used the Adam optimizer with $\beta 1=0.9$, $\beta 2=0.999$. Learning rate was
set at $2e-05$ with linear decay. GPU train batch size was set as $4$. Maximum
sequence length was $350$ tokens. For GNN, the graph hidden channel dimensions
(i.e. output dimension of the first SAGEConv layer) was $1024$, the graph
output dimension (i.e. output dimension of the second SAGEConv layer) was
$512$, and the BiLSTM output dimension was also $512$. Probability for all
dropout layers was $0.1$.

  * •

Device: All experiments were ran on the NVIDIA A100-SXM4-40GB GPU.

  * •

Time taken: For 3 folds over 10 epochs each, the Proposed model took us on
average (over the 5 random seeds) $1hour:48minutes:28seconds$ for bert-base-
cased and $1hour:22minutes:24seconds$ for bert-large-cased to train, validate
and predict. For a single run over 10 epochs to generate our submission, the
code took $28minutes:17seconds$ and $20minutes:52seconds$ to train and predict
for the base and large models respectively.

###  A.3 Qualitative Results

In Table 3, we provide examples where the Proposed versus Baseline model
predicts correctly when the other predicts wrongly. The predictions are
obtained from the CV set of the bert-large-cased model with $seed=916$ and the
first fold.

Index |  Baseline |  Proposed |  Right?
---|---|---|---
0036 .000 11 |  <E>Future sales agreements with suppliers increased during the period, and aggregate contracted sales volumes are now 11.7m tonnes per annum</E>, following <C>new European supply agreements.</C> |  <C>Future sales agreements with suppliers increased during the period, and</C> <E>aggregate contracted sales volumes are now 11.7m tonnes per annum</E>, following new European supply agreements. |  Base-line
0270 .000 09 |  <E> It comes with a £250 free overdraft and requires a £1,000 monthly deposit</E> to <C>avoid a £10 monthly fee.</C> |  <C>It comes with a £250 free overdraft</C> and requires a £1,000 monthly deposit to <E>avoid a £10 monthly fee.</E> |  Base-line
0209 .000 33 |  <C>Fiserv believes that this business combination makes sense from the complementary assets between the two companies, projecting higher revenue growth than</C> <E>it would achieve on its own and costs savings of about $900 million over five years.</E> |  <C>Fiserv believes that this business combination makes sense from the complementary assets between the two companies</C>, <E>projecting higher revenue growth than it would achieve on its own and costs savings of about $900 million over five years.</E> |  Prop-osed
0003 .000 19 |  <E>Additionally, the Congress provided $125 million in the current fiscal year for sustainable landscapes programming</E> to <C>prevent forest loss.</C> |  <E>Additionally, the Congress provided $125 million in the current fiscal year</E> for <C>sustainable landscapes programming to prevent forest loss.</C> |  Prop-osed

Table 3: Predicted Cause-Effect spans for CV set from $seed=916$ on first fold
(i.e. $K0$). _Notes._ Cause and Effect spans highlighted in green and orange
respectively.


i|G|j\rangle=t^{|x^{i}\oplus x^{j}|}$ for a suitable choice of probabilities
$q_{w}$. Indeed, the identity $d(x,y)=d(0^{n},x\oplus y)$ implies that a
matrix element $\langle i|G_{\sigma}|j\rangle$ depends only on $x^{i}\oplus
x^{j}$. By the symmetry, matrix elements $\langle i|G^{\prime}|j\rangle$
depend only on the Hamming weight $h=|x^{i}\oplus x^{j}|$. Therefore it
suffices to compute $\langle i|G^{\prime}|j\rangle$ for the special case when
$x^{i}=0^{n}$ is the all-zero string and $x^{j}$ is any fixed bit string with
the Hamming weight $h$, for example, $x^{j}=1^{h}0^{n-h}$. Then

$\langle i|G^{\prime}|j\rangle=\frac{1}{n!}\sum_{\sigma\in
S_{n}}\;\;\sum_{w\geq d(0^{n},\sigma(1^{h}0^{n-h}))}\;q_{w}.$ (113)

By definition of the distance $d(x,y)$ one gets
$d(0^{n},\sigma(1^{h}0^{n-h}))\leq w$ iff $h\leq w$ and
$\sigma_{1},\ldots,\sigma_{h}\leq w$. The number of such permutations $\sigma$
is ${w\choose h}h!(n-h)!$. Exchanging the sums over $\sigma$ and $w$ in Eq.
(113) one gets

$\langle i|G^{\prime}|j\rangle=\frac{1}{n!}\sum_{w=h}^{n}{w\choose
h}h!(n-h)!\,q_{w}.$ (114)

We shall choose $q_{w}$ as a binomial distribution,

$q_{w}={n\choose w}t^{w}(1-t)^{n-w}.$ (115)

Substituting Eq. (115) into Eq. (114) and introducing a variable $p=w-h$ one
gets

$\langle i|G^{\prime}|j\rangle=\sum_{p=0}^{n-h}{n-h\choose
p}t^{p+h}(1-t)^{n-h-p}=t^{h}.$ (116)

By definition, $h=|x^{i}\oplus x^{j}|$, so that $G^{\prime}=G$ as claimed.
Thus $G$ is indeed a probabilisitic mixture of ultrametric matrices and the
lemma is proved. ∎

∎

##  VI Stabilizer fidelity and Stabilizer extent

In the previous Section we established upper bounds on the approximate
stabilizer rank of a state $\psi$ which depend on the the squared $1$-norm
$\|c\|_{1}^{2}$, where

$|\psi\rangle=\sum_{j}c_{j}|\phi_{j}\rangle,$

is a given stabilizer decomposition. Recall that the stabilizer extent
$\xi(\psi)$ denotes the minimum value of $||c||_{1}^{2}$ over all stabilizer
decompositions of $\psi$. We find that $\xi$ is easier to work with than the
approximate stabilizer rank. For any fixed $n$-qubit state $\psi$, $\xi(\psi)$
can be computed using a simple convex optimization program, although the size
of this computation scales poorly with $n$. In this section we develop tools
that allow us to efficiently compute $\xi(\psi)$ whenever $\psi$ is a tensor
product of $1,2$ and $3$ qubit states. In particular, we prove Proposition 1
which establishes that $\xi$ is multiplicative for tensor products of $1$,
$2$, and $3$-qubit states.

In subsection VI.1 we use standard convex duality to give a characterization
of $\xi$ in terms of the stabilizer fidelity, defined as the maximum overlap
with respect to the set of stabilizer states

$F(\psi):=\mathrm{max}_{\phi\in\mathrm{STAB}_{n}}|\langle\psi|\phi\rangle|^{2}.$
(117)

As a consequence, multiplicativity of $\xi$ is directly related to
multiplicativity of the stabilizer fidelity. In subsection VI.2 we give
sufficient and necessary conditions for multiplicativity of the stabilizer
fidelity. In particular, we define the class of stabilizer-aligned states for
which multiplicativity holds. In subsection VI.3 we investigate the class of
stabilizer-aligned states and prove that all tensor products of $1,2$ and $3$
qubit states are stabilizer-aligned. Finally, in section VI.4 we use these
results to prove Proposition 1.

###  VI.1 Convex duality

Here we show that the optimization of $\xi(\psi)$ can be recast as a dual
convex problem and we prove the following:

######  Theorem 4.

For any $n$-qubit state $\psi$ we have

$\xi(\psi)=\max_{\omega}\frac{|\langle\psi|\omega\rangle|^{2}}{F(\omega)},$
(118)

where the maximum is over all $n$-qubit states $\omega$.

Thus any $n$-qubit state $\omega$ can act as a witness to provide a lower
bound on $\xi$ and, furthermore, there exists at least one optimal witness
state $\omega_{\star}$ which achieves the maximum in Eq. (118). For example,
choosing $\omega=\psi$, we get the lower bound

$\xi(\psi)\geq\frac{1}{F(\psi)}.$ (119)

For Clifford magic states this lower bound is tight as stated in Proposition
2.

###### Proof.

We shall map the problem into the language of convex optimization and use
standard results in that field Boyd and Vandenberghe (2004). Using the
computation basis $\\{|{x}\rangle\\}$ we can decompose any stabilizer state
$|\psi_{j}\rangle=\sum_{{x}}M_{{x},j}|{x}\rangle$. Given a state
$|\psi\rangle=\sum_{{x}}a_{{x}}|{x}\rangle$, the primal optimization problem
can be written as

$\displaystyle\sqrt{\xi(\psi)}$
$\displaystyle=\mathrm{min}_{{c}}f({c})=||{c}||_{1}$ (120) such that
$\displaystyle M{c}-{a}=0$ (121)

This is clearly a convex optimization problem with affine constraints. Because
the coefficient in ${c}$ are complex, rather than real, this is a second order
cone problem Boyd and Vandenberghe (2004). For any convex optimization problem
there exists a dual function

$\displaystyle g(\nu)$
$\displaystyle=\mathrm{inf_{{c}}}\left(||{c}||_{1}+\nu^{T}(M{c}-{a})\right)$
(122) $\displaystyle=\begin{cases}-\nu^{T}{a}&\mbox{ when
}||M^{T}\nu||_{\infty}\leq 1\\\ -\infty&\mbox{ otherwise }\end{cases}$ (123)

where for any value of the dual variables $\nu$ we have
$g(\nu)\leq\sqrt{\xi(\psi)}$. The dual optimization problem is the
maximisation of $g(\nu)$ over $\nu$ to obtain the best lower bound on
$\sqrt{\xi(\psi)}$. We can discount the need for two cases by adding
$||M^{T}\nu||_{\infty}\leq 1$ as a constraint, to obtain the problem

$\displaystyle d^{\star}(\psi)$ $\displaystyle=\mathrm{max}_{\nu}-\nu\cdot{a}$
(124) such that $\displaystyle||M^{T}\nu||_{\infty}\leq 1,$

or more simply

$\displaystyle d^{\star}(\psi)$
$\displaystyle=\mathrm{max}_{\nu}\frac{-\nu\cdot{a}}{||M^{T}\nu||_{\infty}}.$
(125)

Because the primal problem has affine constraints, we have strong duality and
there must exist a $\nu_{\star}$ such that
$g(\nu_{\star})=-\nu_{\star}^{T}{a}=\sqrt{\xi(\psi)}$. Next, we restate this
dual problem in terms of quantum states. For every $\nu$ we can associate a
normalised quantum state

$|\omega_{\nu}\rangle:=\frac{1}{||\nu||_{2}}\sum_{{x}}(-\nu^{*}_{{x}})|{x}\rangle,$
(126)

so that

$\langle\omega_{\nu}|\psi\rangle=\frac{-\nu\cdot{a}}{||\nu||_{2}}.$ (127)

Next we note that

$||M^{T}\nu||_{\infty}=\frac{\mathrm{Max}_{|\phi\rangle\in\mathrm{STAB}}|\langle\omega_{\nu}|\phi\rangle|}{||\nu||_{2}}=\frac{\sqrt{F(\omega_{\nu})}}{||\nu||_{2}}$
(128)

Therefore, the dual problem can also be stated as

$\displaystyle d^{\star}(\psi)$
$\displaystyle=\mathrm{max}_{|\omega_{\nu}\rangle}\frac{\langle\omega_{\nu}|\psi\rangle}{\sqrt{F(\omega_{\nu})}},$
(129)

where the factors $||\nu||_{2}$ have cancelled out. The optimal $\nu_{\star}$
gives the optimal $|\omega_{\star}\rangle$, which completes the proof. ∎

###  VI.2 Stabilizer alignment

Combining Theorems 2 and 1 we get an upper bound
$\chi_{\delta}(\psi)\leq\delta^{-2}F(\psi)^{-1}$ on the approximate stabilizer
rank of any Clifford magic state $\psi$. We shall be interested in the case
when $\psi$ is a tensor product of a large number of few-qubit magic states
such as $T$-type or CCZ-type states. For example, the case $\psi=CCZ^{\otimes
m}$ is relevant to gadget-based simulation of quantum circuits composed of
Clifford gates and $m$ CCZ gates. This motivates the question of whether the
stabilizer fidelity $F(\psi)$ is multiplicative under tensor product, i.e.

$F(\psi\otimes\phi)\stackrel{{\scriptstyle?}}{{=}}F(\psi)F(\phi).$ (130)

Note that $F(\psi\otimes\phi)\geq F(\psi)F(\phi)$ since the set of stabilizer
states is closed under tensor product.

Below we define a set of quantum states $\mathcal{S}$ such that
$F(\phi\otimes\psi)=F(\phi)F(\psi)$ whenever $\phi,\psi\in\mathcal{S}$.
Remarkably, this set is also closed under tensor product, that is
$\phi\otimes\psi\in\mathcal{S}$ whenever $\phi,\psi\in\mathcal{S}$. Moreover,
we show that the stabilizer fidelity is not multiplicative for all states
$\phi\notin{\cal S}$. More precisely, for any $\phi\notin{\cal S}$ there
exists a state $\psi$ such that $F(\phi\otimes\psi)>F(\phi)F(\psi)$. In that
sense, our results provide necessary and sufficient conditions under which the
stabilizer fidelity is multiplicative under tensor product.

To state our results let us generalize the definition of stabilizer fidelity
as follows. For each $n\geq 1$ and $0\leq m\leq n$ define a set $S_{n,m}$
which consists of all stabilizer projectors $\Pi$ on $n$ qubits satisfying
$\mathrm{Tr}[\Pi]=2^{m}$.

######  Definition 9.

For any $n$-qubit state $|\phi\rangle$ define

$F_{m}(\phi)=2^{-m/2}\max_{\Pi\in
S_{n,m}}\langle\phi|\Pi|\phi\rangle.\qquad\qquad m=0,\ldots,n.$

Let us say that $\phi$ is stabilizer-aligned if $F_{m}(\phi)\leq F_{0}(\phi)$
for all $m$.

Note that in the above $F_{0}=F$ is the stabilizer fidelity. Here we
investigate the consequences of stabilizer-alignment. Whether or not a given
state is stabilizer-aligned is discussed in the following subsection.

######  Theorem 5.

Suppose $\phi$ and $\psi$ are stabilizer-aligned. Then $\phi\otimes\psi$ is
stabilizer-aligned and

$F(\phi\otimes\psi)=F(\phi)F(\psi).$

Conversely, suppose $\phi$ is not stabilizer-aligned. Let $\phi^{\star}$ be
the complex conjugate of $\phi$. Then

$F(\phi\otimes\phi^{\star})>F(\phi)F(\phi^{\star}).$

The theorem implies that the stabilizer fidelity is multiplicative for any
stabilizer-aligned states:

######  Corollary 1.

Suppose $\psi_{1},\ldots,\psi_{L}$ are stabilizer-aligned quantum states. Then

$F(\psi_{1}\otimes\psi_{2}\otimes\ldots\otimes\psi_{L})=\prod_{j=1}^{L}F(\psi_{j}).$

We prove Theorem 5 using characterization of entanglement in tripartite
stabilizer states from Ref. Bravyi  _et al._ (2006):

######  Lemma 9 (Bravyi  _et al._ (2006)).

Any pure tripartite stabilizer state can be transformed by local unitary
Clifford operators to a tensor product of states from the set
$\\{|0\rangle,|\Psi^{+}\rangle,|\Psi^{+}_{3}\rangle\\}$ where

$|\Psi^{+}\rangle=\frac{1}{\sqrt{2}}(|00\rangle+|11\rangle)\qquad\qquad|\Psi^{+}_{3}\rangle=\frac{1}{\sqrt{2}}\left(|000\rangle+|111\rangle\right).$

######  Corollary 2 (Bravyi  _et al._ (2006)).

Suppose $\Pi$ be a stabilizer projector describing a bipartite system $AB$.
Then there exists a unitary Clifford operator $U=U_{A}\otimes U_{B}$ and
integers $a,b,c,d\geq 0$ such that

$U\Pi
U^{-1}=\sum_{\alpha=1}^{2^{a}}\sum_{\beta=1}^{2^{b}}\sum_{\gamma=1}^{2^{c}}|\omega_{\alpha\beta\gamma}\rangle\langle\omega_{\alpha\beta\gamma}|,$
(131)

where

$|\omega_{\alpha\beta\gamma}\rangle=2^{-d/2}\sum_{\delta=1}^{2^{d}}|\alpha,\gamma,\delta\rangle\otimes|\beta,\gamma,\delta\rangle.$
(132)

Here $|\alpha,\gamma,\delta\rangle$ and $|\beta,\gamma,\delta\rangle$ are the
computational basis vectors of $A$ and $B$.

###### Proof.

Let us apply Lemma 9 to a tripartite stabilizer state

$|\Psi\rangle=(\Pi\otimes
I)2^{-n/2}\sum_{z\in\\{0,1\\}^{n}}\;|z\rangle_{AB}\otimes|z\rangle_{C},$

where $n=|A|+|B|$ and $C$ is a system of $n$ qubits. The lemma implies that
$\Pi$ is equivalent modulo local Clifford operators to a tensor product of
local stabilizer projectors $|0\rangle\langle 0|$ and $I=|0\rangle\langle
0|+|1\rangle\langle 1|$ as well as bipartite projectors $|00\rangle\langle
00|+|11\rangle\langle 11|$ and $|\Psi^{+}\rangle\langle\Psi^{+}|$ shared
between $A$ and $B$. Let $a$ and $b$ be the number of times $\Pi$ contains the
identity factor on $A$ and $B$ respectively. Let $c$ be the number of times
$\Pi$ contains the projector $|00\rangle\langle 00|+|11\rangle\langle 11|$
shared between $A$ and $B$. Let $d$ be the number of times $\Pi$ contains the
EPR projector $|\Psi^{+}\rangle\langle\Psi^{+}|$. The desired family of states
$\omega_{\alpha\beta\gamma}$ is then obtained by writing each projector $I$
and $|00\rangle\langle 00|+|11\rangle\langle 11|$ as a sum of rank-$1$
projectors onto the computational basis vectors. ∎

###### Proof of Theorem 5.

To prove the first two claims of the theorem it suffices to show that

$F_{m}(\phi\otimes\psi)\leq F_{0}(\phi)F_{0}(\psi).$ (133)

for all $m$. Indeed, combining Eq. (133) and the obvious bound
$F_{0}(\phi)F_{0}(\psi)\leq F_{0}(\phi\otimes\psi)$ shows that
$F_{m}(\phi\otimes\psi)\leq F_{0}(\phi\otimes\psi)$, that is,
$\phi\otimes\psi$ is stabilizer-aligned. Using Eq. (133) for $m=0$ gives
multiplicativity of the stabilizer fidelity
$F_{0}(\phi\otimes\psi)=F_{0}(\phi)F_{0}(\psi)$.

Define a bipartite system $AB$ such that $\phi$ and $\psi$ are states of $A$
and $B$. Let $\Pi$ be a stabilizer projector of rank $2^{m}$ such that

$F_{m}(\phi\otimes\psi)=2^{-m/2}\langle\phi\otimes\psi|\Pi|\phi\otimes\psi\rangle.$

We shall write $\Pi$ as a sum of rank-$1$ stabilizer projectors as stated in
Corollary 2. Since local Clifford unitary operators do not change the
stabilizer fidelity, we shall absorb the unitaries $U_{A}$ and $U_{B}$ into
the states $\phi$ and $\psi$ respectively. Accordingly, below we set $U=I$.
Consider a single term $\omega_{\alpha\beta\gamma}$ in the decomposition of
$\Pi$. Applying the Cauchy-Schwarz inequality one gets

$|\langle\phi\otimes\psi|\omega_{\alpha\beta\gamma}\rangle|^{2}=2^{-d}\left|\sum_{\delta=1}^{2^{d}}\langle\phi|\alpha,\gamma,\delta\rangle\cdot\langle\psi|\beta,\gamma,\delta\rangle\right|^{2}\leq
2^{-d}\langle\phi|\Pi^{A}_{\alpha\gamma}|\phi\rangle\cdot\langle\psi|\Pi^{B}_{\beta\gamma}|\psi\rangle,$
(134)

where we defined stabilizer projectors

$\Pi^{A}_{\alpha,\gamma}=\sum_{\delta=1}^{2^{d}}|\alpha,\gamma,\delta\rangle\langle\alpha,\gamma,\delta|\quad\mbox{and}\quad\Pi^{B}_{\beta,\gamma}=\sum_{\delta=1}^{2^{d}}|\beta,\gamma,\delta\rangle\langle\beta,\gamma,\delta|.$
(135)

By assumption, $\psi$ is stabilizer-aligned. Thus

$\max_{\gamma}\langle\psi|\sum_{\beta=1}^{2^{b}}\Pi^{B}_{\beta\gamma}|\psi\rangle\leq
2^{(b+d)/2}F_{0}(\psi).$ (136)

Here we noted that $\sum_{\beta=1}^{2^{b}}\Pi^{B}_{\beta\gamma}$ is a
projector of rank $2^{b+d}$ for all $\gamma$. Combining Eq. (134,136) gives

$\langle\phi\otimes\psi|\Pi|\phi\otimes\psi\rangle=\sum_{\alpha=1}^{2^{a}}\sum_{\beta=1}^{2^{b}}\sum_{\gamma=1}^{2^{c}}|\langle\phi\otimes\psi|\omega_{\alpha\beta\gamma}\rangle|^{2}\leq
2^{(b-d)/2}F_{0}(\psi)\cdot\langle\phi|\sum_{\alpha=1}^{2^{a}}\sum_{\gamma=1}^{2^{c}}\Pi^{A}_{\alpha,\gamma}|\phi\rangle$
(137)

The assumption that $\phi$ is stabilizer-aligned gives

$\langle\phi|\sum_{\alpha=1}^{2^{a}}\sum_{\gamma=1}^{2^{c}}\Pi^{A}_{\alpha,\gamma}|\phi\rangle\leq
2^{(a+c+d)/2}F_{0}(\phi).$ (138)

Here we noted that
$\sum_{\alpha=1}^{2^{a}}\sum_{\gamma=1}^{2^{c}}\Pi^{A}_{\alpha,\gamma}$ is a
projector of rank $2^{a+c+d}$. Combining Eqs. (137,138) gives

$\langle\phi\otimes\psi|\Pi|\phi\otimes\psi\rangle\leq
2^{(a+b+c)/2}F_{0}(\psi)F_{0}(\phi).$

It remains to notice that $\Pi$ has rank $2^{m}$, where $m=a+b+c$. This
establishes Eq. (133).

We now prove the converse statement from Theorem 5.

######  Lemma 10.

Let $\phi$ be an $n$-qubit state which is not stabilizer-aligned. Then

$F_{0}(\phi\otimes\phi^{\star})>F_{0}(\phi)F_{0}(\phi^{\star}).$

###### Proof.

If $\phi$ is not stabilizer-aligned then we have $F_{m}(\phi)>F_{0}(\phi)$ for
some $m\in\\{1,\ldots,n\\}$. Let $\Pi$ be a stabilizer projector with

$F_{m}(\phi)=\frac{1}{\sqrt{2}^{m}}\langle\phi|\Pi|\phi\rangle.$

Let $C$ be an $n$-qubit Clifford such that

$\Pi=C\left(|0\rangle\langle 0|_{n-m}\otimes I_{m}\right)C^{\dagger}.$

Next consider a system of $2n$ qubits and partition them as
$[2n]=ABA^{\prime}B^{\prime}$ where $|A|=|A^{\prime}|=n-m$ and
$|B|=|B^{\prime}|=m$. Define a $2n$-qubit stabilizer state

$|\theta\rangle=C\otimes\alpha|0\rangle_{A}|\Phi\rangle_{BB^{\prime}}|0\rangle_{A^{\prime}},$

where

$|\Phi\rangle_{BB^{\prime}}=\frac{1}{\sqrt{2}^{m}}\sum_{z\in\\{0,1\\}^{m}}|z\rangle_{B}|z\rangle_{B^{\prime}}.$

Also define a normalized $m$-qubit state

$|\omega\rangle=\frac{1}{2^{m/4}\sqrt{F_{m}(\phi)}}\left(\langle
0|_{n-m}\otimes I_{m}\right)C|\phi\rangle.$ $\displaystyle
F_{0}(\phi\otimes\phi^{\star})$
$\displaystyle\geq\langle\phi\otimes\phi^{\star}|\theta\rangle\langle\theta|\phi\otimes\phi^{\star}\rangle$
(139)
$\displaystyle=\langle\omega\otimes\omega^{\star}|\Phi\rangle\langle\Phi|\omega\otimes\omega^{\star}\rangle
2^{m}(F_{m}(\phi))^{2}$ (140) $\displaystyle=(F_{m}(\phi))^{2}$ (141)
$\displaystyle>F_{0}(\phi)F_{0}(\phi^{\star}).$ (142)

where in the last line we used the fact that
$F_{m}(\phi)>F_{0}(\phi)=F_{0}(\phi^{\star})$. ∎

∎

###  VI.3 Proving and disproving stabilizer alignment

In this section we prove that all states of $n\leq 3$ qubits are stabilizer-
aligned. We also show that typical $n$-qubit states are not stabilizer-aligned
for sufficiently large $n$. An important lemma is the following

######  Lemma 11.

For any quantum state $\psi$ we have $F_{m}(\psi)\leq F_{0}(\psi)$ for
$m=1,2,3$.

It follows immediately that

######  Corollary 3.

All states of $n\leq 3$ qubits are stabilizer-aligned.

Indeed, if we consider $n$-qubit states, it suffices to check that
$F_{m}(\psi)\leq F_{0}(\psi)$ for $m\leq n$.

######  Corollary 4.

If $F_{0}(\psi)\geq 1/4$ then $\psi$ is stabilizer-aligned.

Indeed, if $m\geq 4$ then $F_{m}(\psi)\leq 2^{-m/2}\leq 1/4\leq F_{0}(\psi)$.

Finally, we show that Haar-random $n$-qubit states are not stabilizer-aligned
for sufficiently large $n$.

######  Claim 1.

Let $\psi$ be a Haar-random $n$-qubit state. Then

$\mathrm{Pr}[F_{0}(\psi\otimes\psi^{\star})\neq
F_{0}(\psi)F_{0}(\psi^{\star})]\geq 1-o(1).$

and so for large enough $n$ a typical state $\psi$ is not stabilizer-aligned.

Highly structured states on a large number of qubits may be stabilizer-
aligned, and for instance it is an open question whether or not all Clifford
magic states are stabilizer-aligned.

###### Proof of Lemma 11.

First, we claim that

$F_{m-1}(\psi)\geq
2^{-1/2}\left(1+\left[\frac{2^{m}-1}{4^{m}-1}\right]^{1/2}\right)\cdot
F_{m}(\psi)$ (143)

for all $m\geq 1$. Indeed, consider a fixed $m$ and a rank-$2^{m}$ stabilizer
projector $\Pi\in S_{n,m}$ such that
$F_{m}(\psi)=2^{-m/2}\langle\psi|\Pi|\psi\rangle$. Using the standard
stabilizer formalism one can show that

$U\Pi U^{-1}=I^{\otimes m}\otimes|0\rangle\langle
0|^{\otimes(n-m)}\equiv\Pi^{\prime}$

for some $n$-qubit unitary Clifford operator $U$. Define a state
$|\psi^{\prime}\rangle=U|\psi\rangle$. We have

$\Pi^{\prime}|\psi^{\prime}\rangle=\Gamma^{1/2}|\omega\rangle\otimes|0^{n-m}\rangle$

for some $m$-qubit normalized state $|\omega\rangle$ and
$\Gamma=\langle\psi^{\prime}|\Pi^{\prime}|\psi^{\prime}\rangle=\langle\psi|\Pi|\psi\rangle$.
Since $\omega$ is normalized,

$\sum_{P\neq I}\langle\omega|P|\omega\rangle^{2}=2^{m}-1,$

where the sum runs over all $4^{m}-1$ non-trivial Pauli operators on $m$
qubits. Thus there exists an $m$-qubit Pauli operator $P\neq I$ such that

$\langle\omega|P|\omega\rangle\geq\left(\frac{2^{m}-1}{4^{m}-1}\right)^{1/2}.$
(144)

Define a stabilizer projector

$\Pi^{\prime\prime}=\frac{1}{2}(I+P)\otimes|0\rangle\langle
0|^{\otimes(n-m)}\in S_{n,m-1}.$

Recalling that $\Gamma=\langle\psi|\Pi|\psi\rangle=2^{m/2}F_{m}(\psi)$ we
arrive at

$F_{m-1}(\psi)=F_{m-1}(\psi^{\prime})\geq
2^{-(m-1)/2}\langle\psi^{\prime}|\Pi^{\prime\prime}|\psi^{\prime}\rangle=2^{-(m-1)/2}\frac{\Gamma}{2}(1+\langle\omega|P|\omega\rangle)=2^{-1/2}(1+\langle\omega|P|\omega\rangle)\cdot
F_{m}(\psi).$

Combining this identity and Eq. (144) proves Eq. (143). Applying Eq. (143)
inductively gives

$F_{0}(\psi)\geq 2^{-1/2}(1+\sqrt{1/3})\cdot F_{1}(\psi)\approx 1.115\cdot
F_{1}(\psi),$ (145) $F_{0}(\psi)\geq 2^{-1/2}(1+\sqrt{1/3})\cdot
2^{-1/2}(1+\sqrt{3/15})\cdot F_{2}(\psi)\approx 1.141\cdot F_{2}(\psi),$ (146)
$F_{0}(\psi)\geq 2^{-1/2}(1+\sqrt{1/3})\cdot 2^{-1/2}(1+\sqrt{3/15})\cdot
2^{-1/2}(1+\sqrt{7/63})\cdot F_{3}(\psi)\approx 1.076\cdot F_{3}(\psi).$ (147)

Thus $F_{0}(\psi)\geq F_{m}(\psi)$ for $m=1,2,3$ proving the lemma. ∎

Next, we prove claim 1.

###### Proof.

Let $w$ be any $n$-qubit state. For Haar-random $\psi$ the probability density
function $p(y)$ of $y=|\langle w|\psi\rangle|^{2}$ does not depend on $w$ and
is equal to (equation (9) of Ref. Zyczkowski and Sommers (2000)),

$p(y)=(2^{n}-1)(1-y)^{2^{n}-2}.$

Integrating this we obtain the cumulative distribution function

$\mathrm{Pr}\left[|\langle w|\psi\rangle|^{2}\geq
x\right]=(1-x)^{2^{n}-1}\leq\exp(-x(2^{n}-1)).$

Since an $n$-qubit stabilizer state is specified by $O(n^{2})$ bits the
cardinality of the set $\mathrm{STAB}_{n}$ of $n$-qubit stabilizer states is
$|\mathrm{STAB}_{n}|\leq 2^{O(n^{2})}$. Choosing $x=n^{3}/2^{n}$ and applying
a union bound we get

$\mathrm{Pr}\left[\left(\max_{w\in\mathrm{STAB}_{n}}|\langle\psi|w\rangle|^{2}\right)\geq
n^{3}/2^{n}\right]\leq e^{-\Omega(n^{3})}.$

This says that with probability very close to 1 a random $\psi$ has
$F_{0}(\psi)=F_{0}(\psi^{\star})\leq n^{3}/2^{n}$. Next suppose $\psi$ has
this property. Then

$F_{0}(\psi\otimes\psi^{\star})\geq\left|\frac{1}{\sqrt{2}^{n}}\sum_{z\in\\{0,1\\}^{n}}\langle
z|\psi\rangle\langle z|\psi^{\star}\rangle\right|^{2}=\frac{1}{2^{n}},$

which is strictly greater than $F_{0}(\psi)F_{0}(\psi^{\star})\leq
2^{-2n}(n^{3})^{2}$. ∎

###  VI.4 Multiplicativity of stabilizer extent

This subsection considers tensor products of few-qubit states that involve at
most three qubits each and shows that $\xi$ behaves multiplicatively for such
products, proving Proposition 4. The proof will draw heavily on Theorem 4 and
Corollary 3.

###### Proof of Proposition 4.

By Theorem 4 there exist witness states
$\\{\omega_{\star,1},\omega_{\star,2},\ldots,\omega_{\star,L}\\}$ such that

$\frac{|\langle\psi_{j}|\omega_{\star,j}\rangle|^{2}}{F(\omega_{\star,j})}=\xi(\psi_{j}).$
(148)

We consider the product witness
$|\Omega\rangle=\bigotimes_{j}|\omega_{\star,j}\rangle$ for which

$|\langle\Psi|\Omega\rangle|^{2}=\prod_{j}|\langle\psi_{j}|\omega_{\star,j}\rangle|^{2}.$
(149)

Furthermore, using Corollary 3 and Theorem 5 we get

$F(\Omega)=\prod_{j}F(\omega_{\star,j}).$ (150)

Putting this together yields

$\frac{|\langle\Psi|\Omega\rangle|^{2}}{F(\Omega)}=\prod_{j}\frac{|\langle\psi_{j}|\omega_{\star,j}\rangle|^{2}}{F(\omega_{\star,j})}=\prod_{j=1}^{L}\xi(\psi_{j}).$
(151)

Thus, using $\Omega$ as a witness, we get

$\prod_{j=1}^{L}\xi(\psi_{j})\leq\xi(\Psi).$ (152)

Furthermore, $\xi$ is inherently sub-multiplicative and so we must have
equality. ∎

Now let us see how this can be used to bound the approximate stabilizer rank
of a product state $\alpha^{\otimes n}$ where $\alpha$ is a single-qubit
state. Combining Theorem 4 with Lemma 6 we get

$\chi_{\delta}(\alpha^{\otimes n})\leq\delta^{-1}\xi(\alpha^{\otimes
n})=\delta^{-2}(\xi(\alpha))^{n}.$ (153)

Note that since $\alpha$ is a single-qubit state we can easily compute
$\xi(\alpha)$ by solving a small convex optimization program. In Figure 7 we
plot $\xi(\alpha)$ as a function of the single-qubit state $\alpha$ on the
first octant of the Bloch sphere.

Figure 7: The color indicates the value of $\xi$ for single-qubit states in
the first octant of the Bloch sphere. This function controls the upper bound
on the approximate stabilizer rank as in Eq. (153).

The maximum value plotted in Figure 7 is $\xi(f)=2/(1+1/\sqrt{3})\approx
1.2679$, which is achieved by the so-called face state $|f\rangle$ which lies
in the center of the surface and is defined by

$|f\rangle\langle f|=\frac{1}{2}\left(I+\frac{1}{\sqrt{3}}(X+Y+Z)\right).$

The single-qubit states in Figure 7 which lie in the $x$-$z$ plane are of the
form

$|\theta\rangle=\cos(\theta/2)|0\rangle+\sin(\theta/2)|1\rangle=\left(\cos(\theta/2)-\sin(\theta/2)\right)|0\rangle+\sqrt{2}\sin(\theta/2)|+\rangle$
(154)

for $\theta\in[0,\pi/2]$. In this case, the stabilizer decomposition on the
right hand side achieves the optimal value of $\xi$. We can use this example
to show that in the general case the upper bound on approximate stabilizer
rank given in Theorem 1 is not tight (for $\delta=O(1)$, say). When $\theta$
is close to $0$ it becomes advantageous to expand $\theta^{\otimes n}$ in the
standard $0,1$ basis and truncate amplitudes which are very small. Using this
approach one obtains an approximate stabilizer rank scaling as
$2^{h_{2}(\cos^{2}(\theta/2))}$ where $h_{2}$ is the binary entropy. In Figure
8 we compare the performance of these upper bounds as a function of $\theta$.

Figure 8: The approximate stabilizer rank of $|\theta^{\otimes n}\rangle$ is
upper bounded as $\chi_{\delta}(\theta^{\otimes
n})\leq\delta^{-2}\xi(\theta)^{n}$, where
$\xi(\theta)=(\cos(\theta/2)+\tan(\pi/8)\sin(\theta/2))^{2}$ is attained by
the stabilizer decomposition from Eq. (154). The red line shows the function
$\xi(\theta)$ for $\theta\in[0,\pi/4]$ and the blue line shows the function
$g(\theta)=2^{h_{2}(\cos^{2}(\theta/2))}$ where $h_{2}$ is the binary entropy.
Our upper bound on the approximate stabilizer rank of $\theta^{\otimes n}$
performs better that obtained by a naive expansion in the $0,1$ basis whenever
the red line lies below the blue line.

##  VII Acknowledgements

EC and MH are supported by the EPSRC (Grant No. EP/M024261/1). PC is supported
by the EPSRC (Grant No. EP/L015242/1). The collaboration benefited from
support by the NQIT project partnership fund (Grant No. EP/M013243/1), an
EPSRC IIKE award and the IBM Research Frontiers Institute.

## References

  * Smelyanskiy  _et al._ (2016) M. Smelyanskiy, N. P. Sawaya, and A. Aspuru-Guzik, arXiv preprint arXiv:1601.07195 (2016).
  * Häner and Steiger (2017) T. Häner and D. S. Steiger, in  _Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_ (ACM, 2017) p. 33.
  * Pednault  _et al._ (2017) E. Pednault, J. A. Gunnels, G. Nannicini, L. Horesh, T. Magerlein, E. Solomonik, and R. Wisnieff, arXiv preprint arXiv:1710.05867 (2017).
  * Chen  _et al._ (2018) J. Chen, F. Zhang, M. Chen, C. Huang, M. Newman, and Y. Shi, arXiv preprint arXiv:1805.01450 (2018).
  * Gottesman (1998) D. Gottesman, Physical Review A 57, 127 (1998).
  * Aaronson and Gottesman (2004) S. Aaronson and D. Gottesman, Physical Review A 70, 052328 (2004).
  * Bravyi and Kitaev (2005) S. Bravyi and A. Kitaev, Physical Review A 71, 022316 (2005).
  * Garcia  _et al._ (2012) H. J. Garcia, I. L. Markov, and A. W. Cross, arXiv preprint arXiv:1210.6646 (2012).
  * García  _et al._ (2014) H. J. García, I. L. Markov, and A. W. Cross, Quantum Information & Computation 14, 683 (2014).
  * Bravyi  _et al._ (2016) S. Bravyi, G. Smith, and J. A. Smolin, Physical Review X 6, 021043 (2016).
  * Bravyi and Gosset (2016) S. Bravyi and D. Gosset, Physical Review Letters 116, 250501 (2016).
  * Pashayan  _et al._ (2015) H. Pashayan, J. J. Wallman, and S. D. Bartlett, Physical Review Letters 115, 070501 (2015).
  * Delfosse  _et al._ (2015) N. Delfosse, P. Allard Guerin, J. Bian, and R. Raussendorf, Physical Review X 5, 021003 (2015).
  * Howard and Campbell (2017) M. Howard and E. Campbell, Physical Review Letters 118, 090501 (2017).
  * Bennink  _et al._ (2017) R. S. Bennink, E. M. Ferragut, T. S. Humble, J. A. Laska, J. J. Nutaro, M. G. Pleszkoch, and R. C. Pooser, Physical Review A 95, 062337 (2017).
  * Martínez  _et al._ (1994) S. Martínez, G. Michon, and J. S. Martín, SIAM Journal on Matrix Analysis and Applications 15, 98 (1994).
  * Nabben and Varga (1994) R. Nabben and R. S. Varga, SIAM Journal on Matrix Analysis and Applications 15, 107 (1994).
  * Preskill (2018) J. Preskill, arXiv preprint arXiv:1801.00862 (2018).
  * Aharonov  _et al._ (2017) D. Aharonov, M. Ben-Or, E. Eban, and U. Mahadev, arXiv preprint arXiv:1704.04487 (2017).
  * Morimae and Fitzsimons (2018) T. Morimae and J. F. Fitzsimons, Physical Review Letters 120, 040501 (2018).
  * Jozsa and Strelchuk (2017) R. Jozsa and S. Strelchuk, arXiv preprint arXiv:1705.02817 (2017).
  * McKay  _et al._ (2017) D. C. McKay, C. J. Wood, S. Sheldon, J. M. Chow, and J. M. Gambetta, Physical Review A 96, 022330 (2017).
  * Nielsen and Chuang (2002) M. A. Nielsen and I. Chuang, “Quantum computation and quantum information,” (2002).
  * Bravyi  _et al._ (2006) S. Bravyi, D. Fattal, and D. Gottesman, Journal of Mathematical Physics 47, 062106 (2006).
  * Gottesman and Chuang (1999) D. Gottesman and I. L. Chuang, Nature 402, 390 (1999).
  * Note (1) By simulating a postselective quantum circuit one could solve 3-SAT using a polynomial number of T-gates, see e.g., Ref. Huang  _et al._ (2018).
  * Webb (2016) Z. Webb, Quantum Information & Computaion 16, 1379 (2016).
  * Kueng and Gross (2015) R. Kueng and D. Gross, arXiv preprint arXiv:1510.02767 (2015).
  * Farhi  _et al._ (2014) E. Farhi, J. Goldstone, and S. Gutmann, arXiv preprint arXiv:1412.6062 (2014).
  * Fried  _et al._ (2017) E. S. Fried, N. P. Sawaya, Y. Cao, I. D. Kivlichan, J. Romero, and A. Aspuru-Guzik, arXiv preprint arXiv:1709.03636 (2017).
  * Rötteler (2010) M. Rötteler, in  _Proceedings of the 21st ACM-SIAM Symposium on Discrete Algorithms_ (2010) pp. 448–457.
  * Van den Nest (2011) M. Van den Nest, Quantum Information & Computation 11, 784 (2011).
  * Wolff  _et al._ (2004) U. Wolff, Alpha Collaboration, _et al._ , Computer Physics Communications 156, 143 (2004).
  * van Dam  _et al._ (2006) W. van Dam, S. Hallgren, and L. Ip, SIAM Journal on Computing 36, 763 (2006).
  * Rötteler (2010) M. Rötteler, in  _Proceedings of the Twenty-first Annual ACM-SIAM Symposium on Discrete Algorithms_ , SODA ’10 (Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2010) pp. 448–457.
  * De Raedt  _et al._ (2007) K. De Raedt, K. Michielsen, H. De Raedt, B. Trieu, G. Arnold, M. Richter, T. Lippert, H. Watanabe, and N. Ito, Computer Physics Communications 176, 121 (2007).
  * Markov and Shi (2008) I. L. Markov and Y. Shi, SIAM Journal on Computing 38, 963 (2008).
  * Boixo  _et al._ (2017) S. Boixo, S. V. Isakov, V. N. Smelyanskiy, and H. Neven, arXiv preprint arXiv:1712.05384 (2017).
  * Aaronson and Chen (2016) S. Aaronson and L. Chen, arXiv preprint arXiv:1612.05903 (2016).
  * Li  _et al._ (2018) R. Li, B. Wu, M. Ying, X. Sun, and G. Yang, arXiv preprint arXiv:1804.04797 (2018).
  * Boixo  _et al._ (2018) S. Boixo, S. V. Isakov, V. N. Smelyanskiy, R. Babbush, N. Ding, Z. Jiang, M. J. Bremner, J. M. Martinis, and H. Neven, Nature Physics 14, 595 (2018).
  * Huang  _et al._ (2018) C. Huang, M. Newman, and M. Szegedy, arXiv preprint arXiv:1804.10368 (2018).
  * Fowler  _et al._ (2013) A. G. Fowler, S. J. Devitt, and C. Jones, Scientific reports 3, 1939 (2013).
  * Jones (2013) C. Jones, Physical Review A 87, 022328 (2013).
  * Alon (1990) N. Alon, Graphs and Combinatorics 6, 1 (1990).
  * Van Den Nest (2010) M. Van Den Nest, Quantum Information & Computation 10, 258 (2010).
  * Maslov and Roetteler (2017) D. Maslov and M. Roetteler, arXiv preprint arXiv:1705.09176 (2017).
  * Eldar and Harrow (2017) L. Eldar and A. W. Harrow, in  _Foundations of Computer Science (FOCS), 2017 IEEE 58th Annual Symposium on_ (IEEE, 2017) pp. 427–438.
  * Crosson and Bowen (2017) E. Crosson and J. Bowen, arXiv preprint arXiv:1703.10133 (2017).
  * Bremner  _et al._ (2016) M. J. Bremner, A. Montanaro, and D. J. Shepherd, Physical Review Letters 117, 080501 (2016).
  * Zhu  _et al._ (2016) H. Zhu, R. Kueng, M. Grassl, and D. Gross, arXiv preprint arXiv:1609.08172 (2016).
  * Gross (2006) D. Gross, Journal of mathematical physics 47, 122107 (2006).
  * Harrow (2013) A. W. Harrow, arXiv preprint arXiv:1308.6595 (2013).
  * Åke Björck and Golub (1973) Åke Björck and G. H. Golub, Mathematics of Computation 27, 579 (1973).
  * Hoeffding (1963) W. Hoeffding, Journal of the American Statistical Association 58, 13 (1963).
  * Campbell (2011) E. T. Campbell, Physical Review A 83, 032317 (2011).
  * Boyd and Vandenberghe (2004) S. Boyd and L. Vandenberghe,  _Convex optimization_ (Cambridge university press, 2004).
  * Zyczkowski and Sommers (2000) K. Zyczkowski and H.-J. Sommers, Journal of Physics A: Mathematical and General 33, 2045 (2000).


* [6] H. Godinho and T. C. de Souza Neto.  Pairs of additive forms of degrees $2.3^{\tau}$ and $4.5^{\tau}$.  J. Comb. Number Theory, 3(2):87–102, 2011.
  * [7] H. Godinho and T. C. de Souza Neto.  Pairs of additive forms of degree $p^{\tau}(p-1)$.  Funct. Approx. Comment. Math., 48(part 2):197–211, 2013.
  * [8] H. Godinho, M. P. Knapp, and P. H. A. Rodrigues.  Pairs of additive sextic forms.  J. Number Theory, 133(1):176–194, 2013.
  * [9] H. Godinho and L. Ventura.  Pairs of diagonal forms of degree $3^{\tau}.2$ and Artin’s conjecture.  J. Number Theory, 177:211–247, 2017.
  * [10] C. Kränzlein.  Pairs of additive forms of degree 2 n.  PhD thesis, Universität Stuttgart, 2009.
  * [11] H. B. Mann and J. E. Olson.  Sums of sets in the elementary Abelian group of type $(p,\,p)$.  J. Combinatorial Theory, 2:275–284, 1967.
  * [12] J. E. Olson.  A combinatorial problem on finite Abelian groups. I.  J. Number Theory, 1:8–10, 1969.
  * [13] J. E. Olson.  A combinatorial problem on finite Abelian groups. II.  J. Number Theory, 1:195–199, 1969.
  * [14] S. Poehler.  Two additive quartic forms.  PhD thesis, Universität Stuttgart, 2007.