* [81] A. Knyazev, Locally Optimal Block Preconditioned Conjugate Gradient, https://github.com/lobpcg/blopex, https://www.mathworks.com/matlabcentral/fileexchange/48-locally-optimal-block-preconditioned-conjugate-gradient, Accessed: 2022-07-21 (2022).
  * [82] A. Knyazev, Toward the optimal preconditioned eigensolver: Locally optimal block preconditioned conjugate gradient method, SIAM Journal on Scientific Computing 23 (2) (2001) 517–541.
  * [83] A. V. Knyazev, M. E. Argentati, I. Lashuk, E. Ovtchinnikov, Block locally optimal preconditioned eigenvalue xolvers (BLOPEX) in HYPRE and PETSc, SIAM Journal on Scientific Computing 29 (5) (2007) 2224–2239.
  * [84] J. A. Duersch, M. Shao, C. Yang, M. Gu, A robust and efficient implementation of lobpcg, SIAM Journal on Scientific Computing 40 (5) (2018) C655–C676.
  * [85] M. P. Teter, M. C. Payne, D. C. Allan, Solution of Schrödinger’s equation for large systems, Physical Review B 40 (1989) 12255–12263.
  * [86] A. Mayer, Band structure and transport properties of carbon nanotubes using a local pseudopotential and a transfer-matrix technique, Carbon 42 (10) (2004) 2057–2066.
  * [87] A. S. Banerjee, L. Lin, W. Hu, C. Yang, J. E. Pask, Chebyshev polynomial filtered subspace iteration in the Discontinuous Galerkin method for large-scale electronic structure calculations, arXiv preprint arXiv:1606.03416 (2016).
  * [88] A. S. Banerjee, L. Lin, P. Suryanarayana, C. Yang, J. E. Pask, Two-level chebyshev filter based complementary subspace method: pushing the envelope of large-scale electronic structure calculations, Journal of chemical theory and computation 14 (6) (2018) 2930–2946.
  * [89] Y. Zhou, J. R. Chelikowsky, Y. Saad, Chebyshev-filtered subspace iteration method free of sparse diagonalization for solving the Kohn–Sham equation, Journal of Computational Physics 274 (2014) 770–782.
  * [90] A. A. Laturia, M. L. Van de Put, W. G. Vandenberghe, Generation of empirical pseudopotentials for transport applications and their application to group iv materials, Journal of Applied Physics 128 (3) (2020) 034306.
  * [91] M. L. Van de Put, M. V. Fischetti, W. G. Vandenberghe, Scalable atomistic simulations of quantum electron transport using empirical pseudopotentials, Computer Physics Communications 244 (2019) 156–169.
  * [92] J. R. Chelikowsky, Introductory Quantum Mechanics with MATLAB: For Atoms, Molecules, Clusters, and Nanocrystals, John Wiley & Sons, 2019.
  * [93] A. Rohatgi, Webplotdigitizer: Version 4.5 (2021).
URL https://automeris.io/WebPlotDigitizer

  * [94] A. Mayer, J.-P. Vigneron, Accuracy-control techniques applied to stable transfer-matrix computations, Physical Review E 59 (4) (1999) 4659.
  * [95] J. Pendry, A. MacKinnon, Calculation of photon dispersion relations, Physical Review Letters 69 (19) (1992) 2772.
  * [96] H. Tamura, T. Ando, Conductance fluctuations in quantum wires, Physical Review B 44 (4) (1991) 1792.
  * [97] M. L. Van de Put, A. A. Laturia, M. V. Fischetti, W. G. Vandenberghe, Efficient modeling of electron transport with plane waves, in: 2018 International Conference on Simulation of Semiconductor Processes and Devices (SISPAD), 2018, pp. 71–74.
  * [98] N. Troullier, J. L. Martins, Efficient pseudopotentials for plane-wave calculations, Physical Review B 43 (3) (1991) 1993.
  * [99] W. Kohn, L. J. Sham, Self-consistent equations including exchange and correlation effects, Phys. Rev. 140 (1965) A1133–A1138.
  * [100] J. P. Perdew, Y. Wang, Accurate and simple analytic representation of the electron-gas correlation energy, Physical Review B 45 (1992) 13244–13249.
  * [101] D. Zhang, R. D. James, T. Dumitrica, Electromechanical characterization of carbon nanotubes in torsion via symmetry adapted tight-binding objective molecular dynamics, Physical Review B 80 (11) (2009) 115418.
  * [102] L. Yang, M. Anantram, J. Han, J. Lu, Band-gap change of carbon nanotubes: Effect of small uniaxial and torsional strain, Physical Review B 60 (19) (1999) 13874.
  * [103] L. Yang, J. Han, Electronic structure of deformed carbon nanotubes, Physical review letters 85 (1) (2000) 154.
  * [104] S. Agarwal, A. S. Banerjee, On the solution of the Poisson problem for twisted geometries, (In preparation) (2022).
  * [105] S. Pathrudkar, H. M. Yu, S. Ghosh, A. S. Banerjee, Machine learning based prediction of the electronic structure of quasi-one-dimensional materials under strain, Phys. Rev. B 105 (2022) 195141.
  * [106] Y. Hakobyan, E. Tadmor, R. James, Objective quasicontinuum approach for rod problems, Physical Review B 86 (24) (2012) 245435.
  * [107] C. Le Bris (Ed.), Computational Chemistry, Vol. X of Handbook of Numerical Analysis, North-Holland, 2003.
  * [108] I. Nikiforov, B. Hourahine, B. Aradi, T. Frauenheim, T. Dumitrică, Ewald summation on a helix: A route to self-consistent charge density-functional based tight-binding objective molecular dynamics, The Journal of Chemical Physics 139 (9) (2013) 094110.
  * [109] O. Liebman, S. Agarwal, A. Banerjee, Helical waves for self-consistent first principles calculations of chiral one-dimensional nanomaterials, Bulletin of the American Physical Society (2022).
  * [110] NIST Digital Library of Mathematical Functions, http://dlmf.nist.gov/, Release 1.1.5 of 2022-03-15, f. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds. (2022).
URL http://dlmf.nist.gov/

  * [111] L. J. Laslett, W. Lewish, Evaluation of the zeros of cross-product bessel functions, Mathematics of Computation 16 (78) (1962) 226–232.
  * [112] H. Hellmann, Hans Hellmann: Einführung in Die Quantenchemie: Mit Biografischen Notizen von Hans Hellmann Jr., Springer-Verlag, 2015.
  * [113] R. P. Feynman, Forces in molecules, Physical Review 56 (4) (1939) 340.
  * [114] A. Deaño, D. Huybrechs, A. Iserles, Computing highly oscillatory integrals, SIAM, 2017.
  * [115] Y. Ma, Y. Xu, Computing highly oscillatory integrals, Mathematics of Computation 87 (309) (2018) 309–345.
  * [116] A. Iserles, S. P. Nørsett, Efficient quadrature of highly oscillatory integrals using derivatives, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 461 (2057) (2005) 1383–1399.
  * [117] A. Iserles, S. P. Nørsett, On quadrature methods for highly oscillatory integrals and their implementation, BIT Numerical Mathematics 44 (4) (2004) 755–772.
  * [118] R. Chen, Numerical approximations to integrals with a highly oscillatory bessel kernel, Applied Numerical Mathematics 62 (5) (2012) 636–648.
  * [119] A. Ralston, P. Rabinowitz, A first course in numerical analysis, Courier Corporation, 2001.
  * [120] N. Hale, A. Townsend, Fast and accurate computation of gauss–legendre and gauss–jacobi quadrature nodes and weights, SIAM Journal on Scientific Computing 35 (2) (2013) A652–A674.
  * [121] D. Laurie, Calculation of gauss-kronrod quadrature rules, Mathematics of Computation 66 (219) (1997) 1133–1145.
  * [122] L. F. Shampine, Vectorized adaptive quadrature in matlab, Journal of Computational and Applied Mathematics 211 (2) (2008) 131–140.


# Performance of the wavelength-shifting fiber upgrade for the Mu2e cosmic-ray
veto detector

D. Coveyou  E.C. Dukes  R.C. Group  Y. Oksuzian  S. Roberts  M.
Solt${}^{1},$11footnotetext: Corresponding author.

###### Abstract

The Cosmic-Ray Veto detector for the muon-to-electron conversion experiment at
Fermilab consists of four plastic scintillating counter layers read out by
silicon photo-multipliers through embedded wavelength-shifting fibers. In
order to increase the light yield in the most critical regions of the Cosmic-
Ray Veto detector for improved background rejection, a 1.8 mm diameter fiber
is being used in many of the detector’s critical modules instead of the
previously planned 1.4 mm diameter fiber. This paper reports the testing
procedure and light properties of thirty-four 1.8 mm fiber spools, with
measurements performed using a custom-built scanner. We compare these new
results with previously published data from the 1.4 mm diameter fiber used for
regions of the cosmic-ray veto where the increased light yield is not
required. In addition, measurements of fiber aging were performed.

##  1 Introduction

The muon-to-electron conversion (Mu2e) experiment [1, 2, 3] will search for
charged lepton flavor violation (CLFV) through the neutrinoless conversion of
a muon into an electron in the presence of a nucleus. Muon-to-electron
conversion occurs because of neutrino oscillations, but at an unobservable
rate (approximately forty orders of magnitude below the current experimental
limits). However, many extensions to the Standard Model of particle physics
predict that this process could occur at an observable rate [4]. An
observation of this conversion process would be an unambiguous sign of physics
beyond the Standard Model.

To search for CLFV with the desired sensitivity, Mu2e must keep the total
background to fewer than one event over the full multi-year experimental run.
However, without a cosmic-ray veto system, approximately one background event
will be produced per day at Mu2e from cosmic-ray muons. Thus, for the desired
background suppression, Mu2e needs to reduce the cosmic-ray background by four
orders of magnitude. By surrounding much of the Mu2e apparatus, the Cosmic-Ray
Veto (CRV) is designed to veto cosmic-ray backgrounds with an overall 99.99%
efficiency.

The CRV consists of plastic scintillator counters up to 7 m long containing
embedded wavelength-shifting (WLS) fibers that are read out by silicon photo-
multipliers (SiPMs). The CRV efficiency critically depends on the performance
of the fiber. Therefore, a WLS fiber scanner was designed and fabricated to
measure fiber light yield properties and ensure the CRV’s fiber quality
requirements. More information about the fiber scanner and 1.4 mm production
fiber performance is given in previous publications [5, 6].

In order to improve cosmic ray background rejection efficiency, larger 1.8 mm
diameter fiber will be used for the most critical modules of the CRV. Using
the Fermilab test beam facility, the light-yield improvement of counters
outfitted with 1.8 mm diameter fiber was measured to be 24% [7]. This light
yield increase is consistent with a simple geometric argument that a fiber
should absorb an amount of light proportional to its diameter ($\frac{1.8\
\mathrm{mm}}{1.4\ \mathrm{mm}}\sim 29\%$). To ensure sufficient light yield
for the CRV’s sectors with the largest flux of cosmic ray muons, 14 km of
larger diameter fiber was purchased. This paper reviews important fiber
characteristics, describes the fiber scanner, and presents studies measuring
the 1.8 mm diameter fiber quality, including preliminary fiber aging studies.

##  2 Wavelength-shifting fiber

The CRV uses Kuraray Y11 non-S-type WLS fibers [8] containing fluorescent dye,
K27, that absorbs blue light (375-475 nm) from scintillating counters and
emits light in the green (450-600 nm) spectral region. These fiber signals are
read out by $2.0{\times}2.0$ mm2 (model S13360-2050VE with a pixel size of 50
$\mu$m) Hamamatsu silicon photo-multipliers (SiPMs) [9]. The chosen SiPM
spectral photon detection efficiency response is well-matched to the emission
spectrum for Kuraray Y11 fiber.222 SiPMs in the CRV feature a peak sensitivity
at 450 nm, with 80% (50%) peak sensitivity at 550 nm (620 nm) [10]. Multi-clad
non-S-type fiber was selected for the CRV due to its enhanced light yield and
longer attenuation length [8].

##  3 Fiber scanner

Figure 1: A schematic (left) and picture (right) of the fiber scanner.

A fiber scanner (figure 1) was designed and fabricated333Fiber scanner design
is by D. Shooltz of Shooltz Solutions LLC. to study fiber quality for the CRV
detector. It measures the WLS fiber light yield as a function of wavelength by
varying the distances between a light source and a receiver (either a
spectrophotometer or a photodiode). It consists of a 62 cm diameter take-up
drum holding the optical readout electronics, to which the fiber end is
connected. The first 25 m of fiber is wound onto the drum’s grooved outer rim.
A blue light-emitting diode (LED) source [11] not affixed to the rotating
take-up drum excites the fiber, and the resulting light is captured and
readout by a large-area Hamamatsu S1227-1010BR photodiode [12] or a STS-VIS
USB Ocean Optics spectrophotometer [13], depending on the desired measurement.
The take-up drum stops at regular intervals to measure the response from the
blue LED source at different distances along the fiber, thus providing light
yield measurements as a function of wavelength and distance from the readout
end. Sample light yield responses from the photodiode and spectrometer are
presented in the following.

The fiber is delivered and stored on 90 cm diameter cardboard spools
containing 0.5 km of fiber. The fiber remains on the spool during the
measurement, which takes place in a blue-light filtering dark room to minimize
fiber damage from harmful UV light. The fiber end is epoxied into a ferrule,
polished with a diamond-bit flycutter, then connected to either the photodiode
or spectrometer port affixed on the rotating take-up drum. A stepping motor
drives the take-up drum with controlled gentle acceleration and deceleration.
The drum is moved in 50 cm increments over the first 25 m of fiber for each
spool, and light spectra and intensity measurements are performed at each
point. These measurements are assumed to be sufficient to determine the
quality of the entire spool. The fiber scan procedure is controlled remotely
via WiFi through a web interface on a Raspberry Pi [16] mounted on the take-up
drum. During data collection, the Raspberry Pi controls the stepper motor and
gathers and stores the photodiode or spectrometer data. The scanner is
described in more detail in reference [5].

##  4 Results

###  4.1 1.8 mm production fiber

The CRV detector consists of 5344 counters and 52.3 km of WLS fiber. The
original 1.4 mm diameter fiber was delivered on 104 spools in May 2018 and all
spools were tested with the fiber scanner by the end of June 2018. For
improved background rejection, the fiber is being upgraded to larger 1.8 mm
diameter fiber in the most critical regions of the CRV. The fabrication of the
modules in these regions was scheduled at the end of production so they could
be built with maximum knowledge and fabrication experience. This larger
diameter fiber was delivered on 34 spools in May 2021 and all spools were
tested with the fiber scanner by mid-June 2021.444Fiber diameters were also
measured for all 34 1.8 mm diameter fiber spools as part of the quality
control procedure. A mean with standard error of $1.801\pm 0.002$ mm was
obtained.

The 1.8 mm production spools light yield is shown as a function of the length
of the fiber between the photo-cathode and the LED source in figure 2, with
each curve representing a different fiber spool. This figure also includes a
1.8 mm diameter fiber spool from 2018 that was re-measured in 2021 (referred
to as the benchmark spool). In figure 3 each spools light yield is normalized
to a 1.4 mm diameter fiber spool from the initial 2018 fiber batch, which was
then re-tested in 2021. On average, the production spools light yield response
shows a 60%-70% increase over the original 1.4 mm fiber and a 30% increase
with respect to the benchmark spool. This unexpected increase in light yield
is not observed in measurements of fully constructed di-counters (as reported
section 4.3). The light yield increases can be explained by a combination of
aging effects (discussed in section 4.4) and an increase in fiber diameter.

Figure 2: Light yield as a function of the length of the fiber between the
source and readout for all 1.8 mm diameter spools (left) and their average
(right). The lowest black line on the left is the 1.8 mm benchmark spool.
Measurements are taken in 50 cm increments over the first 25 m of each fiber
spool. All spools achieved the desired light yield performance and showed a
significant improvement over the original 1.4 mm spools. All spools
demonstrate the typical rapid fall in light yield at short distances followed
by a less rapid decline. Figure 3: Light yield as a function of the distance
from the LED source for all thirty-four 1.8 mm spools shown as a ratio to an
original 1.4 mm fiber manufactured in 2017 and re-tested in June 2021. All
tested fibers exceeded the 1.4 mm fiber’s light yield performance. The
measurements of the two spools with the lowest response were confirmed with
repeated measurements.

Short and long attenuation lengths are quantified by separating the data into
two distinct categories by fitting an exponential to the range 1-3 m and 3-25
m from the fiber end, respectively, and extracting the exponential constants.
Distributions of fiber short and long attenuation lengths, as measured from
the readout end, are shown in figure 4.

Figure 4: The distribution of short (left) and long (right) attenuation
lengths for the thirty-four 1.8 mm spools with a bin width of 0.2 m. The short
attenuation length mean and standard deviation are $6.3\pm 0.5$ m, while the
long attenuation length statistics are $10.4\pm 0.2$ m.

The spectrometer enables light yield measurements, separated by wavelength in
bins of 1 nm, at 50 cm increments along the first 25 m of fiber. A single
attenuation length is extracted for each wavelength using an exponential fit
of the spectral light intensity as a function of distance from the readout
end. This single exponential function yields a reasonable fit over the
measured 25 m of fiber when applied to single spectral measurements. Figure 5
shows the attenuation lengths for all production fiber and the benchmark spool
as a function of wavelength. The attenuation length drop around 610 nm is due
to an absorption resonance from the fiber’s polystyrene core while the
attenuation lengths around the optimal SiPM sensitivity range show low
variance. Overall, all new 1.8 mm production fibers displayed higher light
yield responses and longer attenuation lengths than the benchmark spool across
the wavelength ranges.

Figure 5: Attenuation length as a function of fiber emission wavelength for
all production spools and the benchmark spool (lowest curve). These
attenuation lengths are extracted using a single exponential fit to the
spectrometer data taken from the first 25 m of each fiber spool.

###  4.2 S-type versus non-S-type fiber comparison

There are two fiber types considered in this work, S-type and non-S-type,
whose differences arise from the presence or absence of a molecular core
(Polystyrene chain) orientation, respectively. Mu2e uses non-S-type fibers due
to its greater transparency (more than 10% longer attenuation length in
comparison to S-type fiber), but a sample 1.8 mm S-type fiber was purchased
and tested for a hadronic calorimeter (HCal) prototype for the Light Dark
Matter eXperiment (LDMX) [17]. A minor inconvenience of the non-S-type fiber
is a lower flexibility, which requires additional precaution during fiber
handling. Figure 6 shows a comparison of the S-type fiber purchased for the
LDMX HCal prototype with the 1.8 mm non-S-type spools purchased for the Mu2e
CRV. As expected, greater light yield and attenuation-length performance is
observed for the non-S-type fiber.

Figure 6: Comparison of 1.8 mm S-type and non-S-type fiber light yield
performance (left) and attenuation length vs wavelength (right) for three
fibers manufactured in the displayed years. The Mu2e experiment uses non-S-
type fiber while the 1.8 mm LDMX fiber was S-type fiber.

###  4.3 Fiber light yield comparison

As mentioned earlier in the paper, light yield measurements were conducted in
the Fermilab test beam to determine the light yield gain from larger fiber
diameter for fibers of the same age. A 24% light yield gain was measured for
the 1.8 mm fiber over the original 1.4 mm fiber. During the UVA fiber testing
runs, a sample di-counter was created by threading 1.4 mm diameter fiber
through one scintillator side and 1.8 mm through the other side and tested
with UVA’s cosmic test stand [14] to determine the actual light yield increase
expected in the CRV sectors. The larger diameter fiber once again recorded
about a 25% light yield increase over the 1.4 mm fiber which is consistent
with geometrical expectations. In addition, several fully constructed CRV
modules with the new 1.8 mm fiber also show a similar light yield increase.
Both the test-beam data and measurements from the cosmic ray test stand
confirm that larger fiber diameter achieves the expected gain in light yield.
However, the larger than expected light yield increase as measured by the
fiber scanner, beyond what was measured in the fully constructed di-counters,
hinted at potential aging effects in the fiber for wavelengths outside of the
spectral sensitivity range of the SiPMs. This is described in the following
section.

###  4.4 Fiber aging studies

After observing that the benchmark 1.8 mm fiber spool, fabricated in 2015, had
a significant decrease ($\sim 30\%$) in light yield in comparison to the new
spools from 2021, and that all of the production 1.8 mm spools had a larger
than expected increase in light yield over the older 1.4 mm fiber ($\sim
60\%-70\%$), it was decided to investigate potential fiber aging effects.
Specifically, five of the thirty-four new 1.8 mm fiber spools, and the
benchmark 1.8 mm fiber spool, were remeasured after approximately one year
with both the spectrometer and photodiode to provide a direct measurement of
aging effects.555Only five spools were available for this measurement since
the other spools had already been used for CRV production.

The results, shown in figure 7, indicate a few percent decrease in light yield
in the newer spools but no appreciable aging effects in the benchmark spool.
This hints that aging might occur within the first years of a fiber’s
lifetime, but slows down significantly in later years. In addition, the
increase in light yield loss with the distance from the readout end points to
a decrease in attenuation length. From figure 5, there are hints of decreasing
attenuation lengths for these wavelengths above 520 nm; however, there is also
large run-to-run variance for these larger wavelengths which can range between
1-5%. Thus, narrowing this measured aging effect to a specific wavelength is
inconclusive in this work.

Figure 7: The average ratio of new to old photodiode measurements taken
approximately one year apart for the five new 1.8 mm spools (black) and the
benchmark spool (blue). There is some apparent aging for the new spools, while
there is no measurable aging in the benchmark spool. This indicates that there
may be aging effects occurring in the first few years of the fiber’s lifetime.
The shaded region represents the standard deviation of the five measured
spools.

Caution should be exercised in interpreting these results in the context of
potential aging that may occur within fully constructed Mu2e CRV di-counters.
To start, this is not a measurement of the decrease in light yield due to
fiber aging within a fully constructed counter as measured by a SiPM.
Measurements of light yield for several fully constructed CRV modules do not
support rapid fiber aging effects. The difference is that both the
spectrometer and photodiode measure a broad spectrum of wavelengths, while the
SiPM measurements occur in a much narrower region of wavelength as illustrated
in Fig. 8. The measured aging of a fully constructed module in the CRV is
still likely dominated by scintillator aging effects as discussed in the
previous section. In addition, the longest CRV module is 7 m while much of the
measured aging occurs beyond that length scale.

Figure 8: A zoomed in version of figure 5 in which the shaded region
illustrates the wavelengths of interest for the Mu2e CRV. Wavelengths below
490 nm have attenuation lengths much shorter than the length of a Mu2e CRV di-
counter. Above 530 nm, the SiPM quantum efficiency and fiber emission spectrum
is about 20% of its peak and drops rapidly.

It is also possible that, in addition to real aging effects, manufacturing
differences between batches of spools produced several years apart could
contribute to the apparent aging between the benchmark spool and the new
spools. In order to obtain a precise quantitative measure of aging in these
fibers, a more careful dedicated study will need to be performed which is
beyond the scope and interest of the Mu2e CRV. Specifically, testing many
spools from the same manufactured batch multiple times at regular time
intervals (on a month to year time scale) using a dedicated setup in a
temperature-controlled environment would narrow the possibilities discussed
above.

##  5 Conclusion

Light yield and attenuation length measurements were conducted for thirty-four
non-S-type and one S-type of 1.8 mm diameter Kuraray Y11 WLS fiber using a
fiber scanner designed for the Mu2e CRV detector. The fibers were purchased to
enhance the light yield in the CRV’s critical sectors. A significant light
yield increase of about 25% was observed for the 1.8 mm over the 1.4 mm
diameter fibers. Preliminary measurements report less than 5% fiber aging per
year and suggests negligible aging effects for fibers more than five years
old.

## Acknowledgments

The bulk of this research was funded by the DOE HEP UVA base grant under DE-
FOA-0001781. The S-type fiber was funded with grant number 20190875 from The
Crafoord Foundation and purchased for a test-beam prototype for the LDMX
experiment.

This document was prepared by members of the Mu2e Collaboration using the
resources of the Fermi National Accelerator Laboratory (Fermilab), a U.S.
Department of Energy, Office of Science, HEP User Facility. Fermilab is
managed by Fermi Research Alliance, LLC (FRA), acting under Contract No. DE-
AC02-07CH11359. Finally, the authors would like to thank Kuraray and their
staff for providing the QA measurements and allowing us to discuss them in
this paper.

## References

  * [1] L. Bartoszek et al., _Mu2e Technical Design Report_ , (2015) [arXiv:1501.05241].
  * [2] F. Abusalma et al., [Mu2e II], arXiv:1802.02599, 2018.
  * [3] F. Abdi et al. [Mu2e Collaboration], arXiv:2210.11380, 2022.
  * [4] R. H. Bernstein and P. S. Cooper, “Charged Lepton Flavor Violation: An Experimenter’s Guide,” Phys. Rept. 532, 27 (2013).
  * [5] E. C. Dukes, P. J. Farris, R. C. Group, T. Lam, D. Shooltz and Y. Oksuzian, JINST 13, no. 12, P12028 (2018).
  * [6] G. DeZoort, E. C. Dukes, R. C. Group, H. Kessenich, Y. Oksuzian, T. Rase and D. Shooltz, arXiv:1511.06225, 2015.
  * [7] A. Artikov et al., Nucl. Instrum. Meth. A 890, 84 (2018).
  * [8] Kuraray America, Inc., http://kuraraypsf.jp/psf/, [Accessed May 2021].
  * [9] G. Blazey et al., Nucl. Instrum. Meth. A volume 927 (2019), arXiv:1906.07237
  * [10] Hamamatsu Photonics K.K., https://www.hamamatsu.com/eu/en/product/type/S13360-2050VE/index.html, [Accessed July 2021].
  * [11] CREE brand XLamp 3, model 3XL7090RY-L00-D5-14, https://www.cree.com/led-components/, [Accessed May 2021].
  * [12] Hamamatsu Photonics K.K., https://www.hamamatsu.com/us/en/product/type/S1227-1010BR/index.html, [Accessed July 2021].
  * [13] Ocean Insight, https://www.oceaninsight.com/products/spectrometers/microspectrometer/sts-series/sts-vis/, [Accessed July 2021].
  * [14] S. Boi, (2021). “Design and Fabrication of a Novel Large-area, High-efficiency Cosmic Ray Veto Detector for the Mu2e Experiment” [Doctoral Thesis]. University of Virginia.
  * [15] C. Group, P. Farris and Y. Oksuzian, PoS NuFACT2018, 040 (2019).
  * [16] https://oceanoptics.com/product/sts-developers-kit/, [Accessed August 2018].
  * [17] T. Åkesson et al. [LDMX Collaboration], JHEP 04 (2020), 003.


# Variational Quantum Multi-Objective Optimization

Linus Ekstrøm  Honda Research Institute Europe GmbH, Carl-Legien-Str. 30,
63073 Offenbach, Germany  Hao Wang  LIACS/aQa, Leiden University, Niels
Bohrweg 1, 2333 CA Leiden, Netherlands  Sebastian Schmitt  Honda Research
Institute Europe GmbH, Carl-Legien-Str. 30, 63073 Offenbach, Germany

###### Abstract

Solving combinatorial optimization problems using variational quantum
algorithms to be executed on near-term quantum devices has gained a lot of
attraction in recent years. Currently, most works have focused on single-
objective problems. In contrast, many real-world problems need to consider
multiple conflicting objectives simultaneously, which is not well studied
using variation quantum algorithms. In multi-objective optimization, one seeks
the optimal trade-offs among conflicting objectives - the well-known Pareto
set/front. We present a variational quantum multiple-objective optimization
(QMOO) algorithm, which allows us to solve multi-objective optimization
problems using NISQ computers. At the core of the algorithm is a variational
quantum circuit (VQC) tuned to produce a quantum state which is a
superposition of Pareto-optimal solutions, solving the original multi-
objective optimization problem. The VQC achieves this by incorporating all
cost Hamiltonians representing the classical objective functions. We retrieve
a set of solutions from the quantum state prepared by the VQC, and utilize the
widely-applied hypervolume indicator to determine the quality of it as
approximation to the Pareto-front. The variational parameters of the VQC are
tuning by maximizing the hypervolume indicator. As many realistic problems are
integer optimization problems we formulate the whole scheme for qudit quantum
systems. We show the effectiveness of the proposed algorithm on several
benchmark problems with up to five objectives.

##  I Introduction

Over the past few years, quantum computing experiments have indicated that the
threshold of practically useful NISQ quantum computation might be nearing [1,
2]. Also, practical application problems have increasingly moved into the
focus of quantum algorithm research, in particular for the currently available
noisy-intermediate scale quantum (NISQ) devices [3, 4, 5, 6]. A promising
class of algorithms is the variational algorithms such as the quantum
approximate optimization algorithm (QAOA) [7, 8]. However, the exact nature of
the potential quantum advantage for real-world combinatorial optimization
problems remains an open problem [9, 10, 11].

Current versions of variational quantum optimization approaches only deal with
single-objective optimization problems. The final target for the optimized
quantum circuit is to produce a quantum state where the weight, i.e. the
probability to sample the optimal solution of the original problem is as large
as possible. However, many realistic application problems need to consider
more than one objective, such as cost and quality, simultaneously [12, 13].
Since these objectives are mostly in a trade-off relation where increasing the
performance in one objective leads to a degradation in the other objectives,
finding one single optimal solution is not a proper target setting anymore.
Instead, that goal is to find a whole set of Pareto-optimal solutions, which
represent the best possible trade-off solutions. Solving such multi-objective
optimization problems directly with a variational quantum algorithm is
currently not possible.

A recent proposal for solving multi-objective optimization problems with NISQ
devices[14] utilizes a scalarizing approach. The objectives are combined into
several single objectives by weighted sums. Each scalarized objective is then
solved using a regular QAOA to obtain one Pareto-optimal solution. Combining
the solutions of several different scalarizing functions allows to reconstruct
convex parts of the Pareto-front. In contrast, we propose a quantum multi-
objective algorithm where the multiple objectives enter the quantum circuit
and the whole set of Pareto-optimal solutions is extracted from the resulting
one quantum state once the optimal parameters for the quantum circuit are
determined once.

Another recent approach utilized the classical NSGAII multi-objective
algorithm to obtain the solutions to a constraint optimization problem [15].
The quantum state prepared by a variational quantum circuit (VQC) is used to
calculate two objectives which are fed into the classical algorithm. The first
objective is the classical cost function to be minimized while the second
objective is given by the fraction of constraint-fulfilling solutions in the
quantum state. The outcome of this procedure is still only one optimal state
which represents the minimal cost function value in the feasible subspace. In
contrast, our proposed algorithm targets to represent multiple solutions to an
multi-objective optimization problem in a single quantum state.

While the majority of investigations focus on qubit quantum computation,
qubits can lead to significant encoding overhead for integer optimization
problems [5, 16]. In the NISQ-era of quantum computation, such overheads may
be prohibitive to obtaining useful results. It has been suggested some of
these limitations can be handled effectively by utilizing multi-level systems
called qudits. There are some indications that qudit systems may provide some
algorithmic benefits, for example, for circuit simplification [17, 18, 19],
error correction [20, 21] or algorithm speed-up [22]. Qudit-based approaches
have been used for several applications like max-cut, graph-coloring, or
machine learning problems [23, 24, 4, 25, 26, 27]. Qudits also represent a
very promising route to larger Hilbert spaces for actual quantum computing
hardware [28, 29, 30, 31] Therefore, we formulate and evaluate the proposed
algorithm with qudit variables.

The structure of the paper is as follows. After a brief introduction to multi-
objective optimization we repeat the basics of qudit-based formulation of
quantum algorithms. Then we describe the proposed quantum multi-objective
(QMOO) algorithm and introduce several benchmark functions for evaluating its
performance. We present empirical results from simulations of the algorithm.
Finally, we discuss the results, present our conclusions from this work and
give an outlook for future directions of investigations.

##  II Multi-objective optimization

A multi-objective optimization problem (MOP) involves minimizing multiple cost
functions $C_{i}(\mathbf{x})$ simultaneously, i.e.,
$\vec{C}=(C_{1},\ldots,C_{K})$,
$C_{i}:\mathbf{x}\in\mathcal{X}\rightarrow\mathbf{y}\in\mathbb{R}$,
$i\in[1..K]$. For every $\mathbf{y}^{(1)}$ and
$\mathbf{y}^{(2)}\in\mathbb{R}^{K}$, we say $\mathbf{y}^{(1)}$ weakly
dominates $\mathbf{y}^{(2)}$ (denoted by
$\mathbf{y}^{(1)}\preceq\mathbf{y}^{(2)}$) iff. $y^{(1)}_{i}\leq y^{(2)}_{i}$,
$i\in[1..m]$. The Pareto order $\prec$ on $\mathbb{R}^{m}$ is defined:
$\mathbf{y}^{(1)}\prec\mathbf{y}^{(2)}$ iff.
$\mathbf{y}^{(1)}\preceq\mathbf{y}^{(2)}$ and
$\mathbf{y}^{(1)}\neq\mathbf{y}^{(2)}$. A point $\mathbf{x}\in\mathcal{X}$ is
(Pareto-)efficient iff.
$\nexists\mathbf{x}^{\prime}\in\mathcal{X}(\vec{C}(\mathbf{x}^{\prime})\prec\vec{C}(\mathbf{x}))$.
The set of all efficient points in $\mathcal{X}$ is called the _efficient set_
or _Pareto set_. The image of the efficient set under $\vec{C}$ is called the
_Pareto front_.

The arguably most important performance metric in multi-objective optimization
is the well-known hypervolume (HV) indicator [32, 33, 34, 35, 36]. HV of a set
$Y\subseteq\mathbb{R}^{K}$ is defined as the $K$-dimensional Lebesgue measure
$\lambda$ of the subset dominated by $Y$ and bounded above by a reference
point $\vec{r}\in\mathbb{R}^{K}$, i.e.,
$\operatorname{HV}(Y,\vec{r})=\lambda(\\{\vec{y}\in\mathbb{R}^{K}\colon\vec{y}\prec\vec{r}\wedge\exists\vec{p}\in
Y,\vec{p}\prec\vec{y}\\})$. The HV indicator is Pareto-compliant [32, 37, 38,
39]: for two subsets $A$ and $B\subseteq\mathbb{R}^{K}$, $A\prec
B\implies\operatorname{HV}(A,\vec{r})>\operatorname{HV}(B,\vec{r})$, given the
same reference point $\vec{r}$. Due to this property, the HV indicator is
often employed as the performance indicator of multi-objective optimization
algorithms. For continuous problem domains, i.e.,
$\mathcal{X}\subseteq\mathbb{R}^{N}$, maximizing the HV indicator leads to an
asymptotic distribution of points on the Pareto front [40], which depends on
the curvature thereof. To our knowledge, no theoretical results have been
reported for discrete problems, e.g.,
$\mathcal{X}\subseteq\\{0,\dots,d-1\\}^{N}$. In this work, we deal with
unconstrained combinatorial integer optimization problems, namely
$\mathcal{X}=[0..d-1]^{N}$. A schematic illustration of the concepts of multi-
objective optimization is shown in Fig. 1.

Figure 1: Illustration of a multi-objective optimization problem for two-
objectives ($K=2$). The left graph illustrates the search space where each
point represents one potential solution, while the right graph shows the
corresponding objective space. Solutions belonging to the Pareto-optimal set,
and the corresponding points on the Pareto-front are marked in red. The gray
shared area in the objective space indicates the hypervolume of the true
Pareto-front with respect to a reference point (yellow star).

##  III Qudit Representation

The search variable $\mathbf{x}\in\mathcal{X}$ is a vector of $N$ finite-
ranged integers: $\mathbf{x}=(x_{1},\dots,x_{N})^{\top}$, where each variable
takes its value $\\{0,1,\dots,d-1\\}$.

In the quantum formulation these integer variables are represented by a set of
$N$ qudits, i.e. $d$-level quantum systems. For $d=2$ the audits reduce to
commonly investigated qubits. The computational basis is chosen such that each
possible search vector corresponds to one $N$ qudit quantum state,

$\displaystyle\mathbf{x}=(x_{1},\dots,x_{N})^{T}\Leftrightarrow\ket{\mathbf{x}}=\ket{x_{1},\dots,x_{N}}$
(1)

The total dimension of the Hilbert space is $\text{dim}\,\mathcal{H}=d^{N}$
any state vector $\ket{\psi}\in\mathcal{H}$ can be written as

$\displaystyle\ket{\psi}=\sum_{x_{1}=0}^{d-1}\ldots\sum_{x_{N}=0}^{d-1}\alpha_{x_{1}\ldots
x_{N}}\ket{\mathbf{x}}\,,$ (2)

where $\alpha_{x_{1}\ldots x_{n}}$ is the complex amplitude and the states
$\ket{\mathbf{x}}=\ket{x_{1},\ldots,x_{N}}$ form an orthonormal basis, i.e.,
$\braket{\mathbf{x}}{\mathbf{x}^{\prime}}=\delta_{\mathbf{x},\mathbf{x}^{\prime}}=\delta_{x_{1},x_{1}^{\prime}}\dots\delta_{x_{N},x_{N}^{\prime}}$.

In general, one needs $d^{2}-1$ Hermitian basis operators in order to generate
all unitary operations of the $SU(d)$ group, which represent all allowed
actions on a single qudit in the quantum circuit[17]. However, it was
shown[41, 42] that it is sufficient to consider only the three operators
$L_{x}$, $L_{z}$ and $L_{z^{2}}=L_{z}^{2}$ as generators. Here $L_{x}$ and
$L_{z}$ are the $x$\- and $z$-angular momentum operators which act on the
single qudit state in the usual manner,

$\displaystyle L_{x}\ket{x}$
$\displaystyle=\tfrac{1}{2}\big{(}\gamma_{d,x+1}\ket{x+1}+\gamma_{d,x-1}\ket{x-1}\big{)}$
(3) $\displaystyle L_{z}\ket{x}$ $\displaystyle=\tfrac{2x-d+1}{2}\ket{x}$ (4)

where $x\in[0,d-1]$ and $\gamma_{d,x}=\sqrt{(d-x-1)(x+1)}$. The operator
$L_{z^{2}}$ is called squeezing or one-axis twisting operator in the context
of cold-atom systems. These three generators are sufficient to implement
arbitrary unitary operation by (possibly many) repeated finite rotations. The
reason for this is that the iterated commutators of these three operators
generate all $d^{2}-1$ Hermitian basis operators, which are then necessary to
generate all unitary operations of the $SU(d)$ group.

##  IV Quantum Multi-Objective Optimization

Each objective $C_{k}(\textbf{x})$ is represented by a cost Hamiltonian
$H_{k}$, which is diagonal in the measurement basis $\ket{\textbf{x}}$. The
eigenvalue of the basis state is given by the cost function value of that
state,

$\displaystyle H_{k}\ket{\mathbf{x}}$
$\displaystyle=C_{k}(\mathbf{x})\ket{\mathbf{x}}\quad(k=1,\dots,K).$ (5)

Figure 2: Schematic overview of the QMOO algorithm.

We employ a variational quantum circuit to solve many-objective optimization
problems. Analogous to the traditional QAOA, the circuit is built up from $L$
layers, where the functional form of the unitary operation in each layer is
the same across all layers, but each layer has independently adjustable
variational parameters $\bm{\beta}^{L}$ and $\bm{\gamma}^{L}$, i.e.,

$\displaystyle\ket{\psi(\bm{\beta},\bm{\gamma})}$
$\displaystyle=U(\bm{\beta}^{L},\bm{\gamma}^{L})\cdots
U(\bm{\beta}^{1},\bm{\gamma}^{1})\ket{\psi_{0}}\;.$ (6)

$\ket{\psi_{0}}$ denotes a predefined fixed initial state. We choose the
uniform superposition of all possible states,

$\ket{\psi_{0}}=\frac{1}{\sqrt{d^{N}}}\sum_{\mathbf{x}\in\\{0,\cdots,d-1\\}^{N}}\ket{\mathbf{x}},$
(7)

which is an eigenvector to the generalized Pauli $X$-operator [17].
Alternative initial states are also possible, such as one single basis state
or an eigenstate to the $x$ angular momentum operator $L_{x}$. The latter is
typically done for qubit circuits. However, the specific choice of the initial
state should not change the algorithm’s performance qualitatively: the
different initial states are related by unitary transformations, which can be
generated implicitly by the variation quantum circuit. Therefore, the initial
states merely alter the variational optimization problem slightly without
changing its characteristics. We experimentally validated this argument with
various different initial states (not shown in this paper).

Each layer $l$ of the circuit comprises $K$ blocks, where $K$ is the number of
classical objectives. Each block consists of one mixing unitary $U_{M}$ and
one phase unitary $P_{k}$,

$\displaystyle U(\bm{\beta}^{l},\bm{\gamma}^{l})$
$\displaystyle=U_{M}(\bm{\beta}^{l}_{K})P_{K}(\bm{\gamma}^{l}_{K})\cdots
U_{M}(\bm{\beta}^{l}_{1})P_{1}(\bm{\gamma}^{l}_{1})$ (8)

This is in contrast to the regular QAOA where only one phase and one mixing
operators is included, i.e. $K=1$.

While, in principle, the mixing unitaries of each block in Eq. (8) can have a
different structure, we choose them to be the same but allow for different
variational parameters of each mixer. We focus on mixing unitaries which are
generated by the $x$-angular momentum $L_{x}$ and the local squeezing
operator, $L_{z^{2}}=L_{z}^{2}$, which is sufficient to generate all unitaries
as described in Sec. III above. This leads to the mixing unitary of block $k$

$\displaystyle U_{M}(\bm{\beta}^{l}_{k})$
$\displaystyle=e^{-i\beta^{l}_{1,k}\sum_{j=1}^{N}L_{x;j}+\beta^{l}_{2,k}\sum_{j=1}^{N}L_{z^{2};j}}\,,$
(9)

where the sum in the exponent extends of all qudits $j=1,\dots,N$. Note that
for qubits ($d=2$) the local squeezing operator reduces to the identity,
$L_{z^{2}}=\tfrac{1}{4}\mathds{1}$, and consequently we drop that term in that
case.

The phase operators are generated by the $K$ cost Hamiltonians representing
the classical objectives, i.e.,

$\displaystyle P_{k}(\bm{\gamma}^{l}_{k})$
$\displaystyle=e^{-i\gamma^{l}_{k}H_{k}}\,.$ (10)

where $k=1,\dots,K$. The idea behind this choice is that each phase operator
produces a relative phase with the energies with respect to one objective. In
that way, the information of each objective can be encoded in relative phases
of states contributing to the quantum state prepared by the circuit. Given the
individual variation parameters for each objective, the classical optimization
procedure can, therefore, in principle, discriminate between states based on
the set of objective values.

For a given set of variation parameters $(\bm{\beta},\bm{\gamma})$ the quantum
circuit generates a final state as depicted in Eq. (6). In general, this state
is a superposition of all measurement basis states, i.e., a superposition of
all classical solutions,

$\displaystyle\ket{\psi(\bm{\beta},\bm{\gamma})}=\sum_{\text{all
}\textbf{x}}\xi_{\textbf{x}}\ket{\textbf{x}},$ (11)

where $\xi_{\textbf{x}}\in\mathbb{C}$ and
$\sum_{\textbf{x}}|\xi_{\textbf{x}}|^{2}=1$. In standard QAOA, this state is
used to calculate an energy expectation value of the cost Hamiltonian, which
is then minimized by tuning the variational parameters with a classical
optimization routine. The ideal outcome of QAOA is a final quantum state whose
spectral weights concentrate on the global optimal solution. Therefore, the
superposition of Eq. (11) reduces to one classical measurement basis state
$\ket{\mathbf{x}_{\text{opt}}}$ with $|\xi_{\text{opt}}|^{2}\approx 1$.

In contrast, the final quantum states of the proposed multi-objective approach
should represent a whole set of classical solutions as a superposition at
once. A complete set of $N_{S}$ classical solutions is extracted from the one
quantum state by selecting the $N_{S}$ most probable basis states contributing
to the quantum state. Practically, this can be done by implementing a full
state tomography or state discrimination protocol, see, e.g., [43], and
selecting the $N_{S}$ states with largest spectral weights
$|\xi_{\mathbf{x}}|^{2}$. Alternatively, the set of solutions can be extracted
by performing multiple measurement shots and selecting the $N_{S}$ most
frequently measured basis states 111The effect of variations in the selected
solutions induced by a finite number of measurement shots needs to be
investigated..

So the ideal outcome of the proposed QMOO is a final quantum state which
encodes the whole set of Pareto-optimal classical solutions as superposition
with $\sum_{\mathbf{x}\in\mathbb{PS}}|\xi_{\mathbf{x}}|^{2}\approx 1$, and
$\mathbb{PS}$ denotes the efficient set of Pareto-optimal solutions.

Practically, we extract $N_{S}$ different classical solutions in each
iteration of the algorithm,

$\displaystyle\ket{\psi(\bm{\beta},\bm{\gamma})}\Rightarrow\mathbb{S}_{\bm{\beta},\bm{\gamma}}=\\{\mathbf{x}_{i}\\}_{i=1}^{N_{S}}\,.$
(12)

For each of those extracted classical solutions, all classical cost objectives
are calculated,

$\displaystyle\mathbb{C}_{\bm{\beta},\bm{\gamma}}$
$\displaystyle=\\{\\{C_{k}(\mathbf{x}_{i})\\}_{k=1}^{K}\\}_{i=1}^{N_{S}}$ (13)
$\displaystyle=\\{\vec{C}(\mathbf{x}_{i})\\}_{i=1}^{N_{S}}\,.$ (14)

From this set, we extract the non-dominated solutions, which approximate the
Pareto front. The quality of the current Pareto front is then evaluated by
calculating the hypervolume indicator (with a predefined reference point
$\vec{r}$),

$\displaystyle\operatorname{HV}\left(\\{\vec{C}(\mathbf{x}_{i})\\}_{i=1}^{N_{S}};\vec{r}\right).$
(15)

We employ a classical optimization algorithm to optimize the parameters
$(\bm{\beta},\bm{\gamma)}$ of the quantum circuit such that the hypervolume is
maximized. The target quantum state is a superposition of classical solutions
with large weights for non-dominated solutions, which form a good
approximation to the Pareto front, ideally, that classical solution which is
part of the Pareto-optimal set.

The total number of variational parameters of the proposed algorithm is
$N_{\text{params}}=3LK$, where $L$ is the number of layers and $K$ is the
number of objectives. For the case of qubits ($d=2$) without a local squeezing
operator, the number of parameters reduces to $N_{\text{params}}=2LK$.

##  V Benchmark Functions and experimental setup

###  V.1 Benchmark Functions

The performance of the proposed multi-objective quantum optimization algorithm
is evaluated on several benchmark functions. We combine linear and quadratic
cost functions into two-, three-, and five-objective problems where all
objectives are subject to minimization. To define a useful multi-objective
benchmark, the objectives should be in some sort of trade-off relation, where
improving on one objective typically leads to performance degradation in other
objectives.

The first set of benchmark problems is defined by a linear two-objective
function of the form

$\displaystyle\vec{C}_{\mathrm{linear}}(\mathbf{x})$
$\displaystyle=\begin{pmatrix}C_{1}(\mathbf{x})\\\
C_{2}(\mathbf{x})\end{pmatrix}=\begin{pmatrix}\mathbf{c}_{1}\cdot\mathbf{x}\\\
\mathbf{c}_{2}\cdot\mathbf{x}\end{pmatrix}\,.$ (16)

Here, $\mathbf{x}\in[0..d-1]^{N}$ and the cost coefficients are sampled from
$[-1,1]^{N}$ uniformly at random, i.e.,
$\mathbf{c}_{1}\sim\mathcal{U}([-1,1]^{N})$. The trade-off between the
objectives is generated by choosing the second cost coefficients to have
substantial anti-correlations with the first coefficients, i.e.,
$\mathbf{c}_{2}\sim-\frac{1}{2}\mathbf{c}_{1}+\frac{1}{2}\mathcal{U}([-1,1]^{N})$,
which amount to a correlation coefficient of $-\frac{1}{2}$.

We also composed two-, three-, and five-objective benchmark problems using
randomized quadratic cost functions of the form

$\displaystyle C(\mathbf{x})$
$\displaystyle=\mathbf{x}\cdot(\mathbf{J}\mathbf{x})+\mathbf{h}\cdot\mathbf{x}\,,$
(17)

where $\mathbf{J}=\mathbf{J}^{\top}\in\mathbb{R}^{N\times N}$ are symmetric
randomized all-to-all coupling matrices and $\mathbf{h}\in\mathbb{R}^{N}$ are
randomized local fields. To generate trade-offs between the objectives, we are
guided by physical spin models that realize such a cost Hamiltonians.

For a two-objective problem ($K=2$), we chose one objective to favor
ferromagnetic (FM) alignment of the variables and the other to favor anti-
ferromagnetic (AFM) alignment. Consequently, the coupling matrix elements are
chosen from uniform distributions, i.e., $J_{ij}=J_{ji}\sim\mathcal{U}(a,b)$
where for the FM case all matrix elements are negative, i.e., $a<b<0$, and for
the AFM case all are positive, i.e., $b>a>0$. For the cost function to truly
represent spin models, the local fields are chosen as
$\mathbf{h}=\mathbf{g}-d\mathbf{1}\cdot\mathbf{J}$ where $\mathbf{g}$ is a
random vector with $g_{i}\sim\mathcal{U}(-1,1)$, $d$ the local dimension of
the qudits and $\mathbf{1}$ the $N$-dimensional vector of ones. This problem
type is named quadratic FM-AFM problem.

We define another two-objective benchmark as follows: the first objective is
implemented with AFM couplings described above, and the second quadratic
objective is designed to minimize the distance of the search variables to a
random point $\mathbf{x}_{0}\sim\mathcal{U}([0..d-1]^{N})$. For the second
objective, the cost function matrix is just the identity matrix, i.e.,
$J_{ij}=\delta_{ij}$, and the local fields are given by the random point
$\mathbf{h}=-2\mathbf{x}_{0}$. We call this problem quadratic
$\mathbf{x}_{0}$-AFM.

For the three-objective problem ($K=3$), We take the first and second
objectives as AFM and FM-coupled spin models, respectively. The third
objective is given by this quadratic distance to a given random vector
$\textbf{x}_{0}$.

The five objective benchmark problems are constructed as follows. The first
two objectives are chosen to be the AFM and FM spin chains with nearest
neighbor coupling, i.e., $J_{ij}=J_{ji}\sim\delta_{j,i+1}\mathcal{U}(a,b)$
with $b>a>0$ and $a<b<0$ for the AFM and the FM case, respectively. The third
objective is the distance to a given random vector
$\mathbf{x}_{0}\sim\mathcal{U}([0..d-1]^{N})$, i.e., $J_{ij}=\delta_{ij}$ and
$\mathbf{h}=-2\mathbf{x}_{0}$. The fourth (fifth) objective is also realized
as a nearest neighbor spin-chain, but where the first half of the variables,
$x_{q}$ with $q\leq N/2$, have FM (AFM) couplings while the latter half,
$x_{q}$ with $q>N/2$, have opposite, i.e. AFM (FM), couplings.

Figure 3 shows schematic examples of the linear and quadratic objective
functions for two- and three-objective benchmarks for qubits ($d=2$), qutrits
($d=3$) and ququints ($d=5$). Note that we normalized each objective cost
function such that the values are between $C(\textbf{x})\in[0,1]$ for all
functions, numbers of qudits $N$ and dimensionality of the qudits $d$.

Figure 3: Objective space of example instances of several benchmark problems.
The problems are from left to right: linear $K=2$ objective, quadratic FM-AFM
$K=2$, quadratic $\mathbf{x}_{0}$-AFM $K=2$, and quadratic $K=3$ objective
problems. The top row is for $N=12$ qubits ($d=2$), the middle row is $N=8$
qutrits ($d=3$), and the lowest row is $N=5$ ququints ($d=5$). The grey dots
indicate all possible solutions, while the blue dots indicate the true Pareto-
optimal front.

###  V.2 Experimental Setup

The proposed quantum multi-objective optimization algorithm is evaluated on
several instances of the previously introduced benchmark functions. For each
specific randomly generated instance of a benchmark problem, the algorithm is
run 40 times with different random initializations of the parameter values of
the quantum circuits. For each run, all parameters are initialized randomly in
the interval $[-\pi,\pi]$ and $N_{S}=20$ classical solutions are extracted
from the quantum state.

We tested several classical optimization algorithms for determining the
optimized parameters of the quantum circuit. We used COBYLA, L-BFGS-B, and
Powell from the scipy library [45], covariance matrix adaptation evolutionary
strategies [46] taken from [47], and in-house implementations of particle
swarm optimization (PSO) [48], and differential evolution (DE) [49]. The
overall performance of all tested optimization algorithms was roughly
comparable over all settings and benchmarks. An exception was the L-BFGS-B
algorithm, which yields consistently lower hypervolume values, most likely due
to the numerical evaluation of the gradients. The population-based optimizers
like PSO, CMAES, and DE were slightly outperforming Powell and COBYLA in most
instances, but given their large resource requirements (long algorithms run-
times and a large number of circuit evaluations), we focus on Powell and
COBYLA algorithms.

For comparison, we also run several state-of-the-art classical multi-objective
evolutionary optimization algorithms (MOEAs) in the Platypus [50] package.
Specifically, these are NSGA-II, IBEA, and MOEAD. We set the population size
of the MOEAs to 20, comparable to $N_{S}=20$ extracted solutions from the
quantum states in QMOO. We execute each MOEA with 10 independent runs on each
benchmark problem, and we limit the number of iterations thereof to 200.

##  VI Results

Firstly, we show the evolution of the hypervolume as a function of the
classical optimization iterations for a specific instance of the five-
objective benchmark problem in Fig. 4(a). It is evident that the algorithm can
increase the hypervolume over time and output a state containing a
superposition of Pareto efficient solutions.

Figure 4: (a) Hypervolume as a function of the iterations of the QMOO
algorithm for a particular instance of the five objective optimization problem
using $N=12$ qubits ($d=2)$. The thick line represents the median taken over
40 runs with different random initialization of the QMOO parameters, and the
shaded region denotes the $20$th and $80$th quantiles. The classical Powell
optimizer is used to optimize the VQC, consisting of $L=2$ layers. (b)
Progression of the solution set of one run for a specific instance of the two-
objective quadratic $\mathbf{x}_{0}$-AFM problem for $N=12$ qubits. Gray dots
show all possible solutions, while the colored dots show the solutions
extracted from the QMOO at the indicated iterations. The true Pareto front is
shown with black crosses. The classical Powell optimizer is used for
optimizing the VQC, which consists of $L=1$ layer.

To explicitly observe how the approximation to the Pareto front improves over
the iterations, we show an example of progress for a specific optimization run
for a two-objective quadratic problem in Fig. 4(b). The algorithm is randomly
initialized, which leads to solutions in the top right part of the objective
space. With increasing iterations, better parameters of the VQC are found, and
the solution set represented in the quantum state progresses toward the true
Pareto front.

The algorithm’s behavior is further highlighted by considering the number of
non-dominated and true Pareto-optimal solutions contained in the $N_{S}$
solutions extracted from the quantum state in each iteration. This is shown in
Fig. 5 for three different setups. In general, the number of non-dominated
solutions, as well as the number of true Pareto-optimal solutions, increase
over the optimization iterations. Typically, 50% or larger fractions of the
extracted $N_{S}=20$ solutions are true Pareto-optimal solutions. The left
panels show a very good case where many true Pareto-optimal solutions are
generated from the beginning and their number increases rather quickly. In
contrast, the middle panel shows a case where the algorithm does not fully
converge to the true Pareto-front as only about 10%-20% of the extracted
$N_{S}=20$ solutions are on the true Pareto-front. But this run shows the
interesting behavior that the number of true Pareto-optimal solutions actually
declines toward the end of the iterations, while the non-dominated solutions,
as well as the hypervolume, still increase. This reflects the fact the the
number of true Pareto-optimal solutions does not take the spread of the
current non-dominated solutions into account. So, an improved approximation to
the true Pareto front is achieved by having fewer solutions on the actual
front, but instead spreading the non-dominated solutions better across the
objectives. The rightmost panel shows a case typical for lager number of
objectives where the number of non-dominated solutions is very large directly
from the first iteration. This is a general characteristic of higher-
dimensional objective spaces since the probability that randomly sampled
solutions are non-dominated is increasing rapidly with increasing
dimensionality of the objective space $K$. However, since the QMOO algorithm
operates on the hypervolume, it manages to distribute the non-dominated
solutions more favorably and thus increasingly converges to true Pareto-
optimal solutions.

Figure 5: Median number of non-dominated (orange) as well as true Pareto-
optimal (blue) solutions in the extracted $N_{S}=20$ solutions as function of
the iterations for one specific instance of the quadratic $\textbf{x}_{0}$-AFM
problem for $N=12$ qubits (left) and $N=5$ ququints (middle), and the five
objective problem for $N=5$ ququints (right). The shaded regions are the
$20$th and $80$th quantiles over the $40$ runs. The classical Powell optimizer
is used for optimizing the VQCs, which consist of $L=2$ layers.

The effect of the choice of the classical optimizer was also investigated, and
a representative result is shown in Fig 6(a). For both optimizers, the
hypervolume increases with an increasing number of layers, and the Powell
optimizer is slightly outperforming COBYLA. For this specific problem
instance, the QMOO’s performance is comparable to the collection of classical
algorithms.

Figure 6: Hypervolume as a function of VQC layers. (a) Comparison between two
optimizers for the VQC for a specific instance of the $K=3$ objective problem.
(b) Comparison of situations with and without local squeezing in the mixing
operator for an instance of the $K=5$ objective problem. Both results are for
$N=8$ qutrits ($d=3$), and the result of the purely classical optimizes are
shown on the right in both panels.

As discussed earlier, the quantum circuit must include the local squeezing
operator $L_{z^{2}}$ to generate arbitrary unitaries for qudits. To explicitly
demonstrate its effect, we show a comparison of the results with and without
$L_{z^{2}}$ in Fig. 6(b). As observed, the squeezing operator increases the
performance and reduces the results’ variance substantially.

Figure 7: Normalized hypervolume for the QMOO on all benchmark problems for
different setups compared to the results from purely classical multi-objective
solvers for (a) $N=12$ qubits, (b) $N=8$ qutrits, and (c) $N=5$ ququints. The
VQC consisted of $N_{L}=2$ layers and was optimized with the Powell optimizer.
The results are averaged over 4 different problem instances for each problem
type and 40 randomized runs for each specific problem instance. Normalized
hypervolume refers to the fact that for each specific problem instance, we
normalize the hypervolume with the maximally achievable hypervolume calculated
from the exact Pareto front.

We summarize the overall performance of the proposed QMOO algorithm in Fig. 7
where we show the performance for all five types of benchmarks problems for a
setup with $N=12$ qubits ($d=2$), $N=8$ qutrits ($d=3$), and $N=5$ ququints
($d=5$). We compare the achieved hypervolume values to the ones obtained with
purely classical algorithms. While the mean and median performance of the
proposed QMOO approach is typically lower than the classical approaches, for
most problems, the best-performing QMOO instances can achieve the same
performance as the classical algorithms. In particular, for the three- and
five-objective problems for qutrits and ququints the average performance of
the classical and QMOO algorithms is comparable.

##  VII Discussion and Conclusion

The proposed QMOO algorithm utilizes quantum superpositions as a resource to
represent a whole set of classical solutions in one quantum state. We could
show that this allows us to solve multi-objective optimization problems by
optimizing a parametrized quantum circuit such that the hypervolume indicator
for the extracted classical solutions is maximized. Unlike the standard single
objective QAOA, the algorithm does not target the quantum state generated by
the VQC to reduce to a single measurement basis state (i.e., one classical
solution). Rather, the goal of the QMOO is to obtain a quantum state
consisting of a superposition of non-dominated solutions, ideally belonging to
the true Pareto front of the multi-objective optimization problem.

We showed the effectiveness of the QMOO approach by evaluating its performance
over a variety of benchmark problems, circuit structures, and classical
optimizer algorithms. Over the course of the variational optimization
iterations, the hypervolume indicator, the number of non-dominated solutions
extracted from the quantum state, as well as true Pareto-optimal solutions
tend to increase, indicating a quality improvement of the approximation to the
Pareto-front. For the final optimized quantum circuit, typically around 50% or
more of the solutions extracted from the quantum state were true Pareto-
optimal solutions.

The proposed QMOO approach was able to find solutions of similar quality to
the purely classical methods for many cases. But overall, the performance is
not yet fully competitive with state-of-the-art classical multi-objective
algorithms. However, we consider this work a first step toward a truly quantum
algorithm for solving multi-objective optimization problems, which utilizes
quantum resources for computing Pareto-optimal sets. We are aware that this is
merely a first step, as many issues need to be investigated.

We reconfirmed the importance of the local squeezing operator for qudit
circuits and showed that the performance strongly degrades when this operator
is not included. However, it would be interesting to study the effect of
including more generators of the $SU(d)$ into the mixing operators. In a
recent work [25], it was found that the performance increase from more
generators could be matched by fewer generators but more layers in the VQC. It
needs to be evaluated how the number of parameters/layers in comparison to the
number of included generators influences the performance in the proposed QMOO
approach.

In a similar direction, allowing each qubit to have an individual mixing
parameter to increase the flexibility of the quantum circuit was shown to be
beneficial in some applications [51, 52, 53]. However, this comes at the cost
of a greatly increased number of variational parameters and leads to a more
complicated classical optimization problem.

Notably, we did not make any efforts to improve the performance of the
classical optimization algorithms used for the VQC. We took default
configuration and parameters for all algorithms despite it being known that
tuning those can have a rather large impact [54, 55]. Additionally, techniques
commonly used in the literature for multi-objective optimization can also
improve the algorithm’s performance. For example, employing an archive [56,
57] of evaluated solutions and using the previous observations to inform the
future circuits is an obvious and promising route for further performance
improvements.

Random initialization and direct optimization of all parameters of VQCs is not
very efficient. There exist several initialization and parameter setting
schemes [58, 59, 60, 61, 62, 8] which can boost performance quite
substantially. Even though further investigations are necessary, one can
expect that similar schemes will be effective for the proposed QMOO as well.

The presented approach extracts the set of classical solutions based on the
exact spectral weights in the quantum superpositions. On realistic hardware,
this would require performing a full state tomography, which is rather
resource-expensive. In addition, and probably more importantly, the extraction
of the solution set on actual hardware necessarily leads to noise in the
results. At the current point, it is not clear if the proposed algorithm is
robust to noise and still achieves the desired goal in the presence of shot
noise or sampling errors.

In the current formulation, the Pareto-dominance relations between the
extracted classical solutions are imposed indirectly by the evaluation of the
hypervolume indicator. This leaves it to the classical optimizer to generate
the variational parameters, which produce a large fraction of non-dominated
solutions. It would be much more desirable, and most likely also much more
efficient, if the quantum circuit would incorporate unitary operators that are
sensitive to or act as proxies for Pareto-dominance relations between the
contributions of the quantum state directly. However, designing such operators
is not straightforward and is left for future work.

An immediate first step in that direction could be to gain more insight into
the desired final optimal quantum state for a quantum multi-objective circuit.
In that regard, it would be interesting to investigate the entanglement
structure of quantum states that contain non-dominated solutions compared to
more random superpositions. Also, an investigation of entanglement and its
relation to global squeezing operator similar to a recent work [63] would be
insightful. Similarly, being able to construct the Hamiltonian to which the
desired final quantum state would be the ground state could contribute to
elucidating the nature of quantum circuits for quantum multi-objective
optimization. All these issues pose interesting research questions regarding
the nature of the optimal QMOO quantum state.

We only studied rather small problem instances where the total number of
solutions, i.e., the Hilbert space dimensions, was below $10^{4}$.
Additionally, we did not investigate the relation between the characteristics
of the benchmark problems, such as shape and topology of the Pareto-font
(e.g., convexity, continuity) or features of the objective functions (e.g.,
deceptiveness), and the algorithm’s performance. But these aspects are very
well-suited for elucidating the differences between the quantum state encoding
of the efficient solution set and classical approaches and need to be
considered in follow-up works.

In summary, the proposed QMOO algorithm is a very promising novel approach to
solving multi-objective optimization problems natively on quantum hardware. We
are confident that this approach will widen the application domain for
quantum-optimization approaches and hopefully allow for successful solutions
to realistic problems of practical relevance in the future.

##  VIII Acknowledgements

We acknowledge fruitful discussions with T. Bäck, A. Bottarelli, V. Dunjko, P.
Hauke and B. Sendhoff.

L.E. and S.S. acknowledge funding by the European Union under Horizon Europe
Programme. Views and opinions expressed are however those of the author(s)
only and do not necessarily reflect those of the European Union or European
Climate, Infrastructure and Environment Executive Agency (CINEA). Neither the
European Union nor the granting authority can be held responsible for them.
Grant Agreement 101080086 — NeQST.

## References

  * Acharya  _et al._ [2022] R. Acharya, I. Aleiner, R. Allen, T. I. Andersen, M. Ansmann, F. Arute, K. Arya, A. Asfaw, J. Atalaya, R. Babbush, D. Bacon, J. C. Bardin, J. Basso, A. Bengtsson, S. Boixo, G. Bortoli, A. Bourassa, J. Bovaird, L. Brill, M. Broughton, B. B. Buckley, D. A. Buell, T. Burger, B. Burkett, N. Bushnell, Y. Chen, Z. Chen, B. Chiaro, J. Cogan, R. Collins, P. Conner, W. Courtney, A. L. Crook, B. Curtin, D. M. Debroy, A. D. T. Barba, S. Demura, A. Dunsworth, D. Eppens, C. Erickson, L. Faoro, E. Farhi, R. Fatemi, L. F. Burgos, E. Forati, A. G. Fowler, B. Foxen, W. Giang, C. Gidney, D. Gilboa, M. Giustina, A. G. Dau, J. A. Gross, S. Habegger, M. C. Hamilton, M. P. Harrigan, S. D. Harrington, O. Higgott, J. Hilton, M. Hoffmann, S. Hong, T. Huang, A. Huff, W. J. Huggins, L. B. Ioffe, S. V. Isakov, J. Iveland, E. Jeffrey, Z. Jiang, C. Jones, P. Juhas, D. Kafri, K. Kechedzhi, J. Kelly, T. Khattar, M. Khezri, M. Kieferová, S. Kim, A. Kitaev, P. V. Klimov, A. R. Klots, A. N. Korotkov, F. Kostritsa, J. M. Kreikebaum, D. Landhuis, P. Laptev, K.-M. Lau, L. Laws, J. Lee, K. Lee, B. J. Lester, A. Lill, W. Liu, A. Locharla, E. Lucero, F. D. Malone, J. Marshall, O. Martin, J. R. McClean, T. Mccourt, M. McEwen, A. Megrant, B. M. Costa, X. Mi, K. C. Miao, M. Mohseni, S. Montazeri, A. Morvan, E. Mount, W. Mruczkiewicz, O. Naaman, M. Neeley, C. Neill, A. Nersisyan, H. Neven, M. Newman, J. H. Ng, A. Nguyen, M. Nguyen, M. Y. Niu, T. E. O’Brien, A. Opremcak, J. Platt, A. Petukhov, R. Potter, L. P. Pryadko, C. Quintana, P. Roushan, N. C. Rubin, N. Saei, D. Sank, K. Sankaragomathi, K. J. Satzinger, H. F. Schurkus, C. Schuster, M. J. Shearn, A. Shorter, V. Shvarts, J. Skruzny, V. Smelyanskiy, W. C. Smith, G. Sterling, D. Strain, M. Szalay, A. Torres, G. Vidal, B. Villalonga, C. V. Heidweiller, T. White, C. Xing, Z. J. Yao, P. Yeh, J. Yoo, G. Young, A. Zalcman, Y. Zhang, and N. Zhu, arXiv:2207.06431 (2022).
  * Kim  _et al._ [2023] Y. Kim, A. Eddins, S. Anand, K. Wei, E. Berg, S. Rosenblatt, H. Nayfeh, Y. Wu, M. Zaletel, K. Temme, and A. Kandala, Nature 618, 500 (2023).
  * Bharti  _et al._ [2022] K. Bharti, A. Cervera-Lierta, T. H. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. S. Kottmann, T. Menke, W.-K. Mok, S. Sim, L.-C. Kwek, and A. Aspuru-Guzik, Rev. Mod. Phys. 94, 015004 (2022).
  * Deller  _et al._ [2023] Y. Deller, S. Schmitt, M. Lewenstein, S. Lenk, M. Federer, F. Jendrzejewski, P. Hauke, and V. Kasper, Physical Review A 107, 062410 (2023).
  * Yarkoni  _et al._ [2022] S. Yarkoni, E. Raponi, T. Bäck, and S. Schmitt, Reports on Progress in Physics 85, 104001 (2022).
  * Abbas  _et al._ [2023] A. Abbas, A. Ambainis, B. Augustino, A. Bärtschi, H. Buhrman, C. Coffrin, G. Cortiana, V. Dunjko, D. J. Egger, B. G. Elmegreen, N. Franco, F. Fratini, B. Fuller, J. Gacon, C. Gonciulea, S. Gribling, S. Gupta, S. Hadfield, R. Heese, G. Kircher, T. Kleinert, T. Koch, G. Korpas, S. Lenk, J. Marecek, V. Markov, G. Mazzola, S. Mensa, N. Mohseni, G. Nannicini, C. O’Meara, E. P. Tapia, S. Pokutta, M. Proissl, P. Rebentrost, E. Sahin, B. C. B. Symons, S. Tornow, V. Valls, S. Woerner, M. L. Wolf-Bauwens, J. Yard, S. Yarkoni, D. Zechiel, S. Zhuk, and C. Zoufal, arXiv:2312.02279 (2023).
  * Farhi  _et al._ [2014] E. Farhi, J. Goldstone, and S. Gutmann, arXiv:1411.4028 (2014).
  * Blekos  _et al._ [2023a] K. Blekos, D. Brand, A. Ceschini, C.-H. Chou, R.-H. Li, K. Pandya, and A. Summer, arXiv:2306.09198 (2023a).
  * Blekos  _et al._ [2023b] K. Blekos, D. Brand, A. Ceschini, C.-H. Chou, R.-H. Li, K. Pandya, and A. Summer, arXiv:2306.09198 (2023b).
  * Farhi  _et al._ [2015] E. Farhi, J. Goldstone, and S. Gutmann, arXiv:1412.6062 (2015).
  * Farhi and Harrow [2019] E. Farhi and A. W. Harrow, arXiv:1602.07674 (2019).
  * Sharma and Kumar [2022] S. Sharma and V. A. Kumar, Arch Computat Methods Eng 29, 5605 (2022).
  * Pereira  _et al._ [2022] J. L. J. Pereira, G. A. Oliver, M. B. Francisco, S. S. C. Jr, and G. F. Gomes, Arch Computat Methods Eng 29, 2285 (2022).
  * Chiew  _et al._ [2023] S.-H. Chiew, K. Poirier, R. Mishra, U. Bornheimer, E. Munro, S. H. Foon, C. W. Chen, W. S. Lim, and C. W. Nga, arxiv:2308.08245 (2023).
  * Díez-Valle  _et al._ [2023] P. Díez-Valle, J. Luis-Hita, S. Hernández-Santana, F. Martínez-García, Álvaro Díaz-Fernández, E. Andrés, J. J. García-Ripoll, E. Sánchez-Martínez, and D. Porras, arXiv:2302.04196 (2023).
  * Chancellor [2019] N. Chancellor, Quantum Science and Technology 4, 045004 (2019).
  * Wang  _et al._ [2020] Y. Wang, Z. Hu, B. C. Sanders, and S. Kais, Frontiers in Physics 8, 589504 (2020).
  * Lanyon  _et al._ [2009] B. P. Lanyon, M. Barbieri, M. P. Almeida, T. Jennewein, T. C. Ralph, K. J. Resch, G. J. Pryde, J. L. O’Brien, A. Gilchrist, and A. G. White, Nature Physics 5, 134–140 (2009).
  * Gao  _et al._ [2023] X. Gao, P. Appel, N. Friis, M. Ringbauer, and M. Huber, Quantum 7, 1141 (2023).
  * Duclos-Cianci and Poulin [2013] G. Duclos-Cianci and D. Poulin, Phys. Rev. A 87, 062338 (2013).
  * Watson  _et al._ [2015] F. H. E. Watson, H. Anwar, and D. E. Browne, Phys. Rev. A 92, 032309 (2015).
  * Gedik  _et al._ [2015] Z. Gedik, I. Silva, B. Çakmak, G. Karpat, E. L. G. Vidoto, D. O. Soares-Pinto, E. R. deAzevedo, and F. F. Fanchini, Sci Rep 5, 14671 (2015).
  * Bottrill  _et al._ [2023] G. Bottrill, M. Pandey, and O. Di Matteo, in  _2023 IEEE International Conference on Quantum Computing and Engineering (QCE)_, Vol. 01 (2023) pp. 177–183.
  * Bravyi  _et al._ [2022] S. Bravyi, A. Kliesch, R. Koenig, and E. Tang, Quantum 6, 678 (2022).
  * Wach  _et al._ [2023] N. L. Wach, M. S. Rudolph, F. Jendrzejewski, and S. Schmitt, Quantum Mach. Intell. 5, 36 (2023).
  * Roca-Jerat  _et al._ [2023] S. Roca-Jerat, J. Román-Roche, and D. Zueco, arXiv:2308.16230 (2023).
  * Useche  _et al._ [2022] D. H. Useche, A. Giraldo-Carvajal, H. M. Zuluaga-Bucheli, J. A. Jaramillo-Villegas, and F. A. González, Quantum Inf Process 21, 12 (2022).
  * Ringbauer  _et al._ [2022] M. Ringbauer, M. Meth, L. Postler, R. Stricker, R. Blatt, P. Schindler, and T. Monz, Nat. Phys. 18, 1053–1057 (2022).
  * Hrmo  _et al._ [2023] P. Hrmo, B. Wilhelm, L. Gerster, M. W. van Mourik, M. Huber, R. Blatt, P. Schindler, T. Monz, and M. Ringbauer, Nature Communications 14, 2242 (2023).
  * Chi  _et al._ [2022] Y. Chi, J. Huang, Z. Zhang, J. Mao, Z. Zhou, X. Chen, C. Zhai, J. Bao, T. Dai, H. Yuan, M. Zhang, D. Dai, B. Tang, Y. Yang, Z. Li, Y. Ding, L. K. Oxenløwe, M. G. Thompson, J. L. O’Brien, Y. Li, Q. Gong, and J. Wang, Nature Communications 13, 1166 (2022).
  * Blok  _et al._ [2021] M. S. Blok, V. V. Ramasesh, T. Schuster, K. O’Brien, J. M. Kreikebaum, D. Dahlen, A. Morvan, B. Yoshida, N. Y. Yao, and I. Siddiqi, Phys. Rev. X 11, 021010 (2021).
  * Zitzler  _et al._ [2003] E. Zitzler, L. Thiele, M. Laumanns, C. M. Fonseca, and V. G. da Fonseca, IEEE Trans. Evol. Comput. 7, 117 (2003).
  * Coello  _et al._ [2007] C. Coello, G. Lamont, and D. van Veldhuizen,  _Evolutionary Algorithms for Solving Multi-Objective Problems_, Genetic and Evolutionary Computation (Springer US, 2007).
  * Beume  _et al._ [2009] N. Beume, C. M. Fonseca, M. López-Ibáñez, L. Paquete, and J. Vahrenhold, IEEE Trans. Evol. Comput. 13, 1075 (2009).
  * Beume [2009] N. Beume, Evol. Comput. 17, 477 (2009).
  * Guerreiro  _et al._ [2022] A. P. Guerreiro, C. M. Fonseca, and L. Paquete, ACM Comput. Surv. 54, 119:1 (2022).
  * Zitzler  _et al._ [2007] E. Zitzler, D. Brockhoff, and L. Thiele, in  _Evolutionary Multi-Criterion Optimization, 4th International Conference, EMO 2007, Matsushima, Japan, March 5-8, 2007, Proceedings_, Lecture Notes in Computer Science, Vol. 4403, edited by S. Obayashi, K. Deb, C. Poloni, T. Hiroyasu, and T. Murata (Springer, 2007) pp. 862–876.
  * Falcón-Cardona  _et al._ [2021] J. G. Falcón-Cardona, S. Z. Martínez, and A. García-Nájera, in  _GECCO ’21: Genetic and Evolutionary Computation Conference, Lille, France, July 10-14, 2021_, edited by F. Chicano and K. Krawiec (ACM, 2021) pp. 395–402.
  * Falcón-Cardona  _et al._ [2022] J. G. Falcón-Cardona, M. T. M. Emmerich, and C. A. C. Coello, Evol. Comput. 30, 381 (2022).
  * Auger  _et al._ [2009] A. Auger, J. Bader, D. Brockhoff, and E. Zitzler, in  _Proceedings of the Tenth ACM SIGEVO Workshop on Foundations of Genetic Algorithms_, FOGA ’09 (Association for Computing Machinery, New York, NY, USA, 2009) p. 87–102.
  * Kasper  _et al._ [2022] V. Kasper, D. González-Cuadra, A. Hegde, A. Xia, A. Dauphin, F. Huber, E. Tiemann, M. Lewenstein, F. Jendrzejewski, and P. Hauke, Quantum Science and Technology 7, 015008 (2022).
  * Giorda  _et al._ [2003] P. Giorda, P. Zanardi, and S. Lloyd, Physical Review A 68, 062320 (2003).
  * Barnett and Croke [2009] S. M. Barnett and S. Croke, Advances in Optics and Photonics 1, 238 (2009).
  * Note [1] The effect of variations in the selected solutions induced by a finite number of measurement shots needs to be investigated.
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, Nature Methods 17, 261 (2020).
  * Hansen [2006] N. Hansen, The cma evolution strategy: A comparing review, in  _Towards a New Evolutionary Computation: Advances in the Estimation of Distribution Algorithms_, edited by J. A. Lozano, P. Larrañaga, I. Inza, and E. Bengoetxea (Springer Berlin Heidelberg, Berlin, Heidelberg, 2006) pp. 75–102.
  * Shibata [2022] M. Shibata, CMA-ES, https://github.com/CyberAgentAILab/cmaes (2022).
  * Kennedy and Eberhart [1995] J. Kennedy and R. Eberhart, in  _Proceedings of ICNN’95 - International Conference on Neural Networks_, Vol. 4 (1995) pp. 1942–1948 vol.4.
  * Storn and Price [1997] R. Storn and K. Price, Journal of Global Optimization 11, 341–359 (1997).
  * Hadka [2022] D. Hadka, Platypus - Multiobjective Optimization in Python, v1.1.0, https://github.com/Project-Platypus/Platypus (2022).
  * Herrman  _et al._ [2022] R. Herrman, P. Lotshaw, J. Ostrowski, T. Humble, and G. Siopsis, Sci Rep 12, 6781 (2022).
  * Saleem  _et al._ [2021] Z. H. Saleem, T. Tomesh, B. Tariq, and M. Suchara, arXiv:2010.06660 (2021).
  * Shi  _et al._ [2022] K. Shi, R. Herrman, R. Shaydulin, S. Chakrabarti, M. Pistoia, and J. Larson, in  _2022 IEEE/ACM 7th Symposium on Edge Computing (SEC)_ (2022) pp. 414–419.
  * van Rijn  _et al._ [2017] S. van Rijn, H. Wang, B. van Stein, and T. Bäck, in  _Proceedings of the Genetic and Evolutionary Computation Conference_, GECCO ’17 (Association for Computing Machinery, New York, NY, USA, 2017) p. 737–744.
  * López-Ibáñez  _et al._ [2016] M. López-Ibáñez, J. Dubois-Lacoste, L. Pérez Cáceres, M. Birattari, and T. Stützle, Operations Research Perspectives 3, 43 (2016).
  * Schütze and Hernández [2021] O. Schütze and C. Hernández, Archiving in evolutionary multi-objective optimization: A short overview, in  _Archiving Strategies for Evolutionary Multi-objective Optimization Algorithms_ (Springer International Publishing, Cham, 2021) pp. 17–20.
  * Li  _et al._ [2023] M. Li, M. López-Ibáñez, and X. Yao, arXiv:2303.09685 (2023).
  * Egger  _et al._ [2021] D. J. Egger, J. Mareček, and S. Woerner, Quantum 5, 479 (2021).
  * Akshay  _et al._ [2021] V. Akshay, D. Rabinovich, E. Campos, and J. Biamonte, Phys. Rev. A 104, L010401 (2021).
  * Sack  _et al._ [2023] S. H. Sack, R. A. Medina, R. Kueng, and M. Serbyn, Phys. Rev. A 107, 062404 (2023).
  * Jain  _et al._ [2022] N. Jain, B. Coyle, E. Kashefi, and N. Kumar, Quantum 6, 861 (2022).
  * Bravyi  _et al._ [2020] S. Bravyi, A. Kliesch, R. Koenig, and E. Tang, Phys. Rev. Lett. 125, 260505 (2020).
  * Santra  _et al._ [2023] G. C. Santra, F. Jendrzejewski, P. Hauke, and D. J. Egger, arxiv:2205.10383 (2023).