# Localization of Malaria Parasites and White Blood Cells in Thick Blood
Smears

Rose Nakasi
Makerere University
<EMAIL_ADDRESS>
Aminah Zawedde
Makerere University
<EMAIL_ADDRESS>
Ernest Mwebaze
Sunbird AI
<EMAIL_ADDRESS>
Jeremy Francis Tusubira
Makerere University
<EMAIL_ADDRESS>
Gilbert Maiga
Makerere University
<EMAIL_ADDRESS>

###### Abstract

Effectively determining malaria parasitemia is a critical aspect in assisting
clinicians to accurately determine the severity of the disease and provide
quality treatment. Microscopy applied to thick smear blood smears is the _de
facto_ method for malaria parasitemia determination. However, manual
quantification of parasitemia is time consuming, laborious and requires
considerable trained expertise which is particularly inadequate in highly
endemic and low resourced areas. This study presents an end-to-end approach
for localisation and count of malaria parasites and white blood cells (WBCs)
which aid in the effective determination of parasitemia; the quantitative
content of parasites in the blood. On a dataset of slices of images of thick
blood smears, we build models to analyse the obtained digital images. To
improve model performance due to the limited size of the dataset, data
augmentation was applied. Our preliminary results show that our deep learning
approach reliably detects and returns a count of malaria parasites and WBCs
with a high precision and recall. We also evaluate our system against human
experts and results indicate a strong correlation between our deep learning
model counts and the manual expert counts (p=0.998 for parasites, p=0.987 for
WBCs). This approach could potentially be applied to support malaria
parasitemia determination especially in settings that lack sufficient
Microscopists.

##  1 Introduction

Accurately determining malaria parasitemia, the quantitative content of
malaria parasites in the blood, is key in determining measure of response to
treatment, likely cause of disease and severity phases of malaria especially
in non-immune patients [21] [11]. However, most highly endemic regions lack
sufficient skilled experts to efficiently determine parasitemia density of
patients and their manual counts of parasites and WBCs are time consuming and
prone to subjectivity and errors [7].

Different methods have been implemented to detect and count malaria parasites
(trophozoites) and blood cells using image analysis techniques, but most have
majorly been applied on thin blood smears [13] [17] [9] [16] [12]. According
to WHO, thick blood smears are the gold standard for the determination of
malaria parasitemia [19]. Thick blood smears allow the microscopist to examine
a larger number of red cells for the presence of parasites, and low
parasitemia can more readily be identified in thick blood smears as compared
to thin films [1]. Some previous work has used a semi-automated method for
parasitemia determination in thick smears by using an open source java based
image processing program _ImageJ_ (version 1.41) that is based on hand
engineered features used to do manual counts of malaria parasites only [4].

This study explores an end-to-end machine learning based algorithm for
detection and counting of trophozoites and WBCs for improved parasitemia
determination in thick blood smears. The pre-trained deep learning approach
utilised is significantly faster, efficient and requires no additional
operator preprocessing of the images for automated detection and counts. This
is well suited for improving malaria parasitemia quantification in resource-
constrained settings with few skilled lab technologists to interpret
microscopy test results.

##  2 Methods

###  2.1 Data acquisition and preparation

Data was collected from a national referral hospital and annotated by malaria
experts. Annotations are in form of bounding boxes around the parasites in an
image of a blood slide. Proper permissions were obtained from both the
patients and the hospital for this data collection.

A standardized data collection protocol was followed for this data collection.
Thick blood slides where the experts had determined the presence of malaria
parasites were used to generate the training data. From each blood slide, the
microscope was adjusted to obtain different view point images. Each image with
the corresponding set-up parameters including the microscope slide number,
stage micrometer grid readings(x and y), zoom level of the smart phone used
for image collection, objective size of the microscope, and the staining
reagent used were captured. Data collected from the lab was sent to a
centralized server using an ODK form on a smartphone device.

To collect the data in the lab, a 3-D printed microscope attachment mechanism
[14] was attached on an Olympus microscope with objective magnification
(1000×). A Samsung J6 model smart phone with an installed phone zoom at 10 was
attached onto the mounted 3-D printed device. Images captured with the phone
and uploaded to the server, were downloaded and manual data annotations were
performed by microscopists using the annotation tool LabelImg [18].
Annotations were saved in the Pascal VOC format [3]. The xml annotation file
contained coordinates of bounding boxes of trophozoites and WBCs creating a
multi-class object detection task. The annotated images were randomly split
into a 9:1 train and test set. The dataset was encoded as a _TFRecord_ record;
the data input format for a training job optimized for utilising the
Tensorflow Object detection API.

###  2.2 Model training and deployment

In this study, we explored the potential of Faster R-CNN [15] for the task of
localization, and counting of trophozoites and WBCs in thick blood smears. The
model is an improvement to the Region based CNN model (R-CNN) whose detection
procedure first uses Region Proposal Networks (RPNs) to extract features using
selective search purposefully to generate proposals (boxes) which are further
cropped and processed. For our case, Faster R-CNN was implemented with the
fully convolutional Residual Network (ResNet 101) [6] as the backbone network.
Comparison in performance between our proposed approach with single-class
object detections (trophozoites-only) and (WBCs-only) was also conducted.
Further, we compared our proposed model with a baseline Single Shot Multi Box
Detector (SSD) [22] for model performance comparison.

Training and testing was implemented on a 5th Gen Intel core i7 series
processor with 16GB RAM and NVIDIA GTX 1060 GPU. Models pre-trained on the
COCO dataset [10] were adopted in a transfer learning scheme to reduce
training time. Data augmentation using random vertical and horizontal flipping
of images was used to improve sample training size. An image size of 750
$\times$ 750 was configured as input to the Faster R-CNN model in a batch size
of 1 (due to our low resources) with a base learning rate kept very low at
0.0003. For the baseline model, SSD MobileNet, images of 300x300 input size
with a batch size of 12 at a learning rate of 0.0004 were considered for a
light model like SSD MobileNet.

Initially, the learning rate was kept at 0.0003 and 0.0004. However beyond
90,000 steps for each model, the learning rate was decreased to 0.00003 and
0.00004 for Faster R-CNN and SSD respectively. A momentum optimizer of 0.9 at
0.5 IOU threshold was maintained. Though 200,000 time steps were specified,
each model attained optimal performance at different numbers of steps (see
Table 1). Checkpoint files were frozen into a protobuff file to generate a
tensorflow inference model.

##  3 Results and Discussion

All models used were applied on a validation set of 100 images. An example of
model inference detections are shown in Figure [a]. It was observed that our
proposed model performs and matches well with the ground truth as shown in
Figure [b]. Results using the common metrics of mAP (mean Average precision)
[2], precision and recall were used and results are summarised [see Table 1].

Table 1<EMAIL_ADDRESS>precision and recall performance for Faster R-CNN model for detection of malaria parasites and WBC. Algorithm | Steps |<EMAIL_ADDRESS>| Precision | Recall
---|---|---|---|---
Faster R-CNN (trophozoites_only) | 138,700 | 0.5506 | 0.6720 | 0.802
Faster R-CNN (WBC_only) | 42,800 | 0.892 | 0.855 | 0.959
Faster R-CNN (trophozoites + WBC) | 98,700 | 0.6609 |  parasites: 0.686 WBC: 0.805 |  parasites:0.9303 WBC: 0.884
SSD MobileNet (trophozoites + WBC) | 80,000 | 0.6292 |  parasites: 0.760 WBC: 0.89 |  parasites:0.501 WBC: 0.806

As shown in Table 1, Faster R-CNN performs better than the SSD MobileNet
baseline model in terms of mAP and recall except in terms of precision. This
model shows the importance and potential of Regional Proposal Networks for
object localisation.

Furthermore, results show that a multi-class detection task of both malaria
trophozoites produces better results than that from single-class tasks
(trophozoites (only)), but lower than of a single-class of WBCs (0.892). This
could be due to the fact that WBCs are bigger objects and as such are easier
examples to learn as compared to smaller objects (trophozoites only). The
better performance of the multi-class task could be attributed to having
bigger objects complementing small objects (trophozoites). Generally, the
multi-class mAP value of 0.6609 (@0.5IOU) was considered to be good enough
performance given the presence of smaller objects(trophozoites).

Furthermore, we show that our model-generated counts generally correlate well
with manual-experts on a held-out dataset of 8 thick blood smear films with
more than 100 images each (high Spearman’s correlation coefficients for
trophozoites (p=0.998) and WBCs (p=0.987)). In Figure 2, model-generated
parasite counts are slightly higher than expert-manual counts, this implies
effect of false-positives that could be due to low quality annotations thus a
lower precision registered for parasites. In general the precision and recall
values from our proposed model is high compared to the baseline model results.

(a) Groundtruth

(b) Faster R-CNN

Figure 1: Ground truth annotations (a) by lab experts Vs inference by Faster-
RCNN model (b). The blue rectangles indicate WBCs and the green rectangles
indicate malaria trophozoites. Figure 2: Model counts Vs manual expert counts
per film

Following from Table 2, film-wise Parasitemia density based on
Parasites/$\mu$l method per film [20][8] was evaluated from the equation;

$Parasites/\mathrm{\mu l}=\frac{\textrm{ No. of counted
parasites}\times\textrm{No. of counted WBC}}{\textrm{No. of assumed
WBC(8000)}}\\\ $

Table 2: Malaria parasitemia on a unique set of each validation thick blood smear slide with corresponding clinical interpretation.  Slide | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8
---|---|---|---|---|---|---|---|---
Model Parasites/$\mu$l | 2326 | 5740 | 3160 | 4682 | 3011 | 7238 | 1906 | 47982

Malaria parasitemia clinical interpretations [5] were deduced, Slides
1,2,3,4,5,6 and 7 represent parasitemia level above which immune patients will
exhibit malaria symptoms whereas Slide 8 represents a patient with maximum
parasitemia.

##  4 Conclusion and future work

In this study, we provide an end-to-end machine learning approach for the
determination of malaria parasite density on thick blood smears. We address
the diagnostic challenges associated with conventional microscopy for
parasitemia determination. Specifically, we leverage transfer learning and
propose a Faster R-CNN ResNet 101 model pre-trained on COCO dataset for the
task of object detection. We show that our proposed model can accomplish
parasite and WBC localization and count effectively in a fast, accurate and
consistent way for malaria parasitemia determination. In future, we intend to
work towards reducing false positives and false negatives generated by the
model. This could be through enforcing quality data annotation.

## References

  * [1] P. Bejon, L. Andrews, A. Hunt-Cooke, F. Sanderson, S.C. Gilbert, and A.V.S. Hill.  Thick blood film examination for plasmodium falciparum malaria has reduced sensitivity and underestimates parasite density.  Malaria Journal, 5(104), January 2006.
  * [2] J. Cartucho.  map.  https://github.com/Cartucho/mAP, 2019.
  * [3] M. Everingham, S. M. A. Eslami, L. Van Gool, C. K. I. Williams, J. Winn, and A. Zisserman.  The pascal visual object classes challenge: A retrospective.  International Journal of Computer Vision, 111(1):98–136., 2015\.
  * [4] J. Frean.  Reliable enumeration of malaria parasites in thick blood films using digital image analysis.  Malaria journal, 8:218, 2009.
  * [5] L.S Garcia.  Parasitemia determined from conventional light microscopy: clinical correlation and interpretation, 2015.
  * [6] K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
  * [7] A.H. Kilian, W.G. Metzger, E.J. Mutschelknauss, G. Kabagambe, P. Langi, R. Korte, and F. von Sonnenburg.  Reliability of malaria microscopy in epidemiological studies: results of quality control.  https://doi.org/10.1046/j.1365-3156.2000.00509.x, January 2000.
  * [8] A.A. Kloub.  Methods of estimating and counting malaria parasites density.  EC Microbiology: 800-803, 2019.
  * [9] M.T. Le, T.R. Bretschneider, C. Kuss, S. Xiang, and P.R. Preiser.  A novel semi-automatic image processing approach to determine plasmodium falciparum parasitaemia in giemsa-stained thin blood smears.  BMC Cell Biol, 9:15, January 2008.
  * [10] T. Lin, M. Maire, S. Belongie, J. Hays, P. Perona, D. Ramanan, P. Dollar, and C. L. Zitnick.  Microsoft coco: Common objects in context.  In European conference on computer vision, pages740–755, 2014\.
  * [11] W. P. Omeara, F. E. Mckenzie, A.J. Magill, J.R. Forney, B. Permpanich, C. Lucas, R.A Gasser, and C. Wongsrichanalai.  Sources of variability in determing malaria parasite density by microscopy".  The American Journal of Tropical Medicine and Hygiene", 73(3):593–598, 2005.
  * [12] M. Poostchi, I. Ersoy, K. McMenamin, E. Gordon, N. Palaniappan, S. Pierce, R.J. Maude, A. Bansal, P. Srinivasan, L. Miller, K. Palaniappan, G. Thoma, and S. Jaeger.  Malaria parasite detection and cell counting for human and mouse using thin blood smear microscopy. journal of medical imaging.  Journal of medical imaging (Bellingham, Wash.), 5(4), 044506, 2018\.
  * [13] 0\. Proudfoot, N. Drew, A. Scholzen, S. Xiang, and M. Plebanski.  Investigation of a novel approach to scoring giemsa-stained malaria-infected thin blood films.  https://doi.org/10.1186/1475-2875-7-62, January 2007.
  * [14] J.A. Quinn, R. Nakasi, P.K.B. Mugagga, P. Byanyima, W. Lubega, and A. Andama.  Deep convolutional neural metworks for microscopy-based point of care diagnosis.  In preceedings of International Conference on Machine Learning for Health Care, Volume 50., 2016.
  * [15] S. Ren, K. He, R. Girshick, and J. Sun.  Faster r-cnn: Towards real-time object detection with region proposal networks.  IEEE Trans. Pattern Anal. Mach. Intell., 39, 1137–1149., 2015\.
  * [16] S.S. Savkare and S. Narote.  Automatic detection of malaria parasites for estimating.  International Journal of Computer Science and Security, pages 2011–310, 01 2011.
  * [17] S.W. Sio, W. Sun, S. Kumar, W.Z. Bin, S.S. Tan, S.H. Ong, H. Kikuchi, Y. Oshima, and K.S.W. Tan.  Malariacount: an image analysis-based program for the accurate determination of parasitemia.  J Microbiol Methods, 68:11-18, January 2007.
  * [18] D. Tzutalin.  Labelimg.  https://github.com/tzutalin/labelImg, 2015.
  * [19] WHO.  Informal consultation on quality control of malaria microscopy., 2006\.
  * [20] WHO.  Malaria parasite counting., 2009.
  * [21] WHO.  Severe malaria.  Tropical Medicine and International Health,19 (Suppl. 1), 7–131, 2014.
  * [22] W.Liu, D. Anguelov, D. Erhan, C. Szegedy, S. Reed, C. Fu, and A.C. Berg.  Ssd: Single shot multibox detector.  In Proceedings of the European Conference on Computer Vision—ECCV, Amsterdam, The Netherlands, pp. 21–37., 2016.


# Gaussian Boson Sampling with Pseudo-Photon-Number Resolving Detectors and
Quantum Computational Advantage

Yu-Hao Deng  Yi-Chao Gu  Hua-Liang Liu  Si-Qiu Gong  Hao Su  Zhi-Jiong Zhang
Hao-Yang Tang  Meng-Hao Jia  Jia-Min Xu  Ming-Cheng Chen  Jian Qin  Li-Chao
Peng  Jiarong Yan  Yi Hu  Hefei National Laboratory for Physical Sciences at
Microscale and School of Physical Sciences, University of Science and
Technology of China, Hefei, Anhui, 230026, China  CAS Centre for Excellence
and Synergetic Innovation Centre in Quantum Information and Quantum Physics,
University of Science and Technology of China, Shanghai, 201315, China  Hefei
National Laboratory, University of Science and Technology of China, Hefei
230088, China  Jia Huang  Hao Li  State Key Laboratory of Functional Materials
for Informatics, Shanghai Institute of Micro system and Information Technology
(SIMIT), Chinese Academy of Sciences, 865 Changning Road, Shanghai, 200050,
China  Yuxuan Li  Yaojian Chen  Department of Computer Science and Technology
and Beijing National Research Center for Information Science and Technology,
Tsinghua University, Beijing, China  Xiao Jiang  Hefei National Laboratory for
Physical Sciences at Microscale and School of Physical Sciences, University of
Science and Technology of China, Hefei, Anhui, 230026, China  CAS Centre for
Excellence and Synergetic Innovation Centre in Quantum Information and Quantum
Physics, University of Science and Technology of China, Shanghai, 201315,
China  Hefei National Laboratory, University of Science and Technology of
China, Hefei 230088, China  Lin Gan  Guangwen Yang  Department of Computer
Science and Technology and Beijing National Research Center for Information
Science and Technology, Tsinghua University, Beijing, China  Lixing You  State
Key Laboratory of Functional Materials for Informatics, Shanghai Institute of
Micro system and Information Technology (SIMIT), Chinese Academy of Sciences,
865 Changning Road, Shanghai, 200050, China  Li Li  Han-Sen Zhong  Hui Wang
Nai-Le Liu  Hefei National Laboratory for Physical Sciences at Microscale and
School of Physical Sciences, University of Science and Technology of China,
Hefei, Anhui, 230026, China  CAS Centre for Excellence and Synergetic
Innovation Centre in Quantum Information and Quantum Physics, University of
Science and Technology of China, Shanghai, 201315, China  Hefei National
Laboratory, University of Science and Technology of China, Hefei 230088, China
Jelmer J. Renema  Adaptive Quantum Optics Group, Mesa+ Institute for
Nanotechnology, University of Twente, P.O. Box 217, 7500 AE Enschede,
Netherlands  Chao-Yang Lu  Hefei National Laboratory for Physical Sciences at
Microscale and School of Physical Sciences, University of Science and
Technology of China, Hefei, Anhui, 230026, China  CAS Centre for Excellence
and Synergetic Innovation Centre in Quantum Information and Quantum Physics,
University of Science and Technology of China, Shanghai, 201315, China  Hefei
National Laboratory, University of Science and Technology of China, Hefei
230088, China  New Cornerstone Science Laboratory, Shenzhen 518054, China
Jian-Wei Pan  Hefei National Laboratory for Physical Sciences at Microscale
and School of Physical Sciences, University of Science and Technology of
China, Hefei, Anhui, 230026, China  CAS Centre for Excellence and Synergetic
Innovation Centre in Quantum Information and Quantum Physics, University of
Science and Technology of China, Shanghai, 201315, China  Hefei National
Laboratory, University of Science and Technology of China, Hefei 230088, China

###### Abstract

We report new Gaussian boson sampling experiments with pseudo-photon-number-
resolving detection, which register up to 255 photon-click events. We consider
partial photon distinguishability and develop a more complete model for the
characterization of the noisy Gaussian boson sampling. In the quantum
computational advantage regime, we use Bayesian tests and correlation function
analysis to validate the samples against all current classical spoofing
mockups. Estimating with the best classical algorithms to date, generating a
single ideal sample from the same distribution on the supercomputer Frontier
would take $\sim 600$ years using exact methods, whereas our quantum computer,
Jiǔzhāng 3.0, takes only $1.27\text{\,}\mathrm{\SIUnitSymbolMicro s}$ to
produce a sample. Generating the hardest sample from the experiment using an
exact algorithm would take Frontier $\sim$3.1\text{\times}{10}^{10}$$ years.

Quantum computational advantage (QCA) [1, 2, 3, 4, 5] marks an important
milestone in the development of quantum computers. By solving certain quantum
sampling problems [6, 7] that are intractable for classical supercomputers,
QCA experiments [8, 9, 10, 11, 12] have provided strong evidence for the long-
anticipated quantum speed-up theoretically conceived $\sim 40$ years ago,
Similar to Bell tests [13] which were designed to refute the local hidden
variable theories [14], these QCA experiments offered evidence of the
violations of the Extended Church-Turing thesis [3].

Interestingly, these QCA experiments have motivated growing study of faster
classical simulations. For boson sampling, these efforts can be divided into
three approaches. The first approach is reducing the classical simulation
overhead of exactly simulating an ideal implementation of the Gaussian boson
sampling protocol [15, 16, 17]. In this case, one expects that the exponential
gap between the classical and the quantum persists, but it is possible to
bring down the overhead. The leading result is Ref. [17] which uses photon
collisions to achieve better exact classical simulation that could reduce the
computational overhead by orders of magnitudes lower than brute-force
algorithms. The second approach is using imperfections in the experimental
setup [18, 19, 20], such as photon loss and partial photon distinguishability,
to come to a faster simulation of the experimental data produced by an
imperfect experiment. In this case, one can hope for more efficient classical
simulation if the level of imperfections or experimental noise become strong
enough. The final approach is spoofed distributions [17, 21, 22, 23], i.e.,
distributions which are not constructed to be close in variational distance to
the ideal boson sampling distribution, but which are designed to reproduce
some statistical aspect of the experiment. In particular, the recently
proposed more competitive squashed state mockup [12] exhibits better agreement
with the ground truth than the thermal state mockup. Moreover, the treewidth
sampler [22], designed to spoof experiments with restricted circuit-
connectivity, and the Independent Pairs and Singles (IPS) sampler [17],
designed to emulate experiments of limited quantum interference, have not been
reported to be ruled out in previous GBS QCA experiments.

In return, these new challenges from the classical counterparts motivate the
development of higher-fidelity and larger-scale quantum computers, new methods
for validation and characterization [24, 25, 26, 27], as well as better
modeling and understanding of the increasingly complex system [28, 29, 30, 31,
32, 33, 34, 35, 36], which is a fundamental endeavor in its own right [37, 38,
39]. Indeed, only by such a continuous quantum-classical competition can the
QCA milestone be progressively better established.

Figure 1: The experimental setup. 25 stimulated two-mode squeezed state photon
sources are all phase-locked to each other and sent into a 144-mode ultralow-
loss fully-connected optical interferometer. The photons go through 72 units
of fiber loop setups for temporal-spatial demultiplexing and are detected by
144 superconducting nanowire single-photon detectors, which together
constitute the pseudo-photon-number resolving detection scheme. Each fiber
loop setup includes two input modes as represented by distinct colors. Photons
from each mode are temporally demultiplexed by the fiber beam splitters and
the delay lines into four time bins, and each time bin is furthermore split
into two path bins at the last fiber beam splitter. The photons corresponding
to each of the two input modes of the same fiber loop setup can be
distinguished by their parity in time bin with a coincidence event analyzer
(not shown).

In this direction, we report in this Letter a new, higher-efficiency Gaussian
boson sampling (GBS) [40, 41] machine with up to 255 photon clicks in the
output using pseudo-photon-number-resolving detectors (PPNRD), which overcomes
both the most powerful exact classical simulation algorithms and all known
approximate algorithms and spoofing algorithms. The PPNRD scheme overcomes the
shortcoming of threshold detection that it is incapable of resolving photon-
number information, and significantly increases the exact sampling task’s
computational complexity by more than six orders of magnitude. We validate the
samples against all the emerging classical mockups, particularly including the
squashed state [12], the treewidth sampler [22] and the IPS sampler [17]
mentioned above. A new model that includes the partial photon
indistinguishability is used to characterize the system, and exhibits better
agreement with the experiment than previous modelling methods. The
computational complexity of this new GBS device, Jiǔzhāng 3.0, is analyzed and
a new QCA frontier is established.

The first GBS experiments [9, 10] in the QCA regime used threshold detectors
to register the samples. In those cases, there was a possibility of photon
collision, that is, the photons can bunch at the outputs. New classical
algorithms [17] could exploit photon collision to reduce the simulation
overhead. A recent work [12] has reported time-bin-encoded GBS with photon-
number resolved detection, in a fiber loop-based configuration similar to an
earlier single-photon boson sampling experiment [42]. However, the relatively
high photon loss in the fiber loops has allowed only three loops in the
implementation, restricting the depth and universality of the interferometer.
Additionally, the long recovery time of the transition-edge sensors could
render it an unsuitable choice for high-repetition-rate experiments.

A schematic of our GBS experiment is shown in Fig. 1. Transform-limited laser
pulses double pass periodically poled potassium titanyl phosphate (PPKTP)
crystals to create 25 pairs of two-mode squeezed state (TMSS) by a stimulated
emission process [10]. The TMSSs have an average coupling efficiency of 88.4%
and photon indistinguishability of 96.2%, simultaneously. In the experiment,
the laser power is tuned to generate different squeezing parameters ranging
from 1.2 to 1.6. The TMSSs are then fed into an ultralow-loss three-
dimensional interferometer with full connectivity among all 144 modes. The
transmission rate of the interferometer is 97% for each mode, and the average
wave-packet overlap inside the interferometer is above 99.5%. The whole set-up
is actively phase locked within a precision of 15 nm.

For detection, we implement pseudo-photon number resolution using one-to-eight
demultiplexing of the optical modes (see Fig. 1). We conduct detector
tomography [43] of each mode to validate the pseudo-photon-number-resolved
detection (PPNRD) scheme [44]. As shown in Fig. S1(a), the probability of the
nine photon-click number response of the PPNRD agrees well with the
theoretical model. The positive-operator valued measures of the nine photon-
click number measurement can be constructed from the data, which is shown in
Fig. S1(c), together with its Wigner function and fidelity with the
theoretical prediction. While we model the experiment with the PPNRD detection
scheme in the following discussions, we address the question of how well our
detection scheme emulates true photon number-resolving detection in [44].

Figure 2: The experimental photon-click number distribution and the Bayesian
validation results. (a) Photon-click number distribution of this work. Data
from experiments of three pump laser power ranging from 0.72 W and 1.30 W are
displayed. A maximum photon-click number of 129, 203, and 255 are registered
for each of the experiments. (b) Photon-click number distribution of the
experimental results, the ground-truth theory, the squashed state and the
thermal state mockups for the lowest laser intensity configuration. The plot
includes labels indicating the standard deviation of the photon-click number
for each distribution. It is evident that the experimental distribution agrees
with the ground truth best. The error bars are too small to be displayed.
(c-d) Bayesian confidence of the ground-truth theory as a function of
subsystem size, against the thermal state (c) and the squashed state (d)
hypothesis in the high power (1.30 W) experiment. It can be observed that as
the subsystem size grows, the Bayesian confidence is above zero and exhibits a
clearly increasing trend, indicating an exclusion of the two mockups and a
stronger Bayesian confidence for the full system. Bayesian results for photon
click number changing from 20 to 26 are displayed in data points of different
colors and show similar results. (e-f) Bayesian confidence of the ground-truth
theory as a function of the subsystem size for experiments of the three
distinct pump laser power. The results show the Bayesian results against the
thermal state (e) and the squashed state (f) hypothesis, for a fixed photon-
click number (25). For all pump laser powers, the Bayesian $\Delta H$ is above
zero and exhibits a clearly increasing trend with subsystem size, implying
stronger Bayesian confidence for the full system. In all plots error bars
indicate standard error. The statistical fluctuation from photon-click number
probability estimated with Monte-Carlo methods [28] has been included in the
error bars. For results of each mode number, the Bayesian score are obtained
by averaging over an ensemble of randomly selected subsystems for unbiased
benchmarking.  Figure 3: Correlation function analysis. The experiment clearly
exhibits a better agreement with the ground truth than the mockups, with two-
norm distance $D$ and slope $K$ (marked in colored dashed lines on each plot).
(a-b) 3-order correlation function in scatter plot showing the experiment
($D=0.0054,K=1.008$), the greedy sampler (a) ($D=0.0062,K=1.650$), and the TAP
sampler (b) ($D=0.0067,K=1.501$) in comparison versus the ground-truth theory.
For (a) and (b), 10 million samples are used for estimation of the
experimental and mockup’s correlation. All 3-order correlation functions of
144 out of 1152 modes (1 from each 1-to-8 fan-out modes) are displayed. (c-d)
Scatter plot of the 2-order correlation function of the experiment ($D=0.045$,
$K=1.004$), the IPS sampler (c) ($D=0.254$, $K=1.544$), and the treewidth
sampler (d) ($D=0.739$, $K=0.733$), versus the ground-truth theory. In (c-d),
all pairs of 1152 modes are displayed, and 10 million samples are used for
estimation of the experimental correlation whereas exact theoretical value is
used for the mockup samplers. (e-f) The 2-order correlation function of the
experiment ($D=0.040,K=1.006$), the thermal state mockup (e)
($D=0.275,K=2.108$), and the squashed state mockup (f) ($D=0.052,K=1.110$)
versus the ground truth theory, in coarse-grained 144 modes. In (e-f), 10
million samples are used for estimation of the experimental correlation and
the mockup samplers’ correlation. All data are from the highest power (1.30W)
experiment.

The overall linear efficiency of the whole experimental set-up is 43%,
including the quantum light sources, transmission, and detection. This
efficiency is much higher than that in ref. [12], which was 33% although
involving only three loops. The photon-click number distribution under three
laser powers is shown in Fig. 2(a). The maximum photon-click number reaches
255, which is higher than all the previous GBS experiments.

We analyze the obtained GBS samples and validate them against known classical
hypotheses in the QCA regime. The most powerful method to spoof the GBS is to
use classical states which can maximally approximate the experiment’s quantum
light sources which, under photon loss, can be gradually degraded into
squeezed thermal states [18]. Thus, in this work, in addition to the thermal
state hypothesis which has been tested in Jiǔzhāng 1.0 and 2.0, we consider
the more competitive hypothesis, namely the squashed state [12, 18]. The
squashed state is a classical mixture of coherent states and has vacuum
fluctuations in one quadrature and larger fluctuation in the other. It can be
optimized to have the same mean photon number as the input lossy squeezed
state [45], which is the convention we take in this work. Note that there are
other plausible hypotheses such as using coherent light and distinguishable
photons, which are much easier to be ruled out. As shown in [44], the
validation strength against these other hypotheses is typically 2-4 orders of
magnitude larger than the squashed states. Therefore, we will focus on the
discussion of the latter in the main text.

We start by comparing the experimental photon number distribution with the
ground-truth model and the possible mockups in Fig. 2(b). The GBS data set
taken at $0.72\text{\,}\mathrm{W}$ laser power (dot) well overlaps with the
ground truth (red), while obviously deviates from the classical mockups based
on the thermal state (yellow) and the squashed state (blue). These
distributions can be quantitatively distinguished by their standard deviation
as labelled on the plot, where the experiment exhibits the best agreement with
the ground truth.

Then, we continue with the Bayesian test, [46], where two hypothetical
theoretical models are compared against each other based on the likelihood to
generate the experimental samples. We define the Bayesian test score $\Delta
H$ as the difference between the ground-truth model $H_{0}$ and the classical
adversary model hypothesis $H_{1}$ on a set of $n$-photon-click samples:

$\Delta{H}={\frac{1}{N}}\mathrm{ln}\prod_{i=1}^{N}{\frac{{{P}}^{(0)}({\vec{{s}}}_{i}){{P}}^{(1)}(n)}{{P}^{(1)}({\vec{{s}}}_{i}){{P}}^{(0)}(n)}}$
(1)

where $N$ is the number of samples, $\vec{s}_{i}$ represents the $i$th sample,
${P}^{(0/1)}({\vec{{s}}}_{i})$ represents the event probability of sample
$\vec{s}_{i}$ under the $H_{0/1}$ hypothesis , and ${P}^{(0/1)}(n)$ is the
coarse-grained probability for the $n$-photon-click under hypothesis
$H_{0/1}$. When $\Delta H>0$, it is proven that the experimental samples are
more likely from the ground-truth GBS rather than the mockup. A higher
Bayesian test score indicates larger confidence.

For the validation, each of the 144 modes is treated as 1-to-8 fan out,
yielding a total mode number of 1152 for the ground truth. Further, in our
model we consider the noise of the partial photon distinguishability for the
Bayesian test, where the probability of the sample event is calculated by a
modified version of the Torontonian [44]. We find that our model gives a
closer description of the experiment than the previous ones. We use the Sunway
TaihuLight supercomputer to calculate the probability of large photon number
samples in the QCA regime.

The Bayesian test score of the ground truth against the mockups using the
thermal state and the squashed state at $1.3\text{\,}\mathrm{W}$ pump power
are shown in Fig. 2(c)-2(d). Due to the classical computational overhead, we
first start from a subsystem with fewer output bosonic modes, and gradually
increase the subsystem size. We observe not only all the Bayesian scores in
Fig. 2(c)-2(d) are higher than zero, but also they show an evident increasing
trend as the subsystem size ramps up. The positive scores demonstrate that the
experimental GBS samples are more likely generated from the ground-truth
distribution rather than these mockup distributions. More importantly, the
rising trend of the Bayesian score indicates a higher score can be inferred
for the full system. This increasing trend could be understood from that, as
subsystem size increases, more complete information of the full system is
incorporated into the Bayesian test, and therefore results in stronger
validation strength. Therefore, we conclusively infer that for the full-mode
system, though it is not directly computable, stronger Bayesian confidence for
the ground-truth theory is expected.

We continue to investigate the power dependence of the Bayesian test. In Fig.
2(e)-2(f), the Bayesian validation strength is plotted for three different
pump laser powers from $0.72\text{\,}\mathrm{W}$ to $1.3\text{\,}\mathrm{W}$.
The results show that the lower pump power can generate higher Bayesian score.
This is expected because the thermal noise and photon loss become increasingly
sensitive with larger squeezing parameters [18].

Another important tool for characterizing the GBS is the correlation function.
The $k$-order correlation function is recursively defined as

$\kappa(X_{1}X_{2}\cdots X_{k})=E(X_{1}X_{2}\cdots X_{k})-\sum_{p\in
P_{k}}\prod_{b\in p}\kappa[(X_{i})_{i\in b}]$ (2)

where $X_{k}$ represents the experimental measurement operator at the $k$th
output mode, and $P_{k}$ represent all partitions of the set
$\\{1,2,\cdots,k\\}$ excluding the universal set. Based on the correlation
function, we validate our experimental samples against other recently proposed
mockups based on approximating the experimental GBS through low-order marginal
distributions [21].

Reference [21] proposed a greedy algorithm to sample from a distribution that
approximates all first- and second-order correlations of the experiment, which
we call a greedy sampler below. The same work also used the single-mode
marginals and two-mode correlations to configure a Boltzmann machine for
sampling under the Thouless, Anderson, and Palmer (TAP) mean field
approximation. These methods can be generalized to higher orders in principle,
but a full enumeration of all the marginal probabilities of the chosen order
is required, which grow combinatorically, and thus limits it within the low-
order regime. However, these mockups do not capture higher-order correlations,
and can thus be ruled out from this aspect.

In Fig. 3(a)-3(b), we directly compare the third-order correlation functions
between the experimental samples at the QCA regime and the mockup samplers
from the order-2 greedy and the order-2 TAP samplers. The statistics of the
experimental samples agree with the ground truth where they cluster around the
identity line at $45\text{\,}\mathrm{\SIUnitSymbolDegree}$, whereas the mockup
samplers’ third-order correlations show significant divergence from the
ground-truth predictions.

Additionally, Ref. [17] designed the IPS sampler, which was shown to yield
higher heavy output generation (HOG) score than the mockups using the thermal
states and distinguishable photons. We compare the second-order cumulants of
the sampler with the experimental results in Fig. 3(c), where the IPS sampler
shows an evident deviation from the results of the ground truth and the
experimental data.

Figure 4: Quantum computational advantage. The distribution of the estimated
classical simulation time on Frontier with exact methods. We calculate the
paired value $(G,N)$ for each of the samples obtained in the highest power
experiment, and construct the two-dimensional heat map. $G>2$ indicates the
feature of photon number resolving. The hardest sample is marked in star.
Equal simulation time given by equation (3) is shown as the lines. As a
comparison, it takes $1.27\text{\,}\mathrm{\SIUnitSymbolMicro s}$ to produce a
sample in our experiment.

The last mockup to be ruled out is the recently proposed treewidth sampler
[22] which can use local connectivity of the circuits to reduce the overhead
of classical approximate sampling. The sampler could generate higher HOG score
than some GBS experiments, which, however, is due to the limitation of HOG
test itself. In fact, one can show that the treewidth mockup sampler can even
yield a higher HOG score than an ideal ground-truth sampler (see Supplementary
Materials [44]), therefore makes the test unreliable [23] (related arguments
on the limitation of the cross-entropy benchmarking in random circuit sampling
have also been reported [47, 48]). Thus, it is necessary to rule out the
treewidth mockup in other ways. In Fig. 3(d), we show the second-order
correlation of a treewidth mockup (with a propagation length of 65) and the
experimental sampler. The former deviates significantly from the ground-truth
theory, and is thus unambiguously ruled out.

As a complementary test, we also show in 3(e) the second order correlation of
the experiment, the thermal state and the squashed state mockup. The
correlation of the experiment, with two-norm distance $D=0.040$ and slope
$K=1.006$, clearly agrees with the ground truth better the thermal state
mockup ($D=0.275,K=2.108$) and the squashed mockup ($D=0.052,K=1.110$).

Having excluded all currently proposed mockups, we now benchmark the classical
computational cost to simulate a noiseless version of our GBS experiment. Our
pseudo-photon number resolving detection scheme can be modelled by treating
all the 1152 fan-out bins as individual output modes, each with threshold
detection. We can use results from [17] to exploit the fact that different
output bins corresponding to the same optical mode result in repeated entries
in the loop hafnian during the simulation, further reducing the dominant
computational cost to:

$T(\vec{N})=\frac{1}{2}C_{\rm{Frontier}}MN^{3}G^{N/2}$ (3)

where $\vec{N}=\\{n_{1},n_{2},\cdots,n_{144}\\}$ represents the PPNRD sample
with $n_{i}$ being the photon-click number of the $i$th mode, $N$ is the
number of clicked modes, $M$ is the mode number,
$G=\left(\Pi_{i}^{M}(n_{i}+1)\right)^{1/N}$, and $C_{\rm{Frontier}}$ is
estimated based on [12]. The addition of PPNRD therefore substantially
increases the computational complexity of our experiment compared to threshold
detection (which always has $G=2$), due to the strong increase in the number
of registered detection events.

We now estimate the time cost on Frontier, currently the most powerful
supercomputer. For each sample of our high power experiment, we estimate it
would on average take Frontier at least $\sim 600$ years to generate using the
exact methods, while it only takes our machine
$1.27\text{\,}\mathrm{\SIUnitSymbolMicro s}$ to produce a sample, showing an
overwhelming QCA of $1.5\text{\times}{10}^{16}$. Moreover, the hardest sample
from our experiment would take Frontier more than $\sim$
$3.1\text{\times}{10}^{10}$ years to generate using exact algorithm. We show
in Fig. 4 the distribution of the estimated classical overhead of all the
experimental samples. We hope future work will further consider the realistic
experimental imperfection such as photon loss and partial photon
distinguishability for a better benchmark [49].

###### Acknowledgements.

We thank Changhun Oh, Minzhao Liu and Liang Jiang for interesting discussions.
We thank Z.-M. He, R.-X. Wang, X.-Y. Wu, Z.-M. Zhang and S.-T. Zheng for
helpful assistance during the experiment. We thank Y. Liang and Q.-H. Shi for
their help in preparation of the manuscript. This work was supported by the
National Key R&D Program of China (2019YFA0308700), the National Natural
Science Foundation of China, the Chinese Academy of Sciences, the Anhui
Initiative in Quantum Information Technologies, the Science and Technology
Commission of Shanghai Municipality(2019SHZDZX01), the Innovation Program for
Quantum Science and Technology (No. ZD0202010000), the XPLORER PRIZE, the New
Cornerstone Science Foundation, and NWO Veni. Y.-H. D., Y.-C. G., H.-L. L.,
S.-Q. G. and H. S. contributed equally to this work. The dataset of this
experiment can be accessed via [50]. Note added \- After submission, we became
aware of a recent related theoretical proposal for balanced-fan-out PPNRD GBS
[51].

## References

  * Bernstein and Vazirani [1993] E. Bernstein and U. Vazirani, in  _Proceedings of the Twenty-Fifth Annual ACM Symposium on Theory of Computing - STOC ’93_ (ACM Press, San Diego, California, United States, 1993) pp. 11–20.
  * Preskill [2012] J. Preskill, arXiv:1203.5813 (2012).
  * Harrow and Montanaro [2017] A. W. Harrow and A. Montanaro, Nature 549, 203 (2017).
  * Lund  _et al._ [2017] A. P. Lund, M. J. Bremner, and T. C. Ralph, npj Quantum Information 3, 15 (2017).
  * Hangleiter and Eisert [2023] D. Hangleiter and J. Eisert, arXiv:2206.04079 (2023).
  * Aaronson and Arkhipov [2011] S. Aaronson and A. Arkhipov, in  _Proceedings of the 43rd Annual ACM Symposium on Theory of Computing - STOC ’11_ (ACM Press, New York, New York, USA, 2011) pp. 333–342.
  * Boixo  _et al._ [2018] S. Boixo, S. V. Isakov, V. N. Smelyanskiy, R. Babbush, N. Ding, Z. Jiang, M. J. Bremner, J. M. Martinis, and H. Neven, Nature Physics 14, 595 (2018).
  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, _et al._ , Nature 574, 505 (2019).
  * Zhong  _et al._ [2021a] H. S. Zhong, H. Wang, Y. H. Deng, M. C. Chen, L. C. Peng, _et al._ , Science 370, 1460 (2021a).
  * Zhong  _et al._ [2021b] H.-S. Zhong, Y.-H. Deng, J. Qin, H. Wang, M.-C. Chen, _et al._ , Physical Review Letters 127, 180502 (2021b).
  * Wu  _et al._ [2021] Y. Wu, W.-S. Bao, S. Cao, F. Chen, M.-C. Chen, _et al._ , Physical Review Letters 127, 180501 (2021).
  * Madsen  _et al._ [2022] L. S. Madsen, F. Laudenbach, M. F. Askarani, F. Rortais, T. Vincent, J. F. F. Bulmer, _et al._ , Nature 606, 75 (2022).
  * Bell [1964] J. S. Bell, Physics Physique Fizika 1, 195 (1964).
  * Einstein  _et al._ [1935] A. Einstein, B. Podolsky, and N. Rosen, Physical Review 47, 777 (1935).
  * Ágoston Kaposi  _et al._ [2022] Ágoston Kaposi, Z. Kolarovszki, T. Kozsik, Z. Zimborás, and P. Rakyta, arXiv:2109.04528 (2022).
  * Quesada  _et al._ [2022] N. Quesada, R. S. Chadwick, B. A. Bell, J. M. Arrazola, T. Vincent, H. Qi, and R. García-Patrón, PRX Quantum 3, 010306 (2022).
  * Bulmer  _et al._ [2022] J. F. F. Bulmer, B. A. Bell, R. S. Chadwick, A. E. Jones, D. Moise, _et al._ , Science Advances 8, eabl9236 (2022).
  * Qi  _et al._ [2020] H. Qi, D. J. Brod, N. Quesada, and R. García-Patrón, Physical Review Letters 124, 100502 (2020).
  * Renema [2020] J. J. Renema, Physical Review A 101, 063840 (2020).
  * Shi and Byrnes [2022] J. Shi and T. Byrnes, npj Quantum Information 8, 54 (2022).
  * Villalonga  _et al._ [2022] B. Villalonga, M. Y. Niu, L. Li, H. Neven, J. C. Platt, V. N. Smelyanskiy, and S. Boixo, arXiv:2109.11525 (2022).
  * Oh  _et al._ [2022a] C. Oh, Y. Lim, B. Fefferman, and L. Jiang, arXiv:2110.01564 (2022a).
  * Oh  _et al._ [2022b] C. Oh, L. Jiang, and B. Fefferman, arXiv:2210.15021 (2022b).
  * Shchesnovich [2021] V. Shchesnovich, Quantum 5, 423 (2021).
  * Dellios  _et al._ [2022a] A. S. Dellios, M. D. Reid, B. Opanchuk, and P. D. Drummond, arXiv:2211.03480 (2022a).
  * Seron  _et al._ [2022] B. Seron, L. Novo, A. Arkhipov, and N. J. Cerf, arXiv:2212.09643 (2022).
  * Giordani  _et al._ [2023] T. Giordani, V. Mannucci, N. Spagnolo, M. Fumero, A. Rampini, E. Rodolà, and F. Sciarrino, Quantum Science and Technology 8, 015005 (2023).
  * Drummond  _et al._ [2022] P. D. Drummond, B. Opanchuk, A. Dellios, and M. D. Reid, Phys. Rev. A 105, 012427 (2022).
  * Popova and Rubtsov [2021] A. S. Popova and A. N. Rubtsov, arXiv:2106.01445 (2021).
  * Grier  _et al._ [2022] D. Grier, D. J. Brod, J. M. Arrazola, M. B. d. A. Alonso, and N. Quesada, Quantum 6, 863 (2022).
  * Lim and Oh [2022] Y. Lim and C. Oh, arXiv:2211.07184 (2022).
  * Dellios  _et al._ [2022b] A. Dellios, P. D. Drummond, B. Opanchuk, R. Y. Teh, and M. D. Reid, Physics Letters A 429, 127911 (2022b).
  * Iosue  _et al._ [2022] J. T. Iosue, A. Ehrenberg, D. Hangleiter, A. Deshpande, and A. V. Gorshkov, arXiv:2209.06838 (2022).
  * Qiao  _et al._ [2022] Y. Qiao, J. Huh, and F. Grossmann, arXiv:2210.09915 (2022).
  * Liu  _et al._ [2023] M. Liu, C. Oh, J. Liu, L. Jiang, and Y. Alexeev, Complexity of Gaussian boson sampling with tensor networks (2023).
  * Oh  _et al._ [2023a] C. Oh, L. Jiang, and B. Fefferman, arXiv:2301.11532 (2023a).
  * Rahimi-Keshari  _et al._ [2016] S. Rahimi-Keshari, T. C. Ralph, and C. M. Caves, Physical Review X 6, 021039 (2016).
  * Deshpande  _et al._ [2018] A. Deshpande, B. Fefferman, M. C. Tran, M. Foss-Feig, and A. V. Gorshkov, Physical Review Letters 121, 030501 (2018).
  * Chabaud and Walschaers [2023] U. Chabaud and M. Walschaers, Physical Review Letters 130, 090602 (2023).
  * Hamilton  _et al._ [2017] C. S. Hamilton, R. Kruse, L. Sansoni, S. Barkhofen, C. Silberhorn, and I. Jex, Phys. Rev. Lett. 119, 170501 (2017).
  * Quesada  _et al._ [2018] N. Quesada, J. M. Arrazola, and N. Killoran, Phys. Rev. A 98, 062322 (2018).
  * He  _et al._ [2017] Y. He, X. Ding, Z.-E. Su, H.-L. Huang, J. Qin, C. Wang, _et al._ , Physical Review Letters 118, 190501 (2017).
  * Lundeen  _et al._ [2009] J. S. Lundeen, A. Feito, H. Coldenstrodt-Ronge, K. L. Pregnell, Ch. Silberhorn, T. C. Ralph, J. Eisert, M. B. Plenio, and I. A. Walmsley, Nature Physics 5, 27 (2009).
  * [44] See Supplemental Materials for detailed description of 1. Detector tomography of the pseudo-photon-number re- solving detectors, 2. Measurement for the matrix phase, 3. Partial photon indistinguishability model, 4. Bayesian test results against weaker classical hypotheses, 5. On HOG test with the treewidth mockup, which includes Ref. [52, 53].
  * Martínez-Cifuentes  _et al._ [2022] J. Martínez-Cifuentes, K. M. Fonseca-Romero, and N. Quesada, arXiv:2207.10058 (2022).
  * Bentivegna  _et al._ [2014] M. Bentivegna, N. Spagnolo, C. Vitelli, D. J. Brod, A. Crespi, _et al._ , International Journal of Quantum Information 12, 1560028 (2014).
  * Barak  _et al._ [2020] B. Barak, C.-N. Chou, and X. Gao, arXiv preprint arXiv:2005.02421 (2020).
  * Gao  _et al._ [2021] X. Gao, M. Kalinowski, C.-N. Chou, M. D. Lukin, B. Barak, and S. Choi, arXiv:2112.01657 (2021).
  * oh [2] After submission, we became aware of a newly developed faster classical approximate simulation algorithm which exhibits error level of the HOG score and up-to 5-order correlation function comparable to the experiment, by taking advantage of the photon loss [54]. Ref. [54] takes $\sim 10$ minutes to generate the tensor network state corresponding to a specific experiment, while in our experiment it takes $1.27\text{\,}\mathrm{\SIUnitSymbolMicro s}$ to generate a sample and $\sim 0.1$ ms to reconfigure a new experiment (in Jiǔzhāng 2.0, only limited by the speed of the physical elements). In addition, Ref. [54] used a simplified model of ground truth for benchmarking of this work, without taking pseudo-photon detection and partial distinguishability into account. We expect future GBS experiments with higher efficiency, as well as better characterization of the experimental imperfections, in particular the calibration noise, could also help enhance the quantum advantage against the approximate simulation algorithms.
  * JZ [3] Raw data of Jiuzhang 3.0.
  * Bressanini  _et al._ [2023] G. Bressanini, H. Kwon, and M. Kim, arXiv preprint arXiv:2305.00853 (2023).
  * Akhlaghi  _et al._ [2011] M. K. Akhlaghi, A. H. Majedi, and J. S. Lundeen, Optics express 19, 21305 (2011).
  * Christ  _et al._ [2011] A. Christ, K. Laiho, A. Eckstein, K. N. Cassemiro, and C. Silberhorn, New Journal of Physics 13, 033027 (2011).
  * Oh  _et al._ [2023b] C. Oh, M. Liu, Y. Alexeev, B. Fefferman, and L. Jiang, arXiv preprint arXiv:2306.03709 (2023b).