#  Phase conjugation with spatially incoherent light in complex media

YoonSeok Baek<EMAIL_ADDRESS>Laboratoire Kastler Brossel,
ENS–Universite PSL, CNRS, Sorbonne Université, Collège de France, 24 Rue
Lhomond, F-75005 Paris, France  Hilton B. de Aguiar  Laboratoire Kastler
Brossel, ENS–Universite PSL, CNRS, Sorbonne Université, Collège de France, 24
Rue Lhomond, F-75005 Paris, France  Sylvain Gigan  Laboratoire Kastler
Brossel, ENS–Universite PSL, CNRS, Sorbonne Université, Collège de France, 24
Rue Lhomond, F-75005 Paris, France

###### Abstract

Shaping light deep inside complex media, such as biological tissue, is
critical to many research fields. Although the coherent control of scattered
light via wavefront shaping has made significant advances in addressing this
challenge, controlling light over extended or multiple targets without
physical access inside a medium remains elusive. Here we present a phase
conjugation method for spatially incoherent light, which enables the non-
invasive light control based on incoherent emission from multiple target
positions. Our method characterizes the scattering responses of hidden sources
by retrieving mutually incoherent scattered fields from speckle patterns. By
time-reversing scattered fluorescence with digital phase conjugation, we
experimentally demonstrate focusing of light on individual and multiple
targets. We also demonstrate maximum energy delivery to an extended target
through a scattering medium by exploiting transmission eigenchannels. This
paves the way to control light propagation in complex media using incoherent
contrasts mechanisms.

##  I Main

Delivering optical energy and transmitting information through complex media
remains an important challenge in many fields of studies, including optical
manipulation [9], deep-tissue imaging [27, 23] and optogenetics [5, 47]. In
recent years, it has been shown that the coherent control of scattered light
can manipulate spatial, spectral and temporal distributions of light in
scattering media [40, 26, 17, 33, 8]. However, such capabilities are greatly
limited without physical access inside a medium because the scattering
response to target position is difficult to characterize. As a result, non-
invasive light control over extended or multiple targets remains elusive
despite being crucial for real-world applications. Optimizing incident
wavefront based on a feedback signal [20, 11, 3, 4, 24, 31, 37, 38] is mostly
limited to focusing on a single isolated target, and even then it has
limitations that require numerous changes of the wavefront. While time-
reversal or phase conjugation techniques [46, 10, 18, 44, 41, 19, 25, 48, 35,
34, 45, 1] allow for effective light delivery to an optical or virtual source,
they cannot individually control light on multiple targets.

Here we address these challenges by extending phase conjugation to spatially
incoherent light. Our approach utilizes incoherent emission from multiple
targets. We first characterize scattering responses of these hidden sources by
retrieving mutually incoherent fields from spatially modulated speckle
patterns. The retrieved fields are related to the field transmission matrix
[29], and their phase conjugation enables light control over the desired
positions. We demonstrate this experimentally by focusing light on individual
and multiple fluorescent targets through a scattering medium. Finally, we show
that transmission eigenchannels can be identified by decomposing the
incoherent fields and demonstrate maximum energy delivery to a hidden extended
target.

##  II Results

A schematic of the experiment is illustrated in Fig. 1. We consider a scenario
where multiple fluorescent targets are hidden by a scattering medium and act
as guidestars. Fluorescence emitted by these guidestars is scattered,
resulting in an incoherent addition of speckle patterns on the camera. Under
this condition, we aim to deliver light back to each of the guidestars by
time-reversing the scattered fluorescence. Our approach consists of retrieving
multiple incoherent fields that compose the fluorescence, and using them to
generate phase-conjugated beams. To this end, we introduce wavefront
modulation of the scattered fluorescence with a spatial light modulator (SLM).
This modulation induces changes in the measured incoherent speckle patterns,
providing information to retrieve the scattered fields that will later be used
for phase conjugation.

Figure 1: Schematic of phase conjugation with incoherent fluorescence:
Multiple fluorescent guidestars are hidden behind a scattering medium. The
scattered fluorescence, composed of mutually incoherent speckle fields, is
modulated by an SLM. The modulated fluorescence is Fourier transformed by a
lens and results in the incoherently added speckle patterns on the camera.
Then mixed-state phase retrieval recovers a set of scattered fields, whose
phase conjugation enables targeted light control.

To explain the retrieval process in detail, we introduce a partial field
transmission matrix $\bm{T}$ whose input and output are fields at the SLM
plane and guidestar positions, respectively,
$\bm{E}^{guidestar}=\bm{T}\bm{E}^{SLM}$. According to the time-reversal
symmetry, the scattering of fluorescence emitted by $N$ guidestars at the same
wavelength can be expressed by the rows of the transmission matrix,
$\bm{T}=[\bm{t}_{1},\dots,\bm{t}_{N}]^{\top}$, where $\bm{t}_{n}$ represents
the scattered field that corresponds to an individual guidestar. For the sake
of simplicity, we will assume that the brightness of each guidestar is the
same (see Supplementary Information Section 6 for general cases with different
brightness). Then, the incoherent speckle pattern measured by the camera is
expressed as
$\bm{I}_{0}=\textnormal{diag}{\left[(\bm{TF})^{\dagger}\bm{TF}\right]}$, where
$\bm{F}$ represents the discrete Fourier transform between the SLM and camera
plane (see Supplementary Information Section 1 for the derivation). If we
apply $M$ different wavefront modulations using the SLM, and express the $m$th
modulation as a diagonal matrix $\bm{S}_{m}$, the incoherent speckle pattern
after the modulation is

$\bm{I}_{m}=\textnormal{diag}{\left[\left(\bm{TS}_{m}\bm{F}\right)^{\dagger}\bm{TS}_{m}\bm{F}\right]}.$
(1)

To find a set of scattered fields that satisfy Eq. 1, we used an iterative
approach that minimizes the error between the measured and predicted
$\bm{I}_{m}$. We began by making an initial guess of $\bm{t}_{n}$ and applying
the wavefront modulation numerically according to Eq. 1. Then, following the
mixed-state reconstruction [36], we corrected the amplitude parts of the
modulated speckle fields at the camera plane, using the ratio of the measured
to the predicted $\bm{I}_{m}$. Next, we compensated for the wavefront
modulation and updated the guess of $\bm{t}_{n}$. These steps were repeated
for all the wavefront modulations, resulting in a maximum likelihood
estimation of $\bm{T}$ (see Supplementary Information Section 2 for more
information). We note that there is inherent ambiguity in determining
$\bm{T}$, as $\bm{T}$ and its unitary transformation are indistinguishable
based on intensity. This can be confirmed by replacing $\bm{T}$ with $\bm{UT}$
in Eq. 1. For this reason, the scattered fields are retrieved as a mixture of
$\bm{t}_{n}$:

$\bm{H}=\bm{UT},$ (2)

where $\bm{H}$ is a set of retrieved fields $\bm{h}_{n}$,
$\bm{H}=[\bm{h}_{1},...,\bm{h}_{N}]^{\top}$, and $\bm{U}$ is an arbitrary
unitary matrix. We further note that Eq. 2 arises from the nature of mutual
incoherence, rather than from the reconstruction method. Despite this
ambiguity, the retrieved fields offer unique capabilities for phase
conjugation, as we show below.

### Incoherent phase conjugation

The time-reversal of the entire scattered fluorescence will regenerate light
at hidden sources, creating foci on the entire targets. One way to accomplish
this is through the ensemble average of the phase-conjugated scattered fields
$\bm{t}_{n}^{*}$. Alternatively, we chose to use the ensemble average of
$\bm{h}_{n}^{*}$ because it gives the identical phase conjugation result (see
Supplementary Information Section 3). In both cases, the ensemble average
results in the incoherent sum of the phase conjugation of the scattered
fields, which can be realized by shaping either coherent or incoherent light.
We will refer to this approach as “Incoherent phase conjugation” for
simplicity, regardless of the coherence of light used to generate the phase-
conjugated fields (see Discussion for more information).

To demonstrate this incoherent phase conjugation, we introduced several 1
$\mu$m fluorescent beads as guidestars behind the scattering medium. We
retrieved the multiple scattered fields according to the number of the
guidestars. In our experiments, each 1 $\mu$m bead was considered as an
individual guidestar as the speckle grain size at the target plane was
$\sim$0.9 $\mu$m. To implement the incoherent phase conjugation, we used the
SLM and a collimated laser beam to generate $N$ different phase-conjugated
fields of $\bm{h}_{n}^{*}$ over time and measured their time-averaged response
for phase (see Methods and Discussion). To evaluate the performance of the
phase conjugation, we first conducted an experiment with a single fluorescent
bead (Fig. 2a). By phase-conjugating the scattered field, we observed a strong
focus on the bead (Fig. 2b). This is in clear contrast with the random speckle
generated by a beam with a random wavefront (Fig. 2c). The enhancement factor,
defined as the ratio between the optimized focus intensity and mean background
intensity, was $\sim$4,400. We then placed multiple fluorescent beads (Fig.
2d–f) behind the scattering medium. By incoherently phase conjugating the
scattered fields, we successfully generated foci at every guidestar positions
(Fig. 2g–i). Despite the minimal spectral memory effect [43], we were also
able to excite the bead through the scattering medium by shaping the
excitation beam using the retrieved field at the emission wavelength(Fig. S1).

Figure 2: Incoherent phase conjugation for multiple targets: a, The
fluorescence image of an 1 $\mu$m bead taken from the side without a
scattering medium. b, Intensity at the target plane generated by phase
conjugation based on the scattered fluorescence field. (Inset) The phase of
the scattered field shown in the HSV colormap. The central highlighted part is
used for the phase conjugation. c, Intensity at the target plane generated by
a random wavefront. d–f, The fluorescence images of multiple beads hidden
behind the scattering medium. g–i, Intensity at the target plane with
incoherent phase conjugation of scattered fields.

### Selective focusing on individual targets

In order to selectively focus on individual targets, it is necessary to demix
the individual fields $\bm{t}_{n}$ from their mixture $\bm{h}_{n}$. We note
that $\bm{t}_{n}$ is not strictly orthogonal, and thus the orthogonalization
of $\bm{h}_{n}$ can not be a solution. Our solution was to directly invert Eq.
2 by finding $\bm{U}$. To this end, we utilized the memory effect [14, 13],
where neighboring guidestars generate correlated speckle patterns.
Specifically, we iteratively applied a random unitary transformation to the
retrieved fields $\bm{H}$, such that the correlation between the transformed
speckle patterns is maximized (see Supplementary Information Section 5).

Figure 3 shows the experimental result with 5 fluorescent beads. When the
scattered fields $\bm{h}_{n}$ are directly used for phase conjugation, each
phase conjugation generated foci on several guidestars with different
intensities (Fig. 3b). The demixed fields, on the other hand, generated a
focus on a single guidestar, showing that $\bm{t}_{n}$ is successfully
recovered (Fig. 3c). The memory effect range in this experiment, defined as
the full width at half maximum of speckle cross-correlation, was 5 $\mu$m,
which is much smaller than the spatial extent of the guidestars. This result
shows that the selective focusing via demixing is possible as long as a pair
of guidestars lies within the memory effect range.

Figure 3: Selective focusing on individual targets: a, The fluorescence image
of a target comprised of 1 $\mu$m beads. b,c, Intensity at the target plane
after the phase conjugation of individual scattered fields before (b) and
after (c) the demixing process. (Insets) The phase of the scattered fields
used for phase conjugation.

### Targeted energy delivery

Maximum energy delivery through scattering media requires an eigenchannel of
$\bm{T}$ [42, 22]. The transmission eigenchannels correspond to the singular
vectors of $\bm{T}$, and the first singular vector with the largest singular
value, delivers the maximum energy to the target. Although the direct access
to $\bm{T}$ is not always possible, the transmission eigenchannels of $\bm{T}$
can be found using $\bm{H}$. This is because the eigenchannels of $\bm{T}$ and
$\bm{H}$ are identical because $\bm{H}^{\dagger}\bm{H}=\bm{T}^{\dagger}\bm{T}$
according to Eq. 2. Thus, we can deliver the maximum energy to extended
targets using the first singular vector of $\bm{H}$.

To demonstrate the targeted energy delivery, we placed a 5 $\mu$m fluorescent
ink droplet behind the scattering medium (Fig. 4a). Based on the size of the
target, we estimated the number of incoherent fields and retrieved 23
scattered fields. We then performed the singular value decomposition of the
retrieved fields $\bm{H}$. Finally, we injected fields that corresponds to the
singular vectors $\bm{v}_{n}$ and observed the energy delivered to the target.
When a random wavefront is injected to the scattering medium, a speckle
pattern is generated at the target plane (Fig. 4b). In contrast, the singular
vectors produce intensity distributions highly concentrated on the target
(Fig. 4c). By summing the results of all the singular vectors, we confirmed
that the energy is delivered only to the target area (Fig. 4d). The first
singular vector $\bm{v}_{1}$ shows a 174-times increase in the energy on the
target, compared to random wavefronts. The enhancement decays with the
singular vector index [Fig. 4(e)]. We observed that the values are not
perfectly sorted in a descending order, which we believe is due to the
numerical error in the retrieved fields and to the use of the phase-only SLM.

Figure 4: Targeted energy delivery: a, The fluorescence image of an extended
target hidden behind a scattering medium. The image was taken from the side
without the scattering medium. b, Intensity at the target plane when a random
phase pattern is displayed on the SLM. c, Intensity at the target plane using
the first 3 singular vectors of $\bm{H}$. The images are normalized for the
result of the first singular vector. d, The sum of all the results using 23
singular vectors. e, Enhancement of the energy delivered to the target
compared to random realizations. Dash circles in (c) and (d) indicate the
boundary of the target.

##  III Discussion and conclusion

We have presented an approach to control light in scattering media without the
physical access to a target plane by extending phase conjugation to spatially
incoherent light. We have demonstrated focusing and maximum energy delivery to
extended fluorescent targets through a scattering medium. Our method, like
other digital phase conjugation or wavefront shaping techniques, can provide
reliable focusing results with a high contrast provided a sufficiently large
number of incident modes are controlled [26]. Nonetheless, what distinguishes
our approach is its ability to tackle mutually incoherent fields of scattered
fluorescence. Another important aspect is that it does not require the precise
alignment between the camera and SLM [16], since the scattered fields are
retrieved at the SLM plane. Its principle of characterizing the scattering
response is entirely passive, as it does not alter the emission of guidestars,
as opposed to to techniques that modulates the excitation wavefront (e.g.
[4]). It only requires that the emission from the sources be constant on
average over the measurement of $\bm{I}_{m}$ and thus unaffected by the setup
geometry.

To further demonstrate the capability of the proposed method under realistic
conditions, we conducted an additional experiment. Specifically, we employed a
volumetric scattering medium composed of parafilm layers, and we increased the
number of guidestars. More information on the experiment can be found in
Supplementary Information Section 4. The results showed that the proposed
method was able to produce high-contrast foci on 13 and 22 target positions
(Fig. S2), demonstrating the potential of our method for practical
applications.

For the proposed method, it is important to estimate the number of mutually
incoherent waves, $N$. This is because the underestimation of $N$ results in
imperfect reconstruction of the scattered fields. We note that the
overestimation is allowed because it results in redundant reconstruction (see
Fig. S3). Nevertheless, it is recommended to use the smallest possible value
of $N$ for the minimal measurements and computation time. There are different
methods to estimate $N$. The contrast of fluorescence speckle is an useful
indicator for $N$, as it decreases as $\sqrt{N}$ [15]. It is also possible to
find $N$ by analyzing the error in the mixed-state phase retrieval or the
singular value distribution [30] for different values of $N$. In our
demonstrations, we did not consider the spectral degrees of freedom because
narrow spectral responses were measured by using interference filters. If the
detection bandwidth is greater than the spectral memory effect range,
different spectral components should be considered in estimating $N$. This
will allow the proposed technique to handle situations where hidden sources
emit at different wavelengths.

Another important consideration is the number of modulation $M$ required for
the field retrieval, which scales linearly with $N$. In experiments, reliable
phase conjugation results were obtained when $M\geq 6N$ (see Supplementary
Information Section 2). This linearity can be attributed to the multiplexed
information in the intensity of multiple incoherent fields. We emphasize that
$M$ scales with the number of incoherent waves $N$ not with the number of
controlled input modes of the SLM. Recent advances in phase retrieval [7, 6]
show that few measurements are sufficient in retrieving a coherent field
($N=1$). In this regard, we believe that in principle even fewer $M$ may be
used for our method.

In our proof-of-principle experiments, we used the simple algorithms for the
retrieval and demixing of the scattered fields. The performance of the
algorithms can be enhanced by incorporating constraints, convex optimization
[6] or the generalized memory effect [28]. The total measurement time can be
reduced by designing a setup with minimal energy loss, and by using a
sensitive detector, such as an EMCCD.

The current implementation of incoherent phase conjugation is based on the
alternation of multiple coherent fields. This method can accurately generate
the incoherent response of the time-reversed fluorescence, although it
requires some integration time [12]. Importantly, the required integration
time can be reduced to few microseconds by employing fast spatial light
modulators [32, 39]. Alternatively, the incoherent phase conjugation can be
implemented by directly shaping an incoherent source [2]. In this case, the
incoherent response will be obtained at a timescale greater than the coherence
time of the source.

In conclusion, our method enables versatile light control over extended or
multiple targets using incoherent contrast mechanisms. The concept can be
applied to different incoherent emissions, such as spontaneous Raman
scattering [37, 38], and a wide range of photoluminescence [21]. Furthermore,
it enables the passive characterization of a transmission matrix, opening up
the possibility to generalized light control using transmission-matrix-based
operators [8]. We envision that the proposed approach will enable targeted
light delivery through thick biological tissue, facilitating biomedical
applications, such as optogenetic stimulation and phototherapy.

##  IV Methods

### Experimental setup

The experimental setup is shown in Fig. S4. A laser diode ($\lambda$ = 488 nm,
LP488-SF20G, Thorlabs) was used to excite fluorescence. The excitation beam
was delivered to guidestars by a lens (L1, f = 100 mm) and an objective lens
(Plan N 20$\times$ 0.4, Olympus). To monitor the guidestars and phase
conjugation, a dichroic mirror (DMLP490R, Thorlabs), a lens (L2, f = 200 mm),
a bandpass filter (FL532-10, Thorlabs), and a camera (acA5472-17um, Basler)
were placed on the side without a scattering medium. The guidestars were
fluorescent ink mixed with UV glue (NOA 68, Norland), and fluorescent beads
(F8823, Invitrogen) immersed in glycerol. A scattering medium was a 220-grit
ground glass diffuser, placed approximately 170 $\mu$m away from the
guidestars. On the detection side, the scattered fluorescence was collected by
an objective lens (MPlan N 50$\times$ 0.75, Olympus) and two lenses (L3, f =
75 mm; L4, f = 150 mm). An SLM (X10468-04, Hamamtsu) and a linear polarizer
were used to modulate the fluorescence. The modulated fluorescence is Fourier
transformed by lenses (L5, f = 100 mm; L6, f = 200 mm; L7, f = 250 mm) and
then captured by an sCMOS camera (PCO.edge 5.5, PCO) with bandpass filters
(BP2, FL532-3 and FBH520-40, Thorlabs). For the fluorescent beads we used,
these filters transmit 5.35% of the total fluorescence signal. An iris was
placed between L6 and L7 to adjust the speckle grain size at the camera. For
phase conjugation, a laser ($\lambda$ = 532 nm, Compass 215M-50, Coherent) was
collimated using a 5 $\mu$m pinhole and a lens (L8, f = 6 mm). The collimated
laser was then shaped by the SLM to generate a phase-conjugated beam, which
propagated back to the medium. A flip mirror was used to switch between the
fluorescence detection and phase conjugation.

### Phase-conjugated beam generation

A phase-conjugated beam is generated using the collimated laser beam and the
SLM. The collimated beam is shaped to the phase conjugate of a given scattered
field $\bm{E}_{scattered}$ by displaying a phase pattern that corresponds to
$-\textnormal{arg}\left(\bm{E}_{scattered}\right)$ on the SLM. The resultant
phase-conjugated beam propagates back through the scattering medium, retracing
the scattering paths of fluorescence.

## Acknowledgement

This research was supported by H2020 Future and Emerging Technologies
(863203), European Research Council (724473), and Basic Science Research
Program through the National Research Foundation of Korea (NRF) funded by the
Ministry of Education (2022R1A6A3A03072108).

## References

  * [1] D. Aizik, I. Gkioulekas, and A. Levin (2022) Fluorescent wavefront shaping using incoherent iterative phase conjugation.  Optica 9 (7),  pp. 746–754.  Cited by: §I.
  * [2] N. Barré and A. Jesacher (2022) Holographic beam shaping of partially coherent light.  Optics Letters 47 (2),  pp. 425–428.  Cited by: §III.
  * [3] A. Boniface, B. Blochet, J. Dong, and S. Gigan (2019) Noninvasive light focusing in scattering media using speckle variance optimization.  Optica 6 (11),  pp. 1381–1385.  Cited by: §I.
  * [4] A. Boniface, J. Dong, and S. Gigan (2020) Non-invasive focusing and imaging in scattering media with a fluorescence-based transmission matrix.  Nature communications 11 (1),  pp. 1–7.  Cited by: §I, §III.
  * [5] E. S. Boyden, F. Zhang, E. Bamberg, G. Nagel, and K. Deisseroth (2005) Millisecond-timescale, genetically targeted optical control of neural activity.  Nature neuroscience 8 (9),  pp. 1263–1268.  Cited by: §I.
  * [6] E. J. Candes, Y. C. Eldar, T. Strohmer, and V. Voroninski (2015) Phase retrieval via matrix completion.  SIAM review 57 (2),  pp. 225–251.  Cited by: §III, §III.
  * [7] E. J. Candes, X. Li, and M. Soltanolkotabi (2015) Phase retrieval via wirtinger flow: theory and algorithms.  IEEE Transactions on Information Theory 61 (4),  pp. 1985–2007.  Cited by: §III.
  * [8] H. Cao, A. P. Mosk, and S. Rotter (2022) Shaping the propagation of light in complex media.  Nature Physics 18 (9),  pp. 994–1007.  Cited by: §I, §III.
  * [9] T. Čižmár, M. Mazilu, and K. Dholakia (2010) In situ wavefront correction and its application to micromanipulation.  Nature Photonics 4 (6),  pp. 388–394.  Cited by: §I.
  * [10] M. Cui and C. Yang (2010) Implementation of a digital optical phase conjugation system and its application to study the robustness of turbidity suppression by phase conjugation.  Optics express 18 (4),  pp. 3444–3455.  Cited by: §I.
  * [11] A. Daniel, D. Oron, and Y. Silberberg (2019) Light focusing through scattering media via linear fluorescence variance maximization, and its application for fluorescence imaging.  Optics Express 27 (15),  pp. 21778–21786.  Cited by: §I.
  * [12] P. De Santis, F. Gori, G. Guattari, and C. Palma (1986) Synthesis of partially coherent fields.  JOSA A 3 (8),  pp. 1258–1262.  Cited by: §III.
  * [13] S. Feng, C. Kane, P. A. Lee, and A. D. Stone (1988) Correlations and fluctuations of coherent wave transmission through disordered media.  Physical review letters 61 (7),  pp. 834.  Cited by: §II.
  * [14] I. Freund, M. Rosenbluh, and S. Feng (1988) Memory effects in propagation of optical waves through disordered media.  Physical review letters 61 (20),  pp. 2328.  Cited by: §II.
  * [15] J. W. Goodman (2007) Speckle phenomena in optics: theory and applications.  Roberts and Company Publishers.  Cited by: §III.
  * [16] T. R. Hillman, T. Yamauchi, W. Choi, R. R. Dasari, M. S. Feld, Y. Park, and Z. Yaqoob (2013) Digital optical phase conjugation for delivering two-dimensional images through turbid media.  Scientific reports 3 (1),  pp. 1–5.  Cited by: §III.
  * [17] R. Horstmeyer, H. Ruan, and C. Yang (2015) Guidestar-assisted wavefront-shaping methods for focusing light into biological tissue.  Nature photonics 9 (9),  pp. 563–571.  Cited by: §I.
  * [18] C. Hsieh, Y. Pu, R. Grange, and D. Psaltis (2010) Digital phase conjugation of second harmonic radiation emitted by nanoparticles in turbid media.  Optics express 18 (12),  pp. 12283–12290.  Cited by: §I.
  * [19] B. Judkewitz, Y. M. Wang, R. Horstmeyer, A. Mathy, and C. Yang (2013) Speckle-scale focusing in the diffusive regime with time reversal of variance-encoded light (trove).  Nature photonics 7 (4),  pp. 300–305.  Cited by: §I.
  * [20] O. Katz, E. Small, Y. Guan, and Y. Silberberg (2014) Noninvasive nonlinear focusing and imaging through strongly scattering turbid layers.  Optica 1 (3),  pp. 170–174.  Cited by: §I.
  * [21] D. Kim and D. R. Englund (2019) Quantum reference beacon–guided superresolution optical focusing in complex media.  Science 363 (6426),  pp. 528–531.  Cited by: §III.
  * [22] M. Kim, Y. Choi, C. Yoon, W. Choi, J. Kim, Q. Park, and W. Choi (2012) Maximal energy transport through disordered media with the implementation of transmission eigenchannels.  Nature photonics 6 (9),  pp. 581–585.  Cited by: §II.
  * [23] J. Kubby, S. Gigan, and M. Cui (2019) Wavefront shaping for biomedical imaging.  Cambridge University Press.  Cited by: §I.
  * [24] D. Li, S. K. Sahoo, H. Q. Lam, D. Wang, and C. Dang (2020) Non-invasive optical focusing inside strongly scattering media with linear fluorescence.  Applied Physics Letters 116 (24),  pp. 241104.  Cited by: §I.
  * [25] C. Ma, X. Xu, Y. Liu, and L. V. Wang (2014) Time-reversed adapted-perturbation (trap) optical focusing onto dynamic objects inside scattering media.  Nature photonics 8 (12),  pp. 931–936.  Cited by: §I.
  * [26] A. P. Mosk, A. Lagendijk, G. Lerosey, and M. Fink (2012) Controlling waves in space and time for imaging and focusing in complex media.  Nature photonics 6 (5),  pp. 283–292.  Cited by: §I, §III.
  * [27] V. Ntziachristos (2010) Going deeper than microscopy: the optical imaging frontier in biology.  Nature methods 7 (8),  pp. 603–614.  Cited by: §I.
  * [28] G. Osnabrugge, R. Horstmeyer, I. N. Papadopoulos, B. Judkewitz, and I. M. Vellekoop (2017) Generalized optical memory effect.  Optica 4 (8),  pp. 886–892.  Cited by: §III.
  * [29] S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, A. C. Boccara, and S. Gigan (2010) Measuring the transmission matrix in optics: an approach to the study and control of light propagation in disordered media.  Physical review letters 104 (10),  pp. 100601.  Cited by: §I.
  * [30] S. M. Popoff, A. Aubry, G. Lerosey, M. Fink, A. Boccara, and S. Gigan (2011) Exploiting the time-reversal operator for adaptive optics, selective focusing, and scattering pattern analysis.  Physical review letters 107 (26),  pp. 263901.  Cited by: §III.
  * [31] B. Rauer, H. B. de Aguiar, L. Bourdieu, and S. Gigan (2022) Scattering correcting wavefront shaping for three-photon microscopy.  Optics Letters 47 (23),  pp. 6233–6236.  Cited by: §I.
  * [32] B. Rodenburg, M. Mirhosseini, O. S. Magaña-Loaiza, and R. W. Boyd (2014) Experimental generation of an optical field with arbitrary spatial coherence properties.  JOSA B 31 (6),  pp. A51–A55.  Cited by: §III.
  * [33] S. Rotter and S. Gigan (2017) Light fields in complex media: mesoscopic scattering meets wave control.  Reviews of Modern Physics 89 (1),  pp. 015005.  Cited by: §I.
  * [34] H. Ruan, T. Haber, Y. Liu, J. Brake, J. Kim, J. M. Berlin, and C. Yang (2017) Focusing light inside scattering media with magnetic-particle-guided wavefront shaping.  Optica 4 (11),  pp. 1337–1343.  Cited by: §I.
  * [35] H. Ruan, M. Jang, and C. Yang (2015) Optical focusing inside scattering media with time-reversed ultrasound microbubble encoded light.  Nature communications 6 (1),  pp. 1–8.  Cited by: §I.
  * [36] P. Thibault and A. Menzel (2013) Reconstructing state mixtures from diffraction measurements.  Nature 494 (7435),  pp. 68–71.  Cited by: §II.
  * [37] J. V. Thompson, G. A. Throckmorton, B. H. Hokr, and V. V. Yakovlev (2016) Wavefront shaping enhanced raman scattering in a turbid medium.  Optics letters 41 (8),  pp. 1769–1772.  Cited by: §I, §III.
  * [38] B. Tian, B. Rauer, A. Boniface, J. Han, S. Gigan, and H. B. de Aguiar (2022) Non-invasive chemically selective energy delivery and focusing inside a scattering medium guided by raman scattering.  Optics Letters 47 (9),  pp. 2145–2148.  Cited by: §I, §III.
  * [39] O. Tzang, E. Niv, S. Singh, S. Labouesse, G. Myatt, and R. Piestun (2019) Wavefront shaping in complex media with a 350 khz modulator via a 1d-to-2d transform.  Nature Photonics 13 (11),  pp. 788–793.  Cited by: §III.
  * [40] I. M. Vellekoop and A. Mosk (2007) Focusing coherent light through opaque strongly scattering media.  Optics letters 32 (16),  pp. 2309–2311.  Cited by: §I.
  * [41] I. M. Vellekoop, M. Cui, and C. Yang (2012) Digital optical phase conjugation of fluorescence in turbid tissue.  Applied physics letters 101 (8),  pp. 081108.  Cited by: §I.
  * [42] I. M. Vellekoop and A. Mosk (2008) Universal optimal transmission of light through disordered materials.  Physical review letters 101 (12),  pp. 120601.  Cited by: §II.
  * [43] A. G. Vesga, M. Hofer, N. K. Balla, H. B. De Aguiar, M. Guillon, and S. Brasselet (2019) Focusing large spectral bandwidths through scattering media.  Optics express 27 (20),  pp. 28384–28394.  Cited by: §II.
  * [44] X. Xu, H. Liu, and L. V. Wang (2011) Time-reversed ultrasonically encoded optical focusing into scattering media.  Nature photonics 5 (3),  pp. 154–157.  Cited by: §I.
  * [45] J. Yang, L. Li, A. A. Shemetov, S. Lee, Y. Zhao, Y. Liu, Y. Shen, J. Li, Y. Oka, V. V. Verkhusha, et al. (2019) Focusing light inside live tissue using reversibly switchable bacterial phytochrome as a genetically encoded photochromic guide star.  Science advances 5 (12),  pp. eaay1211.  Cited by: §I.
  * [46] Z. Yaqoob, D. Psaltis, M. S. Feld, and C. Yang (2008) Optical phase conjugation for turbidity suppression in biological samples.  Nature photonics 2 (2),  pp. 110–115.  Cited by: §I.
  * [47] J. Yoon, M. Lee, K. Lee, N. Kim, J. M. Kim, J. Park, H. Yu, C. Choi, W. D. Heo, and Y. Park (2015) Optogenetic control of cell signaling pathway through scattering skull using wavefront shaping.  Scientific reports 5 (1),  pp. 1–7.  Cited by: §I.
  * [48] E. H. Zhou, H. Ruan, C. Yang, and B. Judkewitz (2014) Focusing on moving targets through scattering samples.  Optica 1 (4),  pp. 227–232.  Cited by: §I.

Supplementary Information:

##  1Matrix representation of I0 and Im

The intensity speckle measured at the camera plane can be written as

$I_{0}(\bm{k})=\sum_{n=1}^{N}{\left|\mathcal{F}[t_{n}\left(\bm{r}\right)]\right|^{2}}(\bm{k})=\sum_{n=1}^{N}{\left|\tilde{t}_{n}\left(\bm{k}\right)\right|^{2}},$
(S1)

where $t_{n}(\bm{r})$ and $\tilde{t}_{n}(\bm{k})$ is the $n$th scattered field
at the SLM plane and camera plane, respectively, and $\mathcal{F}$ is the
Fourier transform. To express Eq. S1 in a matrix form, we utilize $\bm{T}$ to
express $t_{n}(\bm{r})$ and $\tilde{t}_{n}(\bm{k})$. First, we can express
$t_{n}(\bm{r})$ as the rows of $\bm{T}$, such that $\bm{T}_{ij}=t_{i}(r_{j})$.
Similarly, $\tilde{t}_{n}(\bm{k})$ can be expressed as the rows of $\bm{TF}$,
where $\bm{F}$ is a Fourier transform matrix, such that
$(\bm{TF})_{ij}=\tilde{t}_{i}(k_{j})$. If we consider a matrix
$(\bm{TF})^{\dagger}\bm{TF}$, its elements can be expressed as
$[(\bm{TF})^{\dagger}\bm{TF}]_{ij}=\sum_{n=1}^{N}{\tilde{t}^{*}_{n}(k_{i})\tilde{t}_{n}(k_{j})}$.
Therefore, its main diagonal is given by
$[(\bm{TF})^{\dagger}\bm{TF}]_{ii}=\sum_{n=1}^{N}{\left|\tilde{t}_{n}(k_{i})\right|^{2}}$,
which is equivalent to $I_{0}(k_{i})$. Finally, we can rewrite Eq. S1 as
$\bm{I}_{0}=\textnormal{diag}{\left[(\bm{TF})^{\dagger}\bm{TF}\right]}$ by
expressing the main diagonal of a matrix using the function diag. Equation 1
can be obtained by following a similar process as described above, but by
replacing $\bm{T}$ with $\bm{T}\bm{S}_{m}.$

##  2Mixed-state phase retrieval

Our phase-retrieval algorithm is designed to find $N$ mutually incoherent
fields at the SLM plane from $M$ modulated intensity images. We initialize the
algorithm by letting the scattered fields $\bm{h}_{n}$ as $N$ complex Gaussian
random fields. Then we apply the wavefront modulation $\bm{S}_{m}$ to the
scattered fields $\bm{h}_{n}$. We note that $\bm{S}_{m}$ is a diagonal matrix
whose diagonal elements correspond to the field modulation given by the SLM.
In the experiments, SLM was divided into macro-pixels (composed of 40 $\times$
40 pixels) having random phase values (0 or $\pi$). After the wavefront
modulation $\bm{S}_{m}$, the field at the camera plane is expressed as
$\bm{\tilde{h}}_{n}^{(m)}=\left(\bm{S}_{m}\bm{F}\right)^{\top}\bm{h}_{n}$.
Then we conduct Fourier magnitude projection using an auxiliary function
$\psi$:

$\psi_{n}^{(m)}(\bm{k})=\frac{\left[I_{m}(\bm{k})\right]^{\gamma}}{\sqrt{\sum_{n}{\left|\tilde{h}_{n}^{(m)}(\bm{k})\right|^{2}}}}\tilde{h}_{n}^{(m)}(\bm{k}),$
(S2)

where $\bm{k}$ is a coordinate in the spatial frequency domain,
$\tilde{h}_{n}^{(m)}(\bm{k})$ is the modulated field at the camera plane,
corresponding to $\bm{\tilde{h}}_{n}^{(m)}$, and $\gamma$ is a constant
parameter. Then the fields are updated by compensating the wavefront
modulation:

$\bm{h}_{n}=\left[\left(\bm{S}_{m}\bm{F}\right)^{\top}\right]^{-1}{\bm{\psi}}_{n}^{(m)},$
(S3)

where $\bm{\psi}_{n}^{(m)}$ is the vector representation of
$\psi_{n}^{(m)}(\bm{k})$. The update through Eq. S2–S3 is continued for the
entire $M$ modulations. The whole process is repeated for several times to
obtain consistent $\bm{h}_{n}$ (see Fig. S5). With $\gamma=1/2$, this method
can be interpreted as the maximum likelihood reconstruction and retrieves the
incoherent fields [6]. However, we observed that few initial iterations with
$\gamma=1$ accelerates the convergence greatly, in both numerical simulations
and experiments. This observation agrees with the effect of the nonlinear
modulus constraint that was reported in [9, 8]. It should be noted when
$\gamma$ is greater than 1/2, the brighter area of $I_{m}$ carry more weight
during the iteration. While it has been reported that the larger values of
$\gamma$ lead to faster convergence [9], we expect that the optimal value of
$\gamma$ for the initial iterations would depend on the specific application.
In our experiments, we used $\gamma=1$ for the first 20 iterations and
$\gamma=1/2$ for the rest of the iterations. We stopped the algorithm when
there was no significant change in $\bm{h}_{n}$. Specifically, we evaluated
the convergence by using the mean squared error (MSE) between the current and
previous estimates of $\bm{h}_{n}$:

$\textnormal{MSE}=\frac{1}{NK}\sum_{n}^{N}\sum_{k}^{K}{\left|h_{nk}^{\textnormal{(current)}}-h_{nk}^{\textnormal{(previous)}}\right|}^{2},$
(S4)

where $h_{nk}$ is the $k$th element of $\bm{h}_{n}$. In our experiments, we
stopped the algorithm when the MSE was less than $10^{-5.6}$. See Fig. S6 for
a change of the MSE over iterations. We confirmed numerically that the
algorithm retrieves a correct set of fields for $M\geq 4N$ in the absence of
measurement noise. The size of macro-pixel had almost no effect on the
reconstruction, except when its size is comparable to the SLM. We observed
that the minimum value of $M$ required for correct reconstruction increases
depending on the noise level. For the experimental results shown in the main
text, we used 6–8$N$ modulations.

##  3Incoherent phase conjugation

When all the mutually incoherent components of fluorescence are time-reversed,
the intensity at the nth guidestar can be expressed by multiplying the nth row
of the transmission matrix, $\bm{t}_{n}^{\top}$, and phase-conjugated field
$\bm{t}_{m}^{*}$:

$\sum_{m}\left|{\bm{t}_{n}^{\top}\bm{t}_{m}^{*}}\right|^{2}=\bm{t}_{n}^{\top}\bm{T}^{\dagger}\bm{T}\bm{t}_{n}^{*},$
(S5)

where $\ast$ denotes the complex conjugate. For a medium with negligible
reflection and absorption, $\bm{T}^{\dagger}\bm{T}\approx\mathds{1}$. Thus Eq.
S5 is simplified to $\bm{t}_{n}^{\top}\bm{t}_{n}^{*}$, which remains more or
less constant regardless of n. As a result, the incoherent phase conjugation
of $\bm{t}_{n}$ generates foci on the entire guidestars with roughly the same
intensity. It should be noted that $\bm{T}$ is not unitary in practice, since
we are not able to measure the entire scattered light. Nevertheless, if the
loss of information due to absorption or measurement is minimal, the phase
conjugation can approximate the case of unitary $\bm{T}$ for complex media,
producing high-contrast foci on the entire guidestars [5, 2]. Similarly, the
intensity for the incoherent phase conjugation of $\bm{h}_{n}$ is expressed
as,

$\sum_{m}\left|{\bm{t}_{n}^{\top}\bm{h}_{m}^{*}}\right|^{2}=\bm{t}_{n}^{\top}\bm{H}^{\dagger}\bm{H}\bm{t}_{n}^{*}.$
(S6)

Since $\bm{H}^{\dagger}\bm{H}=\bm{T}^{\dagger}\bm{T}$ according to Eq. 2, Eq.
S6 is identical to Eq. S5. Thus, the incoherent phase conjugation of
$\bm{h}_{n}$ and $\bm{t}_{n}$ are identical.

##  4Additional experiment with volumetric scattering medium

In addition to the incoherent phase conjugation described in Fig. 2, we
conducted an additional experiment with a volumetric scattering medium. In the
experiment, we placed 1 $\mu$m fluorescent beads directly behind three layers
of parafilm, with a total thickness of approximately 360 $\mu$m corresponding
to half of the transport mean free path [1]. We replaced the objective lens
with a long-working-distance objective lens (MPlanFL N 100$\times$ 0.9,
Olympus) and used a 10-nm bandpass filter (FLH05532-10, Thorlabs) which
transmits 14.3% of the fluorescence signal. The rest of the setup was the
same. The experiment was conducted in the same manner as the experiment for
Fig. 2.

The results showed that that the proposed method was able to produce high-
contrast foci on 13 and 22 guidestars (Fig. S2). In Fig. S2d, we observed that
foci were not generated at a few target positions located far from the center.
This is due to the limited field-of-view (FoV) in which the SLM can control
light. By evaluating the speckle pattern generated by a back-propagating beam
with a random wavefront, we confirmed that the size of the controllable FoV
was approximately 60 $\mu$m (Fig. S2e). This is smaller than the field of view
of Fig. S2 (75.6 $\mu$m). We note that the FoV is determined by the scattering
media and the experimental setup, such as the focal length of an objective
lens and the SLM pixel size.

##  5Demixing incoherent fields

The scattered fields $\bm{t}_{n}$ can be recovered by reversing Eq. 2. To this
end, we developed an algorithm inspired by the simulated annealing [3] that
exploits the memory effect between the speckle patterns at the camera plane.
The algorithm starts by applying a random unitary transformation to the
retrieved fields $\bm{H}$ and expressing the result at the camera plane:

$\bm{\overline{H}}={\bm{U}_{j}}^{f_{j}}\bm{\hat{U}}\bm{H}\bm{F},$ (S7)

where $j$ is the iteration index, $\bm{U}_{j}$ is a random unitary matrix,
$f_{j}$ is an exponent that decays to 0 over the iteration, and $\bm{\hat{U}}$
is the current best estimate of $\bm{U}^{-1}$. The decaying behavior of
$f_{j}$ makes $\bm{U}_{j}$ induce smaller changes as the iteration progresses:
$\lim_{j\to\infty}{\bm{U}_{j}}^{f_{j}}=\mathds{1}$. Then the correlation
metric $C$ is calculated at every iteration,

$C=\sum_{n=1}^{N}\Bigl{[}\max\left({|\overline{h}_{n}|^{2}\star|\overline{h}_{(n\bmod
N)+1}|^{2}}\right)\Bigr{]}^{2},$ (S8)

where ${\overline{h}}_{n}(\bm{k})$ is the unitary-transformed speckle field at
the camera plane, corresponding to the $n$th row of $\bm{\overline{H}}$, and
$\star$ is the cross-correlation. If the value of $C$ is greater than the
previous values, $\bm{\hat{U}}$ is updated to
${\bm{U}_{j}}^{f_{j}}\bm{\hat{U}}$. If not, $\bm{\hat{U}}$ is unchanged. At
the end of the iteration, the algorithm results in
$\bm{\hat{U}}\approx\bm{U}^{-1}$, and the scattered fields are recovered by
inverting Eq. 2: $\bm{T}=\bm{\hat{U}}\bm{H}$. We initialized the iteration
with $\bm{\hat{U}}=\mathds{1}$ and $f_{j}=100/(j+1)$, and the maximum
iteration number of $10^{4}$. We note that the experimental noise in $\bm{H}$
can lead to an incorrect demixing. To mitigate this issue, we calculate $C$
only using speckles brighter than the average, $|\overline{h}_{n}|^{2}>\mu$,
where $\mu$ is the mean value of every $|\overline{h}_{n}|^{2}$.

##  6Guidestars of varying brightness

If the brightness of the guidestars is different, we can express the scattered
field corresponding to $n$th guidestar as $b_{n}\bm{t}_{n}$, where $b_{n}$ is
a constant that is the square-root of the energy emitted by the guidestar and
$\bm{t}_{n}$ is a normalized field, such that
$\bm{t}_{n}^{\dagger}\bm{t}_{n}=1$. Then we can rewrite the expressions in the
main text by replacing $\bm{T}$ with $\bm{BT}$, where
$\bm{B}=\textnormal{diag}(b_{1},...,b_{n})$. For example, the unmodulated
speckle image becomes
$\bm{I}_{0}=\textnormal{diag}{\left[(\bm{BTF})^{\dagger}\bm{BTF}\right]}$.
Similarly, Eq. 1 is rewritten as

$\bm{I}_{m}=\textnormal{diag}{\left[\left(\bm{BTS}_{m}\bm{F}\right)^{\dagger}\bm{BTS}_{m}\bm{F}\right]},$
(S9)

and Eq. 2 as

$\bm{H}=\bm{UBT}.$ (S10)

We note that this generalization is related to the interpretation of the
measured image and retrieved fields. Therefore, the mixed-state reconstruction
and the demixing process are unaffected by the brightness of the guidestars.
For incoherent phase conjugation (Supplementary Section 3), foci will be
generated on the entire guidestar with the intensities that correspond to the
brightness of the guidestars. Regarding maximum energy delivery, the
eigenchannel of $\bm{H}$ will preferentially couple light into the brighter
guidestar positions because
$\bm{H}^{\dagger}\bm{H}=\bm{T}^{\dagger}\bm{B}^{2}\bm{T}$, where
$\bm{B}^{2}=\textnormal{diag}(b_{1}^{2},...,b_{n}^{2})$.

##  7Maximum number of guidestars

In principle, the mixed-state algorithm can find solutions regardless of the
number of guidestars $N$, provided that there are enough wavefront
modulations. However, in practice, the number of guidestars could be limited
because the dynamic range of a camera causes quantization of images. For the
phase retrieval of coherent light, it is known that 3-bit representation of
speckle can produce reliable phase retrieval results [4]. Considering that the
mixed-state reconstruction can be thought of as the information multiplexing
of mutually incoherent signals, we would expect the following condition to be
satisfied,

$N\times 2^{3}\leq 2^{\textnormal{(Camera bit depth)}}.$ (S11)

For example, with a 16-bit resolution, the maximum theoretical $N$ will be
8192, although in practice it would be significantly lower. Finally, we note
that even a larger number could be allowed by incorporating the quantization
error into the reconstruction algorithm [7].

Figure S1: Phase conjugation using an excitation beam: a Fluorescence image of
a bead when an excitation beam ($\lambda=475\leavevmode\nobreak\ \mu m$) is
shaped using the phase conjugation pattern of scattered fluorescence
($\lambda=532\leavevmode\nobreak\ \mu m$). b Fluorescence image of the same
bead when the excitation beam is shaped by a random phase pattern. c Intensity
profiles of (a) and (b) along the horizontal line that crosses the center.
Figure S2: Incoherent phase conjugation in volumetric scattering medium: a–b
Fluorescence image of beads hidden behind 3 layers of parafilm. c–d Results of
the incoherent phase conjugation. e Intensity at the target plane with a
random wavefront. Figure S3: Phase conjugation with different numbers of
retrieved fields: a Fluorescence image of hidden targets. b Results of
incoherent phase conjugation when different numbers of scattered fields are
retrieved. c–e Phase conjugation of individual fields when two (c), three (d),
and four (e) scattered fields are retrieved. Figure S4: Experimental setup: L,
lens; Obj., objective lens; BS, beam splitter; BP, banspass filter. Figure S5:
Flowchart of the mixed-state phase retrieval Figure S6: Experimental mean
squared error

## References

  * [1] A. Boniface, B. Blochet, J. Dong, and S. Gigan (2019) Noninvasive light focusing in scattering media using speckle variance optimization.  Optica 6 (11),  pp. 1381–1385.  Cited by: §4.
  * [2] R. Horstmeyer, H. Ruan, and C. Yang (2015) Guidestar-assisted wavefront-shaping methods for focusing light into biological tissue.  Nature photonics 9 (9),  pp. 563–571.  Cited by: §3.
  * [3] S. Kirkpatrick, C. D. Gelatt Jr, and M. P. Vecchi (1983) Optimization by simulated annealing.  science 220 (4598),  pp. 671–680.  Cited by: §5.
  * [4] A. M. S. Maallo, P. F. Almoro, and S. G. Hanson (2010) Quantization analysis of speckle intensity measurements for phase retrieval.  Applied optics 49 (27),  pp. 5087.  Cited by: §7.
  * [5] M. Tanter, J. Thomas, and M. Fink (2000) Time reversal and the inverse filter.  The Journal of the Acoustical Society of America 108 (1),  pp. 223–234.  Cited by: §3.
  * [6] P. Thibault and A. Menzel (2013) Reconstructing state mixtures from diffraction measurements.  Nature 494 (7435),  pp. 68–71.  Cited by: §2.
  * [7] X. Yang, P. C. Konda, S. Xu, L. Bian, R. Horstmeyer, et al. (2021) Quantized fourier ptychography with binary images from spad cameras.  Photonics Research 9 (10),  pp. 1958–1969.  Cited by: §7.
  * [8] F. Zhang, G. Pedrini, and W. Osten (2007) Phase retrieval of arbitrary complex-valued fields through aperture-plane modulation.  Physical Review A 75 (4),  pp. 043805.  Cited by: §2.
  * [9] F. Zhang and J. Rodenburg (2010) Phase retrieval based on wave-front relay and modulation.  Physical Review B 82 (12),  pp. 121104.  Cited by: §2, §2.


# Non-Gaussian state generation certified using the EPR-steering inequality

E. S. Gómez<EMAIL_ADDRESS>G. Cañas  E. Acuña  Center for Optics and
Photonics, MSI-Nucleus on Advanced Optics, Departamento de Física, Universidad
de Concepción, 160-C, Concepción, Chile  W. A. T. Nogueira  Departamento de
Física, ICE, Universidade Federal de Juiz de Fora, Juiz de Fora, CEP
36036-330, Brazil  Universidade Federal de Minas Gerais, Caixa Postal 702,
30123-970, Belo Horizonte, MG, Brazil  G. Lima  Center for Optics and
Photonics, MSI-Nucleus on Advanced Optics, Departamento de Física, Universidad
de Concepción, 160-C, Concepción, Chile

###### Abstract

Due to practical reasons, experimental and theoretical continuous-variable
(CV) quantum information (QI) has been heavily based on Gaussian states.
Nevertheless, many CV-QI protocols require the use of non-Gaussian states and
operations. Here, we show that the Einstein-Podolsky-Rosen steering inequality
can be used to obtain a practical witness for the generation of pure bipartite
non-Gaussian states. While the scenario require pure states, we show its broad
relevance by reporting the experimental observation of the non-Gaussianity of
the CV two-photon state generated in the process of spontaneous parametric
down-conversion (SPDC). The observed non-Gaussianity is due only to the
intrinsic phase-matching conditions of SPDC.

###### pacs:

42.50.Dv

##  I Introduction

Continuous-variable (CV) quantum information (QI) is a research field that has
been increasingly growing in the past few years Braunstein . The need to cover
larger Hilbert spaces is motivated by QI protocols, such as quantum key
distribution (QKD), that has advantages when implemented in higher dimensions
Gisin1 ; Lima13 . The extension of many QI protocols, first proposed
considering discrete quantum systems, to the realm of CV systems has been
heavily based on Gaussian states Eisert1 ; Cirac ; Acin . This is due to the
fact that their covariance matrices are fully determined by the first and
second order moments Paris1 ; Weedbrock , and also because of the practicality
in the creation, manipulation and detection of Gaussian states Vaidman ; CV-
QKD .

Nevertheless, several recent works have shown the relevance of non-Gaussian
states and operations SteveNG . For instance, they are required for quantum
computation with cluster states Lund08_1 , entanglement distillation
Eisert02_1 ; Fiurasek02_1 , quantum error correction Fiurasek09_1  and
loophole-free Bell tests Banaszek99_1 ; Carmichael04_1 . Besides, they provide
advantages to the quantum teleportation, quantum cloning and state estimation
tasks Bonifacio03_1 ; Cerf05_1 ; Genoni09_1 ; Adesso09_1 .

In this work, we show that the Einstein-Podolsky-Rosen (EPR) steering
inequality EPR ; Schro ; Reid ; Steer  can be used to obtain a witness for the
non-Gaussianity of pure bipartite CV quantum states. While the scenario
require pure states, we show its broad relevance by reporting the observation
of the non-Gaussianity of the CV two-photon state generated in the process of
spontaneous parametric down-conversion (SPDC)Mandel ; Monken1 ; Saleh1 ; PR .
SPDC is up to date the most used source for experimental investigations in the
field of quantum information, and our work highlights the simplicity of using
this source for new applications in CV-QI. The generated down-converted
photons are correlated in their transverse momenta and can be used to test the
EPR-paradox BoydEPR . The observed non-Gaussianity is due only to the
intrinsic phase-matching conditions of the SPDC process TboNG ; ExterNF2 ,
thus, highlighting the simplicity of using SPDC sources for new applications
of CV-QI.

##  II The EPR-steering inequality with bipartite Gaussian states

Consider a bipartite system described by a pure CV-state
$\rho_{12}=|\psi\rangle\langle\psi|$, and a generic pair of complementary
noncommuting observables $\hat{u}_{i}$ and $\hat{v}_{i}$, with $i=1,2$ used to
denote the operation at each subsystem. The spectral decomposition of
$\hat{u}_{i}$ and $\hat{v}_{i}$ is a infinite set of continuous variables. In
Ref. Reid  it was introduced the EPR-steering criterion as
$\Delta_{inf}^{2}(\hat{u}_{2})\Delta_{inf}^{2}(\hat{v}_{2})\geq C$, where $C$
is is a value that depends on the chosen observables. The violation of this
inequality implies the implementation of a EPR-paradox BoydEPR ; Steve . The
_inferred_ variances are given by

$\Delta_{inf}^{2}(\hat{u}_{2})=\int du_{1}\,P(u_{1})\Delta^{2}(u_{2}|u_{1}),$
(1)

where $\Delta^{2}(u_{2}|u_{1})$ is the variance of the conditional probability
distribution $P(u_{2}|u_{1})$, and $P(u_{1})$ is the marginal probability
distribution of one party’s outcomes.

Now, let us consider a general continuous variable two mode pure state
described by a Gaussian amplitude
$\mathcal{A}(\mathbf{q}_{1},\mathbf{q}_{2})$, where the vectors
$\mathbf{q}_{1}$ and $\mathbf{q}_{2}$ are the outputs of the observable
$\hat{q}$ on each party

$\mathcal{A}_{G}(\mathbf{q}_{1},\mathbf{q}_{2})\propto\exp\left(-\frac{|\mathbf{q}_{1}+\mathbf{q}_{2}|^{2}}{4\sigma_{+}^{2}}\right)\exp\left(-\frac{|\mathbf{q}_{1}-\mathbf{q}_{2}|^{2}}{4\sigma_{-}^{2}}\right),$
(2)

with $\sigma_{+}$ and $\sigma_{-}$ being the widths of the corresponding
Gaussian functions. Several CV physical systems can be modelled in this way,
for example: an atomic ensemble interacting with an electromagnetic field, two
entangled photons sent through Gaussian channels, photoionization of atoms and
photodissociation of molecules Fedorov2004 , spontaneous emission of a photon
by an atom Eberly2002 ; Eberly2003 ; Fedorov2005 , and multiphoton pair
production Fedorov2006 ; Eisert . Since Gaussian amplitudes have the same
functional form in both transverse directions we can, without loss of
generality, work in one dimension and consider only its scalar form.

From Eq. (1) one can obtain the limit of $C$ for Eq. (2) Eisert . Note that
$\Delta_{inf,G}^{2}(\hat{q}_{2})=\Delta^{2}(q_{2}|q_{1})=\sigma_{+}^{2}\sigma_{-}^{2}/(\sigma_{+}^{2}+\sigma_{-}^{2})$.
One can also find the inferred variance of the complementary observable
$\hat{x}_{2}$, and it is given by
$\Delta_{inf,G}^{2}(\hat{x}_{2})=\Delta^{2}(x_{2}|x_{1})=1/(\sigma_{+}^{2}+\sigma_{-}^{2})$.
Therefore, the EPR-steering inequality reads

$\Delta^{2}(q_{2}|q_{1})\Delta^{2}(x_{2}|x_{1})\\!=\\!\frac{\sigma_{+}^{2}\sigma_{-}^{2}}{\left(\sigma_{+}^{2}+\sigma_{-}^{2}\right)^{2}}\geq\frac{1}{4}.$
(3)

Let us define one parameter $P\equiv\sigma_{+}/\sigma_{-}$ for simplicity. The
Schmidt Number $K_{G}$ for the general Gaussian state of Eq. (2) is given by
Eberly

$K_{G}=\frac{1}{4}\left(\frac{1}{P}+P\right)^{2}.$ (4)

Thus, the EPR-steering inequality can be written as

$\Delta^{2}(q_{2}|q_{1})\Delta^{2}(x_{2}|x_{1})=\frac{1}{4K_{G}}\geq\frac{1}{4}.$
(5)

Note that due to the symmetry of Eq. (2) one has that $K_{G}=K_{G\alpha}\times
K_{G\beta}$, where $K_{Gj}$ represents the Schmidt number in the transverse
direction $j=\alpha,\beta$ Kfact . Thus, one may further simplify the EPR-
steering inequality for Gaussian states to

$W\equiv\Delta^{2}(q_{2}|q_{1})\Delta^{2}(x_{2}|x_{1})=\frac{1}{4K_{G\alpha}^{2}}\geq\frac{1}{4}.$
(6)

Thus, one may test it by performing the measurements in only one transverse
direction. Moreover, it does not depend of the values chosen for $q_{1}$ and
$x_{1}$. When $P=1$, $K_{G}=K_{G_{\alpha}}=K_{G_{\beta}}=1$ and the Gaussian
state is a product state. In Fig. 1(a) [Fig. 1(b)] we show $K_{G_{\alpha}}$
($W$) with a solid blue line, while varying $P$.

Figure 1: (Color Online) (a) Shows the Schmidt number for bipartite Gaussian
states and for the CV spatial state of SPDC [Eq. (7)], while varying $P$ and
considering one transverse direction of $\mathbf{q}_{j}$. (b) The values of
$W$ plotted in terms of $P$.

##  III A witness for the non-Gaussianity of CV quantum states

From the results obtained above it is possible to envisage a simple and
practical way to determine if a certain bipartite pure state is Gaussian or
not. Note from Eq. (6) and Fig. 1(b) that pure bipartite Gaussian states will
always violate the EPR-steering inequality. The reason is that they are in
general entangled states [See Fig. 1(a)]. The only exception is the point
marked with the horizontal dashed line in Fig. 1(b), which represents the
point where the Gaussian state is a product state. In this case, we have that
$W=\frac{1}{4K_{G}}=0.25$, which corresponds to the upper quantum bound for
pure bipartite Gaussian states. Thus, the observation of a value greater than
0.25 for $W$ with pure entangled states can only be achieved while considering
non-Gaussian states.

Figure 2: (Color Online) Experimental setup. (a) State preparation stage. (b)
Setup configuration for measuring $\Delta^{2}(q_{2}|q_{1})$. (c) Setup
configuration for measuring $\Delta^{2}(x_{2}|x_{1})$. See the main text for
details.

While the scenario require pure states, now we show its broad relevance by
reporting the observation of the non-Gaussianity of the CV spatial two-photon
state generated in the process of SPDC Mandel ; Monken1 ; Saleh1 ; PR . When
perfect colinear phase-matching is considered and neglecting effects of
anysotropy, we can write the spatial two-photon state as Monken1 ; Saleh1

$|\psi\rangle_{12}\propto\iint
d\mathbf{q}_{1}d\mathbf{q}_{2}\,\tilde{E}_{p}(\mathbf{q}_{1}+\mathbf{q}_{2})\tilde{G}(\mathbf{q}_{1}-\mathbf{q}_{2})|1\mathbf{q}_{1}\rangle|1\mathbf{q}_{2}\rangle,$
(7)

where $|1\mathbf{q}_{i}\rangle$ represents one photon in mode $i$ ($i=1,2$)
usually called signal or idler, and with the transverse momentum $\mathbf{q}$.
$\tilde{E}_{p}(\mathbf{q})$ is the angular spectrum of the pump beam. Usual
experimental configurations adopt a Gaussian pump beam and in this case
$\tilde{E}_{p}(\mathbf{q})\propto\exp\left[-c^{2}|\mathbf{q}|^{2}/4\right]$.
$c$ represents the beam radius at the crystal plane.
$\tilde{G}(\mathbf{q})=\mbox{\hskip 1.3ptsinc}\left(b|\mathbf{q}|^{2}\right)$
defines the phase-matching conditions of the SPDC process, with $\mbox{\hskip
1.3ptsinc}(\xi)\equiv\frac{\sin(\xi)}{\xi}$. $b$ is defined by
$b\equiv\frac{L}{8k}$, where $L$ is the crystal length, and $k$ the wavenumber
of the down-converted photons. In terms of these definitions, $P$ reads
$\frac{1}{c}\sqrt{\frac{L}{2k}}$. This state can be rewritten in the
complementary transverse position representation as ExterNF2 ; Monken2

$|\psi\rangle_{12}\propto\iint
d\bm{x}_{1}d\bm{x}_{2}\,E_{p}\left(\frac{\bm{x}_{1}+\bm{x}_{2}}{2}\right)G\left(\frac{\bm{x}_{1}-\bm{x}_{2}}{2}\right)|1\bm{x}_{1}\rangle|1\bm{x}_{2}\rangle,$
(8)

where the functions $E_{p}(\bm{x})$ and $G(\bm{x})$ are the Fourier transform
of $\tilde{E}_{p}(\mathbf{q})$ and $\tilde{G}(\mathbf{q})$. Thus,
$E_{p}(\bm{x})\propto\exp\left[-|\bm{x}|^{2}/c^{2}\right]$ and
$G(\bm{x})\propto 1-\frac{2}{\pi}\mbox{\hskip
1.3ptSi}\left(\frac{1}{4b}|\bm{x}|^{2}\right)\equiv\mathrm{sint}\left(\frac{1}{4b}|\bm{x}|^{2}\right)$,
where $\mbox{\hskip 1.3ptSi}(x)\equiv\int_{0}^{x}dt\,\mbox{\hskip
1.3ptsinc}(t)$. Here, the transverse position $\hat{x}_{j}$ and transverse
momentum $\hat{q}_{j}$ are the complementary observables for the EPR-steering
inequality test BoydEPR . Clearly, the CV spatial two-photon state of SPDC is
a non-Gaussian state TboNG .

An experimental observation of the effects that arise form the phase-matching
conditions, namely, the non-Gaussianity of $|\psi\rangle_{12}$, was reported
in Ref. ExterNF2 . They demonstrated how the spatial correlations in the near
field plane of a non-linear crystal changes when the phase matching conditions
varies. Now, we demonstrate experimentally how $W$ can be used to detect the
non-Gaussianity of this state. The demonstration is based on the fact that the
spatial state of SPDC process is pure and entangled Eberly , even considering
a post-selected one transverse direction TboNG ; PR . This means that for any
value of $P$, the Schmidt number $K_{S\alpha}$ is always greater than 1. In
Fig. 1(a) [Fig. 1(b)] we show $K_{S_{\alpha}}$ ($W$) with a dashed red line,
while varying $P$ for the state $|\psi\rangle_{12}$. One can see that for some
values of $P$ ($0.56\leq P\leq 2.58$) the values of $W$ are greater than
$\frac{1}{4}$, thus, witnessing the non-Gaussianity of this state.

##  IV Experiment

The experimental setup is illustrated in Fig. 2. We used a solid-state laser
source, at $355\,nm$, to pump a $\beta$-barium-borate type-II (BBO-II) non-
linear crystal ($L=1.8\,cm$) for the generation of the down-converted photons.
Initially, the Gaussian pump beam had a waist of $c=200\,\mu m$ at the crystal
plane. However, to experimentally observe the dependence of $W$ with $P$ [See
Fig. 1(b)], we generated five more different states by changing the waist of
the beam at the crystal plane. This has been done by using a configurable set
of doublet achromatic lenses placed before the non-linear crystal [See Fig.
2(a)]. The corresponding values of $P$ for the six generated states are shown
in Tab. 1.

Table 1: Corresponding values of $P$ for the generated states. State | $c\,[\mu m]$ | $P$
---|---|---
1 | 200 | 0.1595
2 | 100 | 0.3189
3 | 70 | 0.4556
4 | 45 | 0.7087
5 | 40 | 0.7973
6 | 35 | 0.9112

To guarantee the purity of the measured states, we used spatial and spectral
filters in each measurement apparatus. For instance, interference filters were
used to select degenerated down-converted photons at $710\,nm$ with $5\,nm$ of
bandwidth. This fact implies no entanglement between the frequencies and the
transverse spatial coordinates, and then the reduced spatial state must be
pure. Besides, the usage of spatial filters is explained afterwards.

Furthermore, a polarizer beam splitter (PBS) separates the signal and idler
photon modes. For measuring $\Delta^{2}(q_{2}|q_{1})$ and
$\Delta^{2}(x_{2}|x_{1})$, we performed conditional coincidence measurements
between the idler and signal photons in two different transverse planes (See
Fig. 2): the far- and near-field planes of the non-linear crystal,
respectively BoydEPR .

###  IV.1 Measurement of $\Delta^{2}(q_{2}|q_{1})$

In order to compute $\Delta^{2}(q_{2}|q_{1})$, one shall perform conditional
coincidence measurements at the far-field plane of the non-linear crystal. The
reason is very simple: at this plane, the coincidence rate
$C_{q}(x_{1},x_{2})$ is given by

$C_{q}(x_{1},x_{2})\propto\left|\tilde{E}_{p}\left[\frac{k}{f_{q}}(x_{1}+x_{2})\right]\tilde{G}\left[\frac{k}{f_{q}}(x_{1}-x_{2})\right]\right|^{2},$
(9)

and, since $q_{j}=kx_{j}/f_{q}$, it maps the square of the probability
amplitude of Eq. (7), i.e., the momentum correlation of the photons generated
in the SPDC process.

Figure 2(b) shows our setup configuration. One lens $L_{q}$ with focal
distance $f_{q}=15\,cm$ was placed before the PBS to create the far-field
plane for both signal and idler beams. At these planes, vertical slits were
placed for post-selecting the desirable state LeoGen ; LVNRS09  and performing
the conditional coincidence measurements SteveNG ; BoydEPR ; Yuan . The width
of each slit is $50\,\mu m$. The effect of using a non point-like detector is
that the transmittance function of the slit may broad the far-field
distribution to be measured. However, in the case of 50 $\mu m$ slits and the
experimental configuration adopted, this effect is negligible Kfact . This can
be easily checked through the calculation of the convolution between the
transmittance function of the slit and the predicted far-field distribution.

After the transmission through each slit, the down-converted photons were
collected with a $10$x objective lens and multi-mode fibers. The fibers were
connected to single-photon counting modules, and then a coincidence circuit
(with $4\,ns$ of coincidence window) recorded the data. To perform the
coincidence conditional measurements of $\Delta^{2}(q_{2}|q_{1})$, we scanned
in the horizontal direction one slit (of mode 2) while the other one was fixed
at the center ($q_{1}=0$).

To give an example of the results obtained while scanning the slit at mode-2,
we show in Fig. 3(a) [Fig. 3(b)] the far-field conditional distribution
measured for the second (sixth) state generated. The experimental results are
represented by red points (error bars lie inside the points due to the
observed high rate of coincidence counts) and the black dotted-line is the
theoretical curve for these distributions arising from Eq. (9).

Figure 3: (Color Online) Experimental measurement of the far- and near-field
conditional distributions. In (a) [(c)] we show the results for the second
state at the far-field (near-field) plane. In (b) [(d)] we show the results
for the sixth state at the far-field (near-field) plane. The (red) dots
represent the experimental data and the dotted (black) lines are the
theoretical predictions from Eq. (9) and Eq. (10).

###  IV.2 Measurement of $\Delta^{2}(x_{2}|x_{1})$

To measure $\Delta^{2}(x_{2}|x_{1})$ it is now necessary to measure at the
near-field plane of the non-linear crystal. Again, the reason is very simple:
at this plane, the coincidence distribution is proportional to the square of
the amplitude of Eq. (8), that is

$C_{x}(x_{1},x_{2})\propto\left|E_{p}\left(\frac{x_{1}+x_{2}}{2}\right)G\left(\frac{x_{1}-x_{2}}{2}\right)\right|^{2}.$
(10)

For performing the conditional coincidence measurements at the near-field
plane we used 4x objective lenses to form the image of the center of the BBO-
II onto the transverse plane of the fiber-couplers [see Fig. (2)(c)]. For
doing this, we removed the slits used to scan the coincidence rate in the far-
field plane and the lens $L_{q}$. The multi-mode fibers were replaced with
single-mode fibers whose core diameter were $4.7\,\mu m$. The small size of
the fibers core allows for post-selecting and measuring with high accuracy
$\Delta^{2}(x_{2}|x_{1})$ ExterNF2 ; Eisert . Again, the effect of the
transmittance function of the fiber over the broadening of the near-field
distribution is negligible Kfact . This can be easily checked through the
calculation of the convolution between the transmittance function of the fiber
with the predicted near-field distribution for our experimental
configurations. By scanning transversely the single-mode fiber at mode-2, we
recorded the coincidence conditional distribution at the near-field plane. An
example of the results obtained is shown in Fig. 3(c) [Fig. 3(d)] for the
second (sixth) state generated. The experimental results are represented by
red points. The black dotted-line is the theoretical curve for these
distributions arising from Eq. (8).

We have measured the coincidence conditional distribution of the down-
converted photons by imaging the center of the non-linear crystal to the
transverse plane of the detection system. However, as it has been shown in
Ref. ExterNF2 , the conditional distribution at the near-field plane depends
strongly of which part of the crystal is imaged at the detection system. This
is specially relevant when thicker crystals are considered. In our case, we
have checked that the near-field distribution does not change significantly
while considering different planes of our thin crystal to be imaged at the
detection system. For doing this, we moved the crystal around the longitudinal
central position $z_{0}=0$, imaging 5 different planes $z_{c}$ of the crystal
at the detection plane. For each $z_{c}$, we recorded the conditional
distribution at the near-field plane. The experimental results are shown in
Fig. (4). One can observe that our results are in agreement with the
theoretical prediction, which takes into account our crystal length and our
imaging system. A longitudinal crystal displacement around its center
introduces a phase factor of $\exp[-i|\mathbf{q}_{j}^{2}|z/(2k)]$ onto Eq.
(7). From our results, one can see that there is only a slight narrowing of
the near-field conditional distribution such that this effect does not affect
significantly our test of $W$.

Figure 4: (Color Online) (a) Experimental result and (b) theoretical
prediction for the near-field conditional distribution while moving
longitudinally the non-linear crystal (See the main text for details).

###  IV.3 Testing $W$

Table 2: Results of the conditional variances measured at the near- and far-field planes. $P$ | $\Delta^{2}(x_{2}|0)_{E}\,\left[m^{2}\right]$ | $\Delta^{2}(q_{2}|0)_{E}\,\left[\frac{1}{m^{2}}\right]$ | $W_{\mathrm{E}}$ | $W_{\mathrm{T}}$
---|---|---|---|---
$0.1595$ | $(6.62\pm 0.26)\times 10^{-10}$ | $(3.55\pm 0.14)\times 10^{7}$ | $0.024\pm 0.002$ | $0.033$
$0.3189$ | $(9.17\pm 0.37)\times 10^{-10}$ | $(1.25\pm 0.05)\times 10^{8}$ | $0.115\pm 0.009$ | $0.11$
$0.4556$ | $(7.76\pm 0.31)\times 10^{-10}$ | $(2.60\pm 0.1)\times 10^{8}$ | $0.20\pm 0.02$ | $0.19$
$0.7087$ | $(7.75\pm 0.31)\times 10^{-10}$ | $(4.67\pm 0.19)\times 10^{8}$ | $0.36\pm 0.03$ | $0.34$
$0.7973$ | $(7.23\pm 0.29)\times 10^{-10}$ | $(5.28\pm 0.21)\times 10^{8}$ | $0.38\pm 0.03$ | $0.38$
$0.9112$ | $(8.16\pm 0.33)\times 10^{-10}$ | $(5.18\pm 0.21)\times 10^{8}$ | $0.42\pm 0.03$ | $0.43$

In order to test the $W$ witness, and then certify the non-gaussian feature of
the CV spatial state of SPDC, we compute the variances of the conditional
coincidence measurements at the near- and far-field planes for the six
generated states. The results of all the conditional variances measured are
shown in Tab. 2. The errors of the variances were obtained by minimizing the
squared two-norm of the residuals between the analytical and experimental
results (See Fig. 3).

Figure 5: (Color Online) Measurement of $W$ for the six generated states. The
continuous (black) line corresponds to the theoretical value of the inequality
for the spatial SPDC state. The dashed (blue) line shows the Gaussian bound.

Figure 5 shows the values of the $W$ values for each state generated. The
continuous (black) line indicate the predicted theoretical value of the
inequality when we consider the spatial state of SPDC. The dotted (blue) line
is the upper limit for this inequality for Gaussian states. For the first
three entangled states, the values obtained for this inequality are below the
Gaussian limit as predicted by theory. However, for the last three entangled
states, there is a clear experimental violation of this bound. Due to momentum
conservation, the spatial state of the SPDC process is always entangled Eberly
[See Fig. 1(a)] and, thus, we have a clear experimental demonstration of the
non-Gaussianity of the CV two-photon spatial state of SPDC.

##  V Conclusion

We have introduced a novel application for the EPR-steering inequality by
showing that it can be used for witnessing the non-Gaussianity of CV quantum
states. To demonstrated this we performed an experiment using the CV spatial
state of entangled down-converted photons. Due to the phase-matching
conditions of the SPDC process, the generated is state is naturally a pure
entangled non-Gaussian state. A clear violation of the Gaussian bound of the
EPR-steering inequality has been observed. Since non-Gaussian states are
required for many new protocols of CV-QI, our work highlights the simplicity
and relevance of using SPDC sources for new applications in CV quantum
information processing.

### Acknowledgments

We thank C. H. Monken for discussions of this paper. This work was supported
by Grants FONDECYT 1120067, Milenio P10-030-F and CONICYT FB0824/2008. E. S.
G. and G. C. acknowledge the financial support of CONICYT. W. A. T. N. thanks
CNPq (Brazil) for financial support.

## References

  * (1) S. Braunstein and P. van Lock, Rev. Mod. Phys. 77, 513 (2005).
  * (2) N. J. Cerf, M. Bourennane, A. Karlsson, and N. Gisin, Phys. Rev. Lett. 88, 127902 (2002).
  * (3) S. Etcheverry, G. Cañas, E. S. Gómez, W. A. T. Nogueira, C. Saavedra, G. B. Xavier, and G. Lima, Sci. Rep. 3, 2316 (2013).
  * (4) J. Eisert, and M. B. Plenio, Int. J. Quant. Inf. 1, 479-506 (2003).
  * (5) M. M. Wolf, G. Giedke, and J. I. Cirac, Phys. Rev. Lett. 96, 080502 (2006).
  * (6) M. Navascues, F. Grosshans, and A. Acin, Phys. Rev. Lett. 97, 190502 (2006).
  * (7) D. Buono, G. Nocerino, V. D’Auria, A. Porzio, S. Olivares, and M. G. A. Paris, J. Opt. Soc. Am. B 27, 110 (2010) .
  * (8) C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, Rev. Mod. Phys. 84, 621 (2012).
  * (9) L. Vaidman, Phys. Rev. A 49, 1473 (1994).
  * (10) L. S. Madsen, V. C. Usenko, M. Lassen, R. Filip, and U. L. Andersen, Nat. Comm. 3, 1083 (2012).
  * (11) R. M. Gomes, A. Salles, F. Toscano, P. H. S. Ribeiro, and S. P. Walborn, Pro. Nat. Aca. Sci. 106, 21517 (2009).
  * (12) A. P. Lund, T. C. Ralph, and H. L. Haselgrove, Phys. Rev. Lett. 100, 030503 (2008).
  * (13) J. Eisert, S. Scheel, and M. B. Plenio, Phys. Rev. Lett. 89, 137903 (2002).
  * (14) J. Fiurasek, Phys. Rev. Lett. 89, 137904 (2002).
  * (15) J. Niset, J. Fiurasek, and N. J. Cerf, Phys. Rev. Lett. 102, 120501 (2009).
  * (16) K. Banaszek and K. Wódkiewicz, Phys. Rev. Lett. 82, 2009 (1999).
  * (17) H. Nha and H. J. Carmichael, Phys. Rev. Lett. 93, 020401 (2004).
  * (18) C. Invernizzi, S. Olivares, M. G. A. Paris, and K. Banaszek, Phys. Rev. A 72, 042105 (2005).
  * (19) S. Olivares, M. G. A. Paris, and R. Bonifacio, Phys. Rev. A 67, 032314 (2003).
  * (20) N. J. Cerf, O. Kruger, P. Navez, R. F. Werner, and M. M. Wolf, Phys. Rev. Lett. 95, 070501 (2005).
  * (21) M. G. Genoni, C. Invernizzi, and M. G. A. Paris, Phys. Rev. A , 80, 033842 (2009).
  * (22) G. Adesso, F. Dell’Anno, S. De Siena, F. Illuminati, and L. A. M. Souza, Phys. Rev. A 79, 040305 (2009).
  * (23) A. Einstein, B. Podolsky and N. Rosen, Phys. Rev. 47, 777 (1935).
  * (24) E. Schrödinger, Naturwiss 23, 807 (1935).
  * (25) M. D. Reid, Phys. Rev. A 40, 913 (1989).
  * (26) H. M. Wiseman, S. J. Jones, and A. C. Doherty, Phys. Rev. Lett. 98, 140402 (2007).
  * (27) C. K. Hong and L. Mandel, Phys. Rev. A 31, 2409 (1985).
  * (28) C. H. Monken, P. H. S. Ribeiro and S. Pádua, Phys. Rev. A 57, 3123 (1998).
  * (29) B. E. A. Saleh, A. F. Abouraddy, A. V. Sergienko, and M. C. Teich, Phys. Rev. A 62, 043816 (2004).
  * (30) S. P. Walborn, C. H. Monken, S. Pádua, and P. H. S. Ribeiro, Physics Reports 495, 87 (2010).
  * (31) J. C. Howell, R. S. Bennink, S. J. Bentley, and R. W. Boyd, Phys. Rev. Lett. 92, 210403 (2004).
  * (32) E. S. Gómez, W. A. T. Nogueira, C. H. Monken, and G. Lima, Opt. Express 20, 3753 (2012).
  * (33) H. D.-L. Pires and M. P. van Exter, Phys. Rev. A 80, 053820 (2009).
  * (34) S. P. Walborn, A. Salles, R. M. Gomes, F. Toscano, ans P. H. Souto Ribeiro, Phys. Rev. Lett. 106, 130402 (2011).
  * (35) M. V. Fedorov, M. A. Efremov, A. E. Kazakov, K. W. Chan, C. K. Law, and J. H. Eberly, Phys. Rev. A 69, 052117 (2004).
  * (36) K. W. Chan, C. K. Law, and J. H. Eberly, Phys. Rev. Lett. 88, 100402 (2002).
  * (37) K. W. Chan, C. K. Law, and J. H. Eberly, Phys. Rev. A 68, 022110 (2003).
  * (38) M. V. Fedorov, M. A. Efremov, A. E. Kazakov, K. W. Chan, C. K. Law, and J. H. Eberly, Phys. Rev. A 72, 032110 (2005).
  * (39) M. V. Fedorov, M. A. Efremov, and P. A. Volkov, Opt. Comm. 264, pp. 413-418 (2006).
  * (40) M. Ostermeyer, D. Korn, D. Puhlmann, C. Henkel, and J. Eisert, J. Mod. Opt. 56, 1829-1837 (2009).
  * (41) C. K. Law and J. H. Eberly, Phys. Rev. Lett. 92, 127903 (2004).
  * (42) M. P. van Exter, A. Aiello, S. S. R. Oemrawsingh, G. Nienhuis, and J. P. Woerdman, Phys. Rev. A 74, 012309 (2006).
  * (43) H. D.-L. Pires, C. H. Monken, and M. P. van Exter, Phys. Rev. A 80, 022307 (2009).
  * (44) L. Neves, G. Lima, J. G. Aguirre Gómez, C. H. Monken, C. Saavedra, and S. Pádua, Phys. Rev. Lett. 94, 100501 (2005).
  * (45) G. Lima, A. Vargas, L. Neves, R. Guzmán, and C. Saavedra, Opt. Express 17, 10688 (2009).
  * (46) Wei-Tao Liu, Ping-Xing Chen, Cheng-Zu Li, and Jian-Min Yuan, Phys. Rev. A 79, 061802(R) (2009).