Further author information: (Send correspondence to O.G.)
O.G.: E-mail<EMAIL_ADDRESS>
# High Contrast Imaging at the Photon Noise Limit with WFS-based PSF
Calibration

Olivier Guyon  Subaru Telescope, National Astronomical Observatory of Japan,
National Institutes of Natural Sciences (NINS), 650 North A‘ohōkū Place, Hilo,
HI 96720, United States  Steward Observatory, University of Arizona, Tucson,
AZ 87521, United States  College of Optical Sciences, University of Arizona,
Tucson, AZ 87521, United States  Astrobiology Center of NINS, 2 Chome-21-1,
Osawa, Mitaka, Tokyo, 181-8588, Japan  Barnaby Norris  Sydney Institute for
Astronomy, Institute for Photonics and Optical Science, School of Physics,
University of Sydney, NSW 2006, Australia  Marc-Antoine Martinod  Sydney
Institute for Astronomy, Institute for Photonics and Optical Science, School
of Physics, University of Sydney, NSW 2006, Australia  Kyohoon Ahn  Subaru
Telescope, National Astronomical Observatory of Japan, National Institutes of
Natural Sciences (NINS), 650 North A‘ohōkū Place, Hilo, HI 96720, United
States  Vincent Deo  Subaru Telescope, National Astronomical Observatory of
Japan, National Institutes of Natural Sciences (NINS), 650 North A‘ohōkū
Place, Hilo, HI 96720, United States  Nour Skaf  Subaru Telescope, National
Astronomical Observatory of Japan, National Institutes of Natural Sciences
(NINS), 650 North A‘ohōkū Place, Hilo, HI 96720, United States  Observatoire
de Paris, LESIA, 5 Place Jules Janssen, 92190 Meudon, France  Department of
Physics and Astronomy, University College London, London, United Kingdom
Julien Lozi  Subaru Telescope, National Astronomical Observatory of Japan,
National Institutes of Natural Sciences (NINS), 650 North A‘ohōkū Place, Hilo,
HI 96720, United States  Sébastien Vievard  Subaru Telescope, National
Astronomical Observatory of Japan, National Institutes of Natural Sciences
(NINS), 650 North A‘ohōkū Place, Hilo, HI 96720, United States  Astrobiology
Center of NINS, 2 Chome-21-1, Osawa, Mitaka, Tokyo, 181-8588, Japan
Observatoire de Paris, LESIA, 5 Place Jules Janssen, 92190 Meudon, France
Sebastiaan Y. Haffert  Steward Observatory, University of Arizona, Tucson, AZ
87521, United States  Thayne Currie  Subaru Telescope, National Astronomical
Observatory of Japan, National Institutes of Natural Sciences (NINS), 650
North A‘ohōkū Place, Hilo, HI 96720, United States  Jared R. Males  Steward
Observatory, University of Arizona, Tucson, AZ 87521, United States  Alison
Wong  Sydney Institute for Astronomy, Institute for Photonics and Optical
Science, School of Physics, University of Sydney, NSW 2006, Australia  Peter
Tuthill  Sydney Institute for Astronomy, Institute for Photonics and Optical
Science, School of Physics, University of Sydney, NSW 2006, Australia

###### Abstract

Speckle Noise is the dominant source of error in high contrast imaging with
adaptive optics system. We discuss the potential for wavefront sensing
telemetry to calibrate speckle noise with sufficient precision and accuracy so
that it can be removed in post-processing of science images acquired by high
contrast imaging instruments. In such a self-calibrating system, exoplanet
detection would be limited by photon noise and be significantly more robust
and efficient than in current systems. We show initial laboratory and on-sky
tests, demonstrating over short timescale that residual speckle noise is
indeed calibrated to an accuracy exceeding readout and photon noise in the
high contrast region. We discuss immplications for the design of space and
ground high-contrast imaging systems.

###### keywords:

Adaptive Optics, High Contrast Imaging, Coronagraphy, Exoplanets, Image
Procesing

##  1 Scientific motivation

Direct imaging of exoplanets and circumstellar dust with space and ground
telescopes is challenging due the high flux ration between the central star
and the object of interest (planet or disk), combined with the small angular
separation between the sources. High contrast imaging (HCI) systems are
designed to overcome these challenges by combining optical starlight
suppression with wavefront control. At optical and near-IR wavelengths,
detection with HCI systems is almost always limited by speckle noise: residual
uncontrolled and uncalibrated starlight behind the coronagraph mask imposes a
contrast floor below which detections become too unreliable.

Differential detection techniques can be employed to mitigate speckle noise.
One approach is to use either high spectral reolution signatures in the
exoplanet light [1] or polarized light in starlight by circumstellar dust or
exoplanet atmophseres [2]. The fraction of an exoplanet light containing
polarization of spectral signatures is small, typically up to a few percent,
so these differential imaging techniques, although generally quite reliable,
are inefficient and poorly suited to detect the faintest sources. Another
approach is to separate the speckles from real sources by angular differential
imaging (ADI) [3], which detects the planet image thanks to the known field
rotation, or spectral differential imaging (SDI) [4], which identifies
speckles through their wavelength geometrical scaling. The differential
geometrical effects (rotation or scaling) that ADI and SDI rely on are
proportional to angular separation, so they do not provide significant
contrast gains at the smallest angular separations where some of the most
valuable exoplanets lie. The underlying assumptions of PSF stationarity or
wavevength-scaling are approximations, so ADI and SDI also do not perfectly
remove speckle noise. A third line of research is to use the mutual coherence
between starlight and speckles to identify them, by spatial [5] or temporal
[6] modulation. These promising approaches have recently seen significant
progress, but remain challenging to implement on-sky, as speckles can be mis-
identified as incoherent due to chromaticity and finite temporal resolution.

In this paper, we explore a complementary approach, where wavefront sensing
telemetry could be used to reconstruct the speckle cloud with sufficient
accuracy and precision so that it can be numerically removed down to the
photon noise residual. If successful, this would allow reliable detection of
Earth-like habitable exoplanets orbiting sun-like stars with a 6-m space
telescope without relying on extreme long-term stability of the telescope
optics, or orbiting M-type stars with a ground-based 30-m aperture equiped
with extreme adaptive optics.

To illustrate the scientific potential of high contrast imaging at the photon
noise limit, we consider examples representative of space and ground
telescopes imaging habitable exoplanets:

  * •

An Earth-like planet orbiting a Sun-like star at 10pc distance, observed by a
6-m space telescope

  * •

An Earth-size planet orbiting in the habitable zone of a M4 type star at a 4pc
distance, observed by a 30-m ground-based telescope

Table 1: Observation examples | Space-6m-Earth-G2 | Ground-30m-Earth-M4
---|---|---
Star | G2 at 10pc | M4 at 4pc
Bolometric luminosity [$L_{Sun}$] | 1.000 | 0.0072
Planet orbital radius [au] | 1.0 | 0.085
Maximum angular separation [arcsec] | 0.1 | 0.021
Reflected light planet/star contrast | 1.5e-10 | 2.1e-8
Telescope diameter [m] | 6 | 30
Science spectral bandwidth | 20% | 20%
Central Wavelength | 797 nm (I band) | 1630 nm (H band)
Maximum angular separation [$\lambda$/D] | 4.5 | 1.9
Efficiency | 20 % | 20 %
Total Exposure time | 10 ksec | 10 ksec
Star brightness | $m_{I}=4.04$ | $m_{H}=5.65$
Photon flux in science band (star) | 1.06e9 ph/s | 5.62e9 ph/s
Photon flux in science band (planet) | 0.16 ph/s | 118 ph/s
Background zodi+exozodi [contrast] | 3.1e-10 |
Background starlight [contrast] | 1e-10 | 1e-5
Total background | 4.1e-10 | 1e-5
Background flux in science band | 0.43 ph/s | 56200 ph/s
Photon-noise limited SNR (10 ksec) | 20.8 | 49.7
Exposure time for photon-noise limited SNR=10 | 38 mn | 7 mn
Required PSF calibration accuracy (SNR=10) | 0.15 | 5e-3

Photon-noise limited sensitivity derivations are shown in table 1.

For the space-based observation, the background flux is dominated by the
combined zodiacal and exozodiacal light components, which contribute 2.7
$\times$ more light that the planet. A SNR $=$ 10 detection requires 38 mn
integration. The ground-based telescope benefits from a larger collecting
area, so the planet photon rate is much larger ($>$100 ph/s), but the
background is also significantly higher due to residual atmospheric
turbulence, assumed here to be at the $1e5$ contrast level. A 7 mn integration
allows for SNR=10 detection of the planet.

In order to calibrate speckle noise to 10 $\times$ below the planet flux, the
residual PSF light must be calibrated to 15% accuracy for the space
observation, and to 5e-3 accuracy for the ground observation. We discuss in
this paper steps towards achieving this accuracy, which is especially
challenging for the ground observation example discussed in this section.

##  2 Wavefront Sensor to High Contrast Image Mapping

###  2.1 Self-Calibration Approach

High-contrast imaging (HCI) systems rely on wavefront sensor(s) (WFS) to
measure optical aberrations. The measurements serve as the input to a
wavefront control loop maintaining near-flat wavefront. WFS measurements
usually serve no other purpose. Our calibration approach is to use the WFS
telemetry to reconstruct the point spread function (PSF) with sufficient
accuracy and precision to subtract residual starlight. These additional steps
are shown in Figure 1, colored in red.

Figure 1: WFS-based PSF calibration principle. In a conventional high contrast
imaging system, wavefront sensor (WFS) measurements are used to drive a
wavefront control loop. We discuss in this paper how the WFS telemetry can
also be used to reconstruct the point spread function (PSF) so that speckle
noise can be subtracted in high contrast images (red boxes).

The main challenge to this approach is that it reaquires accurate knowledge of
the relationship between WFS telemetry and PSF. The relationship must be
inferred either by modeling, or learned by analysis of telemetry. We explore
in this paper the latter option.

The WFS-to-PSF relationship is non-linear, so it is not possible to
reconstruct the average PSF from an average of the WFS telemetry. For example,
a rapid zero-mean oscillation of a pupil plane sine wave mode will generate a
pair of focal plane speckles, even through the average wavefront state is free
of aberration. This non-linearity requires the WFS-to-PSF reconstruction to be
performed sufficiently fast to capture wavefront temporal variations. For a
ground-based system, the WFS-to-PSF reconstruction should be performed at
$\approx$ kHz speed, while it can be considerably slower on a space-based
system. Reconstructed intensity PSFs may then summed to temporally match the
science integration and yield a reconstructed PSF that can be subtracted from
the science data. The resulting calibrated science image should then be free
of speckle noise. In this scheme, detection contrat is limited by the photon
noise of the speckle field, as opposed to the speckle field itself.

We note that previous PSF reconstruction efforts have followed a different
approach, where the long exposure PSF is estimated from the temporal variance
of each wavefront mode [7]. This approach does not require each WFS frame to
be processed, but requires an accurate model of the WFS-to-PSF relationship,
and makes several approximations that do not allow for high fidelity PSF
reconstruction for high contrast imaging, namely: WFS linearity, quadratic
relationship between WF and PSF intensity, and statistical independence
between modes. This is not the approach explored in this paper.

###  2.2 Post-processing vs. closed loop control

Since, as shown in Figrue 1, the WFS telemetry serves as the input to the
wavefront control (WFC) loop, we discuss here what is the value of also using
WFS telemetry for PSF reconstruction. We must consider what information is
contained in WFS telemetry that is valuable for PSF reconstruction but cannot
be used for WFC.

There are two fundamental advantages in using WFS telemetry for PSF
reconstruction over WFC :

  * •

Correction Null Space. WFS telemetry contains information about the WF state
that is not correctable by the WFC system. For example, WFS frames may extend
to higher spatial frequencies than correctable by the deformable mirror(s).
WFS frames may also measure amplitude errors. Such errors are in the
correction null space, but affect the PSF.

  * •

Immunity to Latency and Arrow of Time for WF estimation In a WFC loop, each DM
command is based on past WFS measurements. Due to hardware latencies such as
camera exposure time, the current WFS measurement is not accessible (it will
only be available at a later time). Post-processing PSF reconstruction is
immune from such latency.

Figure 2: Use of WFS temetry for control vs. PSF calibration. In wavefront
conrol (top), only past measurements are available to estimate the current WF
state. Due to latency, the current and most recent measurements are not yet
available. Post-processing (bottom) does not suffer from this limitation:
past, current and future measurements are available for PSF estimation.

The arrow of time concept is illustated in Figure 2. In a WFC correction loop
(top), the only WFS measurements available to the control algorithm are in the
past, and the WFS telemetry corresponding to the current (and recent) time is
not yet available. The control command relies on old and noisy estimates,
yielding relatively large WF esimation errors. While predictive control
algorithms may help mitigate this limitation, the correction is fundamentally
blind to recent and current WF changes. In post-processing (bottom), past,
current and future measurements are all available, allowing for a more
accurate WF estimation free of latency. By interpolating/averaging between all
meassurements, WF estimates can also average down measurement noise.

In addition to the null space and arrow of time advantages, post-processing
has access to a large sample set of WFS-to-PSF samples, so an accurate model
of the WFS-to-PSF relationship can be infered from the telemetry.

##  3 On-Sky Experimental Validation of WFS-to-PSF mapping

###  3.1 Solution Uniqueness

The WFS-to-PSF calibration approach requires a unique relationship from WFS
measurement to PSF intensity. If such a relationship exists, and is stable
over time, it can be learned from synchronized samples of WFS and PSF
realizations. The first step in developing a PSF calibration process is
therefore to test is the relationship exists: does WFS telemetry impose a
unique PSF solution ?

Figure 3: Solution uniqueness. In case (A) there is a unique mapping from WFS
to PSF, while in case (B) the mapping is not unique.

Figure 3 illustrates two possible cases, and indicates how statistical tests
can be engineered to test solution uniqueness. WFS measurements are shown on
the left, and corresponding PSF intensity on the right. In case (A), there is
a unique mapping between input and output spaces and WFS-to-PSF calibration is
possible. In case (B) there is no such mapping, and WFS telemetry is
insufficient to constrain the PSF. If a unique mapping exists, clusters of
nearby points in the input space map to clusters on nearby points in the
output space. If no such mapping exists, nearby points in the input space can
map to distant points in the output space.

A statistical test to validate the existence of a WFS-to-PSF calibration
algorithm is to identify nearby WFS realizations and check that the
corresponding output PSFs are similar. We note that this empirical test is not
explicitely relying on a WFS-to-PSF reconstruction algorithm. Yet, a
multidimentional mapping could be constructed by grouping input points in
small clusters and measuring the corresponding output state. Each new
measurement would then be matched to one of the input clusters and the
corresponding solution would be read from a lookup table. In a high dimension
space, this algorithms is not practically feasible and some interpolation is
required.

###  3.2 Experimental Validation of Unique WFS-to-PSF Mapping

The first step toward an on-sky WFS-to-PSF calibration algorithm is to perform
the statistical mapping test described in §3.1. We used the Subaru
Coronagraphic Extreme Adaptive Optics (SCExAO) system [8] on the Subaru
Telescope for this test. The input is a pyramid WFS operating in visible
light, and the output is a focal plane image at 750nm wavelength acquired with
the VAMPIRES instrument [9]. The test was performed in support of the PSF-
sharpening iterative DrWHO algorithm where the AO control loop rewards WFS
frames corresponding to sharp PSF realizations [10, 11].

Figure 4: Solution uniqueness between WFS and PSF: on-sky validation using
SCExAO’s pyramid WFS (top right) and VAMPIRES’s visible light imaging (top
right). Adapted from the DrWHO algorithm on-sky validation [10].

Results are shown in Figure 4. Example WFS and PSF frames are shown at the
top. We first constructed synchronized (WFS, PSF) realizations, and then
compared pairs of such samples. The bottom left plot shows how euclidian
distance between WFS frames compares to euclidian distance between the
corresponding PSF frames. The left side of the diagram corresponds to pairs
with similar WFS realizations, while the bottom part of the diagram shows
similar PSFs. The distribution of points reveals that WFS-similarity implies
PSF-similarity: all points to the left side of the cloud exhibit small PSF
distance. The converse does not hold: PSF-similarity (bottom of the cloud)
does not correspond to WFS-similarity. There are numerous points in the lower
right part of the cloud. This experimental results confirms that there is a
unique mapping from WFS to PSF, but also reveals there is not unique mapping
from PSF to WFS.

This result is consistent with expectations: the WFS provides an unambiguous
measurement of WF state for the modes that are sensed. The focal plane image
does not map unambiguously to WF state: for example, the signature of the
focus WF mode on an otherwise perfect PSF is sign-ambiguous, as the opposite
focus values yield the same image. In this experiment, the DrWHO algorithm
selected high quality PSF realizations for which the focus WF mode was close
to zero. In this second sample (bottom right diagram), WFS similarity is
enforced by PSF similarity.

###  3.3 PSF Reconstruction Validation in Standard Imaging Mode

Figure 5: Prediction of the PSF from the Pyramid WFS data for an on-sky
observation, using a neural network. The predicted PSF image (centre) is
determined entirely from the current WFS image (left), and is seen to closely
match the true PSF measured at that instant (right column). This example shows
a PSF with a large amount of wavefront error (including strong coma) to
provide a clear illustration.

Having verified that a unique mapping exists from WFS telemetry to PSF, we can
attempt to reconstruct the PSF from WFS telemetry. Results of a PSF
reconstruction for on-sky data are shown in Figure 5. Here, a fully-connected
neural network consisting of two 2000-unit layers was trained on 5 minutes of
on-sky data, consisting of synchronised images from the pyramid wavefront
sensor camera and VAMPIRES visible camera[12] (wavelength 750 nm), running at
approximately 500 Hz frame rate. The network used ReLU activation functions
and dropout between each layer as a regularizer, the latter proving to be
crucial for successful reconstruction. While the fully-connected network shown
here provides good results, certain advantages (such as reduced parameter
number and resistance to pupil alignment drift) could be expected from a
convolutional neural network, which is the focus of a current study.

The main features of the PSF speckle cloud are recovered by the
reconstruction. The PSF core shape and bright diffraction features are
accurately reproduced, but fainter speckles are less accurately reconstructed.

##  4 High Contrast Validation

###  4.1 Self-Calibration in a Coronagraphic System

We extend here the WFS-to-PSF calibration to coronagraph systems optimized for
high contrast imaging. When considering HCI systems, we refer to the input WFS
as the bright field (BF) and the output speckle cloud as the drak field (DF).
The self-calibration’s goal is to derive the DH from BF. An ideal system for
self-calibration would have a stable BF-to-DH relationship. This can be
challenging with a configuration such as the one described in §3.3, where the
WFS and PSF use different optical trains and cameras.

We therefore used a configuration where BF and DF are co-located on the same
detector and share the same optical train. The SCExAO near-IR Lyot coronagraph
was configured with speckle control to produce a deep contrast area (DF) over
one side of the focal plane mask, with the opposite side of the focal plane
image remaining relatively bright and serving as the input BF, as shown in
Figure 6 top left image. The experiment is closely related to linear dark
field control (LDFC) [13], where a linear control loop uses the BF as input
for wavefront control. LDFC has been demonstrate to stabilize the DH both in
the laboratory [14] and on-sky [15]. Here, we explore extending LDFC as a DH
calibration algorithm, without assuming linearity.

###  4.2 Laboratory Validation with a Lyot Coronagraph

Results are compiled in Figure 6. The average of all 60000 frames (top left)
shows the BF in the lower half of the image and the DH in the top half. The
intensity variance across the 60000 images is shown for each pixel of the
image at the bottom left, with both readout noise and photon noise variance
terms subtracted to reveal actual speckle intensity variance.

A clustering algorithm is used to derive the mapping between BF and DH, as
detailed in a previous publication [16]. Results are shown in Figure 6,
comparing the full dataset (left column) with one of the clusters serving as
an entry in the BF-to-DH mapping table (center column). The cluster consists
of 128 samples with similar BF realizations: the variance within the BF is
35.7 $\times$ smaller within the cluster than across the full dataset. The
corresponding measured DH variance is 30.7 $\times$ smaller within the cluster
set than across the full input dataset. This last result demonstrates that
images with similar BFs also have similar DHs. Image selection using BF
intensity successfully constrains DH intensity, demonstrating that a BF-to-DH
calibration algorithm can be derived to calibrate residual speckles in high
contrast images.

Figure 6: Self-calibration of high contrast images in a half dark hole
configuration on the SCExAO bench. A near-IR Lyot type coronagraph is used to
block starlight. Focal plane wavefront control creates a deep contrast area
above the coronagraph mask, and a relatively bright area below it. The DH is
calibrated from the input sensing area (BF) below the optical axis. A
clustering algorithm is employed to derive BF-to-DH mapping.

###  4.3 Null Calibration with Photonic Nulling Interferometer: On-Sky
Validation

A photonic nuller is an alternative solution to the high contrast imaging
challenge. Unlike a coronagraph constructed from bulk optics between which
light freely propagates, the photonic nuller couples starlight into a small
number of coherent singlemode waveguides. The waveguides are coherently
combined to produce starlight destructive interference in null output(s).
Bright starlight is directed to bright outputs which measure the intensity in
input waveguides (photometry output(s)) and phase offset between input
waveguides (WFS output(s)). The photonic nuller concept and its implementation
are discussed in publications from the GLINT instrument team[17, 18].

The photonic nuller is well-suited for self calibration :

  * •

Starlight is coupled in a small number of coherent waveguides. At each
wavelength, light into the photonic device input is fully described by phase
and amplitude (and possibly polarization), so the number of dimension in the
input is a few times the number of waveguides

  * •

The relationship between input variables (phase and amplitude of each
waveguide) and output intensities is entirely established within the photonic
chip so it is significantly more stable than an optical train of optical
components subject to relative misalignments.

Figure 7: Null calibration with the GLINT photonic nuller: Laboratory (left)
and on-sky (right) demonstration. Each of the 8 vertical panel shows the GLINT
detector ouput, consisting of 16 horizontal spectra (one per output of the
photonic chip) ranging from 1340 nm (right edge) to 1690 nm (left edge).
Images are grouped by pair, with the average image on the left and standard
deviation on the right. The leftmost pair is the average (a) and standard
deviation (b) of the full set of 12,796 consecutive images. The average (c)
and standard deviation (d) of one of the clusters that define the BF-to-DH is
shown. Blue boxes indicate the BF used for selection, and the red boxes show
the DH signals. The right side of the figure shows the same data for on-sky
observations.

We validated the BF-to-DH calibration approach on the GLINT instrument[18]
installed on the Subaru Telescope. Figure 7 shows laboratory (left) and on-sky
(right) results. Wavefront errors were added with the system’s deformable
mirror in the lab experiment. All datasets were aquired at 1.4 kHz frame rate.

Both experiments demonstrate that the DH signal is well constrained by the BF
state. In the on-sky result, the standard deviation of the DH signal (d) after
calibration is below the readout and photon noise of the DH signal, indicating
that the self-calibration is able to estimate the residual starlight to an
accuracy below photon and readout noise. Details of the experiments are
provided in a previous publication [16]. The experiment validates uniqueness
of the null solution for a given BF measurement, with no evidence of a
measurement null space which would induce a variation in the null ouputs
without a corresponding signature in the BF.

##  5 Conclusion and Perspectives

Self-calibration of high contrast imaging data using WFS telemetry is a
promising solution to the current speckle noise limitation, potentially
providing photon-noise limited detection limits on space and ground high
contrast imaging systems. Preliminary on-sky experiments described in this
paper are encouraging, and demonstrate that speckle noise is well constrained
by WFS telemetry. The experiments were however limited in duration and did not
validate that the WFS-to-PSF relationship is stable over long periods of time,
which is critical to subtract speckle noise without removing planet signal. If
the calibration holds, then the WFS-to-PSF mapping can be recorded on a
calibration star (without a planet) and used to remove speckle noise on the
science target. This demanding stability requirement will most likely be
achieved by minimizing the time-variable non-common path aberrations (NCPA)
between the input WFS and output science image.

Our findings suggest that future space-based high contrast imaging systems
could be designed with less demanding telescope wavefront stability, as
speckle noise can be accurately calibrated and removed from the science
images. In this regime, speckle noise at or below the natural zodi+exozodi
background at $\approx$ 3e9 contrast will not affect detetion limits, and
stronger residual starlight can be compensated for by longer integration time
to average photon noise. The coronagraph throughput and bandpass would be
optimized for sensitivity in the photon-noise regime instead of reaching
contrast levels below the exoplanet level.

For both space and ground systems, accurate self-calibration will require a
stable or well-calibrated relationship between WFS and science data. The
photonic nulling chip approach appears particularlty promising, as (1) the WFS
and nulling functions are imprinted in the same physically small device,
leaving little room for the WFS-to-DF relationship to change, (2) the WF state
is discretized in a small number of coherent waveguides and (3) the output
channels can easily be dispersed in wavelength.

We have relied on an explicit BF-to-DH mapping in this study. Our experiments
built the mapping by clustering input samples and averaging the corresponding
output samples. This allows for inspection of the output samples corresponding
to the samples of the input cluster, so that the mapping uniqueness can be
verified and quantified. This ”lookup table” approach does not scale
efficiently to high dimensions, for which more efficient reconstruction
algorithms are needed. Neural network approaches have recently been shown to
address this challenge [19, 20], and perform well on the same photonic nulling
data as presented in §4.3. Validation over larger datasets and longer time
spans will be required to confirm algorithms are suitable for the high
contrast imaging application discussed in this paper.

###### Acknowledgements.

This work was supported by NASA grants #80NSSC19K0336 and #80NSSC19K0121. This
work is based on data collected at Subaru Telescope, which is operated by the
National Astronomical Observatory of Japan. The authors wish to recognize and
acknowledge the very significant cultural role and reverence that the summit
of Maunakea has always had within the Hawaiian community. We are most
fortunate to have the opportunity to conduct observations from this mountain.
The authors also wish to acknowledge the critical importance of the current
and recent Subaru Observatory daycrew, technicians, telescope operators,
computer support, and office staff employees. Their expertise, ingenuity, and
dedication is indispensable to the continued successful operation of these
observatories. The development of SCExAO was supported by the Japan Society
for the Promotion of Science (Grant-in-Aid for Research #23340051, #26220704,
#23103002, #19H00703 & #19H00695), the Astrobiology Center of the National
Institutes of Natural Sciences, Japan, the Mt Cuba Foundation and the
director’s contingency fund at Subaru Telescope. KA acknowledges support from
the Heising-Simons foundation.

## References

  * [1] Wang, J., Mawet, D., Ruane, G., Hu, R., and Benneke, B., “Observing exoplanets with high dispersion coronagraphy. i. the scientific potential of current and next-generation large ground and space telescopes,” The Astronomical Journal 153, 183 (mar 2017).
  * [2] Hunziker, S., Schmid, H. M., Mouillet, D., Milli, J., Zurlo, A., Delorme, P., Abe, L., Avenhaus, H., Baruffolo, A., Bazzon, A., Boccaletti, A., Baudoz, P., Beuzit, J. L., Carbillet, M., Chauvin, G., Claudi, R., Costille, A., Daban, J. B., Desidera, S., Dohlen, K., Dominik, C., Downing, M., Engler, N., Feldt, M., Fusco, T., Ginski, C., Gisler, D., Girard, J. H., Gratton, R., Henning, T., Hubin, N., Kasper, M., Keller, C. U., Langlois, M., Lagadec, E., Martinez, P., Maire, A. L., Menard, F., Meyer, M. R., Pavlov, A., Pragt, J., Puget, P., Quanz, S. P., Rickman, E., Roelfsema, R., Salasnich, B., Sauvage, J. F., Siebenmorgen, R., Sissa, E., Snik, F., Suarez, M., Szulágyi, J., Thalmann, C., Turatto, M., Udry, S., van Holstein, R. G., Vigan, A., and Wildi, F., “RefPlanets: Search for reflected light from extrasolar planets with SPHERE/ZIMPOL,” A&A 634, A69 (Feb. 2020).
  * [3] Marois, C., Lafrenière, D., Doyon, R., Macintosh, B., and Nadeau, D., “Angular Differential Imaging: A Powerful High-Contrast Imaging Technique,” ApJ 641, 556–564 (Apr. 2006).
  * [4] Racine, R., Walker, G. A. H., Nadeau, D., Doyon, R., and Marois, C., “Speckle Noise and the Detection of Faint Companions,” PASP 111, 587–594 (May 1999).
  * [5] Baudoz, P., Boccaletti, A., Baudrand, J., and Rouan, D., “The Self-Coherent Camera: a new tool for planet detection,” in [IAU Colloq. 200: Direct Imaging of Exoplanets: Science & Techniques ], Aime, C. and Vakili, F., eds., 553–558 (Jan. 2006).
  * [6] Guyon, O., “Imaging Faint Sources within a Speckle Halo with Synchronous Interferometric Speckle Subtraction,” ApJ 615, 562–572 (Nov. 2004).
  * [7] Gendron, E., Clénet, Y., Fusco, T., and Rousset, G., “New algorithms for adaptive optics point-spread function reconstruction,” A&A 457, 359–363 (Oct. 2006).
  * [8] Lozi, J., Guyon, O., Jovanovic, N., Goebel, S., Pathak, P., Skaf, N., Sahoo, A., Norris, B., Martinache, F., N’Diaye, M., Mazin, B., Walter, A. B., Tuthill, P., Kudo, T., Kawahara, H., Kotani, T., Ireland, M., Cvetojevic, N., Huby, E., Lacour, S., Vievard, S., Groff, T. D., Chilcote, J. K., Kasdin, J., Knight, J., Snik, F., Doelman, D., Minowa, Y., Clergeon, C., Takato, N., Tamura, M., Currie, T., Takami, H., and Hayashi, M., “SCExAO, an instrument with a dual purpose: perform cutting-edge science and develop new technologies,” in [Adaptive Optics Systems VI ], Proc. SPIE 10703, 1070359 (July 2018).
  * [9] Norris, B. R. M., Tuthill, P., Jovanovic, N., Lozi, J., Guyon, O., Cvetojevic, N., and Martinache, F., “Diffraction-limited polarimetric imaging of protoplanetary disks and mass-loss shells with VAMPIRES,” Proc. SPIE 11203, 112030S (Jan. 2020).
  * [10] Skaf, N., Guyon, O., Gendron, É., Ahn, K., Bertrou-Cantou, A., Boccaletti, A., Cranney, J., Currie, T., Deo, V., Edwards, B., Ferreira, F., Gratadour, D., Lozi, J., Norris, B., Sevin, A., Vidal, F., and Vievard, S., “On-sky validation of image-based adaptive optics wavefront sensor referencing,” A&A 659, A170 (Mar. 2022).
  * [11] Skaf, N., Guyon, O., Boccaletti, A., Deo, V., Lozi, J., Vievard, S., Currie, T., Norris, B., Ahn, K., Gendron, E., Vidal, F., Ferreira, F., Sevin, A., and Bertrou-Cantou, A., “High contrast imaging wavefront sensor referencing from coronagraphic images,” Proc. SPIE 11823, 387 – 401, SPIE (2021).
  * [12] Norris, B., Schworer, G., Tuthill, P., Jovanovic, N., Guyon, O., Stewart, P., and Martinache, F., “The VAMPIRES instrument: imaging the innermost regions of protoplanetary discs with polarimetric interferometry,” Monthly Notices of the RAS 447, 2894–2906 (Mar. 2015).
  * [13] Miller, K. L., Bos, S. P., Lozi, J., Guyon, O., Doelman, D. S., Vievard, S., Sahoo, A., Deo, V., Jovanovic, N., Martinache, F., Snik, F., and Currie, T., “Spatial linear dark field control on Subaru/SCExAO. Maintaining high contrast with a vAPP coronagraph,” Astronomy and Astrophysics 646, A145 (Feb. 2021).
  * [14] Currie, T., Pluzhnik, E., Guyon, O., Belikov, R., Miller, K., Bos, S., Males, J., Sirbu, D., Bond, C., Frazin, R., Groff, T., Kern, B., Lozi, J., Mazin, B. A., Nemati, B., Norris, B., Subedi, H., and Will, S., “Laboratory Demonstration of Spatial Linear Dark Field Control For Imaging Extrasolar Planets in Reflected Light,” PASP 132, 104502 (Oct. 2020).
  * [15] Bos, S. P., Miller, K. L., Lozi, J., Guyon, O., Doelman, D. S., Vievard, S., Sahoo, A., Deo, V., Jovanovic, N., Martinache, F., Currie, T., and Snik, F., “First on-sky demonstration of spatial Linear Dark Field Control with the vector-Apodizing Phase Plate at Subaru/SCExAO,” A&A 653, A42 (Sept. 2021).
  * [16] Guyon, O., Norris, B., Martinod, M.-A., Ahn, K., Tuthill, P., Males, J., Wong, A., Skaf, N., Currie, T., Miller, K., Bos, S., Lozi, J., Deo, V., Vievard, S., Belikov, R., van Gorkom, K., Haffert, S., Mazin, B., Bottom, M., Frazin, R., Rodack, A., Groff, T., Jovanovic, N., and Martinache, F., “High contrast imaging at the photon noise limit with self-calibrating WFS/C systems,” Proc. SPIE 11823, 1182318 (Sept. 2021).
  * [17] Norris, B. R. M., Cvetojevic, N., Lagadec, T., Jovanovic, N., Gross, S., Arriola, A., Gretzinger, T., Martinod, M.-A., Guyon, O., Lozi, J., Withford, M. J., Lawrence, J. S., and Tuthill, P., “First on-sky demonstration of an integrated-photonic nulling interferometer: the GLINT instrument,” Monthly Notices of the RAS 491, 4180–4193 (Jan. 2020).
  * [18] Martinod, M.-A., Norris, B., Tuthill, P., Lagadec, T., Jovanovic, N., Cvetojevic, N., Gross, S., Arriola, A., Gretzinger, T., Withford, M. J., Guyon, O., Lozi, J., Vievard, S., Deo, V., Lawrence, J. S., and Leon-Saval, S., “Scalable photonic-based nulling interferometry with the dispersed multi-baseline GLINT instrument,” Nature Communications 12, 2465 (Jan. 2021).
  * [19] Norris, B., Wei, J., Betters, C., Leon-Saval, S., Xin, Y., Lin, J., Kim, Y. J., Sallum, S., Lozi, J., Vievard, S., Guyon, O., Gatkine, P., Jovanovic, N., Mawet, D., and Fitzgerald, M., “Demonstration of a photonic-lantern focal-plane wavefront sensor using fiber mode conversion and deep learning,” Proc. SPIE 12188 (2022).
  * [20] Norris, B., Martinod, M.-A., Tuthill, P., Gross, S., Cvetojevic, N., Jovanovic, N., Lagadec, T., Klinner-teo, T., Guyon, O., Lozi, J., Deo, V., Vievard, S., Arriola, A., Gretzinger, T., Lawrence, J., and Withford, M., “Optimal self-calibration and fringe tracking in photonic nulling interferometers using machine learning,” Proc. SPIE 12188 (2022).


# Robustness of optimal probabilistic storage and retrieval of unitary
channels to noise

Jaroslav Pavličko and Mário Ziman  RCQI, Institute of Physics, Slovak Academy
of Sciences, Dúbravská cesta 9, 84511 Bratislava, Slovakia

###### Abstract

We investigate robustness of probabilistic storage and retrieval device
optimized for phase gates to noise. We use noisy input composed of convex
combination of unitary channel with either depolarizing or dephasing channel.
We find out that the resistance to dephasing noise is higher than to
depolarization. Interestingly, for the depolarisation the retrieval reduces
the degree of noise. We also examine the possible realizations showing that
their performance is different when the noise is present.

###### pacs:

03.67.-a, 03.67.Ac

##  I Introduction

Quantum information theory and quantum computers themselves, started to draw
more attention since the discovery of the first algorithms with potential
quantum speed-ups such as factorization or quantum search shor ; grover , or
seminal work on quantum cryptography QKD . The theory itself is bounded by no-
go theorems noclon ; nobroadcast ; nohiding ; nodelete  \- theorems that put
limitations on what the quantum processing of information can achieve. There
exist two different conceptual ways how to circumvent the no-go theorems:
either give up the exactness, or the determinism of the performance. In other
words one is either satisfied only with approximate solutions or one keeps the
requirement of exact solutions, with the caveat of only a certain probability
of successfully reaching them.

At the heart of a classical computer science lies the concept of processor.
Therefore, one envisions a similar device also for a quantum computer - a so-
called quantum processor nielsen1 . However, due to the no-programming theorem
nielsen1 , one can never construct a truly universal quantum processor capable
of executing all transformations with absolute precision. One way of dealing
with the no-programming theorem is to define processors that are able to
implement the exact transformation however only with some probability
probProc1 ; probProc2  (the other way is to use approximate processors
appProc3 ; appProc4 ).

In a sense complementary problem to quantum programming is known as quantum
learning. The aim is to design protocols to store and retrieve the action of
an unknown quantum transformations. Informally, the goal of the storage phase
is to exploit the transformation device given number of times $N$ and imprint
its transformation into a quantum state of some memory system in a way that
during the retrieval phase this state enables us (using the quantum processor)
to implement the stored transformation given number of times $M$. There are
several variations of this problem and few of them has been already considered
quantLearning ; storingQuantDyn ; STORunitary ; STORPhaseGate . In this work,
we investigate the noise robustness of optimal probabilistic storage and
retrieval device (PSAR) for qubit phase gates introduced in STORPhaseGate .

Following STORPhaseGate  consider a one-parametric set of unitary gates
$U_{\varphi}=|0\rangle\langle 0|+e^{i\varphi}|1\rangle\langle 1|$ acting on a
qubit Hilbert space $\mathcal{H}_{2}$. We assume that during the storing phase
(see Figure 1)we can access the unknown unitary gate $N$ times and produce a
state $\omega_{U}$ (memory). This state is used in the retrieval phase to
restore the action of the gate with a probability of success $p_{\rm
success}=N/(N+1)$ independent of $U_{\varphi}$.

In this work we aim to analyze the performance of the PSAR devices introduced
in STORPhaseGate  in cases when during the storing phase the ideal unitary
gates are replaced by their noisy versions. In particular, we consider two
different noise models: depolarizing and phase-damping ones. The retrieved
channels will be compared with the stored noisy gates, but also with the
noiseless unitary gates. The paper is organized as follows: in Section II we
describe the mathematical framework of quantum networks used to formulate and
solve the problem. Section III provides mathematical formulation of the
problem and introduces the noise models we investigate. In Section IV we
discuss in details the case of $2\to 1$ PSAR and Section V contains
generalizations for $N\to 1$ case. The implementations are studied in Section
VI and Section VII summarized the conclusions of the performed analysis.

Figure 1: Schematic image of device PSAR optimized for implementation of
channel ${\cal{U}}_{\varphi}$ for a general case when we have access to
unitary channel $N$ times with the input state $|\psi\rangle$ from equation
(2). At the output of register ${\cal{H}}_{D}$, we expect to retrieve the
unitary channel in case of successful implementation. The state $|\xi\rangle$
is the one we desire.

##  II Quantum Networks

Mathematical formalism used in this paper stems from the framework of quantum
networks QN ; higherOrderQT ; architecture ; QN2 . Let us denote by
${\cal{L}}(\mathcal{H}_{a})\equiv{\cal{L}}_{a}$ the set of linear operators on
a finite-dimensional Hilbert space $\mathcal{H}_{a}$ and by
${\cal{L}}(\mathcal{H}_{a},\mathcal{H}_{b})\equiv{\cal{L}}_{ab}$ the set of
linear operators from $\mathcal{H}_{a}$ to $\mathcal{H}_{b}$. Quantum
operation ${\cal{O}}:{\cal{L}}_{a}\to{\cal{L}}_{b}$ is a completely positive
trace non-increasing linear map. Trace-preserving quantum operation, thus
mapping quantum states into quantum states, is called quantum channel.

There exists a one-to-one correspondence (Choi-Jamiolkowski isomorphism)
between linear operators $A\in{\cal{L}}(\mathcal{H}_{a},\mathcal{H}_{b})$ and
vectors $|A\rangle\\!\rangle\in\mathcal{H}_{a}\otimes\mathcal{H}_{b}$:

$\displaystyle A=\sum A_{mn}|m\rangle\langle
n|\quad\leftrightarrow\quad|A\rangle\\!\rangle=\sum_{mn}A_{mn}|m\rangle|n\rangle$

where $A_{mn}=\langle m|A|n\rangle$, $\\{|m\rangle\\}$ and $\\{|n\rangle\\}$
are orthonormal basis of $\mathcal{H}_{a}$ and $\mathcal{H}_{b}$,
respectively. Similarly, every quantum operation ${\cal{O}}$ is associated
with a positive operator (Choi operator)

$\displaystyle
O_{ba}=({\cal{O}}\otimes{\cal{I}}_{a})({|I\rangle\\!\rangle\langle\\!\langle
I|}),$

where ${\cal{I}}_{a}$ is an identity map on ${\cal{L}}(\mathcal{H}_{a})$ and
$I\in{\cal{L}}_{aa}$ is the identity operator, i.e.
$|I\rangle\\!\rangle=\sum_{n}|n\rangle\otimes|n\rangle$.

Consider a composition of two operations
${\cal{O}}:{\cal{L}}_{a}\to{\cal{L}}_{b})$ and
$\widetilde{\cal{O}}:{\cal{L}}_{b}\to{\cal{L}}_{c}$. Then, the Choi operator
of the composition $\widetilde{\cal{O}}\circ\cal{O}$ is expressed via link
product denoted by $\star$ and defined as follows:

$\displaystyle[O\star\widetilde{O}]_{ac}=\operatorname{tr}_{b}\left[(\mathbb{I}_{c}\otimes
O_{ab}^{T_{b}})(\mathbb{I}_{a}\otimes\widetilde{O}_{bc})\right],$

where $T_{b}$ denotes partial transposition on space $\mathcal{H}_{b}$ and
$\mathbb{I}_{a}$ stands for the identity operator on $\mathcal{H}_{a}$.

Figure 2: Depiction of a quantum network formed by a concatenation of $N$
quantum operations $\\{O_{i}\\},i=1,\cdots,N$.

Quantum network ${\cal{R}}$ (also known as quantum comb) is a concatenation of
quantum operations ${\cal{O}}_{1},\cdots,{\cal{O}}_{N}$ (see Fig. (2)), where
some outputs of the preceding operation are connected with some inputs of the
subsequent one. The connectivity structure is included in the definition of
each particular link product in the whole sequence of compositions. As a
result after linking operations some of the outputs and inputs remain open,
thus, allowing the quantum network to accept quantum channels and operations
of a suitable input-output type at its inputs. Overall the quantum network
${\cal{R}}$ transforms $N-1$ operations into one quantum operation. The
associated Choi operator of quantum network ${\cal{R}}$ is an operator
$R^{N}=O_{1}\star O_{2}\star\cdots\star O_{N}$ acting on the Hilbert space
${\cal{H}_{\rm in}}\otimes{\cal{H}_{\rm out}}$, where ${\cal{H}}_{\rm
in}=\otimes_{i=0}^{N-1}{\cal{H}}_{2i}$ represents the causally ordered input
spaces and ${\cal{H}}_{\rm out}=\otimes_{i=0}^{N-1}{\cal{H}}_{2i+1}$ stands
for the output spaces. Choi operators of deterministic (trace-preserving)
quantum networks obey the recursive normalization conditions:

$\displaystyle\operatorname{tr}_{2k-1}[R^{k}]=\mathbb{I}_{2k-2}\otimes
R^{k-1},$

for $k=1,\cdots,N$. For probabilistic quantum network ${\cal{S}}$ there always
exists a deterministic quantum network ${\cal{R}}$ such that $S\leq R$.
Collection of probabilistic quantum networks
${\cal{S}}_{1},\cdots,{\cal{S}}_{m}$ summing up to a deterministic quantum
network $\sum_{x}{\cal{S}}_{x}={\cal{R}}$ form a generalized quantum
instrument. Generalized quantum instrument fulfills a similar role for quantum
networks as does quantum instrument for quantum channels. The index $x$
represents the corresponding classical outcome.

##  III PSAR with noise.

The optimal probabilistic storage and retrieval (PSAR) for phase gates
introduced in STORPhaseGate  consists of the probe state
$|\Psi\rangle_{AA^{\prime}}$, where $A$ labels all inputs of $N$ phase gates
and $A^{\prime}$ represent the ancillary part of the memory, and the retrieval
operation ${\cal R}_{s}$ using the memory $M$ to implement the phase gate
action on the unknown input system $C$ that is transformed into the system $D$
(see Fig. 1).

In particular, during the storing this state is transformed into

$|\Psi_{\varphi}\rangle=U_{\varphi}\otimes\cdots\otimes U_{\varphi}\otimes
I_{A^{\prime}}|\Psi\rangle_{AA^{\prime}}\,,$ (1)

where $U_{\varphi}=|0\rangle\langle 0|+e^{i\varphi}|1\rangle\langle 1|$ is the
phase gate. It was derived in STORPhaseGate  that optimal probe state takes
the form

$|\Psi\rangle_{AA^{\prime}}=\frac{1}{\sqrt{N+1}}\bigoplus_{j}|\Pi_{j}\rangle\\!\rangle_{AA^{\prime}}\,,$
(2)

where $\Pi_{j}$ are projectors onto Hilbert subspaces ${\cal H}_{j}$. These
subpsaces follows from the irreducible representations of $U(1)$

$U_{\varphi}^{\otimes N}=\bigoplus_{j=0}^{N}e^{ij\varphi}\otimes I_{m_{j}}\,,$
(3)

where $I_{m_{j}}$ stands for identity operator in the multiplicity subspace
${\cal H}_{m_{j}}$, thus, inducing the decomposition ${\cal
H}_{A}=\sum_{j}{\cal H}_{j}\otimes{\cal H}_{m_{j}}$ with one-dimensional
${\cal H}_{j}$. Consequently,

$|\Psi_{\varphi}\rangle=\frac{1}{\sqrt{N+1}}\bigoplus_{j}e^{ij\varphi}|\Pi_{j}\rangle\\!\rangle\,.$
(4)

We may set $\Pi_{j}=|\overline{j}\rangle\langle\overline{j}|$ with
$|\overline{j}\rangle=|0^{\otimes(N-j)}\rangle\otimes|1^{\otimes
j}\rangle\in{\cal H}_{2}^{\otimes N}$. Then
$|\Psi\rangle_{AA^{\prime}}=(1/\sqrt{N+1})\sum_{j}|\overline{j}\overline{j}\rangle$.

Let us denote the memory state as
$\Psi_{\varphi}=|\Psi_{\varphi}\rangle\langle\Psi_{\varphi}|$. It was shown in
STORPhaseGate  that using the Choi operator

$R_{s}=\bigoplus_{J=0}^{N-1}\sum_{k,k^{\prime}=0}^{1}|J+k,J+k\rangle_{M}\langle
J+k^{\prime},J+k^{\prime}|\otimes|kk\rangle_{CD}\langle k^{\prime}k^{\prime}|$

associated with the retrieval operation ${\cal R}_{s}$ the action of the
retrieval is evaluated as follows

$\displaystyle R_{s}\star\Psi_{\varphi}$ $\displaystyle=$ $\displaystyle{\rm
tr}_{M}[R_{s}(\Psi_{\varphi}^{T}\otimes
I_{CD})]=\langle\Psi_{\varphi}^{*}|R_{s}|\Psi_{\varphi}^{*}\rangle$
$\displaystyle=$
$\displaystyle\frac{N}{N+1}|U_{\varphi}\rangle\\!\rangle\langle\\!\langle
U_{\varphi}|\,,$

where
$|\Psi_{\varphi}^{*}\rangle=(1/\sqrt{N+1})\oplus_{j}e^{-ij\varphi}|\Pi_{j}\rangle\\!\rangle$
and $N/(N+1)$ is the success probability, because

$R_{s}\star\Psi_{\varphi}\star|\xi\rangle\langle\xi|=\frac{N}{N+1}U_{\varphi}|\xi\rangle\langle\xi|U_{\varphi}^{\dagger}\,.$
(5)

There are several variations how the noise can enter the design of the storing
and retrieval procedures. In what follows we will assume the performance of
the black box PSAR introduced in STORPhaseGate  is unchanged and only the
phase gates we are aiming to store are noisy. That is, instead of ${\cal
U}_{\varphi}(\cdot)=U_{\varphi}\cdot U_{\varphi}^{\dagger}$ the ”phase gates”
implement the following transformation ${\cal E}_{\varphi}=q{\cal
U}_{\varphi}+(1-q){\cal N}$, where ${\cal N}$ is the noise. We will consider
two types of noise:

  * •

Depolarisation. In this case we set ${\cal N}={\cal C}_{I/2}$, where ${\cal
C}_{I/2}$ stands for the completely depolarizing noise (also known as white
noise) transforming any quantum state $\varrho$ into the complete mixture,
i.e. ${\cal{C}}_{I/2}(\varrho)=I/2$. The phase gate implements the channel

$\displaystyle{\cal{E}}_{\varphi}=q{\cal{U}}_{\varphi}+(1-q){\cal{C}}_{I/2},$

being a convex combination ($q\in[0,1]$) of the desired unitary channel and a
completely depolarizing noise ${\cal C}_{I/2}$ with the Choi operator
$C=\frac{1}{2}(|00\rangle\langle 00|+|01\rangle\langle 01|+|10\rangle\langle
10|+|11\rangle\langle 11|)=(I\otimes I)/2$.

  * •

Dephasing. In this case we set ${\cal N}={\cal P}$, where ${\cal P}$ stands
for the completely dephasing channel for all states $\varrho$ diminishing off-
diagonal terms in the computational basis identified with eigenvectors of
$\sigma_{z}$ operator, i.e. ${\cal
P}(\varrho)=\frac{1}{2}(\varrho+\sigma_{z}\varrho\sigma_{z})={\rm
diag}[\varrho]=\langle 0|\varrho|0\rangle|0\rangle\langle 0|+\langle
1|\varrho|1\rangle|1\rangle\langle 1|$. That is, phase gates are implementing
convex combination

$\displaystyle{\cal{F}}_{\varphi}=q{\cal{U}}_{\varphi}+(1-q){\cal{P}},$

where $q\in[0,1]$ and the Choi operator of the completely dephasing channel
reads $P=|00\rangle\langle 00|+|11\rangle\langle 11|$.

##  IV Case-study: $2\rightarrow 1$ PSAR with noise.

In this section we will investigate in details the situation when the unknown
noisy phase gate device is used twice in the storing phase. The situation is
depicted in Fig. 3. For brevity we introduce the following notation for the
involved Hilbert spaces ${\cal{H}}_{13}={\cal{H}}_{A}$,
${\cal{H}}_{24}={\cal{H}}_{B}$, ${\cal{H}}_{0}={\cal{H}}_{C}$ and
${\cal{H}}_{5}={\cal{H}}_{D}$. Choi operator for two uses of the considered
channel is given as follows:

$\displaystyle E_{\varphi,AB}^{\otimes 2}=q^{2}|U_{\varphi}^{\otimes
2}\rangle\\!\rangle_{AB}\langle\\!\langle U^{\otimes
2}_{\varphi}|+(1-q^{2})N_{12}\otimes N_{34}$
$\displaystyle+q(1-q)\left(|U_{\varphi}\rangle\\!\rangle_{12}\langle\\!\langle
U_{\varphi}|\otimes
N_{34}+N_{12}\otimes|U_{\varphi}\rangle\\!\rangle_{34}\langle\\!\langle
U_{\varphi}|\right).$

Unitary phase gate part can be written as follows

$\displaystyle U_{\varphi}^{\otimes 2}$ $\displaystyle=|00\rangle\langle
00|+e^{i\varphi}(|01\rangle\langle 01|+|10\rangle\langle
10|)+e^{i2\varphi}|11\rangle\langle 11|$
$\displaystyle=\bigoplus_{j=0}^{2}e^{ij\varphi}\otimes
I_{m_{j}}=\sum_{k=0}^{2}e^{ik\varphi}|\overline{k}\rangle\langle\overline{k}|+e^{i\varphi}|\overline{3}\rangle\langle\overline{3}|$
$\displaystyle\,,$

where we introduced the vectors $|\overline{0}\rangle\equiv|00\rangle$,
$|\overline{1}\rangle\equiv|01\rangle$,
$|\overline{2}\rangle\equiv|11\rangle$, $|\overline{3}\rangle\equiv|10\rangle$
and $I_{m_{j}}$ denotes the identity operator on multiplicity spaces. Using
this notation the probe state $|\Psi\rangle$ reads

$\displaystyle\begin{split}|\Psi\rangle_{AA^{\prime}}&=\bigoplus_{j=0}^{2}\frac{1}{\sqrt{3}}|\Pi_{j}\rangle\\!\rangle_{AA^{\prime}}=\frac{1}{\sqrt{3}}(|\overline{0}\overline{0}\rangle+|\overline{1}\overline{1}\rangle+|\overline{2}\overline{2}\rangle)_{AA^{\prime}}\,.\end{split}$

Figure 3: Schematic image of device PSAR implementing channel
${\cal{E}}_{\varphi}$ twice with the input state $|\psi\rangle$ from equation
(2). At the output of register ${\cal{H}}_{5}$, we expect to retrieve the
unitary channel, possibly with some noise, in case of successful
implementation.

###  IV.1 Depolarization noise.

For the depolarizing noisy phase gates the storing results in the state

$\displaystyle\Psi_{\varphi}$ $\displaystyle=E_{\varphi,AB}^{\otimes
2}\star|\Psi\rangle_{AA^{\prime}}\langle\Psi|$ (6)
$\displaystyle=\operatorname{tr}_{A}[(E_{\varphi,AB}^{\otimes 2}\otimes
I_{A^{\prime}})(|\Psi\rangle_{AA^{\prime}}\langle\Psi|^{T_{A}}\otimes I_{B})]$
$\displaystyle=q^{2}\varrho^{U,U}+q(1-q)\left(\varrho^{C,U}+\varrho^{U,C}\right)+(1-q)^{2}\varrho^{C,C}\,,$

where

$\displaystyle\varrho^{U,U}$ $\displaystyle=$
$\displaystyle\operatorname{tr}_{A}\bigg{[}\left(|U_{\varphi}^{\otimes
2}\rangle\\!\rangle\langle\\!\langle U_{\varphi}^{\otimes 2}|_{AB}\otimes
I_{A^{\prime}}\right)\left(|\Psi\rangle\langle\Psi|_{AA^{\prime}}^{T_{A}}\otimes
I_{B}\right)\bigg{]}$ (7) $\displaystyle=$
$\displaystyle\bigoplus_{j,k=0}^{2}\frac{1}{3}e^{i(j-k)\varphi}|\overline{j}\overline{j}\rangle\langle\overline{k}\overline{k}|_{BA^{\prime}}\,,$
$\displaystyle\varrho^{C,C}$ $\displaystyle=$
$\displaystyle\operatorname{tr}_{A}\Bigg{[}\left(\frac{1}{2}I_{12}\otimes\frac{1}{2}I_{34}\otimes
I_{A^{\prime}}\right)(|\Psi\rangle_{AA^{\prime}}\langle\Psi|^{T_{A}}\otimes
I_{B})\Bigg{]}$ (8) $\displaystyle=$
$\displaystyle\frac{1}{4}\times\frac{1}{3}\left(I_{B}\otimes\sum_{k=0}^{2}|\overline{k}\rangle\langle\overline{k}|_{A^{\prime}}\right),$
$\displaystyle\varrho^{U,C}$ $\displaystyle=$
$\displaystyle\frac{1}{6}(\Pi_{\overline{00}}+\Pi_{\overline{11}}+\Pi_{\overline{22}}+\Pi_{\overline{01}}+\Pi_{\overline{10}}+\Pi_{\overline{32}})+$
(9)
$\displaystyle+\frac{1}{6}e^{i\varphi}(|\overline{32}\rangle\langle\overline{01}|+|\overline{22}\rangle\langle\overline{11}|)+c.c.$
$\displaystyle\varrho^{C,U}$ $\displaystyle=$
$\displaystyle\frac{1}{6}(\Pi_{\overline{00}}+\Pi_{\overline{11}}+\Pi_{\overline{22}}+\Pi_{\overline{12}}+\Pi_{\overline{21}}+\Pi_{\overline{30}})+$
(10)
$\displaystyle+\frac{1}{6}e^{i\varphi}(|\overline{11}\rangle\langle\overline{00}|+|\overline{21}\rangle\langle\overline{30}|)+c.c.\,,$

where we used the notation
$\Pi_{\overline{jj}}=|\overline{jj}\rangle\langle\overline{jj}|$. The
evaluation of the retrieving instrument

$R_{s}=\bigoplus_{J=0}^{1}\sum_{k,k^{\prime}=0}^{1}|J+k,J+k\rangle_{M}\langle
J+k^{\prime},J+k^{\prime}|\otimes|kk\rangle_{CD}\langle
k^{\prime}k^{\prime}|\,.$

acting on the stored state (Eq. (6)) reduces to its evaluation for individual
terms (Eqs. (7),(8),(9),(10)) gives

$\displaystyle R_{s}\star\varrho^{U,U}$ $\displaystyle=$
$\displaystyle\frac{2}{3}|U_{\varphi}\rangle\\!\rangle\langle\\!\langle
U_{\varphi}|\,,$ $\displaystyle R_{s}\star\varrho^{U,C}$ $\displaystyle=$
$\displaystyle\frac{1}{6}[|U_{\varphi}\rangle\\!\rangle\langle\\!\langle
U_{\varphi}|+(|00\rangle\langle 00|+|11\rangle\langle 11|)]\,,$ $\displaystyle
R_{s}\star\varrho^{C,U}$ $\displaystyle=$ $\displaystyle
R_{s}\star\varrho^{U,C}\,,$ $\displaystyle R_{s}\star\varrho^{C,C}$
$\displaystyle=$ $\displaystyle\frac{1}{6}(|00\rangle\langle
00|+|11\rangle\langle 11|)\,.$

It follows the retrieval results in the transformation

$\displaystyle
R_{s}\star\Psi_{\varphi}=\operatorname{tr}_{M}\left[R_{s,MCD}\left(\varrho_{M}^{T}\otimes
I_{CD}\right)\right]$
$\displaystyle=\frac{2}{3}\left[\frac{q(1+q)}{2}|U_{\varphi}\rangle\\!\rangle\langle\\!\langle
U_{\varphi}|+\frac{1-q^{2}}{4}(|00\rangle\langle 00|+|11\rangle\langle
11|)\right]$
$\displaystyle=\frac{2}{3}\frac{(1+q)^{2}}{4}\left(\frac{2q}{1+q}|U_{\varphi}\rangle\\!\rangle\langle\\!\langle
U_{\varphi}|+\frac{1-q}{1+q}P\right)\,,$

being a convex mixture of the desired phase gate and the phase damping channel
${\cal P}$. The success probability equals $(1+q)^{2}/6$, thus, it depends on
noise parameter $q$. Let us recall that the retrieved channel possesses
qualitatively different noise if compared to the originally stored channel.

###  IV.2 Dephasing Noise.

We will follow the same calculation as for the $2\rightarrow 1$ case of
depolarization. Choi operator in this case has the following form:

$\displaystyle F_{\varphi}^{\otimes 2}$
$\displaystyle=q^{2}|U_{\varphi}\rangle\\!\rangle_{12}\langle\\!\langle
U_{\varphi}|\otimes|U_{\varphi}\rangle\\!\rangle_{34}\langle\\!\langle
U_{\varphi}|+(1-q)^{2}P_{12}\otimes P_{34}$
$\displaystyle+q(1-q)\left(|U_{\varphi}\rangle\\!\rangle_{12}\langle\\!\langle
U_{\varphi}|\otimes
P_{34}+P_{12}\otimes|U_{\varphi}\rangle\\!\rangle_{34}\langle\\!\langle
U_{\varphi}|\right),$

where $P=\sum_{j=0}^{1}|jj\rangle\langle jj|$. After the storage, the memory
is found in the state

$\displaystyle\Psi_{\varphi}$ $\displaystyle=$ $\displaystyle
F_{\varphi,AB}^{\otimes 2}\star|\Psi\rangle_{AA^{\prime}}\langle\Psi|$
$\displaystyle=$ $\displaystyle
q^{2}\varrho^{U,U}+q(1-q)(\varrho^{P,U}+\varrho^{U,P})+(1-q)^{2}\varrho^{P,P}\,,$

where state $\varrho^{U,U}$ is the same as in equation (7) and

$\displaystyle\varrho^{U,P}$ $\displaystyle=$
$\displaystyle\frac{1}{3}(\Pi_{\overline{00}}+\Pi_{\overline{11}}+\Pi_{\overline{22}}+e^{-i\varphi}|\overline{11}\rangle\langle\overline{22}|+c.c.)\,,$
$\displaystyle\varrho^{P,U}$ $\displaystyle=$
$\displaystyle\frac{1}{3}(\Pi_{\overline{00}}+\Pi_{\overline{11}}+\Pi_{\overline{22}}+e^{-i\varphi}|\overline{00}\rangle\langle\overline{11}|+c.c.)\,,$
$\displaystyle\varrho^{P,P}$ $\displaystyle=$
$\displaystyle\frac{1}{3}(\Pi_{\overline{00}}+\Pi_{\overline{11}}+\Pi_{\overline{22}})\,.$

Using the phase gate PSAR retrieving operation we obtain

$\displaystyle
R_{s}\star\varrho^{U,P}=R_{s}\star\varrho^{P,U}=\frac{1}{3}(|U_{\varphi}\rangle\\!\rangle\langle\\!\langle
U_{\varphi}|+P)_{CD}\,,$ $\displaystyle
R_{s}\star\varrho^{P,P}=\frac{2}{3}P_{CD}\,,\quad
R_{s}\star\varrho^{U,U}=\frac{2}{3}|U_{\varphi}\rangle\\!\rangle\langle\\!\langle
U_{\varphi}|\,.$

Putting it all together we find

$\displaystyle R_{s}\star\Psi_{\varphi}$ $\displaystyle=$
$\displaystyle\operatorname{tr}_{M}\left[R_{s,MCD}\left(\Psi_{\varphi,M}^{T_{M}}\otimes
I_{CD}\right)\right]$ $\displaystyle=$
$\displaystyle\frac{2}{3}\Bigg{[}q|U_{\varphi}\rangle\\!\rangle_{CD}\langle\\!\langle
U_{\varphi}|+\left(1-q\right)P_{CD}\Bigg{]}\,.$

The success probability equals $2/3$ and is independent of noise parameter
$q$. Moreover, the retrieved channel coincides with the noisy phase gate, i.e.
the dephasing noisy phase gates are optimally stored and retrieved by PSAR
protocol.

##  V $N\to 1$ PSAR with noisy phase gates

Let us start the general $N\to 1$ analysis with the depolarization noise.
Performing a direct calculation of the retrieval for PSAR derived in
STORPhaseGate  we find out that for arbitrary $N$ if we leave $\frac{1}{N+1}$
in front of all the parentheses then the factors next to
$|U_{\varphi}\rangle\\!\rangle\langle\\!\langle U_{\varphi}|$ follow the
pattern:

$\displaystyle\begin{split}&1\frac{N}{2^{0}}q^{N}(1-q)^{0}+N\frac{N-1}{2^{1}}q^{N-1}(1-q)^{1}+\\\
&\dots+N\frac{1}{2^{N-1}}q^{1}(1-q)^{N-1}+1\frac{0}{2^{N}}q^{0}(1-q)^{N}\,.\end{split}$
(13)

Similarly the factors in front of $|00\rangle\langle 00|+|11\rangle\langle
11|$ are of the form:

$\displaystyle\begin{split}&1\frac{0}{2^{0}}q^{N}(1-q)^{0}+N\frac{1}{2^{1}}q^{N-1}(1-q)^{1}+\\\
&\dots+N\frac{N-1}{2^{N-1}}q^{1}(1-q)^{N-1}+1\frac{N}{2^{N}}q^{0}(1-q)^{N}\,.\end{split}$
(14)

In these expressions we can identify binomial distribution and the whole
action of the retrieval process results in the retrieved transformation

$\displaystyle\begin{split}&R^{s}\star\Psi_{\varphi}=\frac{1}{N+1}\sum_{k=0}^{N}\binom{N}{k}q^{N-k}(1-q)^{k}\\\
&\Bigg{[}\frac{N-k}{2^{k}}|U_{\varphi}\rangle\\!\rangle\langle\\!\langle
U_{\varphi}|+\frac{k}{2^{k}}(|00\rangle\langle 00|+|11\rangle\langle
11|)\Bigg{]}\\\
&=\frac{N(1+q)^{N}}{2^{N}(N+1)}\Bigg{[}\frac{2q}{1-q}|U_{\varphi}\rangle\\!\rangle\langle\\!\langle
U_{\varphi}|\\\ &\qquad+\frac{1-q}{1+q}(|00\rangle\langle
00|+|11\rangle\langle 11|)\Bigg{]}.\end{split}$

Using this result we can formulate the following theorem.

######  Theorem 1.

Implementation of the optimal phase gate $N\to 1$ probabilistic storing and
retrieval device on noisy phase gates with the depolarization noise, i.e.
${\cal E}_{\varphi}=q{\cal U}_{\varphi}+(1-q){\cal C}_{I/2}$ implements the
noisy channel

${\cal E^{\prime}_{\varphi}}=\frac{2q}{1+q}{\cal
U}_{\varphi}+\frac{1-q}{1+q}{\cal P}\,,$ (15)

where ${\cal P}$ denotes the complete phase damping transformation, i.e. full
diagonalization in the computational basis. The probability of success is
given by the formula

$p_{\rm success}=\frac{N}{N+1}\left(\frac{1+q}{2}\right)^{N}\,.$ (16)

We can see that with the increasing number of uses of noisy phase gates, the
probability of successful retrieval is vanishing as $N$ goes to infinity,
however, the quality of the retrieval is exactly the same.

The same analysis for the case of dephasing noise is more straighforward. The
calculation gives that the retrieved channel is independent of the number of
uses $N$, thus, the optimal PSAR for phase gates does the same job also for
their noisy version if the noise is modeled by the discussed dephasing.

######  Theorem 2.

Implementation of the optimal phase gate $N\to 1$ probabilistic storing and
retrieval device on noisy phase gates with the dephasing noise, i.e. ${\cal
F}_{\varphi}=q{\cal U}_{\varphi}+(1-q){\cal P}$ implements the noisy channel

${\cal F^{\prime}_{\varphi}}=q{\cal U}_{\varphi}+(1-q){\cal P}\,.$ (17)

The probability of success is given by the formula $p_{\rm success}=N/(N+1)$.

###  V.1 Comparison

Fig. 4 illustrates the dependence of success probability of the retrieval on
the noise parameter $q$ for both depolarizing (solid lines) and dephasing
(dashed lines) noisy phase gates for cases $N=1$, $N=3$, $N=7$, and $N=15$.
The probability of success for dephasing noise is higher. In case of
depolarizing channel the success probability goes to $0$ together with $q$
(increasing noise), but there is always some interval of high $q$ for which
the success probability is improving if compared with one use success rate
$1/2$. Surprisingly, for larger degrees of depolarizing noise, more uses do
not lead to improving the success probability.

Figure 4: (Color online) The dependence of the success probability $p_{\rm
succ}$ on the noise parameter $q$ for different number $N$ of uses of the
channel. The solid lines represent the cases of depolarizing noise and dashed
lines illustrate the cases of phase damping noise.

In Fig. 5 we compare the degree of noise of the retrieved transformation
$q^{\prime}$ with the original degree of noise $q$ of phase gates for both
types of noise: depolarizing channel (solid lines) and dephasing noise (dashed
lines). For both cases the relation $q\to q^{\prime}$ is independent on the
number of uses $N$. In case of dephasing $q^{\prime}=q$, thus, the noise is
the same. However, for case of depolarizing noise, $q^{\prime}=2q/(1+q)\geq
q$. Therefore, we may conclude that PSAR reduces the depolarizing noise.
Unfortunately, this noise reduction does not improve with the number of uses.

Figure 5: (Color online) Comparison of noise degrees $q$ (before the storing)
and $q^{\prime}$ (after the retrieval) for depolarizing noise (solid line) and
phase damping noise (dashed line).

##  VI Implementations

In this part we will study two different implementations of optimal PSAR
device for phase gate: first one due to Vidal-Masanes-Cirac from Ref.
storingQuantDyn  limited to $N=2^{k}-1$ uses and the one reported in Ref.
STORPhaseGate  minimizing the size of the storage register to single qudit.
Both are equivalent if noiseless phase gates are considered, however, in the
presence of noise they may lead to different results.

###  VI.1 Vidal-Masanes-Cirac

First of the implementations originally proposed by Vidal, Masanes and Cirac
in Ref. storingQuantDyn  is depicted in the figure 6. The $N=1$ storage
consists of the application of the noisy phase gate on the state
$|+\rangle=(|0\rangle+|1\rangle)/\sqrt{2}$ to obtain the state

$\Psi_{\varphi}={\cal{E}}_{\varphi}(|+\rangle\langle+|)={\cal{F}}_{\varphi}(|+\rangle\langle+|)=qU_{\varphi}{|\xi\rangle\langle\xi|}U_{\varphi}^{\dagger}+\frac{1-q}{2}I\,.$

Let us stress this state is the same for both considered noises. The retrieval
step is composed of the application of the controled NOT gate $U_{\rm CNOT}$
on the unknown state $|\xi\rangle=a|0\rangle+b|1\rangle$ (the control system)
and the stored state $\Psi_{\varphi}$ (the target system), i.e. the final
state reads

$\displaystyle\Xi_{\varphi}^{\prime}$ $\displaystyle=$ $\displaystyle U_{\rm
CNOT}[{|\xi\rangle\langle\xi|}\otimes\Psi_{\varphi}]$ $\displaystyle=$
$\displaystyle\frac{1}{2}[qU_{\varphi}{|\xi\rangle\langle\xi|}U_{\varphi}^{\dagger}+(1-q){\rm
diag}(\xi)]\otimes{|0\rangle\langle 0|}$
$\displaystyle+\frac{1}{2}[qU_{-\varphi}{|\xi\rangle\langle\xi|}U_{-\varphi}^{\dagger}+(1-q){\rm
diag}(\xi)]\otimes{|1\rangle\langle 1|}\,,$

where ${\rm diag}(\xi)=\langle 0|\xi|0\rangle{|0\rangle\langle 0|}+\langle
1|\xi|1\rangle{|1\rangle\langle 1|}$. Measuring the memory register in the
computational basis the outcome value 0 identifies the successful realization
of the noisy phase gate

$|\xi\rangle\langle\xi|\mapsto
qU_{\varphi}{|\xi\rangle\langle\xi|}U_{\varphi}^{\dagger}+(1-q){\rm
diag}(\xi)\,,$

with success probability $p_{\rm success}=1/2$. The case of failure can be
corrected by storing and retrieving sequence of two phase gates

$\displaystyle\Psi^{(2)}_{\varphi}$ $\displaystyle=$
$\displaystyle{\cal{E}}_{\varphi}({\cal{E}}_{\varphi}({|+\rangle\langle+|}))={\cal{F}}_{\varphi}({\cal{F}}_{\varphi}({|+\rangle\langle+|}))$
$\displaystyle=$ $\displaystyle
q^{2}U_{2\varphi}{|\xi\rangle\langle\xi|}U_{2\varphi}^{\dagger}+\frac{1-q^{2}}{2}I\,.$

Applying $U_{\rm CNOT}$ as before on the reused state $\xi_{\rm
fail}^{(1)}=qU_{-\varphi}|\xi\rangle\langle\xi|U_{-\varphi}^{\dagger}+(1-q){\rm
diag}(\xi)$ that left after failure and $\Psi^{(2)}_{\varphi}$ results in

$\displaystyle\Xi^{\prime(2)}$ $\displaystyle=$ $\displaystyle U_{\rm
CNOT}(\xi_{\rm fail}^{(1)}\otimes\Psi^{(2)}_{\varphi})$ $\displaystyle=$
$\displaystyle\frac{1}{4}[q^{2}U_{\varphi}{|\xi\rangle\langle\xi|}U_{\varphi}^{\dagger}+(1-q^{2}){\rm
diag}(\xi)]\otimes{|0\rangle\langle 0|}$
$\displaystyle+\frac{1}{4}[q^{2}U_{-3\varphi}{|\xi\rangle\langle\xi|}U_{-3\varphi}^{\dagger}+(1-q^{2}){\rm
diag}(\xi)]\otimes{|1\rangle\langle 1|}\,,$

with success probability $1/4$. We can iteratively continue to correct the
failures by storing and retrieving sequence of $2^{k}$ phase gates in the
$k$th correction. In particular,

$\displaystyle\Xi^{\prime(k)}$ $\displaystyle=$ $\displaystyle U_{\rm
CNOT}(\xi_{\rm fail}^{(k-1)}\otimes\Psi^{(k)}_{\varphi})$ $\displaystyle=$
$\displaystyle\frac{1}{2^{k}}[q^{N}{\cal U}_{\varphi}(\xi)+(1-q^{N}){\rm
diag}(\xi)]\otimes{|0\rangle\langle 0|}$ $\displaystyle+$
$\displaystyle\frac{1}{2^{k}}[q^{N}{\cal U}_{-N\varphi}(\xi)+(1-q^{N}){\rm
diag}(\xi)]\otimes{|1\rangle\langle 1|}\,,$

where $N=2^{k}-1$ equals to the total number of uses of noisy phase gate.
Summing up the success probabilities we obtain for the total success rate

$p_{\rm
success}=\frac{1}{2}+\frac{1}{4}+\cdots\frac{1}{2^{k}}=N/(N+1)=1-\frac{1}{2^{k}}\,,$

but let us stress that for each $k$ the implemented channel is different.
Therefore the above success probability cannot be associated with
implementation of particular noisy channel. The larger the $k$ the smaller the
success probability and the noisier the retrieved phase gate.

Figure 6: Vidal-Masanes-Cirac realization scheme for arbitrary $N=2^{k}-1$
times of applying the channel ${\cal{E}}_{\varphi}$ $N$ times with $k$ being
the number of qubits.

###  VI.2 Virtual Qudit

The implementation of PSAR for arbitrary $N$ introduced in Ref. STORPhaseGate
is removing the use of ancilla $A^{\prime}$. Instead of
$|\Psi\rangle=\frac{1}{\sqrt{N+1}}\sum_{j}|\overline{jj}\rangle$ this
construction uses the state
$|\Omega\rangle=\frac{1}{\sqrt{N+1}}\sum_{j}|\overline{j}\rangle$, where
$|\overline{j}\rangle=|0^{\otimes(N-j)}\rangle\otimes|1^{\otimes j}\rangle$
are states of $N$ qubits system (see Fig. 7). Each of the qubits is
transformed by the phase gate to store the action in the state

$\Omega_{\varphi}={\cal E}_{\varphi}\otimes\cdots\otimes{\cal
E}_{\varphi}(|\Omega\rangle\langle\Omega|)\,.$

Without noise the states $\Omega_{\varphi}$ belong to $d=(N+1)$-dimensional
subspace. Therefore, the retrieval operation is considered only on this
virtual qudit subspace ${\cal H}_{d}\subset{\cal H}_{2}^{\otimes N}\equiv{\cal
H}_{A}$. However, in our case whole Hilbert space matters. Following the
original implementation from Ref. STORPhaseGate  we introduce the conditional
shift as follows

$\displaystyle{\cal{C}}_{\ominus}(|c\rangle\otimes|\overline{t}\rangle)$
$\displaystyle=$ $\displaystyle|c\rangle\otimes|\overline{t\ominus c}\rangle$
(21)
$\displaystyle{\cal{C}}_{\ominus}(|c\rangle\otimes|\overline{t}_{\perp}\rangle)$
$\displaystyle=$
$\displaystyle|c\rangle\otimes|\overline{t}_{\perp}\rangle\,,$

for all vector states $|\overline{t}_{\perp}\rangle$ orthogonal to all vectors
$|\overline{j}\rangle$. In other words the conditional shift acts as identity
on the subspace orthogonal to virtual qudit system.

Figure 7: Implementation of noisy channel ${\cal{E}}_{\varphi}^{\otimes N}$
using virtual qudit and shift-down operator ${\cal{C}}_{\ominus}$ defined in
equation (21).

####  VI.2.1 Dephasing noise.

Let us start with the dephasing noise and perform the calculation explicitly
for the case of $N=2$. The memory state equals

$\displaystyle\Omega_{\varphi}^{(2)}$ $\displaystyle=$
$\displaystyle{\cal{F}}_{\varphi}^{\otimes 2}({|\Omega\rangle\langle\Omega|})$
$\displaystyle=$ $\displaystyle q^{2}{\cal{U}}_{\varphi}^{\otimes
2}({|\Omega\rangle\langle\Omega|})+(1-q)^{2}{\cal{P}}^{\otimes
2}({|\Omega\rangle\langle\Omega|})$
$\displaystyle+q(1-q)\big{[}({\cal{U}}_{\varphi}\otimes{\cal{P}})+({\cal{P}}\otimes{\cal{U}}_{\varphi})\big{]}({|\Omega\rangle\langle\Omega|})\,,$

where

$\displaystyle\begin{split}&{\cal{U}}_{\varphi}^{\otimes
2}({\Omega})=\frac{1}{3}[\Pi_{\overline{012}}+e^{-i\varphi}(X_{\overline{01}}+X_{\overline{12}})+e^{-2i\varphi}X_{\overline{02}}+cc]\\\
&({\cal{U}}_{\varphi}\otimes{\cal{P}})(\Omega)=\frac{1}{3}[\Pi_{\overline{012}}+e^{-i\varphi}X_{\overline{12}}+e^{i\varphi}X_{\overline{21}}],\\\
&({\cal{P}}\otimes{\cal{U}}_{\varphi})({\Omega})=\frac{1}{3}[\Pi_{\overline{012}}+e^{-i\varphi}X_{\overline{01}}+e^{i\varphi}X_{\overline{10}}],\\\
&({\cal{P}}\otimes{\cal{P}})(\Omega)=\frac{1}{3}\Pi_{\overline{012}}\,,\end{split}$
(23)

where we used the notation $\Omega={|\Omega\rangle\langle\Omega|}$,
$\Pi_{\overline{012}}={|\overline{0}\rangle\langle\overline{0}|}+{|\overline{1}\rangle\langle\overline{1}|}+{|\overline{2}\rangle\langle\overline{2}|}$
and $X_{\overline{jk}}={|\overline{j}\rangle\langle\overline{k}|}$. Retrieval
then results in the state

$\displaystyle{\cal{C}}_{\ominus}[\xi\otimes\Omega^{(2)}_{\varphi}]=\frac{1}{3}[q^{2}U_{\varphi}\xi
U_{\varphi}^{\dagger}+(1-q^{2}){\rm diag}(\xi)]\otimes\Pi_{\overline{01}}+$
$\displaystyle+\frac{1}{3}[q^{2}U_{-2\varphi}\xi
U_{-2\varphi}^{\dagger}+(1-q^{2}){\rm
diag}(\xi)]\otimes{|\overline{2}\rangle\langle\overline{2}|}\,,$

where
$\Pi_{\overline{01}}={|\overline{0}\rangle\langle\overline{0}|}+{|\overline{1}\rangle\langle\overline{1}|}$
and ${\rm diag}(\xi)=\langle 0|\xi|0\rangle{|0\rangle\langle 0|}+\langle
1|\xi|1\rangle{|1\rangle\langle 1|}$. It follows the retrieval is successful
with probability $p_{\rm success}=2/3$, when the channel $\xi\rightarrow
q^{2}U_{\varphi}\xi U_{\varphi}^{\dagger}+(1-q^{2}){\rm diag}(\xi)$ is
retrieved. The noise of the retrieved channel is higher.

The calculation for more general case of $N\to 1$ storage and retrieval
protocol gives

$\displaystyle{\cal{C}}_{\ominus}[\xi\otimes\Omega_{\varphi}^{(N)}]=$
$\displaystyle=\frac{1}{N+1}[q^{N}U_{\varphi}\xi
U_{\varphi}^{\dagger}+(1-q^{N}){\rm
diag}(\xi)]\otimes\Pi_{\overline{01...(N-1)}}$
$\displaystyle+\frac{1}{N+1}[q^{N}U_{-N\varphi}\xi
U_{-N\varphi}^{\dagger}+(1-q^{N}){\rm
diag}(\xi)]\otimes{|\overline{N}\rangle\langle\overline{N}|}\,,$

where
$\Pi_{\overline{01...(N-1)}}=\sum_{j=0}^{N-1}{|\overline{j}\rangle\langle\overline{j}|}$.
As $N$ is increasing the probability of ”success” is increasing, however, also
the noise of the retrieved channel is increasing and in the limit completely
diminishes the dependence on $\varphi$ and converges to purely dephasing noise
${\cal P}$.

####  VI.2.2 Depolarizing noise.

It is more involved to get the general formula for the general case of
depolarizing noise. In what follows we will explicitly investigate the case of
$2\to 1$ PSAR of depolarized phase gates and illustrate the behavior. After
applying the phase gate twice its action is stored in the state

$\displaystyle\Omega_{\varphi}^{(2)}$ $\displaystyle=$
$\displaystyle{\cal{E}}_{\varphi}^{\otimes
2}({|\Omega\rangle\langle\Omega|})=q^{2}{\cal{U}}_{\varphi}^{\otimes
2}(\Omega)+(1-q)^{2}{\cal{C}}_{I/2}^{\otimes 2}(\Omega)$
$\displaystyle+q(1-q)[{\cal{U}}_{\varphi}\otimes{\cal{C}}_{{I}/2}+{\cal{C}}_{{I}/2}\otimes{\cal{U}}_{\varphi}](\Omega)\,,$

where

$\displaystyle{\cal{U}}_{\varphi}^{\otimes
2}(\Omega)=\frac{1}{3}[\Pi_{\overline{012}}+e^{i\varphi}(X_{\overline{10}}+X_{\overline{21}})+e^{2i\varphi}X_{\overline{20}}+cc]\,,$
$\displaystyle{\cal{C}}_{{I}/2}^{\otimes
2}(\Omega)=\frac{1}{4}[{|0\rangle\langle 0|}+{|1\rangle\langle
1|}+{|2\rangle\langle 2|}+{|3\rangle\langle 3|}]\,,$
$\displaystyle({\cal{U}}_{\varphi}\otimes{\cal{C}}_{{I}/2})(\Omega)=\frac{1}{6}[2\Pi_{\overline{01}}+\Pi_{\overline{23}}+e^{i\varphi}(X_{\overline{30}}+X_{\overline{21}})+cc]\,,$
$\displaystyle({\cal{C}}_{{I}/2}\otimes{\cal{U}}_{\varphi})(\Omega)=\frac{1}{6}[2\Pi_{\overline{12}}+\Pi_{\overline{03}}+e^{i\varphi}(X_{\overline{10}}+X_{\overline{23}})+cc]\,,$

where we used the same notation as before. The retrieval outputs

$\displaystyle{\cal{C}}_{\ominus}[\xi\otimes\Omega_{\varphi}^{(2)}]=$
$\displaystyle=\frac{1}{6}q(1-q)\left[\xi_{00}^{2}{|0\rangle\langle
0|}\otimes{|\overline{1}\rangle\langle\overline{1}|}+\xi_{11}^{2}{|1\rangle\langle
1|}\otimes{|\overline{0}\rangle\langle\overline{0}|}\right]$
$\displaystyle+\left[\frac{1}{6}q(q+1)U_{\varphi}\xi
U_{\varphi}^{\dagger}+\left(\frac{1}{4}-\frac{1}{12}q^{2}-\frac{1}{6}q\right){\rm
diag}(\xi)\right]\otimes\Pi_{\overline{01}}$
$\displaystyle+\left[\frac{1}{3}q^{2}U_{-2\varphi}\xi
U_{-2\varphi}^{\dagger}+\frac{1}{4}(1-q^{2}){\rm
diag}(\xi)\right]\otimes{|\overline{2}\rangle\langle\overline{2}|}$
$\displaystyle+\left[\frac{1}{4}-\frac{1}{12}q^{2}-\frac{1}{6}q\right]\xi\otimes{|\overline{3}\rangle\langle\overline{3}|}\,.$
(24)

Measurement associated with the projection $\Pi_{\overline{01}}$ corresponds
to successful measurement while measuring $\Pi_{\overline{23}}$ corresponds to
a failure. The success probability equals $p_{\rm success}=(3+q)/6$, thus,
becomes dependent on the initial noise, and the retrieved channel is a mixture
of the desired phase gate and the diagonalisation in the computational basis
(complete decoherence). The retrieved channel is again more noise than the
original although it must be stressed it is from different family of channels.
That is for this implementation the PSAR device does no reduce the noise. Let
us also note that in the case of outcome
${|\overline{3}\rangle\langle\overline{3}|}$ the qubit state is unaffected.

##  VII Conclusions

We addressed the question of noise robustness of particular example of higher-
order quantum information processing task for storage and retrieval (quantum
learning) of quantum processes. Surprisingly we found that for the case of
depolarisation noise the PSAR device can reduce the noise level both
quantitatively and qualitatively. This feature is not generic, but when it
happens PSAR may be used to improve the performance of noisy quantum phase
gate and probabilistically “distill” the noiseless evolution. In particular,
the white noise is not only reduced, but it also changes to dephasing noise.
We also observed that different implementations of PSAR have different
performances for noisy phase gates, because their actions on the subspace
unused in the noiseless case are different.

###### Acknowledgements.

We acknowledge the support by the projects OPTIQUTE (APVV-18-0518), DESCOM
(VEGA 2/0183/21). MZ acknowledges the support of the John Templeton Foundation
under the project ID JTF-61466 (QISS). The opinions expressed in this
publication are those of the authors and do not necessarily reflect the views
of the John Templeton Foundation.

## References

  * (1) P. Shor, SIAM J. Computing 26 (1997), 1484-1509.
  * (2) L. K. Grover, Proceedings, 28th Annual ACM Symposium on the Theory of Computing (STOC), May 1996, pages 212-219
  * (3) C. Bennett and G. Brassard, Theoretical Computer Science - TCS 560, 175-179 (1984)
  * (4) W.K. Wootters and W.H. Zurek, Nature 299, 802 (1982)
  * (5) H. Barnum, C. M. Carlton, C. A. Fuchs, R. Jozsa and B. Schumacher, Phys. Rev. Lett. 76, 2818-2821 (1996)
  * (6) S. L. Braunstein and A. K. Pati, Phys. Rev. Lett. 98, 080502 (2007)
  * (7) A. Pati and S. Braunstein, Nature 404, 164-165 (2000)
  * (8) M. A. Nielsen and Isaac L. Chuang, Phys. Rev. Lett. 79, 321 (1997)
  * (9) M. Hillery, V. Bužek and M. Ziman, Phys. Rev. A 65 (2001)
  * (10) V. Bužek, M. Ziman and M. Hillery, Fortschritte der Physik 52, 1056-1063 (2004)
  * (11) M. Hillery, M. Ziman and V. Bužek, Phys. Rev. A 69, 042311 (2004)
  * (12) A. M. Kubicki, C. Palazuelos, D. Perez-Garcia, Phys. Rev. Lett. 122, 080505 (2019)
  * (13) A. Bisio, G. Chiribella, G. M. D’Ariano, S. Facchini and P. Perinotti, Phys. Rev. A 81, 032324 (2010)
  * (14) G. Vidal, L. Masanes and J. I. Cirac, Phys. Rev. Lett. 88, 047905 (2002)
  * (15) M. Sedlák, A. Bisio and M. Ziman, Phys. Rev. Lett. 122, 170502 (2019)
  * (16) M. Sedlák and M. Ziman, Phys. Rev. A 102, 032618 (2020)
  * (17) G. Chiribella, G. M. D’Ariano and P. Perinotti, Phys. Rev. A 80, 022339 (2009)
  * (18) G. Chiribella, G. M. D’Ariano and P. Perinotti, Phys. Rev. Lett. 101, 060401 (2008)
  * (19) A. Bisio and P. Perinotti, Proc. R. Soc. A 475 20180706 (2019)
  * (20) A. Bisio, G. Chiribella, G. M. D’Ariano and P. Perinotti, Acta Physica Slovaca 61, 273-390 (2011)