# The Details of Limb Brightening Reveal the Structure of the Base of the Jet
in M 87 for the First Time

Brian Punsly  1415 Granvia Altamira, Palos Verdes Estates CA, USA 90274
ICRANet, Piazza della Repubblica 10 Pescara 65100, Italy  ICRA, Physics
Department, University La Sapienza, Roma, Italy

(Received …; Revised …; Accepted …)

###### Abstract

It has become commonplace is astronomy to describe the transverse coarse
structure of jets in loosely defined terms such as “sheath” and “spine” based
on discussions of parsec scale properties. But, the applicability, dimension
and prominence of these features on sub-lt-yr scales has previously been
unconstrained by observation. The first direct evidence of jet structure near
the source in M 87 is extreme limb brightening (a double-rail morphology), 0.3
- 0.6 mas from the source, that is prominent in observations with high
resolution and sensitivity. Intensity cross-cuts of these images provide three
strong, interdependent constraints on the geometry responsible for the double-
rail morphology: the rail to rail separation, the peak to trough intensity
ratio and the rail widths. Analyzing these constraints indicates that half or
more of the jet volume resides in a thick-walled, tubular, mildly
relativistic, protonic jet only $\sim 0.25$ lt-yr (or $\sim 300$ M, where M is
the central black hole mass in geometrized units) from the source. By
contrast, the Event Horizon Telescope Collaboration interprets their
observations with the aid of general relativistic magnetohydrodynamic
simulations that produce an invisible (by construction) jet with a surrounding
luminous, thin sheath. Yet, it is shown that synthetic images of simulated
jets are center brightened 0.3 - 0.6 mas from the source. This serious
disconnection with observation occurs in a region previously claimed in the
literature to be well represented by the simulations. The limb brightening
analysis motivates a discussion of possible simulation modifications to
improve conformance with observations.

black hole physics — galaxies: jets—galaxies: active — accretion, accretion
disks

††journal: ApJ

##  1 Introduction

The powerful jet in the nearby galaxy M 87 ($\approx 16.8$ Mpc distant) is the
best candidate for observing jet launching and the jet base in detail,
therefore a prime target of $\sim$1.3mm VLBI (Very Long Baseline
Interferometry), the Event Horizon Telescope (EHT)(Doeleman et al., 2012). The
EHT has not yet detected a jet, only an annulus of emission (EHT Collaboration
et al., 2019). The observations are very challenging and are interpreted with
the aid of a library of general relativistic magnetohydrodynamic (GRMHD)
numerical simulations (Porth et al., 2019). There is little observational
information on the jet configuration just above the EHT annulus and therefore
little to restrict the jet produced in the GRMHD simulations. To improve this
situation, this study explores the extreme limb brightening in the highest
sensitivity 3.5mm and 7mm VLBI experiments in order to quantify the dimensions
and composition of the transverse structure $\sim 0.25-0.45$ lt-yrs downstream
of the EHT annulus.

The highest signal to noise and sensitivity 43(86) GHz VLBI image is
reproduced in the top (bottom) left hand panel of Figure 1. The details of the
43(86) GHz image and the observations have been discussed previously in Walker
et al. (2018); Punsly (2021) and (Hada et al., 2016; Kim et al., 2018). The
images are formatted in a manner that accentuates the jet’s transverse
structure. The contour levels are chosen to highlight the jet at the expense
of not capturing a much brighter core and an unresolved region of the jet. The
images are rotated so the jet axis, z, is horizontal, i.e. a rotation of
PA=$-23^{\circ}$ (Hada et al., 2016). The elliptical Gaussian restoring beam
(the blue ellipses) have full width at half maximum, FWHM, 0.16 mas x 0.21
mas, PA=$0^{\circ}$ and 0.116 mas x 0.307 mas, PA=$-9^{\circ}$ at 43 and 86
GHz, respectively. These are also rotated by PA=$-23^{\circ}$ in Figure 1. The
86 GHz image in Figure 1 is the Kim et al. (2018) reconstruction of the radio
source that was observed and reconstructed originally in (Hada et al., 2016).
Assuming intrinsic bilateral symmetry, it was deduced from the 43 GHz image
that the associated Doppler boosting (with a line of sight (LOS) of
$18^{\circ}$ to the jet axis EHT Collaboration et al. (2019); Chael et al.
(2019)) corresponds to mildly relativistic speeds for the bulk velocity of the
jet, accelerating from $\sim 0.27$c (z= 0.4 mas) to $\sim 0.38$c (z=0.65
mas)(Punsly, 2021). The bright ridges (“limb brightening” in the sky plane)
are the main topic of this study. An intensity cross-cut at z = 0.5 mas is
reproduced in the right hand panel of Figure 1 (Punsly and Chen, 2021). Notice
that the peak to trough intensity ratios are very large $\sim 6-10$ at 43 GHz.
The figure also compares the width of the ridges (black double Gaussian fit)
to the Gaussian restoring beam projected along the cross-cut (dashed green).
There are 6 cross-cuts for a total of 12 ridges in (Punsly and Chen, 2021). In
that paper the restoring beam was de-convolved from the black Gaussian fits.
The 12 de-convolved ridge FWHM were on average 42% of the jet radius, $R$.
Ostensibly, this indicates intrinsically wide ridges. Wide intensity ridges
spatially juxtaposed to a central deep trough are two behaviors that seem to
be at odds with each other. The intrinsic configuration is obfuscated by the
finite resolution of VLBI and a nearly polar LOS. The main goal of this study
is to understand the intrinsic source geometry that is compatible with the
detailed structure of the intensity cross-cuts.

Figure 1: A high sensitivity 43 (86) GHz VLBI jet image is shown on the top
(bottom) left. The limbs of the approaching jet in the 43 GHz image are very
bright relative to a faint axial trough. The large degree of limb brightening
is shown in the cross-cut on the top right. The bottom right hand frame
compares intensity cross-cuts of a synthetic image in Figure 4 (dashed curves)
from the most promising candidate simulation (“MAD: a/M=0.94”) to cross-cuts
of the VLBI image at 86 GHz (solid curves).

##  2 The Extreme Limb Brightening and Wide Jet Ridges

Limb brightening is a phenomenon associated with astrophysical objects such as
enhanced ovals and ridges of emission in planetary nebulae and supernovae
remnants (Masson, 1990; Corradi et al., 2003; Jones et al., 1998). These
objects are relatively nearby and their dimension are easily resolved. With
very high resolution, the peak to trough ratio diverges as $1/\sqrt{W}$ for a
narrow ridge width, $W$ (Masson, 1990). But the limited resolution of VLBI is
relatively insensitive to jet wall variations if $W$ is less than half of the
beam FWHM (see top right of Figure 2). The other obfuscating circumstance is
the nearly polar LOS and gradients in the emissivity along the jet. The LOS
crosses the near side of the jet at low emissivity and exits the far side at
high emissivity. In a cross-cut, the intensity trough includes contributions
from LOS exit points from the jet wall that are closer to the nucleus than any
other LOS and this is the region of highest emissivity (x=0 in bottom panel of
Figure 2). This makes the very extreme limb brightening even more difficult to
explain.

Figure 2: The top left panel is an experiment to explore ratio of the
intensity peak to trough as a function of the geometrical model of the jet.
The top right experiment explores the effect of wall thickness on the cross-
cut. The nomenclature, $0.05<W/R<0.25$, is an abbreviation for a cone plus
cylinder with $W/R=0.25$ at the junction of the cylinder and cone and tapers
to $W/R=0.05$ when the LOS exits the cone wall at x=0. In order to clarify
this, an example of a cone with tapered walls plus cylinder is shown in detail
in the bottom panel. The image is a snapshot from a 3-D model in Paraview
3.3.0. The blue, (x, y), plane is foliated by the set of all LOS that
constitute a cross-cut. The intersection of the plane with the tubular jet is
in dark olive green. The key feature for extreme limb brightening is that the
intersection is narrow at the exit of the LOS plane from the tubular wall on
the right hand side. The dashed black rectangle shows a cutout of the region
of maximum limb brightening that is magnified at the bottom. This region has
the longest LOS through the tubular jet wall (red line) as indicated by the
white grid overlayed on the blue LOS plane.

The intensity cross-cuts provide three strong, interdependent constraints, the
peak to peak separation, the peak to trough ratio and the ridge widths. Limb
brightening was investigated with simplified axisymmetric geometrical models
of emissivity that are not intended to replicate the actual jet. They are
investigative tools to understand which geometrical properties can contribute
to both the deep trough and the wide ridges with a minimal number of degrees
of freedom. For $z<1$ mas, it was previously found that $\alpha\approx 0.8$
from 43 GHz to 86 GHz, where the flux density, $S(\nu)\sim\nu^{-\alpha}$ Hada
et al. (2016). Thus, the jet appears to be optically thin in the frequency
range of interest. Therefore, absorption is ignored in the equation of
radiative transfer. Each axisymmetric model is used to generate cross-cuts by
numerically approximating the following convolution integral,

$I(x_{o},\,z_{o})=\int{\int{\int{N\epsilon(x,y,z)}\frac{1}{2\pi\sqrt{\mid\mathbf{a}\mid\mid\mathbf{b}\mid}}\rm{e}^{{-\frac{\mid\mathbf{r}\cdot\mathbf{a}\mid^{2}}{2(\mathbf{a}\cdot\mathbf{a})^{2}}}-\frac{\mid\mathbf{r}\cdot\mathbf{b}\mid^{2}}{2(\mathbf{b}\cdot\mathbf{b})^{2}}}}}dy\,dx\,dz\;,\quad\epsilon(x,y,z)\approx\epsilon(y)=-0.781y+1\;.$
(1)

where, the restoring beam is defined by an elliptical Gaussian with major
(minor) axis $\mathbf{a}$ ($\mathbf{b}$) of length of the standard deviation
$\sigma$= FWHM/2.35, the displacement from $(x_{o},z_{o})$ in the sky (or x-z)
plane is $\mathbf{r}$ and the emissivity is $\epsilon$. $N$ is the only free
parameter once a geometrical model is chosen and sets the normalization to
agree with the peak intensity (note 0.1 mas $\approx 2.5\times 10^{16}$ cm).
The intensity, $I(x_{o},\,z_{o})$, is evaluated in a cross-section at
$z=z_{o}$ at all points $x=x_{o}$ along the transverse slice. $y$ is
orthogonal to the sky plane with LOS = $18^{\circ}$. Since
$dy/dz^{\prime}=\cos{18^{\circ}}\approx 1$ ($z^{\prime}$ is measured along the
symmetry axis of the jet), from Table 1, assuming a cross-cut at $z\sim
0.4$mas, we can estimate the emissivity in the axisymmetric approximation. The
emissivity, $\epsilon(x,y,z)\approx\epsilon(y)$ is given in Equation (1),
where the LOS intersects the jet wall, and is 0 elsewhere. The function is
referenced to the point where the outer boundary of the slice through the jet
wall by the LOS plane is tangent to the y axis as indicated in the bottom
panel of Figure 2. At this point, $\epsilon(y)=1$ and $y=0$. This is
emissivity in the observer’s frame of reference, so it includes effects like
Doppler boosting. $\epsilon(x,y,z)$ in Equation (1) is suitable for the
purposes of this study because in the model analysis, the peak to trough
ratios and the ridge width are relatively insensitive to the emissivity
profile (compared to the geometry) unless the gradients increase by at least
an order of magnitude from those indicated by Table 1.

Figure 2 highlights the results of an investigation of the geometrical
properties that can produce deep surface brightness troughs between adjacent,
wide intensity peaks. The top left panel investigates the effect of geometry
on the depth of the trough. The top right panel explores the relationship
between widening the jet walls, the intensity ridge width and the trough
depth. The 43 GHz VLBI cross-cut in the top left panel was chosen for the
investigation of the trough depth for two reasons. First, it is the most
symmetric cross-cut from Punsly and Chen (2021) and the models are
axisymmetric. Secondly, it is the only cross-cut from Punsly and Chen (2021)
in which models were found that achieve the trough depth relative to the peak.
Note that the trough in the top right hand panel of the Figure 1 is deeper.
These deeper troughs that exist downstream of z = 0.405 mas cannot be achieved
with these axisymmetric models. The same situation is true with axisymmetric
models of extreme limb brightening in planetary nebulae. It has been proposed
that fortuitous asymmetric structure is required for extreme limb brightening
in planetary nebulae (Masson, 1990). This might be the case for M87 as well.
The extra degrees of freedom implicit in this scenario are not desirable to
introduce here as the effort is focused on understanding the basic geometry
that is responsible for a persistent prominent phenomenon of deep troughs and
wide ridges. Speculation on the asymmetries required in the jet wall will be
left for future research.

The investigation began with the study of simple cones and cylinders. The wall
widths of the cones and cylinders were varied and the opening angle of the
cone was varied as well. The axisymmetric emissivity function was also allowed
to vary. No model could even produce the shallowest troughs in the 43 GHz
cross-cuts from Punsly and Chen (2021), i.e.that which is shown in the top
left frame of Figure 2. The parameter $W/R$ in of itself cannot produce these
deep troughs. The general trends are summarized in the top left panel of
Figure 2. Choosing a value smaller than $W/R=0.1$ does not deepen the trough
very much due to the relatively low resolution of the restoring beam and it
hinders the ability to create a wide ridge as indicated in the top right hand
panel of Figure 2. The models require another degree of freedom in order to
succeed. The model that produces the deepest troughs in the top left panel of
Figure 2 has, $W/R$, that decreases linearly towards a conical apex as in the
bottom panel of Figure 2 (“tapered cone”, hereafter).

The 3-D spatial relationships are rather complex and the bottom panel of
Figure 2 is provided as a visual aid. The dark, olive green region is where
the LOS, along a cross-cut, intersects the wall of the tubular jet. The LOS in
the middle of the cross-cut at x=0 slices through two regions that are very
far apart. Since, M 87 has a deep trough and the emissivity is highest close
to the nucleus, the tapering of the cone wall helps achieve a deep trough by
decreasing the area of the intersection of the blue LOS plane with the jet
wall on the far right, where the emissivity is high. Table 1 indicates a
nearly cylindrical jet for $z=0.4-0.65$ mas, so our fiducial geometry is a
tapered cone at the base and a cylinder at larger z. A $28^{\circ}$ cone
intrinsic opening angle is chosen (Walker et al., 2018). In summary, the
additional required parameter that creates the deep troughs is the amount of
wall thickness tapering in the conical region. The variable, $W/R$, at the
junction of the cylinder and the cone cannot create the required trough depth.
However, the conical wall taper can produce troughs representative of those in
the 86 GHz cross-cuts for a wide range, $W/R\leq 0.35$. If $W/R$ tapers by a
factor of $\sim 4$, the troughs of the 86 GHz VLBI cross-cuts can be achieved
as in the top right hand frame of Figure 2. The method is limited by the
relatively low resolution of VLBI and cannot produce the deeper troughs seen
at 43 GHz (like the one in the top right hand panel of Figure 1) regardless of
the rate of the linear tapering of the wall thickness.

Table 1: Jet Properties of the Image in Figure 1 (Punsly and Chen, 2021) $z$ |  $z^{\prime}\equiv$ the De-projected | Distance | Flux Density
---|---|---|---
| Axial | Between | of 0.16 mas Wide
| Displacement | Intensity Peaks | Cross-Section
(mas) | (mas) | (mas) | (mJy)
0.405 | 1.31 | $0.514\pm 0.040$ | $66.4\pm 6.7$
0.45 | 1.46 | $0.555\pm 0.040$ | $52.2\pm 5.3$
0.5 | 1.62 | $0.542\pm 0.040$ | $44.6\pm 4.6$
0.55 | 1.78 | $0.548\pm 0.040$ | $33.5\pm 3.5$
0.6 | 1.94 | $0.518\pm 0.040$ | $32.2\pm 3.4$
0.65 | 2.11 | $0.523\pm 0.040$ | $33.2\pm 3.4$

The top left panel of Figure 2 indicates a degeneracy in parameter space: by
just making the walls very thin in the tapered region of the cone, one can get
troughs of depths similar to the z = 0.405 mas cross-cut, regardless of W/R.
However, the wider ridges adjacent to these troughs are a definite challenge
to comprehend. The top right panel of Figure 2, studies the effect of the wall
thickness of the model on the observed ridge width. The models under
consideration are axisymmetric, so the study utilizes the most symmetric cross
section of the 86 GHz image in which the ridges are resolved (z = 0.55 mas)
for analysis in the upper right hand panel of Figure 2. Since the restoring
beam is relatively large, the convolved ridge thickness increases slowly with
increasing $W$ and the trough gets shallower (see the top right of Figure 2).
The best description of the cross-cuts has $W/R=0.35$ at the junction of the
cylinder and cone and tapers to $W/R=0.08$ when the LOS exits the cone wall at
x=0 (abbreviated as $0.08<W/R<0.35$, hereafter). But, the ridges are not as
wide as observed. If one chooses $W/R=0.40$, the trough depth requirement will
be breached at the expense of a little more observed ridge width. This is the
essence of the conflict of wide ridges adjacent to deep troughs, the
fundamental conundrum presented by the observations. These simplest models are
close to observation, but not sufficient.

To gain additional intensity ridge width, without raising the trough
substantially, a sheath (i.e., a dissipative shear layer between the external
medium and the thick-walled tubular jet) was wrapped around the tubular jet.
It has the same shape as the tubular jet, the inner profile of the sheath is
the same as the outer profile of the tubular jet. It is therefore a wider
version of the geometry in Figure 2, with narrow walls and a taper abbreviated
as $0.03<W/R<0.1$. This is successful as a perturbation because the geometric
source for the weak sheath is wider by construction, so it is more resolved by
the convolving beam. This accentuates the difference between long and short
path lengths through the jet wall, i.e. it deepens the trough. This also
explains the deeper troughs as the physical jet expands in the bottom right of
Figure 1 and Figure 3 which shows the resulting fits. Since the models are
axisymmetric by construction, they cannot reproduce a cross-cut with ridges of
different peak intensity. Hence, we introduce the cross-cut at z = 0.35 mas in
the left hand panel of Figure 3, even though the jet is too narrow to be fully
resolved by the VLBI experiment. The most symmetric cross-section that is
transversely resolved is at z = 0.55 mas and the corresponding fit is shown on
the right hand side of Figure 3. Notice that the sheath is merely a small
perturbation of the thick walled jet in the top panels. However, the fit is
not unique. If one chooses a smaller value of the wall thickness in the
tubular jet, say $W/R=0.1$ for the sake of argument, then the contribution of
the sheath must increase considerably in order to fill-in the wide wings of
the ridges (see the bottom panel of Figure 3). It is no longer just a
perturbation. In order to achieve the required limb width there is a gap of
$\sim 0.105$ mas between the intensity peak of the sheath and intensity peak
of the wall. Since the sheath is of comparable surface brightness to the jet
wall, the model of the emission is the aggregation of the two tapered cones
plus cylinders and the gap. One can compute the total thickness of the
emissivity region to get $W_{\rm{total}}/R=0.37$. A wide emission region
surrounding a very low luminosity core seems inescapable regardless of the
partitioning into different labels. The gap is not impossible, but an adjacent
sheath is more compatible with an enveloping shear layer. The choice here is
to make the jet wall as close to the complete solution as possible in order to
minimize the contribution of the outer sheath, i.e., make it a perturbation
representative of an adjacent shear layer. These results do not prove an exact
model this was never the intent of the simplified axisymmetric models), but
they indicate that a very wide structure of a tubular wall and surrounding
sheath are needed to create the wide ridges. An approximate estimate of the
ratio of the tubular jet volume to the total jet volume is $\geq
1-(1-W/R)^{2}=0.58$ at z = 0.35 mas and z = 0.55 mas.

Figure 3: The top panels of the figure show the fit to two intensity cross-
cuts. There are two components: the main body of the jet, a tapered cone plus
cylinder ($0.08<W/R<0.35$) and a sheath, an outer tapered cone plus cylinder
($0.03<W/R<0.1$). The cross-cut on the top left is the most symmetric cross-
cut from the 86 GHz VLBI image and is the most compatible with an assumed
axisymmetric model. The lack of symmetry between the north and south ridges in
the cross-cut on the top right renders a single axisymmetric model imperfect.
The bottom panel is an alternate method of implementing the model. In this
model, the jet wall thickness is narrow ($0.03<W/R<0.1$). This fit has
comparable surface brightness contributions from the putative jet and the
sheath.

It is worth noting that one can also consider the cross-cuts of the super-
resolved version of the 86 GHz image in the bottom left frame of Figure 1 that
was previously published (Hada et al., 2016; Punsly et al., 2018). The image
was restored with a circular beam of 0.11 mas instead of the 0.116 mas x 0.307
mas, PA=$-9^{\circ}$ restoring beam of Figure 1. The axis that is parallel to
the major axis of the synthesized beam is heavily super-resolved (by $\approx
30\%$). However, there are no baselines that are long enough to resolve
objects as small as 0.11 mas along the major axis of the synthesized beam.
Thus, results based on an individual feature of the super-resolved image are
not as reliable as those based on the source reconstruction in Figure 1. The
results are only suggestive taken as a whole and certainly not very reliable
taken individually. The relevant cross-cuts of the super-resolved image from
z= 0.3, 0.4 and 0.5 mas were previously published (Punsly et al., 2018). The
main advantage of a high resolution image is the possibility of achieving a
stronger connection between the FWHM of the de-convolved Gaussian fits of the
peaks, $\rm{FWHM}_{\rm{deconvolved}}$, and the intrinsic width of the jet
walls. Thus motivated, the six intensity peaks from the three cross-cuts were
fit with Gaussian functions then de-convolved from the circular restoring
beam. $R$ is defined as one-half the distance between the outer de-convolved
half maxima of the two ridges. It is found that the mean value of
$\rm{FWHM}_{\rm{deconvolved}}/R=0.41$ and the median value is
$\rm{FWHM}_{\rm{deconvolved}}/R=0.35$. This is consistent with the conclusions
presented above.

##  3 Simulated Jets Contrasted with the Sub-mas Jet in M 87

It was previously noted that synthetic images of numerically simulated jets
(typical of those in the EHT Collaboration library) were much narrower (about
a factor of $0.4-0.7$) compared to 86 GHz VLBI observations of jets for $z\leq
0.25$ mas (Punsly, 2019; Cruz-Osorio et al., 2021). However, the apparent non-
existent limb brightening for $z<0.6$ mas could not be quantified in an easy
to grasp format since the raw data of the simulations were not made public
(Punsly, 2019; Chael et al., 2019; Mościbrodzka et al., 2016). Fortunately,
publicly shared synthetic image (Flexible Image Transport System, FITS) files
generated from general relativistic simulations were recently created with the
same restoring beam as the 86 GHz VLBI observation in Figure 1 (Cruz-Osorio et
al., 2021).

The lower right panel of Figure 1 shows intensity cross-cuts extracted from
the synthetic image FITS file compared to cross-cuts the same distance down
the jet axis for the 86 GHz VLBI observation. The simulated jets have a peak
in the middle and the VLBI jet has a deep trough in the middle. The simulated
jet also fades too rapidly with z, both properties appear in other published
numerically simulated synthetic images (Mościbrodzka et al., 2016; Chael et
al., 2019; Punsly, 2019; Fromm et al., 2022). The EHT Collaboration have
assessed their large library of simulations and have determined that only
their simulations of magnetically arrested accretion (MAD) can explain the
coarse properties of the annulus and produce a jet with enough power (EHT
Collaboration et al., 2019, 2021). It has been shown that high spin black
holes, $0.5<a/M<1$ (where $Ma$ is the angular momentum of the black hole in
geometrized units: black holes are only defined for $-1<a/M<1$) produce the
widest MAD jets and are therefore the least nonconforming with the observed
jet width, noted above (Cruz-Osorio et al., 2021; Fromm et al., 2022). The
MAD, a/M=0.94 simulations have a brighter jet than the MAD, a/M=0.5
simulations (Cruz-Osorio et al., 2021). Considering the rapid fading of the
simulated jet (in the lower right panel of Figure 1) compared to the VLBI jet,
the a/M=0.94 simulations are closer to reality. Thus, the simulation presented
in Figures 1 and 4 was chosen for maximum conformance to observation.

The highest emissivity in the simulated jets occurs in a narrow, surrounding
“sheath” adjacent to an invisible (by construction) jet. Edge brightening
(only the edge of the jet is visible) is maximized by construction, but limb
brightening on the sky plane is never pronounced in synthetic images
(Mościbrodzka et al., 2016; Chael et al., 2019; Cruz-Osorio et al., 2021;
Fromm et al., 2022; Yang et al., 2022). To understand this apparent
contradiction, consider the left hand side of Figure 4 which shows the source
function before convolution with the restoring beam. We superimposed ellipses
representing the restoring beam FWHM at the locations of the peak surface
brightness (in Figure 1) for three cross-cuts. Notice that these peaks are not
due to LOS geometry and limb brightening. They simply represent local
emissivity maxima on a background of low emissivity (almost voids in
brightness). These bright areas arise from a low filling factor set of high
dissipation regions (perhaps shocks). If the tubular jet had thicker walls and
a higher median emissivity, the jet walls would provide a much larger Gaussian
weighted emissivity column within the restoring beam. As such, the bright
areas would not dominate and limb brightening would be conspicuous on the sky
plane. As expected, a parameter study of emissivity functions and locations of
boundaries between the sheath and the invisible leptonic jet never changes
this wafer thin appearance of the sheath (Fromm et al., 2022). Thus, near the
jet base, the limb brightening on the sky plane will never be strong in this
family of simulations. The absence of limb brightening in the sly plane in the
synthetic images is more evidence supporting the conclusion of the analysis in
the previous section, the jet walls in M 87 are very thick.

Figure 4: The source function (left) for the synthetic image on the right
provides some clues to the lack of limb brightening on the sky plane. The
dashed black ellipses are the FWHM of the restoring beam at the peak intensity
of the cross-cuts in Figure 1. The locations of the cross-cuts are indicated
in the image on the right (Cruz-Osorio et al., 2021). The bright peaks are
just local bright spots on an otherwise mostly faint, hollow, conical jet.
There is no significant geometrical limb brightening on the cross-cuts.

##  4 Conclusion

The analysis of limb brightening determined that the source of the observed
jet base emission in M 87 is a tubular jet with more than half of the total
jet volume. If the jet were leptonic and magnetically dominated, being beyond
the light cylinder (Punsly and Chen, 2021), the particles would flow at highly
relativistic velocities (McKinney, 2006), contrary to the mild Doppler
boosting inferred from the top left image in Figure 1. This is consistent with
the unique magnetohydrodynamic solution consistent with that image. It was
determined that the only outflows that satisfy the mildly relativistic bulk
flow velocity profile, the flux density distribution along the jet as well as
the conservation laws of mass, energy, and angular momentum are tubular jets
of protonic plasma driven by a large scale magnetic field with a significant
Poynting flux (Punsly and Chen, 2021). The energy flux at this epoch is
$\gtrsim 10^{42}$ erg/sec and the wind transports $\gtrsim 10^{-4}$ solar
masses per year (Punsly and Chen, 2021).

The results of this study can be used as guidelines to modify the numerical
simulations when considered from a comprehensive perspective. There are three
issues that are inter-related and need to be considered simultaneously,

  1. 1.

At a de-projected axial distance of $\sim 100$M – 200M from the nucleus, the
synthetic images of simulated high black hole spin, MAD, jets are
approximately $\sim 0.4-0.7$ as wide as the jet in M87 (Punsly, 2019; Chael et
al., 2019; Cruz-Osorio et al., 2021). This implies that the luminous jet
streamlines in simulations need to diverge more rapidly as the jet propagates
away from the launch site in order for jet width to be consistent with the
observed transverse structure.111A new simulation claims to have resolved this
issue (Yang et al., 2022). However, since the jet is very center brightened,
they compare the FWHM of the of the central ridge in the simulations to the
ridge to ridge separation in the “double rail” images. This is not an
equitable method of comparison. If the methods of Punsly (2019) are applied to
these same synthetic images, no width difference is seen at $z<0.45$ mas from
the synthetic images of Chael et al. (2019). So this issue seems not to be
resolved.

  2. 2.

300 M - 550 M downstream of the nucleus, the jet is extremely limb brightened
and the simulated jets are center brightened.

  3. 3.

The extreme limb brightening and wide ridges can only be made by a very thick-
walled tubular jet with over half the total jet volume.

A conversion of 3.8 $\mu$as = 1M and $M\approx 9.6\times 10^{14}$ cm was
assumed in the above (EHT Collaboration et al., 2019).

One possible resolution is to extend the emissivity into the jet to fix points
2) and 3). This is a scientific and computational issue because the mass flux
is set by an invoked mass floor and numerical diffusion in the jet (Chael et
al., 2019). Also, this would make the jet appear narrower, therefore more
nonconforming to point 1). Furthermore, it would likely create a bulk flow
velocity that exceeds the mildly relativistic estimates based on Figure 1
(McKinney, 2006; Punsly, 2021). Alternatively, one might want to extend the
emission region outward into the protonic region. The only collimated protonic
outflows in the EHT Collaboration library of simulations are the “funnel wall
jets” created as the accreting gas rams into a centrifugal barrier and is
accelerated by pressure in the accretion disk’s corona (Porth et al., 2019;
Hawley and Krolik, 2006). This is a very thin outflow (perhaps much of what we
see in Figure 4). There is also the possibility of a magneto-centrifugally
driven wind from the disk. But there are few organized large scale polodoidal
fields in the disk in these simulations. It is mainly chaotic (turbulent) with
a net poloidal excess when averaged over time and azimuth, these are not the
properties that lead to a collimated jet (probably the reason for its absence
in synthetic images) (Tchekhovskoy, 2015; Blandford and Payne, 1982; Qian et
al., 2018). Another possibility is that there is a huge amount of entrainment
in this first 0.25 lt-yrs of jet propagation ($\sim 10^{-4}$ solar masses per
year) that is not captured by the simulations. Not only would there need to be
large amounts of entrainment, but it would need to migrate 1/3 of the distance
to the jet center. For geometrical reasons (larger lever arm), these are the
streamlines with most of the created jet Poynting flux power in the
simulations (Krolik et al., 2005), so the interior smaller volume would carry
the minority of the jet power. This would contradict the notion that most of
the jet power is hidden in an invisible core or spine (Ghisellini et al.,
2005). This scenario does not remedy the non-conformance in point 1).

This study benefitted from the comments of a helpful referee. I am indebted to
Sina Chen for the intensity cross-cuts. Christian Fromm generously provided
the image FITs files used in this paper from Cruz-Osorio et al.(2021). Alan
Marscher and Sam Gralla provided comments on earlier versions of this
manuscript that greatly improved the science.

## References

  * Blandford and Payne (1982) Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883
  * Chael et al. (2019) Chael, A., Narayan, R., Johnson, M. D. 2019, MNRAS 486, 2873
  * Corradi et al. (2003) Corradi, R. L. M.; Schönberner, D.; Steffen, M.; Perinotto, M. 2003 MNRAS 340, 417
  * Cruz-Osorio et al. (2021) Cruz-Osorio, A.; Fromm, C. M.; Mizuno, Y. et al. 2021 Nature Astronomy NatAs.tmp..218C doi:10.1038/s41550-021-01506-w
  * Doeleman et al. (2012) Doeleman, S., Fish, V., Schenck, D. et al. 2012, Science 338 355
  * EHT Collaboration et al. (2019) Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2019, ApJL, 875, L5
  * EHT Collaboration et al. (2021) Event Horizon Telescope Collaboration, Akiyama, K., Algaba, J. C., et al. 2021, ApJL, 910, L13
  * Fromm et al. (2022) Fromm, C. M.; Cruz-Osorio, A.; Mizuno, Y. 2022, submitted to Astron. & Astro. 660 107
  * Gebhardt et al. (2011) Gebhardt, K., Adams, J., Richstone, D., et al. 2011, ApJ, 729, 119
  * Ghisellini et al. (2005) Ghisellini, G., Tavecchio, F., & Chiaberge, M. 2005 Astron. & Astro. 432 401
  * Hada et al. (2016) Hada, K., Kino, M., Doi, A., et al. 2016, ApJ, 817, 131
  * Hawley and Krolik (2006) Hawley J., Krolik J., 2006, ApJ, 641, 103
  * Jones et al. (1998) Jones, T., Rudnick, L. Jun, B.-L. et al 1998 PASP 110 125
  * Kim et al. (2018) Kim, J.-Y., Krichbaum, T., Lu, R.-S., Ros, E. et al. 2018 A&A 616 188
  * Krolik et al. (2005) Krolik, J.H., Hawley, J.F. & Hirose, S., 2005 ApJ 622 1008.
  * Masson (1990) Masson, C. 1990, ApJ 348, 580
  * McKinney (2006) Mckinney, J., 2006 MNRAS 368 1561
  * Mościbrodzka et al. (2016) Mościbrodzka, M., Falcke, H., Shiokawa, H. 2016 A & A 586 38
  * Porth et al. (2019) Porth, O., Chaterjee, K., Narayan, R. et al. 2019 ApJS 243 26
  * Punsly et al. (2018) Punsly, B.; Hardcastle, M.; Hada, K. 2018, A & A 614 104
  * Punsly (2019) Punsly, B. 2019, ApJL 879, 11
  * Punsly (2021) Punsly, B. 2021, ApJ 918, 4
  * Punsly and Chen (2021) Punsly, B. and Chen, S. 2021, ApJL 921 L38.
  * Qian et al. (2018) Qian, Q. Fendt, C. Vourellis, C., 2018 ApJ 859 28.
  * Tchekhovskoy (2015) Tchekhovskoy, A. 2015 The Formation and Disruption of Black Hole Jets, ASSL, 414 ISBN 978-3-319-10355-6. Springer International Publishing Switzerland, p. 45
  * Walker et al. (2018) Walker R. C., Hardee P., Davies, F., Ly, C., Junor, W., 2018, ApJ 855 128
  * Yang et al. (2022) Yang, H., Yuan, F., Li, H. et al, 2022, https://arxiv.org/abs/2206.05661


An Analysis of the Quantum Approximation Optimisation Algorithm

Behzad Mansouri*

QuantGates Ltd, London, EC2A 4NE, UK

Abstract

> This article consists of a short introduction to the quantum approximation
> optimisation algorithm (QAOA). The mathematical structure of the QAOA, as
> well as its basic properties, are described. The implementation of the QAOA
> on MaxCut problems, quadratic unconstrained binary optimisation problems
> (QUBOs), and Ising-type Hamiltonians is considered in detail.

††*<EMAIL_ADDRESS>
##  1 The Quantum Approximate Optimisation Algorithm

The quantum approximation optimisation algorithm (QAOA) was first constructed
and implemented on MaxCut problems by Farhi et al. [1]. The QAOA was
specifically designed to run on circuit-based quantum computers. The quality
of the approximation provided by the algorithm directly depends on an integer
$p$, and with $p$ growing slowly with $n$ (the input size, $n$ bits) we hope
the algorithm will be powerful enough to find solutions beyond the
capabilities of classical algorithms. In [2], the QAOA was applied to E3LIN2
(MAX-3XOR) and initially beat the state-of-the-art classical algorithm.
Shortly after, a better classical algorithm was constructed that outperformed
the QAOA in this case [3]. Besides these empirical results, [4] argued that
the QAOA can exhibit a form of quantum supremacy in such a way that even the
lowest depth of the algorithm cannot be simulated efficiently using any
classical computer, based on, of course, reasonable assumptions in complexity
theory. Regarding the implementation of the QAOA for a variety of
combinatorial optimisation problems, Hadfield et al. [5] provided a framework
for designing QAOA circuits for these problems, including graph colouring, the
travelling salesman problem (TSP), and single machine scheduling (SMS).
Moreover, the QAOA has been extended by introducing the quantum alternating
operator ansatz [6] to encompass a more general class of quantum states by
allowing alternation between broader classes of operators in which it can
widen the applicability of the QAOA. Herein, we reviewed how the QAOA works
and its implementation on MaxCut problems, and finally, we studied the
implementation of this promising algorithm on quadratic unconstrained binary
optimisation problems (QUBOs) and Ising-type Hamiltonians, in which the QUBOs
embrace many combinatorial optimisation problems, showing that it has
practical applications in every industry.

The QAOA is mainly designed to tackle constraint satisfaction problems (CSPs).
A large number of combinatorial optimisation problems can be recast in the
form of constrained satisfaction problems, the latter of which can be defined
with $n$ bits and $m$ constraints.

The objective function of a CSP can be written as:

$C(z)=\sum\limits_{a=1}^{m}{{C_{a}}}(z)$ (1)

where $z={z_{1}}{z_{2}}...{z_{n}}$ is the bit string and ${C_{a}}(z)=1$ if $z$
satisfies the constraint $a$, and 0 otherwise [1]. By promoting $n$ bits to
$n$ qubits, we can construct a ${2^{n}}$ dimensional Hilbert space with a
computational basis $\left|z\right\rangle$. The constraint functions in the
above can be written in diagonal matrix forms. By promoting $(1)$ as an
operator (diagonal form) acting on Hilbert space vectors, we have:

$C\left|z\right\rangle=\sum\limits_{a=1}^{m}{{C_{a}}}\left|z\right\rangle=\sum\limits_{a=1}^{m}{{P_{a}}}\left|z\right\rangle=P(z)\left|z\right\rangle$
(2)

where ${{P_{a}}}$ is an eigenvalue corresponding to the constraint operator
${{C_{a}}}$. Here, the largest value in the eigenvalue set, $P(z^{\prime})$,
is the maximum of the objective function $C$. When considering a general state
in our Hilbert space as $\left|\psi\right\rangle$, one can evaluate the
expectation value of the operator $C$ as below:

$\begin{array}[]{*{20}{l}}{\left\langle
C\right\rangle=\langle\psi|C\left|\psi\right\rangle=(\sum\limits_{z}{{b_{z}}^{*}}\langle
z|)(\sum\limits_{z}{{b_{z}}}C\left|z\right\rangle)}\\\
{=\sum\limits_{z}{{{\left|{{b_{z}}}\right|}^{2}}P(z)\leq}\sum\limits_{z}{{{\left|{{b_{z}}}\right|}^{2}}P(z^{\prime})}=P(z^{\prime})}\end{array}$
(3)

The maximum of $\left\langle C\right\rangle$ is equal to $P(z^{\prime})$, the
eigenvalue associated with the eigenvector $\left|{z^{\prime}}\right\rangle$.
This setup allows us to create our algorithm. The operator $C$ with its
eigenvalues will $mark$ the computational basis into the two search problem
categories of $solutions$ and $not$-$solutions$, and the problem of finding
the solutions will be addressed by maximising of the expectation value of the
objective function $C$ in a $p$-layered variational ansatz. To achieve this
aim, we review Farhi’s prescription below.

First, a unitary operator is constructed as:

$U(C,\gamma)={e^{-i\gamma C}}=\prod\limits_{a=1}^{m}{{e^{-i\gamma{C_{a}}}}}$
(4)

As $[{C_{a}},{C_{b}}]=0$, it can be written as the product of exponentials.

This quantum gate can be applied on a general state $\left|\psi\right\rangle$;
then, we get:

$\begin{array}[]{l}U(C,\gamma)\left|\psi\right\rangle=\sum\limits_{z}{{e^{-i\gamma\sum\limits_{a=1}^{m}{{C_{a}}}}}{a_{z}}}\left|z\right\rangle\\\
=\sum\limits_{z}{(\prod\limits_{a=1}^{m}{U({C_{a}}},\gamma){a_{z}}}\left|z\right\rangle\end{array}$
(5)

The diagonal constraint operators allow to write $U(C,\gamma)$ as the product
of exponentials, and this makes the implementation of the gate via universal
gates simpler. Note that here, for the eigenbasis
$\left|{z^{\prime}}\right\rangle$ satisfied by constraints ${C_{a}}$, a phase
change (${e^{-i\gamma}}$) is received, and this can make the
${\left|{\left.{a^{\prime}}\right|}\right.^{2}}$ large enough, and obviously,
by the corresponding measurement, we can obtain the state with the higher
probability. One can implement such quantum gates with the help of controlled
phase gates and one ancilla qubit.

As an explicit example, consider
$\left|\psi\right\rangle=\left|{{y_{1}}}\right\rangle\left|{{y_{2}}}\right\rangle$,
where
$\left|{{y_{1}}}\right\rangle={\alpha_{0}}\left|0\right\rangle+{\alpha_{1}}\left|1\right\rangle$
and
$\left|{{y_{2}}}\right\rangle={\beta_{0}}\left|0\right\rangle+{\beta_{1}}\left|1\right\rangle$,
and two constraints as ${{C_{1}}=1}$ if
${\left|{{y_{1}}}\right\rangle=\left|{{y_{2}}}\right\rangle=\left|1\right\rangle}$
and ${{C_{2}}=1}$ if $\left|{{y_{1}}}\right\rangle=\left|1\right\rangle$

First, we directly applied $U(C,\gamma)$ to our two qubit systems.

$\begin{array}[]{l}U(C,\gamma)\left|{{y_{1}}{y_{2}}}\right\rangle={e^{-i\gamma{C_{2}}}}{e^{-i\gamma{C_{1}}}}({\alpha_{0}}{\beta_{0}}\left|{00}\right\rangle+{\alpha_{0}}{\beta_{1}}\left|{01}\right\rangle+{\alpha_{1}}{\beta_{0}}\left|{10}\right\rangle+{\alpha_{1}}{\beta_{1}}\left|{11}\right\rangle)\\\
=({\alpha_{0}}{\beta_{0}}\left|{00}\right\rangle+{\alpha_{0}}{\beta_{1}}\left|{01}\right\rangle+{e^{-i\gamma}}{\alpha_{1}}{\beta_{0}}\left|{10}\right\rangle+{e^{-2i\gamma}}{\alpha_{1}}{\beta_{1}}\left|{11}\right\rangle)\end{array}$
(6)

Now, with help of two controlled phase gates and one ancilla qubit, we
obtained an equivalent result.

Figure 1: Implementation using two controlled-phase gates and one ancilla
qubit.

Here, we define our phase gate:
$R(\gamma)=\left[{\begin{array}[]{*{20}{c}}1&0\\\
0&{{e^{-i\gamma}}}\end{array}}\right]$

The ancilla qubit was prepared in $\left|1\right\rangle$, resulting in:

$\left|{{y_{1}}{y_{2}}}\right\rangle\left|1\right\rangle={\alpha_{0}}{\beta_{0}}\left|{00}\right\rangle\left|1\right\rangle+{\alpha_{0}}{\beta_{1}}\left|{01}\right\rangle\left|1\right\rangle+{\alpha_{1}}{\beta_{0}}\left|{10}\right\rangle\left|1\right\rangle+{\alpha_{1}}{\beta_{1}}\left|{11}\right\rangle\left|1\right\rangle$
(7)

Consider ${R_{1}}(\gamma)$ and ${R_{2}}(\gamma)$ in their explicit matrix
form:

${R_{1}}(\gamma)=\left[{\begin{array}[]{*{20}{c}}I&0&0&0\\\ 0&I&0&0\\\
0&0&I&0\\\
0&0&0&{R(\gamma)}\end{array}}\right]\begin{array}[]{*{20}{c}}{}\hfil&{}\hfil\end{array}{R_{2}}(\gamma)=\left[{\begin{array}[]{*{20}{c}}I&0&0&0\\\
0&I&0&0\\\ 0&0&{R(\gamma)}&0\\\ 0&0&0&{R(\gamma)}\end{array}}\right]$ (8)

Here, the ${R_{1}}(\gamma)$ gate controls two qubits,
$\left|{{y_{1}}}\right\rangle$ and $\left|{{y_{2}}}\right\rangle$, and acts on
our third qubit (auxiliary qubit) as its target. Based on our constraints
here, the gate ${R_{1}}(\gamma)$ only acts $R(\gamma)$ on ancilla when we have
$\left|{{y_{1}}}\right\rangle=\left|{{y_{2}}}\right\rangle=\left|1\right\rangle$
and we get

${R_{1}}(\gamma)\left|{{y_{1}}{y_{2}}}\right\rangle\left|1\right\rangle={\alpha_{0}}{\beta_{0}}\left|{00}\right\rangle\left|1\right\rangle+{\alpha_{0}}{\beta_{1}}\left|{01}\right\rangle\left|1\right\rangle+{\alpha_{1}}{\beta_{0}}\left|{10}\right\rangle\left|1\right\rangle+{\alpha_{1}}{\beta_{1}}\left|{11}\right\rangle\otimes{R}(\gamma)\left|1\right\rangle$
(9)

After that, ${R_{2}}(\gamma)$ gate comes into play in which it controls only
one qubit $\left|{{y_{1}}}\right\rangle$ and acts on the third qubit as its
target, resulting in:

$\begin{array}[]{l}{R_{2}}(\gamma){R_{1}}(\gamma)\left|{{y_{1}}{y_{2}}}\right\rangle\left|1\right\rangle\\\
={\alpha_{0}}{\beta_{0}}\left|{00}\right\rangle\left|1\right\rangle+{\alpha_{0}}{\beta_{1}}\left|{01}\right\rangle\left|1\right\rangle+{\alpha_{1}}{\beta_{0}}\left|{10}\right\rangle\otimes{R}(\gamma)\left|1\right\rangle+{\alpha_{1}}{\beta_{1}}\left|{11}\right\rangle\otimes{R}(\gamma){R}(\gamma)\left|1\right\rangle\\\
=({\alpha_{0}}{\beta_{0}}\left|{00}\right\rangle+{\alpha_{0}}{\beta_{1}}\left|{01}\right\rangle+{e^{-i\gamma}}{\alpha_{1}}{\beta_{0}}\left|{10}\right\rangle+{e^{-2i\gamma}}{\alpha_{1}}{\beta_{1}}\left|{11}\right\rangle)\left|1\right\rangle\end{array}$
(10)

Note that here, based on the simple nature of our explicit constraints, the
circuit implementation is straightforward. In general, more complicated
constraints require us to apply more complicated gates.

Note ${e^{-ik\gamma}}{a_{z}}\left|z\right\rangle$ is equal to
${a_{z}}\left|z\right\rangle$ up to the phase factor ${e^{-ik\gamma}}$, and
the output of the measurement for these two states are the same. In fact, the
unitary operator $U(C,\gamma)$ is diagonal in the computational basis and does
not push the probability around. The origin of this diagonality comes from
selecting a specific type of optimisation problem (i.e., CSP) in which the
computational basis is taken to be the eigenbasis of the constraint operators.
In order to save and see the effect of this phase change related to
$U(C,\gamma)$ in our measurement, we needed to introduce a rotation operator
that takes $\left|0\right\rangle$ to a combination of $\left|0\right\rangle$
and $\left|1\right\rangle$; inevitably, this operator will increase the
degrees of freedom by one at the lowest depth. A convenient choice of such
operator can be given by:

$U(B,\beta)={e^{-i\beta
B}}=\prod\limits_{i=1}^{n}{{e^{-i\beta{\sigma_{i}}^{x}}}}$ (11)

With operator $B$ defined as

$B=\sum\limits_{i=1}^{n}{{\sigma_{i}}^{x}}={\sigma_{1}}^{x}\otimes{I^{\otimes(n-1)}}+{I_{1}}\otimes{\sigma_{2}}^{x}\otimes{I^{\otimes(n-2)}}+...+{I^{\otimes(n-1)}}\otimes{\sigma_{n}}^{x}$
(12)

Figure 2: $U(B,\beta)$ implementation with depth 1

Again, as $[{\sigma_{i}}^{x}\otimes{I_{j}},{I_{i}}\otimes{\sigma_{j}}^{x}]=0$,
$U(B,\beta)$ can be written as the product of exponentials.

Let us see the effect of $U(B,\beta)$ on an explicit example of one qubit.
Consider a qubit in its superposition,
$\left|y\right\rangle=a\left|0\right\rangle+b\left|1\right\rangle$, and
suppose $U(C,\gamma)$ delivers a phase factor ${e^{-i\gamma}}$ to base
$\left|0\right\rangle$. Now, applying ${e^{-i\beta{\sigma^{x}}}}$ and
performing measurements, we can obtain a probability of
${a^{2}}{\cos^{2}}(\beta)+{b^{2}}{\sin^{2}}(\beta)+ab\sin(2\beta)\sin(\gamma)$
to get $\left|0\right\rangle$, and we can see the contribution of the phase,
related to $U(C,\gamma)$, as desired.

By acting the operator $U(B,{\beta_{\alpha}})U(C,{\gamma_{\alpha}})$ $p$ times
with different ${\gamma_{\alpha}}$ and ${\beta_{\alpha}}$, we arrive at the
below $p$-layered variational ansatz:

$\left|{\vec{\gamma},\vec{\beta}}\right\rangle=U(B,{\beta_{p}})U(C,{\gamma_{p}})...U(B,{\beta_{1}})U(C,{\gamma_{1}})\left|s\right\rangle$
(13)

Where
$\left|s\right\rangle=\frac{1}{{{2^{n/2}}}}\sum\limits_{z}{\left|z\right\rangle}$
is the uniform superposition state. Now consider the expectation of C in this
quantum state

$F(\vec{\gamma},\vec{\beta})=\left\langle{\vec{\gamma},\vec{\beta}}\right|C\left|{\vec{\gamma},\vec{\beta}}\right\rangle$
(14)

Figure 3: The QAOA circuit diagram.

$\begin{array}[]{l}\left\langle{\vec{\gamma},\vec{\beta}}\right|C\left|{\vec{\gamma},\vec{\beta}}\right\rangle=\left\langle
s\right|{U^{\dagger}}(C,{\gamma_{1}})...{U^{\dagger}}(B,{\beta_{p}}).C.U(B,{\beta_{p}})...U(C,{\gamma_{1}})\left|s\right\rangle\\\
=\frac{1}{{{2^{n}}}}{\left|{{f_{1}}(\vec{\gamma},\vec{\beta})}\right|^{2}}P(\left|{00...0}\right\rangle)+...+\frac{1}{{{2^{n}}}}{\left|{{f_{n}}(\vec{\gamma},\vec{\beta})}\right|^{2}}P(\left|{11...1}\right\rangle)\end{array}$
(15)

As one can see, here, the eigenvalues of the objective function $C$ act as
markers (eigenvalues of the solutions will be higher in numerical value), and
by just demanding the maximum of the expectation value, the coefficients of
the solutions ${\left|{{f_{i}}(\vec{\gamma},\vec{\beta})}\right|^{2}}$ will be
maximised for a given set of parameters ${\vec{\gamma},\vec{\beta}}$. Now, one
can call the quantum computer to create the
$\left|{\vec{\gamma},\vec{\beta}}\right\rangle$ and to obtain a string $z$.

The reason for the $p$-times iterations is that there is no promise that,
after
$\left|{{\gamma_{1}},{\beta_{1}}}\right\rangle=U(B,{\beta_{1}})U(C,{\gamma_{1}})\left|s\right\rangle$,
those coefficients related to the solutions in the above become maximised and
the quantum state $\left|{{\gamma_{1}},{\beta_{1}}}\right\rangle$ becomes
close enough to the solution $\left|{z^{\prime}}\right\rangle$. Inevitably, we
need to increase the degrees of freedom in the problem by $2p$ and hope that
$p$ grows slower than $n$, thus becoming able to maximise the expectation
value.

The workflow of the QAOA started with determining the $2p$ angles
$(\vec{\gamma},\vec{\beta})$ that maximise the expectation value $F$; then, we
constructed our quantum state $\left|{\vec{\gamma},\vec{\beta}}\right\rangle$
on the quantum computer and measured it on the computational basis to obtain
the value $\left|z\right\rangle$. Finally, we calculated $C(z)$ on a classical
computer. By repeating these steps, we obtained a distribution on $C(z)$,
whose mean is $F(\vec{\gamma},\vec{\beta})$.

One superior way for the optimisation of the $2p$ parameters
${\vec{\gamma},\vec{\beta}}$ is utilising the Gibbs objective function defined
by $-\log\left\langle{{e^{-\eta C}}}\right\rangle$ [7].

##  2 QAOA Implementation on MaxCut problems

The MaxCut problem, one of the 21 NP-complete problems identified by Karp, has
many applications in computer science and statistical physics. The MaxCut
problem is in the APX class, telling us that, unless P = NP, no polynomial
time approximation scheme is available for it, arbitrarily near to the optimal
solution.

Herein, we first set the stage for the MaxCut problem, and then, by the
explicit representation of its constraints as operators, we reviewed the
construction of its related quantum gates in the QAOA algorithm.

The MaxCut problem can be defined by $n$ vertices and $m$ edges. By
introducing a two-valued function that acts on the vertices,
$k(i)\in\\{+1,-1\\}$, a cut can be defined if $k(i)k(j)=-1$ for two vertices
of an edge.

Figure 4: The MaxCut problem on the butterfly and Moser spindle graphs.

Thus, each vertex can be assigned by either $+1$ or $-1$. The aim is to find a
configuration of assignments, leading to a maximum number of cuts in the
graph.

In order to solve such a problem on a quantum computer, we promote that the
$n$ vertices be $n$ qubits. Thus, all different ${2^{n}}$ assignments for the
vertices will correspond to the ${2^{n}}$ dimensional Hilbert space generated
by $n$ qubits.

Now, one can easily define the objective function in the form of a constrained
satisfaction problem (CSP).

$\begin{array}[]{l}C=\sum\limits_{<ij>}^{m}{{C_{<ij>}}}\\\
{C_{<ij>}}=\frac{1}{2}(+1-k(i)k(j))\end{array}$ (16)

By promoting the objective function as an operator, we can get:

$C\left|z\right\rangle=\sum\limits_{<ij>}^{m}{{C_{<ij>}}}\left|z\right\rangle=P(z)\left|z\right\rangle$
(17)

where ${C_{<ij>}}=\frac{1}{2}(+I-\sigma_{i}^{z}\otimes\sigma_{j}^{z})$.

Here, if $\left\langle{ij}\right\rangle$ is a cut, the related qubits will
land onto different bases, for example,
$\left|{{z_{i}}}\right\rangle=\left|0\right\rangle$ and
$\left|{{z_{j}}}\right\rangle=\left|1\right\rangle$, or vice versa. Now,
consider
${\sigma_{i}}^{z}\otimes{\sigma_{j}}^{z}\left|{01}\right\rangle=-\left|{01}\right\rangle$,
then we can get
${C_{<ij>}}\left|{{z_{i}}{z_{j}}}\right\rangle=\left|{{z_{i}}{z_{j}}}\right\rangle$.

Now, we are ready to construct $U(C,\gamma)$ for the MaxCut so that:

$U(C,\gamma)={e^{-i\gamma\sum\limits_{<ij>}{{C_{<ij>}}}}}=\prod\limits_{<ij>}{{e^{-i\frac{\gamma}{2}(+I-{\sigma_{i}}^{z}\otimes{\sigma_{j}}^{z})}}}$
(18)

Explicitly, on
$\left|\psi\right\rangle=\left|{{y_{1}}}\right\rangle\left|{{y_{2}}}\right\rangle$.

${e^{-i\frac{\gamma}{2}I\otimes
I}}{e^{+i\frac{\gamma}{2}{\sigma_{1}}^{z}\otimes{\sigma_{2}}^{z}}}\left|{{y_{1}}{y_{2}}}\right\rangle={\alpha_{0}}{\beta_{0}}\left|{00}\right\rangle+{\alpha_{0}}{\beta_{1}}{e^{-i\gamma}}\left|{01}\right\rangle+{\alpha_{1}}{\beta_{0}}{e^{-i\gamma}}\left|{10}\right\rangle+{\alpha_{1}}{\beta_{1}}\left|{11}\right\rangle$
(19)

The circuit for the above gate can be implemented easily with the help of one
controlled-phase gate and two ordinary phase gates with the below arrangement:

$\begin{array}[]{l}CR(-2\gamma)\left|{{y_{1}}{y_{2}}}\right\rangle={\alpha_{0}}{\beta_{0}}\left|{00}\right\rangle+{\alpha_{0}}{\beta_{1}}\left|{01}\right\rangle+{\alpha_{1}}{\beta_{0}}\left|{10}\right\rangle+{\alpha_{1}}{\beta_{1}}{e^{+i2\gamma}}\left|{11}\right\rangle\\\
(R(\gamma)\otimes
R(\gamma))CR(-2\gamma)\left|{{y_{1}}{y_{2}}}\right\rangle={\alpha_{0}}{\beta_{0}}\left|{00}\right\rangle+{\alpha_{0}}{\beta_{1}}{e^{-i\gamma}}\left|{01}\right\rangle+{\alpha_{1}}{\beta_{0}}{e^{-i\gamma}}\left|{10}\right\rangle+{\alpha_{1}}{\beta_{1}}\left|{11}\right\rangle\end{array}$
(20)

Now, the angle-dependent quantum state can be considered.

$\left|{\vec{\gamma},\vec{\beta}}\right\rangle=U(B,{\beta_{p}})U(C,{\gamma_{p}})...U(B,{\beta_{1}})U(C,{\gamma_{1}})\left|s\right\rangle$
(21)

In order to determine the angles
$\vec{\gamma}=({\gamma_{1}},...,{\gamma_{p}})$ and
$\vec{\beta}=({\beta_{1}},...,{\beta_{p}})$, we consider the expectation value

$F(\vec{\gamma},\vec{\beta})=\left\langle{\vec{\gamma},\vec{\beta}}\right|C\left|{\vec{\gamma},\vec{\beta}}\right\rangle=\sum\limits_{<ij>}{\left\langle
s\right|{U^{\dagger}}(C,{\gamma_{1}})...{U^{\dagger}}(B,{\beta_{p}}){C_{<ij>}}U(B,{\beta_{p}})...U(C,{\gamma_{1}})\left|s\right\rangle}$
(22)

Consider the operator associated with the edge $\left\langle{ij}\right\rangle$
for the case of $p=1$:

${U^{\dagger}}(C,{\gamma_{1}}){U^{\dagger}}(B,{\beta_{1}}){C_{<ij>}}U(B,{\beta_{1}})U(C,{\gamma_{1}})$
(23)

The factors in ${U}(B,{\beta_{1}})$ which do not involve qubits i or j will
commute with ${C_{<ij>}}$ and we have

${U^{\dagger}}(C,{\gamma_{1}}){e^{+i{\beta_{1}}(\sigma_{i}^{x}+\sigma_{j}^{x})}}{C_{<ij>}}{e^{-i{\beta_{1}}(\sigma_{i}^{x}+\sigma_{j}^{x})}}U(C,{\gamma_{1}})$
(24)

Moreover, the factors in $U(C,{\gamma_{1}})$ that do not involve qubits $i$ or
$j$ commute with ${e^{-i{\beta_{1}}(\sigma_{i}^{x}+\sigma_{j}^{x})}}$ and
obviously with ${C_{<ij>}}$, and then cancel out.

Thus, the operator in (23) only involves the edge
$\left\langle{ij}\right\rangle$ and the edges adjacent to
$\left\langle{ij}\right\rangle$, as well as the qubits on those edges.

Now, consider $p=2$, and we can obtain the operator:

${U^{\dagger}}(C,{\gamma_{1}}){U^{\dagger}}(B,{\beta_{1}}){U^{\dagger}}(C,{\gamma_{2}}){U^{\dagger}}(B,{\beta_{2}}){C_{<ij>}}U(B,{\beta_{2}})U(C,{\gamma_{2}})U(B,{\beta_{1}})U(C,{\gamma_{1}})$
(25)

Here, the factors in $U(B,{\beta_{2}})$ that do not involve qubits $i$ and $j$
commute with ${C_{<ij>}}$, and we can obtain:

${U^{\dagger}}(C,{\gamma_{1}}){U^{\dagger}}(B,{\beta_{1}}){U^{\dagger}}(C,{\gamma_{2}}){e^{+i{\beta_{2}}(\sigma_{i}^{x}+\sigma_{j}^{x})}}{C_{<ij>}}{e^{-i{\beta_{2}}(\sigma_{i}^{x}+\sigma_{j}^{x})}}U(C,{\gamma_{2}})U(B,{\beta_{1}})U(C,{\gamma_{1}})$
(26)

Now, the factors in $U(C,{\gamma_{2}})$ that will survive will be:

$\prod\limits_{<ik>,<jl>}{{e^{-i\frac{{{\gamma_{2}}}}{2}(-\sigma_{j}^{z}\sigma_{l}^{z}+I)}}{e^{-i\frac{{{\gamma_{2}}}}{2}(-\sigma_{i}^{z}\sigma_{k}^{z}+I)}}{e^{-i\frac{{{\gamma_{2}}}}{2}(-\sigma_{i}^{z}\sigma_{j}^{z}+I)}}}$
(27)

Obviously, the factors in $U(B,{\beta_{1}})$ that will survive will be
${e^{-i{\beta_{1}}(\sigma_{i}^{x}+\sigma_{j}^{x}+\sigma_{k}^{x}+\sigma_{l}^{x})}}$.

Finally, the factors in $U(C,{\gamma_{1}})$ that will survive will be:

$\prod\limits_{<ik>,<jl>,<km>,<\ln>}{{e^{-i\frac{{{\gamma_{1}}}}{2}(-\sigma_{l}^{z}\sigma_{n}^{z}+I)}}{e^{-i\frac{{{\gamma_{1}}}}{2}(-\sigma_{m}^{z}\sigma_{k}^{z}+I)}}{e^{-i\frac{{{\gamma_{1}}}}{2}(-\sigma_{j}^{z}\sigma_{l}^{z}+I)}}{e^{-i\frac{{{\gamma_{1}}}}{2}(-\sigma_{i}^{z}\sigma_{k}^{z}+I)}}{e^{-i\frac{{{\gamma_{1}}}}{2}(-\sigma_{i}^{z}\sigma_{j}^{z}+I)}}}$
(28)

Thus, for $p=2$, the operator only involves edges that are, at most, two steps
away from $\left\langle{ij}\right\rangle$.

Consequently, for any $p$, our operator only involves edges that are, at most,
$p$ steps away from $\left\langle{ij}\right\rangle$ [1].

Therefore, each term in Equation (22) depends only on the subgraph involving
qubits $i$ and $j$ and those qubits at a distance less than or equal to $p$
steps away from $\left\langle{ij}\right\rangle$.

Now, let us define the restricted objective function ${C_{S}}$ that lives on a
subgraph $S$ as below:

${C_{S}}=\sum\limits_{<ij>\in S}{{C_{<ij>}}}$ (29)

and the associated operators:

$\begin{array}[]{l}U({C_{S}},\gamma)={e^{-i\gamma{C_{S}}}}\\\
U({B_{S}},\beta)={e^{-i\beta{B_{S}}}}\end{array}$ (30)

Where ${B_{S}}=\sum\limits_{i\in S}{\sigma_{i}^{x}}$.

Each edge $\left\langle{ij}\right\rangle$, with its related subgraph,
generates the below contribution to $F(\vec{\gamma},\vec{\beta})$:

$\left\langle{s,S(i,j)}\right|{U^{\dagger}}({C_{S(i,j)}},{\gamma_{1}})...{U^{\dagger}}({B_{S(i,j)}},{\beta_{p}}){C_{<ij>}}U({B_{S(i,j)}},{\beta_{p}})...U({C_{S(i,j)}},{\gamma_{1}})\left|{s,S(i,j)}\right\rangle$
(31)

It is easy to see that the contribution to the expectation value for different
edges with isomorphic subgraphs are same. Thus, we can rewrite the expectation
value as a sum over subgraph types.

Therefore, we can define:

$\begin{array}[]{l}{f_{S}}(\vec{\gamma},\vec{\beta})=\langle
s,S(i,j)|{U^{\dagger}}({C_{S(i,j)}},{\gamma_{1}})...{U^{\dagger}}({B_{S(i,j)}},{\beta_{p}}){C_{<ij>}}U({B_{S(i,j)}},{\beta_{p}})\\\
...U({C_{S(i,j)}},{\gamma_{1}})\left|{s,S(i,j)}\right\rangle\end{array}$ (32)

And rewrite F as:

$F(\vec{\gamma},\vec{\beta})=\sum\limits_{S}{{m_{S}}}{f_{S}}(\vec{\gamma},\vec{\beta})$
(33)

where ${{m_{S}}}$ is the number of occurrences of the subgraph S [1].

Regarding to the expectation value for MaxCut problem, we have:

$\frac{m}{2}+\frac{{-1}}{2}\sum\limits_{<ij>}^{m}{Tr[{M_{d}}.\prod\limits_{<ij>}^{m}{{e^{+i\frac{\gamma}{2}(+I-\sigma_{i}^{z}\sigma_{j}^{z})}}}.{e^{+i\beta(\sigma_{i}^{x}+\sigma_{j}^{x})}}.(\sigma_{i}^{z}\sigma_{j}^{z}).{e^{-i\beta(\sigma_{i}^{x}+\sigma_{j}^{x})}}.\prod\limits_{<ij>}^{m}{{e^{-i\frac{\gamma}{2}(+I-\sigma_{i}^{z}\sigma_{j}^{z})}}}]}$
(34)

Where ${M_{d}}$ is the initial density matrix.

By considering $U(B,\beta)=\prod\limits_{i}^{n}{({\mathop{\rm
Cos}\nolimits}(}\beta)I-i{\mathop{\rm Sin}\nolimits}(\beta)\sigma_{i}^{x})$,
we can simplify the above trace as:

$\begin{array}[]{l}{\mathop{\rm Cos}\nolimits}(2\beta){\mathop{\rm
Sin}\nolimits}(2\beta)Tr[{M_{d}}.\prod\limits_{<ij>}^{m}{{e^{+i\frac{\gamma}{2}(+I-\sigma_{i}^{z}\sigma_{j}^{z})}}}.(\sigma_{i}^{z}\sigma_{j}^{y}).\prod\limits_{<ij>}^{m}{{e^{-i\frac{\gamma}{2}(+I-\sigma_{i}^{z}\sigma_{j}^{z})}}}]\\\
+{\mathop{\rm Cos}\nolimits}(2\beta){\mathop{\rm
Sin}\nolimits}(2\beta)Tr[{M_{d}}.\prod\limits_{<ij>}^{m}{{e^{+i\frac{\gamma}{2}(+I-\sigma_{i}^{z}\sigma_{j}^{z})}}}.(\sigma_{i}^{y}\sigma_{j}^{z}).\prod\limits_{<ij>}^{m}{{e^{-i\frac{\gamma}{2}(+I-\sigma_{i}^{z}\sigma_{j}^{z})}}}]\\\
+{\mathop{\rm
Sin}\nolimits}{(2\beta)^{2}}Tr[{M_{d}}.\prod\limits_{<ij>}^{m}{{e^{+i\frac{\gamma}{2}(+I-\sigma_{i}^{z}\sigma_{j}^{z})}}}.(\sigma_{i}^{y}\sigma_{j}^{y}).\prod\limits_{<ij>}^{m}{{e^{-i\frac{\gamma}{2}(+I-\sigma_{i}^{z}\sigma_{j}^{z})}}}]\end{array}$
(35)

Moreover, one can derive an analytical expression for the expectation value by
further decomposing and simplifying these traces.

##  3 QAOA Implementation on QUBOs and Ising Type Hamiltonians

One of the main applications of quantum computers is solving NP-hard
combinatorial optimisation problems efficiently. With classical computation
paradigms, in most of the cases, all we have are heuristic algorithms with
which to tackle these computationally intensive problems. Therefore, it worth
seeking quantum algorithms with the ability to find higher-quality solutions
for such high-valued problems. A wide range of real-world optimisation
problems can be recast into higher order binary optimisation (HOBO) and
quadratic unconstrained binary optimisation problems (QUBOs) [8, 9]. In the
case of HOBO, the natural prescription is to transform a higher-order problem
into a quadratic one (${\rm{quadratisation}}$)[10, 11]. Moreover, for solving
very large versions of this problem, it has been shown that decomposing large
QUBO instances into mini-QUBOs and merging the sub-solutions can provide an
effective solution [12, 13, 14]. QUBOs can be connected to Ising Hamiltonian
via linear transformations. Here, we apply the QAOA to a Ising Hamiltonian.
This will pave the way for solving QUBOs on gate-based quantum computers.

The objective function in QUBO formulation can be written as:

${H_{QUBO}}=\sum\limits_{i<j}{{J_{ij}}}{x_{i}}{x_{j}}+\sum\limits_{i}{{J_{i}}}{x_{i}}=\sum\limits_{i\leq
j}{{x_{i}}}{J_{ij}}{x_{j}}=\left\langle{\vec{x},J\vec{x}}\right\rangle$ (36)

The variables ${x_{i}}$ live in the Boolean space $\\{0,1\\}$. Usually, the
${{J_{ij}}}$ matrix is written in symmetric or in upper triangle form.

In the Ising model, the objective function can be written as:

${H_{Ising}}=\sum\limits_{i<j}{{J_{ij}}}{s_{i}}{s_{j}}+\sum\limits_{i}{{J_{i}}}{s_{i}}$
(37)

Here the variables ${s_{i}}$ live in the Ising space $\\{-1,1\\}$.

Translation between Ising and QUBO representations can be done with simple
mapping, $s=2x-1$, where 1 is the vector all of whose components are 1.

With the above linear transformation, one can establish a simple relationship
between the two objective functions of QUBO and Ising.

$\left\langle{\vec{x},J\vec{x}}\right\rangle=\frac{1}{4}(\left\langle{\vec{1},J\vec{1}}\right\rangle+\left\langle{\vec{s},J\vec{s}}\right\rangle+2\left\langle{\vec{1},J\vec{s}}\right\rangle)$
(38)

Now, by promoting binary variables ${s_{i}}$ to Pauli matrices
${\sigma_{i}}^{z}$, we can achieve a quantum Ising Hamiltonian as:

$H=\sum\limits_{i<j}{{J_{ij}}}{\sigma_{i}}^{z}{\sigma_{j}}^{z}+\sum\limits_{i=1}^{n}{{J_{i}}}{\sigma_{i}}^{z}$
(39)

where it is understood that:

${\sigma_{i}}^{z}{\sigma_{j}}^{z}={{\rm I}_{1}}\otimes...\otimes{{\rm
I}_{i-1}}\otimes{\sigma_{i}}^{z}\otimes{{\rm I}_{i+1}}\otimes...\otimes{{\rm
I}_{j-1}}\otimes{\sigma_{j}}^{z}\otimes{{\rm I}_{j+1}}\otimes...\otimes{{\rm
I}_{n}}$ (40)

We can construct $U(H,\gamma)$, so that:

$U(H,\gamma)={e^{-i\gamma
H}}=\prod\limits_{i<j}{{e^{-i\gamma{J_{ij}}{\sigma_{i}}^{z}{\sigma_{j}}^{z}}}}\prod\limits_{i}{{e^{-i\gamma{J_{i}}{\sigma_{i}}^{z}}}}$
(41)

The quantum gate $U(B,\beta)$ should be the same as (11). Now, we can
construct the variational ansatz as:

$\left|{\vec{\gamma},\vec{\beta}}\right\rangle=U(B,{\beta_{p}})U(H,{\gamma_{p}})...U(B,{\beta_{1}})U(H,{\gamma_{1}}){\left|+\right\rangle^{\otimes
n}}$ (42)

Where ${\left|+\right\rangle^{\otimes n}}$ is the uniform superposition, can
be constructed via applying the Hadamard gates
${H^{1}}\otimes...\otimes{H^{n}}$ to ${\left|0\right\rangle^{\otimes n}}$ and
note $\vec{\gamma}=({\gamma_{1}},...,{\gamma_{p}})$ and
$\vec{\beta}=({\beta_{1}},...,{\beta_{p}})$.

In order to run the quantum algorithm, ${\vec{\gamma}}$ and ${\vec{\beta}}$
need to be determined by minimisation (or maximisation) of the expectation
value of the Hamiltonian as below:

$H(\vec{\gamma},\vec{\beta})=\left\langle{\vec{\gamma},\vec{\beta}}\right|H\left|{\vec{\gamma},\vec{\beta}}\right\rangle$

For the simplest case, $p=1$, we arrive at:

$\left\langle H\right\rangle=\left\langle
s\right|{U^{\dagger}}(H,{\gamma_{1}}){U^{\dagger}}(B,{\beta_{1}})(\sum\limits_{i<j}{{J_{ij}}}{\sigma_{i}}^{z}{\sigma_{j}}^{z}+\sum\limits_{i=1}^{n}{{J_{i}}}{\sigma_{i}}^{z})U(B,{\beta_{1}})U(H,{\gamma_{1}})\left|s\right\rangle$
(43)

For $p=1$ and for a triangle-free graph, one can check that the below factors
will survive:

$\begin{array}[]{l}\left\langle
H\right\rangle=\left\langle{{H_{1}}}\right\rangle+\left\langle{{H_{2}}}\right\rangle\\\
=\sum\limits_{i}{{}^{n}\left\langle+\right|}\prod\limits_{(i,j)\in
E}{{e^{+i\gamma{J_{ij}}{\sigma_{i}}^{z}{\sigma_{j}}^{z}}}}{e^{+i\gamma{J_{i}}{\sigma_{i}}^{z}}}{e^{+i\beta{\sigma_{i}}^{x}}}({J_{i}}{\sigma_{i}}^{z}){e^{-i\beta{\sigma_{i}}^{x}}}{e^{-i\gamma{J_{i}}{\sigma_{i}}^{z}}}\prod\limits_{(i,j)\in
E}{{e^{-i\gamma{J_{ij}}{\sigma_{i}}^{z}{\sigma_{j}}^{z}}}}{\left|+\right\rangle^{n}}\\\
+\sum\limits_{(i,j)}{{}^{n}\left\langle+\right|}(\prod\limits_{(j,l)\in
E}{{e^{+i\gamma{J_{jl}}{\sigma_{j}}^{z}{\sigma_{l}}^{z}}})}(\prod\limits_{(i,k)\in
E}{{e^{+i\gamma{J_{ik}}{\sigma_{i}}^{z}{\sigma_{k}}^{z}}}})({e^{+i\gamma{J_{ij}}{\sigma_{i}}^{z}{\sigma_{j}}^{z}}})(\prod\limits_{i,j}{{e^{+i\gamma{J_{i}}{\sigma_{i}}^{z}}})(}(\prod\limits_{i,j}{{e^{+i\beta{\sigma_{i}}^{x}}})}({J_{ij}}{\sigma_{i}}^{z}{\sigma_{j}}^{z})\\\
(\prod\limits_{i,j}{{e^{-i\beta{\sigma_{i}}^{x}}})}(\prod\limits_{i,j}{{e^{-i\gamma{J_{i}}{\sigma_{i}}^{z}}})(}{e^{-i\gamma{J_{ij}}{\sigma_{i}}^{z}{\sigma_{j}}^{z}}})(\prod\limits_{(i,k)\in
E}{{e^{-i\gamma{J_{ik}}{\sigma_{i}}^{z}{\sigma_{k}}^{z}}}})(\prod\limits_{(j,l)\in
E}{{e^{-i\gamma{J_{jl}}{\sigma_{j}}^{z}{\sigma_{l}}^{z}}})}{\left|+\right\rangle^{n}}\end{array}$
(44)

The expectation value can be achieved analytically as:

$\begin{array}[]{l}\left\langle
H\right\rangle=\sum\limits_{i}{{J_{i}}}\sin(2\beta)\sin(2\gamma{J_{i}})\prod\limits_{(i,j)\in
E}{\cos(2\gamma{J_{ij}})}\\\
+\sum\limits_{(i,j)}{\frac{1}{2}{{\rm{J}}_{ij}}\sin(\beta)\cos(\beta)(2\cos(2\beta)\sin(2\gamma{{\rm{J}}_{ij}})(}2\cos(2\gamma{{\rm{J}}_{i}})\prod\limits_{(i,k)}{\cos(2\gamma{{\rm{J}}_{ik}})}\\\
+2\cos(2\gamma{{\rm{J}}_{\rm{j}}})\prod\limits_{(j,l)}{\cos(2\gamma{{\rm{J}}_{jl}})})\\\
+4\sin(2\beta)\sin(2\gamma{{\rm{J}}_{i}})\sin(2\gamma{{\rm{J}}_{j}})\prod\limits_{(i,k)}{\cos(2\gamma{J_{ik}})}\prod\limits_{(j,l)}{\cos(2\gamma{{\rm{J}}_{jl}})})\end{array}$
(45)

##  4 Quantum Algorithm Implementations

Here, we review the implementation of the QAOA on some explicit examples
(MaxCut on butterfly and Moser spindle graphs). In order to find the optimal
angles $\vec{\gamma}$ and $\vec{\beta}$, especially for cases of $p>1$,
different prescriptions have been offered, such as gradient descent and basin-
hopping [15, 16]. Such treatments on the QAOA need to be automated to have
ready-to-run QAOA jobs at peak performance.

Figure 5: QAOA circuit implementation for MaxCut problem on butterfly graph

Returning to the MaxCut problem on a butterfly graph, for $p=1$, the
expectation value is:

$F(\gamma,\beta)=\left\langle
s\right|{U^{\dagger}}(C,\gamma){U^{\dagger}}(B,\beta).C.U(B,\beta)U(C,\gamma)\left|s\right\rangle$
(46)

Based on the simplification procedure discussed in the previous section, we
can write $F(\gamma,\beta)=2{f_{A}}+4{f_{B}}$, where ${f_{A}}$ and ${f_{B}}$
are the expectation values for one edge in the sub-graphs $\\{(0,1),(3,4)\\}$
and $\\{(0,2),(1,2),(3,2),(4,2)\\}$, respectively.

$\begin{array}[]{l}{f_{A}}=\frac{1}{2}(1-\left\langle{{+^{3}}}\right|U_{12}^{\dagger}(\gamma)U_{02}^{\dagger}(\gamma)U_{01}^{\dagger}(\gamma)U_{1}^{\dagger}(\beta)U_{0}^{\dagger}(\beta)({\sigma_{0}}^{z}\otimes{\sigma_{1}}^{z})\\\
{U_{0}}(\beta){U_{1}}(\beta){U_{01}}(\gamma){U_{02}}(\gamma){U_{12}}(\gamma)\left|{{+^{3}}}\right\rangle)\end{array}$
(47)

${f_{A}}=\frac{1}{4}(\sin(2\beta)\sin(\gamma)(\cos(2\beta-\gamma)+3\cos(2\beta+\gamma))+2)$
(48)

$\begin{array}[]{l}{f_{B}}=\frac{1}{2}(1-\left\langle{{+^{5}}}\right|U_{24}^{\dagger}(\gamma)U_{23}^{\dagger}(\gamma)U_{12}^{\dagger}(\gamma)U_{02}^{\dagger}(\gamma)U_{01}^{\dagger}(\gamma)U_{2}^{\dagger}(\beta)U_{0}^{\dagger}(\beta)({\sigma_{0}}^{z}\otimes{\sigma_{2}}^{z})\\\
{U_{0}}(\beta){U_{2}}(\beta){U_{01}}(\gamma){U_{02}}(\gamma){U_{12}}(\gamma){U_{23}}(\gamma){U_{24}}(\gamma)\left|{{+^{5}}}\right\rangle)\end{array}$
(49)

${f_{B}}=\frac{1}{8}(\sin(2\beta)\sin(2\gamma)(\cos(2(\beta+\gamma))+3\cos(2\beta))+4)$
(50)

The six gates of form $U({C_{ij}},\gamma)$ can be written as below using phase
shift gates $R(\gamma)$ and controlled phase shift gates $CR(-2\gamma)$.

$U({C_{ij}},\gamma)={e^{-i\frac{\gamma}{2}(-\sigma_{i}^{z}\otimes\sigma_{j}^{z}+I)}}=({R_{i}}(\gamma)\otimes{R_{j}}(\gamma)).{C_{ij}}R(-2\gamma)$
(51)

Classical Preprocessing to Determine $\gamma$ and $\beta$

    #In Python, set up your Wolfram language session:
    from wolframclient.language import wlexpr
    from wolframclient.evaluation import WolframLanguageSession
    session = WolframLanguageSession()
    #Evaluate any Wolfram Language code from Python:
    Angles=session.evaluate(wlexpr(’NArgMax[1/ 8 (21 +
    3 Cos[4 \[Beta]] + 4 Cos[4 \[Beta] - 2 \[Gamma]] +
    2 Cos[2 \[Gamma]] + Cos[4 \[Gamma]] -  Cos[4 (\[Beta]
    + \[Gamma])] -  6 Cos[2 (2 \[Beta] + \[Gamma])]),
    {\[Gamma], \[Beta]} \[Element]
    DiscretizeRegion[Disk[{0, 0}, {2 Pi, Pi}, {0, Pi/2}]]]’))

    Gamma=Angles[0]
    Beta=Angles[1]

    # importing Qiskit
    from qiskit import QuantumRegister, ClassicalRegister, QuantumCircuit
    from qiskit import IBMQ, execute
    from qiskit.tools.monitor import job_monitor
    from qiskit.visualization import plot_histogram

Figure 6: Running the QAOA on ibmq-16-melbourne device

    # prepare the butterfly graph [17]
    import numpy as np
    m = 5
    Nodes = np.arange(0,m,1)
    Edges =[(0,1),(0,2),(1,2),(3,2),(3,4),(4,2)]
    # prepare the quantum and classical resisters
    QAlgo = QuantumCircuit(len(Nodes), len(Nodes))
    # apply the Hadamard gates to all qubits to have uniform superposition state
    QAlgo.h(range(len(Nodes)))
    QAlgo.barrier()
    # to implement U(C,gamma), apply below equivalent circuit [17]
    # note here u1 is the phase shift gate
    for edge in Edges:
        i = edge[0]
        j = edge[1]
        QAlgo.cu1(-2*Gamma, i, j)
        QAlgo.u1(Gamma, i)
        QAlgo.u1(Gamma, j)

    # to implement U(B,beta), apply R_x(2*beta) to all qubits
    QAlgo.barrier()
    QAlgo.rx(2*Beta, range(len(Nodes)))
    # Measure the result in the computational basis
    QAlgo.barrier()
    QAlgo.measure(range(len(Nodes)),range(len(Nodes)))

    # Use the ibmq_16_melbourne device
    from qiskit.tools.monitor import job_monitor
    provider = IBMQ.load_account()
    backend = provider.get_backend(’ibmq_16_melbourne’)
    shots = 2048
    job = execute(QAlgo, backend=backend, shots=shots)
    job_monitor(job)
    results = job.result()
    plot_histogram(results.get_counts(),bar_labels = False,figsize = (12,10))

For the Moser spindle graph, one can directly compute the single-layer
expectation value and the related angles from the below in Wolfram language:

    MaxCut[x_, y_] := {NE = {Range[x], y}, nodes = NE[[1]],
    edges = NE[[2]], cons = Range[Length[edges]],
    Sig[i_, j_] := ReplacePart[Table[IdentityMatrix[2],
    {k, Length[NE[[1]]]}], {i -> PauliMatrix[3], j -> PauliMatrix[3]}],
    Do[cons[[i]] = Sig[edges[[i]][[1]], edges[[i]][[2]]],
    {i, Length[edges]}],consyy = Range[Length[edges]],Sigyy[i_, j_] :=
    ReplacePart[Table[IdentityMatrix[2], {k, Length[NE[[1]]]}],
    {i -> PauliMatrix[2], j -> PauliMatrix[2]}],
    Do[consyy[[i]] = Sigyy[edges[[i]][[1]], edges[[i]][[2]]],
    {i, Length[edges]}],consizyj = Range[Length[edges]],Sigizyj[i_, j_] :=
    ReplacePart[Table[IdentityMatrix[2], {k, Length[NE[[1]]]}], {i ->
    PauliMatrix[3], j -> PauliMatrix[2]}],
    Do[consizyj[[i]] = Sigizyj[edges[[i]][[1]], edges[[i]][[2]]],
    {i, Length[edges]}],consiyzj = Range[Length[edges]],
    Sigiyzj[i_, j_] := ReplacePart[Table[IdentityMatrix[2],
    {k, Length[NE[[1]]]}], {i -> PauliMatrix[2], j -> PauliMatrix[3]}],
    Do[consiyzj[[i]] = Sigiyzj[edges[[i]][[1]], edges[[i]][[2]]],
    {i, Length[edges]}],Conjugate[\[Gamma]] ^:= \[Gamma],
    Mapi[i_] := Function[k, Exp[-I \[Gamma] (1/2) k]] /@
    Diagonal[(-KroneckerProduct @@ cons[[i]] +
    IdentityMatrix[2^Length[nodes]])],
    UCgamma = DiagonalMatrix[Product[Mapi[i], {i, Length[edges]}]],
    DaggerUCgamma = ConjugateTranspose[UCgamma],
    reducedCijyy[i_] := KroneckerProduct @@ consyy[[i]],
    reducedCizyj[i_] := KroneckerProduct @@ consizyj[[i]],
    reducedCiyzj[i_] := KroneckerProduct @@ consiyzj[[i]],
    Do[reducedCijyy[i]; reducedCizyj[i]; reducedCiyzj[i],
    {i, Length[edges]}],Conjugate[\[Beta]] ^:= \[Beta],
    DenMat = Table[1/2^Length[nodes], {i, 1, 2^Length[nodes]},
    {j, 1, 2^Length[nodes]}],
    TRACE[i_] := Cos[2 \[Beta]] Sin[2 \[Beta]] Diagonal[
    DenMat.DaggerUCgamma.reducedCizyj[i]].Diagonal[UCgamma] +
    Cos[2 \[Beta]] Sin[2 \[Beta]] Diagonal[DenMat.DaggerUCgamma.
    reducedCiyzj[i]].Diagonal[UCgamma] +
    Sin[2 \[Beta]]^2 Diagonal[DenMat.DaggerUCgamma.
    reducedCijyy[i]].Diagonal[UCgamma],
    ParallelDo[TRACE[k], {k, Length[NE[[2]]]}],
    Expectval = FullSimplify[(1/2 Length[NE[[2]]] -
    1/2 ParallelSum[TRACE[i], {i, Length[NE[[2]]]}])],
    Print["The Expectation Value: ", Expectval],
    Angles = NArgMax[Expectval, {\[Gamma], \[Beta]} \[Element]
    Disk[{0, 0}, {2 Pi, Pi}]]}
    Ψ

    MaxCut[7,{{1,2},{1,3},{2,3},{4,3},{2,4},{4,5},{6,5},{5,7},{6,7},{1,6},{1,7}}];

In this example, the single layer expectation value is:

$\frac{1}{{16}}\left({88+8\cos{{\left[\gamma\right]}^{2}}\left({9+2\cos\left[\gamma\right]}\right)\sin\left[{4\beta}\right]\sin\left[\gamma\right]-8\left({2+\cos\left[\gamma\right]}\right)\sin{{\left[{2\beta}\right]}^{2}}\sin{{\left[{2\gamma}\right]}^{2}}}\right)$
(52)

## Acknowledgement

The author would like to thank Nike Dattani and Edward Farhi for their helpful
comments.

## References

  * [1] E. Farhi, J. Goldstone, and S. Gutmann, _A Quantum Approximate Optimization Algorithm_ , _arXiv Prepr. arXiv1411.4028, pages 1-16, 2014._
  * [2] E. Farhi, J. Goldstone, and S. Gutmann, _A Quantum Approximate Optimization Algorithm Applied to a Bounded Occurrence Constraint Problem_ , _arXiv Prepr. arXiv1411.4028, pages 1-13, 2014._
  * [3] B. Barak, A. Moitra, R. O’Donnell, P. Raghavendra, O. Regev et al, _Beating the random assignment on constraint satisfaction problems of bounded degree_ , _arXiv:1505.03424, 2015._
  * [4] E. Farhi, _Quantum Supremacy through the Quantum Approximate Optimization Algorithm_ , _arXiv Prepr.arXiv1602.07674, pages 1-22, 2016._
  * [5] S. Hadfield, Z. Wang, E.G. Rieffel, B. O’Gorman, D. Venturelli, and R. Biswas, _Quantum approximate optimization with hard and soft constraints_ , _in Proceedings of the Second International Workshop on Post Moores Era Supercomputing. ACM, pp. 15-21, 2017._
  * [6] S. Hadfield, Z. Wang, B. O’Gorman, E.G. Rieffel, D. Venturelli, and R. Biswas, _From the Quantum Approximate Optimization Algorithm to a Quantum Alternating Operator Ansatz_ , _arXiv Prepr. arXiv1709.03489, pages 1-46, 2017._
  * [7] L. Li, M. Fan, M. Coram, P. Riley and S. Leichenauer, _Quantum Optimization with a Novel Gibbs Objective Function and Ansatz Architecture Search_ , _Phys. Rev. Research 2, 023-074 , 2020._
  * [8] G. Kochenberger, J.-K. Hao, F. Glover, M. Lewis, Z. Lü, H. Wang, and Y. Wang, _The unconstrained binary quadratic programming problem: A survey_ , _Journal of Combinatorial Optimization, vol. 28, no. 1, pp. 58-81, 2014._
  * [9] A. Lucas, _Ising formulations of many NP problems_ , _ArXiv preprint arXiv:1302.5843, 2013._
  * [10] N. Dattani, _Quadratization in discrete optimization and quantum mechanics_ , _arXiv:1901.04405_ [quant-ph].
  * [11] N. Dattani, H. Chau, _All 4-variable functions can be perfectly quadratized with only 1 auxiliary variable_ , _arXiv:1910.13583_ [cs.DM].
  * [12] F. Glover, Z. Lü, and J.-K. Hao, _Diversification-driven tabu search for unconstrained binary quadratic problems_ , _4OR, vol. 8, no. 3, pp. 239-253, 2010._
  * [13] Y.Wang, Z.Lü, F.Glover, and J.-K. Hao, _Path relinking for unconstrained binary quadratic programming_ , _European Journal of Operational Research, vol. 223, no. 3, pp. 595-604, 2012._
  * [14] Y. Wang, Z. Lü, F. Glover, and J.-K. Hao, _A multilevel algorithm for large unconstrained binary quadratic optimization_ , _Integration of AI and OR Techniques in Constraint Programming for Combinatorial Optimization Problems, Springer, pp. 395-408. 2012._
  * [15] Z. Wang, S. Hadfield, Z. Jiang, and E. G. Rieffel, _Quantum approximate optimization algorithm for MaxCut: A fermionic view_ , _Physical Review A, vol. 97, p. 022304, 2018._
  * [16] Z. Wang, N. C. Rubin, J. M. Dominy, and E. G. Rieffel, _XY mixers: Analytical and numerical results for the quantum alternating operator ansatz_ , _Physical Review A, vol. 101, no. 1, p. 012320, 2020._
  * [17] https://qiskit.org/textbook/ch-applications/qaoa.html