Maximum information states for coherent scattering measurements
MMaximum information states for coherent scattering measurements
Dorian Bouchet, ∗ Stefan Rotter, and Allard P. Mosk Debye Institute for Nanomaterials Science, Utrecht University,P.O. Box 80000, 3508 TA Utrecht, the Netherlands Institute for Theoretical Physics, Vienna University of Technology (TU Wien), A-1040 Vienna, Austria
The use of coherent light for precision measurements has been a key driving force for numerousresearch directions, ranging from biomedical optics [1, 2] to semiconductor manufacturing [3]. Re-cent work demonstrates that the precision of such measurements can be significantly improved bytailoring the spatial profile of light fields used for estimating an observable system parameter [4–10]. These advances naturally raise the intriguing question of which states of light can provide theultimate measurement precision [11]. Here, we introduce a general approach to determine the opti-mal coherent states of light for estimating any given parameter, regardless of the complexity of thesystem. Our analysis reveals that the light fields delivering the ultimate measurement precision areeigenstates of a Hermitian operator which quantifies the Fisher information based on the system’sscattering matrix [12, 13]. To illustrate this concept, we experimentally show that these maximuminformation states can probe the phase or the position of an object that is hidden by a disorderedmedium with a precision improved by an order of magnitude as compared to unoptimized states.Our results enable optimally precise measurements in arbitrarily complex systems, thus establishinga new benchmark for metrology and imaging applications [3, 14].
No physical observable can be determined with ab-solute certainty. Instead, the noise inherent in anymeasurement process sets a fundamental limit on theprecision that a physical observable can be estimatedwith [11, 15]. Whenever light or other kinds of elec-tromagnetic radiation are involved in a measurement, anecessary condition to reach this ultimate precision is tooptimize the spatial distribution of the radiation field inthe measured system [14]. To achieve this goal, a cru-cial task is to identify the spatial pattern that shouldbe imprinted on the incoming field in order to get themaximum information out of it. First progress in this di-rection has recently been made using wavefront shapingtechniques and metasurfaces to precisely estimate lateraldisplacements [4, 8], fluorophore positions [5, 6], spectralshifts [7] or phase variations [9].A central challenge that remains unresolved, however,is to identify a unifying approach to reach the ultimateprecision limit that is applicable even to complex scat-tering systems. Earlier work suggests that such an ap-proach should be connected to the concept of Fisherinformation [5, 10, 11, 15], which quantifies the amount ofinformation relevant to the estimation of a given parame-ter from measured data. However, for the generic case ofa complex medium, the Fisher information is intrinsicallylinked to the microstructure of the medium [10], whichis not only overwhelmingly complex in realistic systemsbut also typically unknown.Here, we overcome this difficulty by expressing theFisher information in terms of a Hermitian operator thatdepends on the system’s optical scattering matrix. Basedon this idea, we introduce and experimentally demon- ∗ [email protected] strate a direct approach to generate optimal coherentstates of light for parameter estimation, regardless of thecomplexity of the system. Such light states are shownto be specifically tailored not only with respect to thespecific observable of interest, but also with respect tothe position of the observer. By unambiguously iden-tifying these optimal light states, we establish a newgeneral benchmark for metrology and imaging applica-tions [3, 14]. Furthermore, in the ideal case for which alloptical modes supported by the system are accessible tothe observer, our analysis reveals that maximum informa-tion states are, at the same time, the optimal states foroptical micro-manipulation [16, 17], thereby uncoveringa fundamental relationship between information theoryand measurement backaction.To set up this approach we recall that a measurementscheme is optimal when the measurements, the estima-tion function, and the choice of the incident state are alloptimal concurrently [11]. To realize this situation forcoherent states of light, we start with a general modelof scattering measurements on a complex medium pa-rameterized by a scalar parameter θ (Fig. 1). This pa-rameter can be a global parameter characterizing the en-tire scattering medium. It can also be a local parameterof limited spatial extent such as the phase or the posi-tion of a small phase object hidden behind a scatteringmaterial, as in our experiments. The medium is illumi-nated from the far field by an incident coherent state | E in (cid:105) , which can be spatially modulated using wavefrontshaping techniques [18]. Quadrature components of theresulting outgoing state | E out (cid:105) are then measured in thefar field at N independent sampling points using a homo-dyne detection scheme. This scheme is known to operateat the optimal noise limit: for a strong reference beam,the variance of the measured quadrature components of a r X i v : . [ phy s i c s . op ti c s ] J u l Fig. 1 | Principle of an optimal coherent scatteringmeasurement.
A scattering medium is characterized by anunknown parameter θ . This parameter is estimated by illu-minating the medium with coherent light and by measuringthe outgoing field state via a homodyne detection scheme.While plane-wave illumination may lead to imprecise estima-tions (top panel), the optimal incident state can be generatedusing wavefront shaping techniques to perform optimal esti-mations (bottom panel). the field is limited by Heisenberg’s uncertainty principle.Noise fluctuations in the measured data fundamentallylimit the achievable precision on the determination of θ .This limit is mathematically expressed by the Cramér-Rao inequality, which sets a lower bound on the varianceof unbiased estimators of θ . In general, the Cramér-Raobound is given by the reciprocal of the Fisher information J ( θ ) = E ([ ∂ θ ln p ( X ; θ )] ) where X is a N -dimensionalrandom variable representing the data, p ( X ; θ ) is a jointprobability density function and E denotes the expecta-tion operator acting over noise fluctuations [15]. In thecase of scattering measurements performed with a detec-tion scheme optimized for the estimation of θ , a strongreference beam must be shaped so that its phase at the k -th sampling point is φ k = arg( ∂ θ E out k ) , where E out k de-notes the outgoing field evaluated at the k -th samplingpoint (see Supplementary Information S1.1). Then, thestatistical independence of the samples allows us to de-rive the following simplified expression for the Fisher in-formation: J ( θ ) = 4 N (cid:88) k =1 (cid:12)(cid:12) ∂ θ E out k (cid:12)(cid:12) . (1)We can compare this expression of the classical Fisher in-formation to its quantum counterpart, which sets a lowerbound on the variance of unbiased estimators withoutany specification on the detection scheme [19]. Calculat- ing the quantum Fisher information for coherent statesand assuming statistical independence of the samples, theclassical expression given by Eq. (1) is shown to reachthe quantum limit (see Supplementary Information S1.2),thereby ensuring that the homodyne detection schemeconsidered here is optimal for the estimation of θ .We now investigate the optimization of the incidentstate. An incident state | E in (cid:105) is considered as optimalfor the Fisher information it generates when the corre-sponding outgoing state | E out (cid:105) maximizes Eq. (1) for agiven number of incident photons. In order to identifythis state, the pivotal quantity is the scattering matrix S of the medium [12, 13], which relates incoming andoutgoing asymptotic states via | E out (cid:105) = S | E in (cid:105) . Usingthis relationship, the Fisher information can be written inthe quadratic form J ( θ ) = 4 (cid:104) E in | F θ | E in (cid:105) , where we usedbra-ket notations for the complex inner product. We des-ignate the term F θ in the center of this expression as the Fisher information operator that takes on the remarkablysimple form (see Supplementary Information S1.3) F θ = ( ∂ θ S ) † ∂ θ S , (2)where † stands for the conjugate transpose. Since this op-erator F θ is Hermitian already by construction, the inci-dent state that maximizes the Fisher information is givenby the eigenstate associated with the largest eigenvalueof F θ . Furthermore, in the limit of small parameter vari-ations, we obtained a closed-form expression of the min-imum variance unbiased estimator (see SupplementaryInformation S1.4), which is the optimal estimation func-tion. Importantly, the variance of this estimator alwaysreaches the Cramér-Rao bound, which confirms that theFisher information is a relevant quantity to compare theestimation precision achievable with different light states.In our experiments, we first choose, as the observ-able parameter θ that we aim to estimate, the phaseshift ϕ generated by a small cross (total length µm)displayed by a spatial light modulator (SLM), as rep-resented in Fig. 2a (see also Supplementary Informa-tion S2). This target is hidden . mm behind a groundglass diffuser (scattering angle ◦ ). We illuminate thisscattering medium using a continuous-wave laser (wave-length nm) via a × objective (numerical aperture . ). The phase of the incident field is controlled using anSLM located in the conjugate plane of the surface of thediffuser. The horizontal component of the field is mea-sured in reflection using an off-axis reference arm and acharge-coupled device (CCD) camera, which images thesurface of the diffuser. Both quadrature components ofthe complex field are measured in a single shot using dig-ital off-axis holography. Using a tilted beam instead ofthe optimally-shaped reference field does not impact theshape of the incident field optimized for the estimation of ϕ , but only leads to a reduction in the Fisher informationby a factor of two (see Supplementary Information S3.1).The experimental procedure starts by measuring threereflection matrices for different values of ϕ , which allowsus to access both the reflection matrix r and its deriva-tive ∂ ϕ r for incident states and outgoing states(see Methods and Supplementary Information S4). Eventhough the measured reflection matrix is large, it consti-tutes only a sub-part of the full S -matrix of the medium.This is not a limitation of our approach, which also ap-plies to non-unitary and non-square matrices. Definingthe operator f ϕ = ( ∂ ϕ r ) † ∂ ϕ r , the eigenvector associatedwith the largest eigenvalue of f ϕ is the optimal incidentstate based on the available knowledge. Illuminating themedium with this state using the input SLM, we mea-sured the spatial distributions of the outgoing signal in-tensity (Fig. 2b) and of the Fisher information per unitarea (Fig. 2c), respectively normalized by the average sig-nal intensity and by the average Fisher information perunit area under plane-wave illumination. Averaging overthe field of view, the maximum information state gen-erates a -fold enhancement of the Fisher informationalong with a -fold intensity enhancement, as comparedto the average values measured under plane-wave illumi-nation.We then checked explicitly in which way the maximuminformation state is shaped in the near-field of the cross-shaped phase perturbation. To this end, we measuredthe single-pixel sensitivity, which we define as being thetotal Fisher information in the outgoing state for phasevariations on each individual pixel on the hidden SLM.These measurements are performed by successively vary-ing the phase shift induced by individual pixels sequen-tially instead of varying the phase shift induced by thecross-shaped target as a whole (see Methods). We findthat the maximum information state is primarily sensi-tive to a few pixels in the center of the cross (Fig. 2d),which confirms their economical wave function design.We emphasize that optimal states are not conceived toreveal the shape of the target, but to estimate the phaseshift it induces. This does not require the intensity tobe uniformly distributed on the area defining the target,explaining why its shape is not always fully revealed.An important characteristic of maximum informationstates is their specificity with respect to both the positionof the observer and to the observable of interest. To ex-plicitly probe the influence of the observer, we repeat theexperiment from above using, however, a reduced field ofview in the detection plane (Fig. 2e–h and Extended DataFig. 1). Remarkably, the maximum information statereadjusts to this new observer by redirecting the spatialdistributions of its outgoing signal intensity (Fig. 2f) andof its Fisher information per unit area (Fig. 2g) straightto the selected observer area. This confirms the essentialrole played by the set of optical modes incorporated bythe experimentalist into the definition of maximum infor-mation states. While the single-pixel sensitivity (Fig. 2h)is similar to what is observed when the whole field of view is taken into account (compare with Fig. 2d), it is heremore uniformly distributed inside the cross-shaped area.Indeed, such a redistribution in the plane of the hiddenSLM is associated with smaller diffraction angles, thusallowing for a reduced spatial extent of the intensity dis-tribution in the plane of the diffuser. To demonstratealso the specificity of maximum information states withrespect to the observable of interest, we repeat the ex-perimental procedure by choosing the horizontal position x of a circular phase object (radius µm) as being theobservable parameter of interest (Fig. 2i–l and ExtendedData Fig. 2). The single-pixel sensitivity is now local-ized on the right or left edge of the object, which atteststhat maximum information states selectively focus lightwaves onto those specific areas of the object that havethe most pronounced dependence on the observable ofinterest. Interestingly, the maximum information statetypically focuses on a single edge rather than on bothedges simultaneously, which is to be expected as the mir-ror symmetry in the system is broken by the diffuser.To exemplify the remarkable advantages these featuresprovide for measurements with very few photons wherethe estimation precision is limited by shot noise, we re-duce the incident photon flux by placing a neutral den-sity filter (fractional transmittance . × − ) in theoptical path. Under these conditions, illuminating themedium with a plane wave does not allow reliable es-timations for the phase shift induced by target. Thisis confirmed by calculating the precision limit σ crb forthe plane waves used to construct the reflectionmatrix (see Supplementary Information S3). For theseplane waves, the median of the distribution is . rad,with a minimum value of . rad (Fig. 3a). In contrast,the precision limit associated with the maximum infor-mation state equals . rad, an entire order of magni-tude smaller than the minimum value measured for planewaves. While reaching this precision limit would requireamplitude and phase modulation of the incident field,we only modulate its phase by the input SLM in theexperiments. We find that the precision limit associatedwith phase-only modulation of the maximum informationstate equals . rad, a value that is only slightly largerthan that for a joint amplitude and phase modulation.This observation corroborates that our approach is alsovery robust with respect to small errors or imperfectionsin the preparation of the incident state.In order to demonstrate the practical implementationof the estimation process with a maximum informationstate, we use this state to illuminate the medium andperform a sequence of measurements. The observable pa-rameter we estimate is the phase shift ϕ induced by thecross-shaped target, which is varied every measure-ments between − . rad and . rad in a step-like man-ner. From the knowledge of the expected outgoing stateand its derivative with respect to ϕ , we can construct theminimum variance unbiased estimator of ϕ , applicable to Fig. 2 | Examination of maximum information states. a , Sketch of the experiment: the observer (left) is a camerawith a field of view covering µm , separated by a diffuser (middle) from a cross-shaped target object (right) that inducesa phase shift ϕ as our observable parameter of interest. b , c , Measured spatial distributions of the intensity and of the Fisherinformation per unit area for the optimal state with maximum overall information content, respectively normalized by theaverage signal intensity and by the average Fisher information per unit area under plane-wave illumination. d , Single-pixelsensitivity measured by shifting the phase of each pixel in the target area of the hidden SLM for the optimal state, normalizedby the average single-pixel sensitivity under plane-wave illumination. The position of the cross-shaped target is delimited bywhite dashed lines. e – h , Analogous to a – d when the field of view of the camera covers a reduced area of µm , as delimitedby white dashed lines in f and g . The maximum information state fully adjusts to the changes in the observer by deliveringthe Fisher information here primarily to the limited field of view. i – l , Analogous to a – d when the observable parameter is thehorizontal position x of a circular phase object. Also here the maximum information state adapts to the change in observableby redirecting its incoming intensity to those pixels on the very right or left edge of the target object that are most stronglyaffected by a lateral displacement ∆ x . small parameter variations (see Supplementary Informa-tion S3.3). We can then estimate the value of ϕ frommeasurements of the outgoing state (see Extended DataFig. 3). Even though only approximately , inci-dent photons are probing the medium per measurement,each step can be clearly resolved (Fig. 3b,c). The ob-served standard error on the estimates is . rad, whichcorresponds to a transverse displacement of the target of . nm (see also Extended Data Fig. 4 for estimations oflateral displacements). Importantly, the observed stan-dard error on the estimates almost reaches the precisionlimit predicted by the Cramér-Rao inequality. This con-firms that shot noise is the dominant source of noisein our experiment, and demonstrates that the precisionlimit is correctly predicted from reflection matrix mea-surements. We applied the same procedure for measure-ments performed by illuminating the medium with the best plane wave used to construct the reflection matrix.The precision limit is then much larger than the step size,thereby prohibiting a clear detection of the phase steps(Fig. 3d,e).The Fisher information operator F θ employed in ouranalysis not only constitutes an operational tool for theidentification of maximum information states but it isalso deeply connected to fundamental concepts in optics.Provided that reciprocity holds in terms of the transposi-tion symmetry of the scattering matrix ( S T = S ), we canshow (see Supplementary Information S1.5) that the iter-ative phase conjugation of a small perturbation ∆ θ con-verges towards the largest eigenvalues of F θ . This insightnot only suggests a potentially useful approach to iden-tify maximum information states but also provides a newunderstanding of existing focusing procedures based ontime-reversed adapted perturbation [20–22]. Moreover, Fig. 3 | Demonstration of optimal estimations. a , His-togram of precision limits for the plane waves used toconstruct the reflection matrix (the observable parameter ishere the phase shift ϕ of the cross-shaped target shown inFig. 2a). As compared to the best plane wave, the preci-sion limit is improved by an order of magnitude with maxi-mum information states (AP, amplitude and phase modula-tion; PO, phase-only modulation). b , Estimated phase shiftof the cross-shaped target as a function of measurement indexfor measurements performed by illuminating the medium withthe maximum information state. c , Histograms of estimatedangles for a positive phase shift ( ∆ ϕ + = +0 . rad) and a neg-ative phase shift ( ∆ ϕ − = − . rad) applied by the hiddenSLM. The length of error bars equals σ crb . d , e , Analogousto b and c for measurements performed by illuminating themedium with the best plane wave. in the ideal case of a unitary S -matrix ( S † = S − ) whereall optical modes supported by the system are accessible,we obtain the identity F θ = Q θ , with Q θ = − iS − ∂ θ S being the generalized Wigner-Smith operator. This oper-ator was recently introduced to design optimal light fieldsfor optical micro-manipulation in complex media [16, 17].The simple relation between F θ and Q θ suggests a newinterpretation of Q θ as the operator representing themeasurement backaction on the conjugate quantity to θ . The eigenstates of Q θ (also called principal modes )have the remarkable property of being insensitive withrespect to small variations in θ except for a global phasefactor [23, 24]. Likewise, the Fisher information of max-imum information states is exclusively enclosed in vari-ations of the global phase of the outgoing state, ratherthan in the state’s intensity variations or speckle decor- relation. Nevertheless, this property strictly holds onlyfor a unitary S -matrix and deviations from this propertyare observed in our experiments (see Supplementary In-formation S4.3). Finally, it is interesting to discuss thespecial case for which the observable of interest is thedielectric constant (cid:15) of a target. Indeed, provided that S † = S − , eigenstates of Q (cid:15) maximise the integratedintensity inside the target [17], which implies that maxi-mum information states are then, at the same time, thelight states that maximize power delivery to the target.To summarize, we introduce a general frameworkto perform coherent scattering measurements in com-plex environments that reach the fundamental precisionbound set by the quantum fluctuations of coherent states.This enables us to experimentally generate, for the firsttime, these optimal states of light that maximize theFisher information for any observable of interest and forany set of optical modes controlled by the observer. Theability to produce optimal light states by far-field wave-front shaping opens up new perspectives to enhance theperformance of imaging techniques [14] and, by simulta-neously engineering the Fisher information operator it-self, could also be used to improve the sensitivity of exist-ing nano-structured sensing devices [25, 26]. We empha-size that our results are generally applicable to any para-metric dependence of a wave field and can thus be trans-posed to other types of waves such as in acoustics [27] orin the microwave regime [28]. Finally, it can be expectedthat the classical bound derived in this work will be sur-passed using quantum metrology protocols [11, 29, 30].In this context, our results provide a general benchmarkto assess the performance of quantum states optimizedfor parameter estimation, and suggest a new path to-wards the identification of optimally-sensitive quantumstates of light using scattering matrices of complex sys-tems. Acknowledgements.
The authors thank M. van Beur-den, J. Bosch, S. Faez, M. Horodynski, M. Kühmayer, P.Pai, and J. Seifert for insightful discussions, and P. Jur-rius, D. Killian and C. de Kok for technical support. Thiswork was supported by the Netherlands Organization forScientific Research NWO (Vici 68047618 and PerspectiveP16-08) and by the Austrian Science Fund (FWF) underproject number P32300 (WAVELAND). [1] Y. Park, C. Depeursinge, and G. Popescu, Nat. Photonics , 578 (2018).[2] R. W. Taylor and V. Sandoghdar, Nano Lett. , 4827(2019).[3] W. Osten and N. Reingand, Optical Imaging and Metrol-ogy: Advanced Technologies (John Wiley & Sons, 2012).[4] E. G. van Putten, A. Lagendijk, and A. P. Mosk, Opt.Lett. , 1070 (2012).[5] Y. Shechtman, S. J. Sahl, A. S. Backer, and W. E. Mo- erner, Phys. Rev. Lett. , 133902 (2014).[6] F. Balzarotti, Y. Eilers, K. C. Gwosch, A. H. Gynnå,V. Westphal, F. D. Stefani, J. Elf, and S. W. Hell, Science , 606 (2017).[7] P. Ambichl, W. Xiong, Y. Bromberg, B. Redding, H. Cao,and S. Rotter, Phys. Rev. X , 041053 (2017).[8] G. H. Yuan and N. I. Zheludev, Science , 771 (2019).[9] T. Juffmann, A. de los Ríos Sommer, and S. Gigan, Opt.Commun. , 124484 (2020).[10] D. Bouchet, R. Carminati, and A. P. Mosk, Phys. Rev.Lett. , 133903 (2020).[11] V. Giovannetti, S. Lloyd, and L. Maccone, Nat. Photon-ics , 222 (2011).[12] S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, A. C.Boccara, and S. Gigan, Phys. Rev. Lett. , 100601(2010).[13] S. Rotter and S. Gigan, Rev. Mod. Phys. , 015005(2017).[14] H. H. Barrett and K. J. Myers, Foundations of ImageScience (John Wiley & Sons, 2013).[15] S. Kay,
Fundamentals of Statistical Processing, VolumeI: Estimation Theory (Prentice Hall, 1993).[16] P. Ambichl, A. Brandstötter, J. Böhm, M. Kühmayer,U. Kuhl, and S. Rotter, Phys. Rev. Lett. , 033903(2017).[17] M. Horodynski, M. Kühmayer, A. Brandstötter, K. Pich-ler, Y. V. Fyodorov, U. Kuhl, and S. Rotter, Nat. Pho-tonics , 149 (2020).[18] A. P. Mosk, A. Lagendijk, G. Lerosey, and M. Fink, Nat.Photonics , 283 (2012).[19] S. L. Braunstein, C. M. Caves, and G. J. Milburn, Ann.Phys. (N. Y.) , 135 (1996).[20] E. H. Zhou, H. Ruan, C. Yang, and B. Judkewitz, Optica , 227 (2014).[21] C. Ma, X. Xu, Y. Liu, and L. V. Wang, Nat. Photonics , 931 (2014).[22] H. Ruan, T. Haber, Y. Liu, J. Brake, J. Kim, J. M.Berlin, and C. Yang, Optica , 1337 (2017).[23] J. Carpenter, B. J. Eggleton, and J. Schröder, Nat. Pho-tonics , 751 (2015).[24] W. Xiong, P. Ambichl, Y. Bromberg, B. Redding, S. Rot-ter, and H. Cao, Phys. Rev. Lett. , 053901 (2016).[25] H. Hodaei, A. U. Hassan, S. Wittek, H. Garcia-Gracia,R. El-Ganainy, D. N. Christodoulides, and M. Kha-javikhan, Nature , 187 (2017).[26] W. Chen, S. Kaya Özdemir, G. Zhao, J. Wiersig, andL. Yang, Nature , 192 (2017).[27] B. Gérardin, J. Laurent, A. Derode, C. Prada, andA. Aubry, Phys. Rev. Lett. , 173901 (2014).[28] Z. Shi and A. Z. Genack, Phys. Rev. Lett. , 043901(2012).[29] V. Giovannetti, S. Lloyd, and L. Maccone, Phys. Rev.Lett. , 010401 (2006).[30] L. J. Fiderer, J. M. E. Fraïsse, and D. Braun, Phys. Rev.Lett. , 250502 (2019).[31] P. Pai, J. Bosch, and A. P. Mosk, OSA Continuum ,637 (2020).[32] M. Takeda, H. Ina, and S. Kobayashi, J. Opt. Soc. Am. , 156 (1982).[33] E. Cuche, P. Marquet, and C. Depeursinge, Appl. Opt. , 4070 (2000). METHODS
Construction of the reflection matrix and its deriva-tive.
We construct densely-sampled reflection matricesrelating incident field states to reflected ones [31]. Reflec-tion matrix measurements are performed with no densityfilter in the signal path. To illuminate the medium, wevary the incidence angle of a Gaussian beam character-ized by a full-width at half maximum of µm, whichis times larger than the field of view (such Gaussianbeams are refereed to as plane waves in the manuscript).We probe M = 2437 different incidence angles, coveringa numerical aperture of . . The sampling is performedusing a triangular lattice (in Fourier space), with a lat-tice constant taken to be the smallest angular separationat which the complex inner product of nearest-neighborfields drops to zero.For each incident angle, we record the reflected fieldusing digital off-axis holography [32, 33]. The methodrelies on a reference beam which is tilted by an anglewith respect to the reflected signal beam. The complexfield is then numerically reconstructed from the mea-sured holograms by digitally filtering the zero-order com-ponent. The reflected field is sampled using a triangu-lar lattice, with a lattice constant that we determine byfinding the distance between the maximum and the firstminimum of the autocorrelation of a random field. Weuse N = 2465 different sampling points, covering an areaof µm on the surface of the diffuser. The reflectionmatrix r is therefore constructed column by column and,as a result, we obtain a × matrix.In the initial experiment, the parameter of interest isthe phase shift ϕ generated by a cross-shaped target ob-ject. The reflection matrix r is then measured at an angle ϕ = ϕ by setting the phase shift induced by all pixels ofthe hidden SLM to a given value. Note that ϕ can be setto zero without loss of generality. In order to estimatethe derivative of the reflection matrix with respect to ϕ , we measure two other reflection matrices r ( ϕ − ∆ ϕ ) and r ( ϕ + ∆ ϕ ) , where ∆ ϕ = 0 . rad. We can thenestimate ∂ ϕ r by applying the centered finite differencescheme ∂ ϕ r (cid:39) [ r ( ϕ + ∆ ϕ ) − r ( ϕ − ∆ ϕ )] / (2∆ ϕ ) .In another experiment, the parameter of interest is cho-sen to be the lateral position x of a circular object. In thiscase, a super-Gaussian function of order is displayed bythe hidden SLM, with a full-with at half maximum equalto µm. The phase difference between the object andthe background is set to π/ rad and the reflection ma-trix r is then measured at a position x = x . In orderto estimate the derivative of the reflection matrix withrespect to x , we measure two other reflection matrices r ( x − ∆ x ) and r ( x + ∆ x ) , where ∆ x = 5 µm. We canthen estimate ∂ x r by applying the centered finite differ-ence scheme ∂ ϕ r (cid:39) [ r ( x + ∆ x ) − r ( x − ∆ x )] / (2∆ x ) . Measurement of the single-pixel sensitivity.
In orderto measure the single-pixel sensitivity, we successivelyvary the phase shift ϕ j induced by each individual pixel j for pixels covering an area of µm on the sur-face of the hidden SLM. This allows us to access thederivative of the outgoing field with respect to ϕ j usinga centered difference scheme, and to calculate the asso-ciated Fisher information using Eq. (1). For each indi-vidual pixel, we perform an averaging of the derivativeof the outgoing field over independent measurements.Once the Fisher information is calculated, we also sub-tract a residual noise floor that we estimate by takingdifferent measurements of the same outgoing state. Weperform the same analysis by illuminating the mediumwith the maximum information state and with differentplane waves, and we normalize the values obtained forthe maximum information state using the average valueobtained with different plane waves. Note that a highsingle-pixel sensitivity cannot be achieved without a highintensity inside the pixel area. Thus, mapping the single-pixel sensitivity in the plane of the hidden SLM also pro-vides us with an indirect way to approximate the inten-sity distribution in the plane of the hidden SLM. Monitoring of the global phase drift.
Due to the sig-nificant acquisition time ( min in total), the globalphase of the measured outgoing field slowly drifts in timebecause of the imperfect thermal stability of the exper-imental setup. During the acquisition, we continuously monitor this drift by regularly measuring a known out-going field as a phase reference [31]. We use differentphase-reference fields depending on the incident poweron the sample.When no density filters are placed in the optical path,the phase-reference field is generated by illuminating themedium with a given plane wave, with a slight angle sothat no reflection from the back-focal plane of the objec-tive can be observed. We then calculate how the globalphase of this field changes over time by using a complexinner product of the phase-reference field measured at agiven time with the phase-reference field measured at thebeginning of the acquisition.When any density filter is placed in the optical path,we first calculate a truncated reflection matrix r (cid:48) , whichdoes not include the few columns for which reflectionfrom the back-focal plane of the objective can be ob-served. The phase-reference field is generated by illumi-nating the medium using the right-singular vector of r (cid:48) associated with its largest singular value. By doing so,we maximize the signal-to-noise ratio of phase-referencemeasurements. We then calculate how the global phaseof this field changes over time by using a complex innerproduct of the phase-reference field measured at a giventime with the field predicted from the knowledge of thereflection matrix.Finally, we perform linear interpolations for estimatingthe global phase drift at any time during the acquisition,and we subsequently apply the appropriate phase correc-tion to any measured data. Extended Data Fig. 1. | Measured intensity and Fisher information distributions for maximum informationstates associated with different detection areas. a – d , Measured spatial distributions of the intensity for optimal states,normalized by the average signal intensity under plane-wave illumination. The observable parameter is the phase shift inducedby the cross-shaped target, and optimal states are defined here with respect to a reduced field of view of the camera that coversan area of µm , as delimited by white dashed lines. e – h , Analogous to a – d for the measured spatial distribution of theFisher information per unit area. Remarkably, the maximum information state always delivers the Fisher information to thedesignated observer window. Extended Data Fig. 2. | Single-pixel sensitivity for maximum information states associated with differentobservable parameters. a – c , Single-pixel sensitivity measured by shifting the phase of each pixel in the target area of thehidden SLM for the optimal state, normalized by the average single-pixel sensitivity under plane-wave illumination. The objectdisplayed on the hidden SLM is a circular phase object, whose position is delimited by a white dashed circle. The field ofview of the detection camera covers here an area of µm . The incident states used to illuminate the scattering mediumare the maximum information states relative to a phase shift ( a ), a horizontal shift ( b ) and a vertical shift ( c ) of the object. d – f , Analogous to a – c when the field of view of the camera covers a reduced area of µm . In all cases, the maximuminformation state directs the incoming intensity to those parts on the hidden SLM that are most affected by the change inthe observable parameter. Interestingly, when the target parameter is either a horizontal or a vertical shift of the object, themaximum information state typically focuses on a single edge rather than on both edges simultaneously, which is to be expectedas the mirror symmetry in the system is broken by the diffuser. Extended Data Fig. 3. | Predicted and measured signal intensity distribution at low photon counts. a , Predicteddistribution of the signal intensity for the optimal incident state expressed in analog-to-digital units (ADU). This distributionis calculated from the measured reflection matrix, considering that the neutral density filter ND6 (fractional transmittance . × − ) is placed in the optical path. b , Measured distribution of the signal intensity for the optimal incident state whenND6 is in the signal path. Such measurements are shot-noise limited, and the observed signal-to-noise per pixel is largely smallerthan unity. Thus, the measured distribution of the signal intensity appears as a random noise, which has been low-pass filteredby the data analysis procedure used to digitally reconstruct complex fields from off-axis intensity measurements. Despite thislow signal-to-noise per pixel, such data allow to correctly estimate the phase shift induced by the hidden target when usingthe minimum variance unbiased estimator. This can be achieved since only a single parameter (the phase shift induced by thetarget) needs to be estimated from a large number of independent sampling points. The Fisher information associated witheach pixel of the detection camera effectively adds up, resulting in a total Fisher information that is sufficient to resolve thephase steps induced on the hidden SLM. c , d , Analogous to a and b for the best plane wave used to construct the reflectionmatrix. Extended Data Fig. 4. | Estimations of lateral displacements at low photon counts.
Analogous to Fig. 3 when theobservable parameter is the horizontal position x of the circular phase object shown in Fig. 2i. a , Histogram of precision limitsfor the plane waves used to construct the reflection matrix and for the maximum information states (AP, amplitude andphase modulation; PO, phase-only modulation). b , Estimated lateral displacements of the circular phase object as a functionof measurement index for measurements performed by illuminating the medium with the maximum information state. Thecalculated precision limit equals . µm, and the observed standard error on the estimates is . µm. c , Histograms of estimatedangles for a positive lateral displacement ( ∆ x + = +2 . µm) and a negative lateral displacement ( ∆ x − = − . µm) applied bythe hidden SLM. The length of error bars equals σ crb . d , e , Analogous to b and c for measurements performed by illuminatingthe medium with the best plane wave. Maximum information states for coherent scattering measurementsSupplementary information
Dorian Bouchet, Stefan Rotter, and Allard P. Mosk Nanophotonics, Debye Institute for Nanomaterials Science,Utrecht University, P.O. Box 80000, 3508 TA Utrecht, the Netherlands Institute for Theoretical Physics, Vienna University of Technology (TU Wien), A-1040 Vienna, Austria
SUPPLEMENTARY SECTION S1 – FISHER INFORMATION FOR THE OPTIMAL DETECTIONSCHEME
In this section, we calculate the classical Fisher information for the optimal homodyne detection scheme, which weshow to be equal to the quantum Fisher information for coherent states and statistically independent samples. Wethen define the Fisher information operator from the scattering matrix of the system, and we provide the expression ofthe minimum variance unbiased estimator. Finally, we show that maximum information states can also be iterativelyidentified by performing time-reversal of a small perturbation using phase conjugation.
S1.1 – Classical Fisher information
Let us describe the measured data by a N -dimensional random variable X and a joint probability density function p ( X ; θ ) parameterized by an arbitrary parameter θ . In general, the Fisher information is expressed by [1] J ( θ ) = E (cid:0) [ ∂ θ ln p ( X ; θ )] (cid:1) , (S1)where E denotes the expectation operator acting over noise fluctuations. In the case of scattering measurementsperformed using a homodyne detection scheme with a strong reference beam, readout noise can be neglected. Thesignal measured at the k -th sampling point can then be modeled by a Poisson distribution of expectation I tot k = | E out k + E ref k | , where E out k and E ref k are the complex values of the outgoing signal field and of the reference field,respectively. Considering N statistically independent sampling points, the joint probability density function p ( X ; θ ) is therefore expressed by p ( X ; θ ) = N (cid:89) k =1 e −| E out k + E ref k | | E out k + E ref k | X k X k ! . (S2)Injecting this expression in Eq. (S1), we obtain J ( θ ) = N (cid:88) k =1 (cid:2) ∂ θ (cid:0) | E out k + E ref k | (cid:1)(cid:3) | E out k + E ref k | . (S3)For a strong reference field which does not depend on θ (i.e., for | E ref k | (cid:29) | E out k | and ∂ θ E ref k = 0 ), this expressionsimplifies to J ( θ ) = N (cid:88) k =1 (cid:16) ∂ θ (cid:104) E out k (cid:0) E ref k (cid:1) ∗ + E ref k ( E out k ) ∗ (cid:105)(cid:17) | E ref k | . (S4)We now introduce E out k = Q out k + iP out k and E ref k = | E ref k | e iφ k . Equation (S4) becomes J ( θ ) = 4 N (cid:88) k =1 (cid:2) ( ∂ θ Q out k ) cos φ k + ( ∂ θ P out k ) sin φ k (cid:3) . (S5)3For any integer m , the phase angles of the reference field which maximize Eq. (S5) are given by φ max k = arg( ∂ θ E out k ) + mπ , and those which minimize Eq. (S5) are given by φ min k = arg( ∂ θ E out k ) + ( m + 1 / π . Choosing φ min k as phase anglesof the reference field yields J ( θ ) = 0 . In contrast, choosing φ max k as phase angles of the reference field yields J ( θ ) = 4 N (cid:88) k =1 | ∂ θ E out k | . (S6)This expression can be identified as the Fisher information for a random vector composed of N complex randomvariables whose real and imaginary parts are independent normally distributed random variables with variance σ =1 / . S1.2 – Quantum Fisher information
The classical Fisher information J ( θ ) sets a lower bound on the variance of unbiased estimators of θ for a definitemeasurement scheme (i.e. the homodyne scheme in our case). A more general lower bound exists, which applies toany quantum measurement described by a positive-operator-valued measure (POVM). This lower bound is given bythe reciprocal of the quantum Fisher information, which is expressed by [2] I ( θ ) = Tr( ρ out L ) , (S7)where ρ out is a density operator describing the quantum state of the system and L out is the symmetrized logarithmicderivative of ρ out with respect to θ defined as follows: ρ out L out + L out ρ out = 2 ∂ θ ρ out . (S8)Considering that ρ out describes a N -modes coherent state composed of simply-separable pure states, Eq. (S7) simplifiesto [3, 4] I ( θ ) = 4 N (cid:88) k =1 (cid:0) (cid:104) ∂ θ α k | ∂ θ α k (cid:105) − |(cid:104) ∂ θ α k | α k (cid:105)| (cid:1) , (S9)where | α k (cid:105) is the single-mode coherent state associated with the k -th mode and | ∂ θ α k (cid:105) is its derivative with respectto θ . We can represent | α k (cid:105) in the basis of Fock states | n (cid:105) labeled by the occupation number n , which reads [5] | α k (cid:105) = e −| α k | / ∞ (cid:88) n =0 α nk ( n !) / | n (cid:105) , (S10)where α k denotes the eigenvalue of the annihilation operator. Taking the derivative of this expression with respect to θ leads to | ∂ θ α k (cid:105) = ∂ θ α k (cid:34) −| α k | e −| α k | / ∞ (cid:88) n =0 α nk ( n !) / | n (cid:105) + e −| α k | / ∞ (cid:88) n =1 nα n − k ( n !) / | n (cid:105) (cid:35) . (S11)We can now use Eqs. (S10) and (S11) to calculate the two terms (cid:104) ∂ θ α k | ∂ θ α k (cid:105) and |(cid:104) ∂ θ α k | α k (cid:105)| that appear in Eq. (S9).Let us first calculate (cid:104) ∂ θ α k | ∂ θ α k (cid:105) from Eq. (S11): using the orthonormality of Fock states, we obtain (cid:104) ∂ θ α k | ∂ θ α k (cid:105) = | ∂ θ α k | (cid:34) | α k | e −| α k | ∞ (cid:88) n =0 | α k | n n ! − | α k | α k e −| α k | ∞ (cid:88) n =1 | α k | n − ( n − −| α k | α ∗ k e −| α k | ∞ (cid:88) n =1 | α k | n − ( n − e −| α k | ∞ (cid:88) n =1 n | α k | n − ( n − (cid:35) . (S12)This expression can be simplified by using the following properties of exponential series: ∞ (cid:88) n =0 | α k | n n ! = e | α k | , (S13)4 ∞ (cid:88) n =0 ( n + 1) | α k | n n ! = (1 + | α k | ) e | α k | . (S14)Injecting Eqs. (S13) and (S14) into Eq. (S12) leads to (cid:104) ∂ θ α k | ∂ θ α k (cid:105) = | ∂ θ α k | (2 | α k | − | α k | α k − | α k | α ∗ k + 1) . (S15)Let us now calculate (cid:104) ∂ θ α k | α k (cid:105) from Eqs. (S10) and (S11): using again the orthonormality of Fock states, we obtain (cid:104) ∂ θ α k | α k (cid:105) = ∂ θ α ∗ k (cid:34) −| α k | e −| α k | ∞ (cid:88) n =0 | α k | n n ! + α k e −| α k | ∞ (cid:88) n =1 | α k | n − ( n − (cid:35) . (S16)Injecting Eq. (S13) into Eq. (S16) leads to |(cid:104) ∂ θ α k | α k (cid:105)| = | ∂ θ α k | (2 | α k | − | α k | α k − | α k | α ∗ k ) . (S17)Finally, injecting Eqs. (S15) and (S17) into Eq. (S9) results in I ( θ ) = 4 N (cid:88) k =1 | ∂ θ α k | . (S18)We can use the same quantum-mechanical formalism to derive the classical Fisher information for the homodynedetection scheme. To this end, we consider the coherent state | α k , r k (cid:105) = | α k + r k (cid:105) which represents the superpositionof the coherent state | α k (cid:105) with a reference coherence state | r k (cid:105) . The probability to observe n k photons in the state | α k , r k (cid:105) is P k = |(cid:104) n k | α k + r k (cid:105)| , which leads to a Poisson distribution of expectation value | α k + r k | . In the sameway that we obtained Eq. (S6) from Eq. (S1) with classical fields, we now obtain J ( θ ) = 4 N (cid:88) k =1 | ∂ θ α k | . (S19)The classical Fisher information given by Eq. (S19) equals the quantum Fisher information given by Eq. (S18), therebyensuring that the homodyne detection scheme that we considered is optimal for the estimation of θ . S1.3 – Fisher information operator
In the formalism of the S -matrix, the outgoing field state is expressed by | E out (cid:105) = S | E in (cid:105) , where S denotes thescattering matrix (or S -matrix) and | E in (cid:105) is the incident field state. We can write E out k as a projection of the outgoingstate on the state associated with the k -th sampling point, which we write as E out k = (cid:104) k | S | E in (cid:105) in bra-ket notations,leading to J ( θ ) = 4 N (cid:88) k =1 |(cid:104) k | ∂ θ S | E in (cid:105)| . (S20)This expression can be expanded into J ( θ ) = 4 N (cid:88) k =1 (cid:104) E in | ( ∂ θ S ) † | k (cid:105)(cid:104) k | ∂ θ S | E in (cid:105) . (S21)Using the completeness relation (cid:80) k | k (cid:105)(cid:104) k | = I N where I N is the N -dimensional identity matrix, we finally obtain J ( θ ) = 4 (cid:104) E in | ( ∂ θ S ) † ∂ θ S | E in (cid:105) . (S22)In this expression, we can identify the operator F θ = ( ∂ θ S ) † ∂ θ S , which we refer to as Fisher information operator .We obtain the quadratic form J ( θ ) = 4 (cid:104) E in | F θ | E in (cid:105) , which is the expression of the Fisher information given in themanuscript. Finally, for a unitary scattering matrix ( S † = S − ), the operator F θ is expressed by F θ = ( − iS − ∂ θ S ) , (S23)5where we used the identity ∂ θ S − = − S − ( ∂ θ S ) S − . Introducing the generalized Wigner-Smith operator [6], whichis defined by Q θ = − iS − ∂ θ S , we obtain the identity F θ = Q θ , (S24)which implies that F θ and Q θ share the same eigenstates. Writing the eigenvalue equation Q θ |E in j (cid:105) = Θ j |E in j (cid:105) where Θ j is the j -th eigenvalue of Q θ and |E in j (cid:105) is the associated eigenstate, the outgoing field state satisfies | ∂ θ E out j (cid:105) = i Θ j |E out j (cid:105) . For a scattering matrix evaluated at θ and a small parameter variation ∆ θ = θ − θ , this results in |E out j ( θ ) (cid:105) (cid:39) e i Θ j ∆ θ |E out j ( θ ) (cid:105) . Thus, in the limit of a unitary scattering matrix, maximum information states areinsensitive with respect to small variations in θ except for a global phase factor. S1.4 – Minimum variance unbiased estimator
For small parameter variations around a given parameter value noted θ , the measured data can be described bythe following linear model: X k = I tot k + ( ∂ θ I tot k )( θ − θ ) + W k , (S25)where X k represents the intensity data measured by the camera at the k -th sampling point, where I tot k and ∂ θ I tot k areevaluated at θ , and where W k are N independent and normally distributed random variables with mean zero andvariance I tot k . The normal distribution is indeed a good approximation of the Poisson distribution for large expectationvalues. In the case of linear models, general expressions exist for the minimum variance unbiased estimator, whichdepends on the noise statistics [1]. For the linear model expressed by Eq. (S25), the minimum variance unbiasedestimator reads ˆ θ ( X ) − θ = 1 J ( θ ) N (cid:88) k =1 ( ∂ θ I tot k )( X k − I tot k ) I tot k . (S26)Writing I tot k = | E out k + E ref k | and recalling that | E ref k | (cid:29) | E out k | and ∂ θ E ref k = 0 , this expression simplifies to ˆ θ ( X ) − θ = 1 J ( θ ) N (cid:88) k =1 (cid:16) ∂ θ (cid:104) E out k (cid:0) E ref k (cid:1) ∗ + E ref k ( E out k ) ∗ (cid:105)(cid:17) ( X k − | E out k + E ref k | ) | E ref k | . (S27)Introducing E ref k = | E ref k | e iφ k and using the phase angles φ max k = arg( ∂ θ E out k ) which maximize the Fisher informationleads to ˆ θ ( X ) − θ = 2 J ( θ ) N (cid:88) k =1 | ∂ θ E out k | ( X k − | E out k + E ref k | ) | E ref k | . (S28)In this expression, J ( θ ) is obtained from Eq. (S6), E out k is obtained from E out k = (cid:104) k | S | E in (cid:105) and ∂ θ E out k is obtainedfrom ∂ θ E out k = (cid:104) k | ∂ θ S | E in (cid:105) , with a scattering matrix S evaluated at θ . S1.5 – Finding maximum information states using phase conjugation
Time reversal using phase conjugation is a well-known technique that allows for instance to focus waves [7–9] orto identify open channels [10] in multiple scattering media. In some implementations, light waves are focused into ascattering medium by performing time-reversal of a perturbation using phase conjugation [11–13]. Here, we show thatan iterative procedure based on such techniques actually converges toward maximum information states. Light wavesare then focused onto those specific areas of the object that are most affected by the perturbation, in such a way thatthe Fisher information available to the observer is maximized. In order to demonstrate the link between maximuminformation states and time-reversed adapted perturbation, we must restrict the analysis to a N × N scattering matrixthat satisfies the reciprocity relation S T = S (or equivalently S † = S ∗ ). The procedure then relies on an iterativeapproach used to calculate the incident state | E in( n ) (cid:105) at the n -th iteration from measurements based on an incident state | E in( n − (cid:105) . The medium is illuminated with the incident state | E in( n − (cid:105) and, before a perturbation ∆ θ has occurred,6the measured outgoing state is | E out( n − ( θ ) (cid:105) = S ( θ ) | E in( n − (cid:105) . For the same incident state, the outgoing state measuredafter the perturbation has occurred is | E out( n − ( θ + ∆ θ ) (cid:105) = S ( θ + ∆ θ ) | E in( n − (cid:105) . The new incident state | E in( n ) (cid:105) is thenobtained by phase-conjugating the difference between these two outgoing states. In the limit of a small perturbation ∆ θ , this leads to | E in( n ) (cid:105) = A − / n ) ∂ θ S ∗ | E in( n − (cid:105) ∗ , (S29)where A − / n ) is a (real-valued) normalization coefficient. Using the relation S ∗ = S † , we can write ∂ θ S ∗ = ∂ θ S † .Then, starting from a random incident state | E in(0) (cid:105) and after n successive iterations, we obtain | E in( n ) (cid:105) = B − / n ) (cid:2) ( ∂ θ S ) † ∂ θ S (cid:3) n/ | E in(0) (cid:105) if n is even, (S30) | E in( n ) (cid:105) = B − / n ) (cid:2) ( ∂ θ S ) † ∂ θ S (cid:3) n/ − ∂ θ S † | E in(0) (cid:105) ∗ if n is odd. (S31)Choosing the normalization condition (cid:104) E in( n ) | E in( n ) (cid:105) = 1 , the (real-valued) normalization coefficient B ( n ) is expressed by B ( n ) = (cid:104) E in(0) | (cid:2) ( ∂ θ S ) † ∂ θ S (cid:3) n | E in(0) (cid:105) . (S32)When the medium is illuminated with the incident state | E in( n ) (cid:105) , the Fisher information associated with the resultingoutgoing state reads J ( n ) ( θ ) = 4 (cid:104) E in( n ) | ( ∂ θ S ) † ∂ θ S | E in( n ) (cid:105) . (S33)Using either Eq. (S30) if n is even or Eq. (S31) if n is odd, Eq. (S33) becomes J ( n ) ( θ ) = 4 (cid:104) E in(0) | [( ∂ θ S ) † ∂ θ S ] n +1 | E in(0) (cid:105)(cid:104) E in(0) | [( ∂ θ S ) † ∂ θ S ] n | E in(0) (cid:105) . (S34)Any incident state can be decomposed in the orthonormal basis formed by the eigenstates of F θ . In this basis, theinitial state | E in(0) (cid:105) is expressed by | E in(0) (cid:105) = N (cid:88) j =1 γ j |E in j (cid:105) , (S35)where |E in j (cid:105) is the j -th eigenstate of F θ and γ j = (cid:104)E in j | E in(0) (cid:105) . From Eq. (S35) and using the eigenvalue equation F θ |E in j (cid:105) = Λ j |E in j (cid:105) , Eq. (S34) becomes J ( n ) ( θ ) = 4 (cid:80) Nj =1 | γ j | Λ n +1 j (cid:80) Nj =1 | γ j | Λ nj . (S36)This expression can also be written in the following form: J ( n ) ( θ ) = 4Λ max × (cid:80) Nj =1 | γ j | (Λ j / Λ max ) n +1 (cid:80) Nj =1 | γ j | (Λ j / Λ max ) n , (S37)where Λ max is the largest eigenvalue of F θ . Assuming that | γ j | (cid:54) = 0 for the associated eigenstate, we end up with lim n → + ∞ J ( n ) = 4Λ max . (S38)Thus, performing time-reversal of a perturbation using phase conjugation allows one to iteratively identify the max-imum information state relative to the perturbation, provided that the scattering matrix describing the medium isa square matrix satisfying the reciprocity condition S T = S , and that the overlap between the initial state and themaximum information state is not equal to zero.7 Supplementary Figure S1 | Fisher information maximization using time-reversed adapted perturbation.
Fisherinformation associated with light states identified using time-reversed adapted perturbation for
10 000 random initial states,normalized by the Fisher information of the maximum information state. A horizontal line goes through each box at the medianvalue, edges of the boxes represent lower and upper quartiles, and whiskers represent minimum and maximum values. Thebox and whiskers for iterations appear as a single horizontal line; indeed, random initial states generate a Fisher informationthat is much smaller than the maximum information state. The solid red line represents the Fisher information for the largesteigenvalue of f ϕ , and the dashed red lines represent the Fisher information for the second, third, fourth and fifth largesteigenvalues. It clearly appears that the eigenstates of F θ constitute a relevant basis to analyze this iterative procedure. Indeed,the number of iterations needed for the procedure to converge depends on both the values of γ j and the eigenvaluespectrum of F θ , according to Eq. (S36). In order to illustrate this feature, we consider the operator f ϕ = ( ∂ ϕ r ) † ∂ ϕ r measured by shifting the phase ϕ induced by the cross-shaped object. We numerically generate
10 000 random initialstates and we use Eq. (S36) to calculate the ratio J ( n ) ( ϕ ) / J max ( ϕ ) , where J max ( ϕ ) is the Fisher information associatedwith the maximum information state. After the first iteration, the median value of the ratio J ( n ) ( ϕ ) / J max ( ϕ ) is equalto . (Fig. S1), and we observe that iterations are generally required to identify a state that reaches ofthe optimal Fisher information. Moreover, some initial states show a small overlap with the maximum informationstate; for such states, the ratio J ( n ) ( ϕ ) / J max ( ϕ ) after iterations is close to . , corresponding to the second largesteigenvalue of f ϕ . SUPPLEMENTARY SECTION S2 – OPTICAL SETUP
The optical setup is represented in Fig. S2. The light source is a continuous wave solid-state laser (Coherent OBIS532-120 LS FP) emitting at nm. The laser light is coupled to a single-mode polarization-maintaining fiber andout-coupled using a collimator (Schäfter-Kirchhoff, 60FC-L-4-M75-01). The beam is separated into a signal path anda reference path using a 50:50 beamsplitter. In the signal path, the light beam passes a 50:50 beamsplitter and isthen modulated by the input SLM (Holoeye Pluto VIS). It can then pass different neutral density filters which aremechanically placed or removed. The incident power on the objective is: . µW when all density filters are removedfrom the optical path; nW when a neutral density filter of optical density (ND2) is placed in the optical path(measured fractional transmittance T (cid:48) = 0 . × − ); pW when a neutral density filter of optical density (ND6)is placed in the optical path (measured fractional transmittance T = 0 . × − ).The surface of the SLM is imaged onto a ground glass diffuser using a 4f system composed by a mm lens anda × objective (Nikon 50X CFI60 TU Plan Epi ELWD, 0.6 NA). A 90:10 beamsplitter is located between the lensand the objective. The diffuser is made by polishing a microscope coverslip. The resulting scattering angle, definedfrom the full width at half maximum of the transmitted intensity distribution, is approximately ◦ . The diffuser ismounted on the windows of the hidden SLM (Holoeye Pluto BB) at a distance of . mm. The surface of the diffuseris then imaged using a 4f system composed by the objective and a mm lens using a CCD camera (AVT StingrayF145-B) with an exposure time of µs. Before reaching the camera, the beam passes a polarizer to ensure that onlythe horizontal component of the field is measured, and also passes a 90:10 beamsplitter to recombine the referencepath with the signal path. Camera acquisition is triggered by the input SLM in order to limit phase noise due to theflicker of the SLM.8 Supplementary Figure S2 | Optical setup.
The phases on the input SLM are modulated to reproduce the maximuminformation state which reaches optimal sensitivity in its output with respect to any specified parameter characterizing theobject displayed on the hidden SLM, such as phase variations or lateral displacements. Neutral density filters are removed forreflection matrix measurements. BS, beamsplitter; ND, neutral density filters; Obj, objective; NA, numerical aperture; Pol,linear polarizer; L and L , lenses with focal length mm. SUPPLEMENTARY SECTION S3 – FISHER INFORMATION IN THE EXPERIMENTS
In this section, we give the expression of the Fisher information for our experiments, we verify the noise statistics ofthe quadrature components of the field, and we provide the expression of the minimum variance unbiased estimatorused in the experiments.
S3.1 – Expression of the Fisher information in the experiments
In our proof-of-principle experiments, the reference beam is not optimally shaped but consists of a tilted planewave. In this case, the Fisher information is obtained by averaging Eq. (S5) over the phase angle of the referencebeam: J ( θ ) = 2 π N (cid:88) k =1 (cid:90) π d φ (cid:2) ( ∂ θ Q out k ) cos φ + ( ∂ θ P out k ) sin φ (cid:3) . (S39)We obtain the following simplified expression: J ( θ ) = 2 N (cid:88) k =1 | ∂ θ E out k | . (S40)Note that this expression is identical to the expression of the optimal Fisher information given in Eq. (S6), except fora factor of two. Hence, the incoming state which maximizes Eq. (S40) also maximizes Eq. (S6).In the experiments, we must also take into account two additional considerations. First, the value of the N = 2465 sampling points are actually estimated from N cam = 53 044 statistically independent camera pixels (oversampling isrequired in off-axis holography). The oversampling ratio β = N cam /N must therefore be included as a prefactor inEq. (S40). Moreover, while we measured the reflection matrix without any neutral density filter, we consistently usethe strongest neutral density filter ND6 (fractional transmittance T = 0 . × − ) to compute all values of the Fisherinformation. Therefore, the expression of the Fisher information relevant to our experiments is J ( θ ) = T σ (cid:104) E in | ( ∂ θ r ) † ∂ θ r | E in (cid:105) , (S41)where we noted σ = 1 / (2 β ) in order to highlight that this expression can be identified as the Fisher information fora random vector composed of N complex random variables whose real and imaginary parts are independent normallydistributed random variables with variance σ . The precision limit σ crb , which bounds the standard deviation of anyestimator of θ in our experiments, is then expressed from Eq. (S41) by σ crb = J − / .9 S3.2 – Gaussian distribution of the quadrature components of the field
We can verify that the measured quadrature components of the field follow a Gaussian distribution when theneutral density filter ND6 is in the signal path. To this end, we illuminate the medium with the optimal incident staterelative to the estimation of ϕ and we perform successive measurements of the outgoing field. We then computethe p-value for each sampling point according to the Shapiro-Wilk test. The uniform distribution of p-value (Fig. S3a)confirms that the measured quadrature components of the field follow a Gaussian distribution. We also constructthe histogram of the sample variance (Fig. S3b), with a mean value of . ADU. This value is in good agreementwith the theoretical value given by / (2 β ) = 0 . . We consistently choose the mean value experimentallydetermined ( σ = 0 . ADU) to calculate the Fisher information in Eq. (S41).
Supplementary Figure S3 | Statistics of the quadrature components of the field. a , P-value distribution for the sampling points. The p-value for each point is calculated using Shapiro-Wilk test statistics from different realizationof the random noise. b , Distribution of the sample variance for the sampling points. S3.3 – Minimum variance unbiased estimator in the experiments
For small variations of θ around θ , measured data can be described by the following linear model: Z k = E out k + ( ∂ θ E out k )( θ − θ ) + W k , (S42)where Z k represents here the complex field retrieved from the measured data using off-axis holography and evaluatedat the k -th sampling point, where E out k and ∂ θ E out k are evaluated at θ , and where W k are N independent complexrandom variables whose real and imaginary parts are independent normally distributed random variables with meanzero and variance σ . For this linear model, the minimum variance unbiased estimator reads [1] ˆ θ ( Z ) − θ = Re [ (cid:104) ∂ θ E out | Z (cid:105) − (cid:104) ∂ θ E out | E out (cid:105) ] (cid:104) ∂ θ E out | ∂ θ E out (cid:105) . (S43)All estimations presented in the manuscript are performed from data measured with the neutral density filter ND6(characterized by a fractional transmittance T ) placed in the signal optical path. Moreover, due to measurementnoise during reflection matrix measurements, the intensity and the Fisher information predicted from the knowledgeof the reflection matrix differ from direct measurements by a factor η i and η f , respectively (see SupplementaryInformation S3). The estimator of θ associated with our experimental setup is thus obtained from Eq. (S43) by taking | E out (cid:105) = ( η i T ) / r | ˜ E in (cid:105) and | ∂ θ E out (cid:105) = ( η f T ) / ∂ θ r | ˜ E in (cid:105) . We obtain the following expression: ˆ θ ( Z ) − θ = Re (cid:104) ( T η f ) − / (cid:104) ˜ E in | ( ∂ θ r ) † | Z (cid:105) − ( η i /η f ) / (cid:104) ˜ E in | ( ∂ θ r ) † r | ˜ E in (cid:105) (cid:105) (cid:104) ˜ E in | ( ∂ θ r ) † ∂ θ r | ˜ E in (cid:105) . (S44)Evaluating this estimator using measured data Z directly yields the estimates shown in the manuscript. Importantly,all experimental parameters involved in this expression (i.e. T , η i and η f ) are characterized using independentmeasurements. Note that, whereas estimates shown in the manuscript are almost unbiased, small biases often appearwhen applying Eq. (S44) to experimental data. Such biases, which are usually smaller than σ crb , can be explainedby the influence of measurement noise in the reflection matrices and by an imperfect correction of the global phasevariations induced by thermal instabilities of the setup.0 SUPPLEMENTARY SECTION S4 – CHARACTERIZATION OF THE OUTGOING FIELD AND ITSDERIVATIVE
In this section, we show that reflection matrix measurements allow to faithfully predict the outgoing field andits derivatives with respect to the phase and to the position of hidden objects, and we characterize the correlationbetween the outgoing field distribution and its derivative.
S4.1 – Predicting the outgoing field
In general, illuminating the scattering medium with an arbitrary state | E in (cid:105) requires amplitude and phase mod-ulation of the incident field state. However, in the experiments, only the phase of the field can be modulated bythe input SLM. We must therefore determine the expression of the phase-only modulated state | ˜ E in (cid:105) that is ex-perimentally used to illuminate the medium instead of | E in (cid:105) . Let us define the SLM pattern | E slm (cid:105) = M | E in (cid:105) ,where M is a transformation matrix mapping incident states (expressed in a basis of plane waves) to SLM pat-terns (expressed in a basis of SLM pixels). The phase-only modulated state can be numerically approximated by | ˜ E in (cid:105) = argmin( (cid:107)| ˜ E SLM (cid:105) − M | E in (cid:105)(cid:107) , | E in (cid:105) ) , where | ˜ E slm (cid:105) is the SLM pattern that has the same phase as | E slm (cid:105) butwith uniform amplitude. Thus, predicted outgoing states are given by | E outap (cid:105) = r | E in (cid:105) for an amplitude and phasemodulation, and by | E outpo (cid:105) = r | ˜ E in (cid:105) for a phase-only modulation.We can test this procedure by measuring the outgoing field when the medium is illuminated using optimal incidentstates. This characterization is experimentally performed by averaging the outgoing field over measurements, withthe neutral density filter ND2 placed in the signal path (measured fractional transmittance T (cid:48) = 0 . × − ) to avoidsaturation of the camera. The complex correlation coefficients between measured and predicted fields are . − . i (for the maximum information state relative to the estimation of the phase shift of the cross-shaped object) and .
95 + 0 . i (for the maximum information state relative to the estimation of the lateral shift of the circular object). Supplementary Figure S4 | Predicted and measured intensity distributions for maximum information states.a , b , Predicted spatial distribution of the normalized intensity for amplitude and phase (AP) modulation and phase-only (PO)modulation, when the observable parameter is the phase shift of the cross-shaped object. c , Measured spatial distributionof the normalized intensity. d – f Analogous to a – c when the observable parameter is the lateral displacement of the circularobject. η i (cid:39) . when compared to the predicted intensity. S4.2 – Predicting the derivative of the outgoing field
In order to faithfully assess the Fisher information in the experiments, we must ensure that ∂ θ r is correctly estimatedusing the finite-difference scheme ∂ θ r (cid:39) [ r ( θ + ∆ θ ) − r ( θ − ∆ θ )] / (2∆ θ ) . To this end, the change in the measuredoutgoing states generated by parameter variations of ± ∆ θ must be larger than the level of noise in the measurements.In Fig. S5a, we consider the case in which the observable of interest is the phase shift induced by the cross-shapeobject, and we show the distribution of the estimated Fisher information for each incident plane wave used to generatethe reflection matrix (blue histogram). For comparison purpose, we show the distribution of noise estimates (redhistogram), obtained for each plane wave by taking two identical measurements of the outgoing state for ϕ = ϕ .We clearly observe that the signal is significantly larger than the noise, with a signal mean value of . rad − anda noise mean value of . rad − . Similarly, in Fig. S5b, we consider the case in which the observable of interest isthe lateral displacement of the circular object, and we show the distribution of the estimated Fisher information foreach incident plane wave used to generate the reflection matrix along with the distribution of noise estimates. Inthis case, the signal is also larger than the noise, with a signal mean value of . µm − and a noise mean value of . µm − . Supplementary Figure S5 | Histograms of Fisher information and measurement noise. a , Histogram of Fisherinformation for the plane waves used to construct the reflection matrix, along with the histogram of associated measurementnoise, when the observable parameter is the phase shift of the cross-shaped object. b , Analogous to a when the observableparameter is the lateral displacement of the circular object. We can verify that this signal-to-noise ratio is sufficient to faithfully estimate ∂ ϕ r and ∂ x r by measuring the derivativeof the outgoing field when the optimal incident state is used to illuminate the medium, and by comparing it to thederivative of the field predicted using the derivative of the reflection matrix. This characterization is experimentallyperformed by averaging the derivative of the outgoing field over measurements, with ND2 in the signal path toavoid saturation of the camera. The complex correlation coefficients between measured and predicted derivative ofthe field are . − . i (for the maximum information state relative to the estimation of the phase shift of thecross-shaped object) and .
99 + 0 . i (for the maximum information state relative to the estimation of the lateralshift of the circular object). These values show that we can correctly predict the derivative of the outgoing field fromreflection matrix measurements. Comparing the spatial distributions of the measured Fisher information per unitarea to the predicted ones, we can see that these distributions are very similar (Fig. S6). However, the total measuredFisher information is lower by a factor of the order of η f (cid:39) . when compared to the predicted Fisher information.This difference is largely explained by the observed difference in the predicted and measured outgoing intensity (forwhich a factor η i (cid:39) . was measured).2 Supplementary Figure S6 | Predicted and measured distributions of the Fisher information per unit area formaximum information states. a , b , Predicted spatial distribution of the normalized Fisher information per unit area foramplitude and phase (AP) modulation and phase-only (PO) modulation, when the observable parameter is the phase shift ofthe cross-shaped object. c , Measured spatial distribution of the normalized Fisher information per unit area. d – f Analogousto a – c when the observable parameter is the lateral displacement of the circular object. S4.3 – Correlation between the outgoing field state and its derivative
In order to characterize the correlation between | E out (cid:105) and | ∂ θ E out (cid:105) , we calculate the following complex correlationcoefficient: C θ = (cid:104) E out | ∂ θ E out (cid:105)(cid:107) E out (cid:107) · (cid:107) ∂ θ E out (cid:107) . (S45)This correlation coefficient is a relevant quantify to assess whether the Fisher information is enclosed in variations ofthe global phase of the outgoing state as expected for a unitary S -matrix ( C θ (cid:39) ± i ), or whether it rather enclosed inthe state’s intensity variations ( C θ (cid:39) ± ) or speckle decorrelation ( C θ (cid:39) ).We first calculate the correlation coefficient for each plane wave used to construct the reflection matrix, whenthe observable parameter is the phase shift ϕ induced by the cross-shaped object (Fig. S7a) and when it is thelateral position x of the circular object (Fig. S7b). On average, the measure correlation coefficients are equal to C ϕ = − .
01 + 0 . i and C x = 0 .
00 + 0 . i , respectively. For the maximum information states, the measuredcorrelation coefficients reach C ϕ = 0 .
01 + 0 . i and C x = 0 .
18 + 0 . i , respectively. Thus, in both experiments, weobserve that the average correlation coefficient for plane wave is close to zero, while the correlation coefficient forthe maximum information state is close to the imaginary unit. This observation is likely to reflect the invarianceproperty of maximum information states in the limit of a unitary S -matrix, in the same way as the principal modes ofa multimode fiber are invariant (to first order) in their output profile with respect to parameter variations except for aglobal phase shift [14–16]. However, in our case, the measured reflection matrices are not unitary, thereby explainingthe deviations from this property that are observed in the experiments. [1] Kay, S. Fundamentals of Statistical Processing, Volume I: Estimation Theory (Prentice Hall, 1993). Supplementary Figure S7 | Correlation coefficient between the outgoing field and its derivative. a , Complexcorrelation coefficient between the outgoing field and its derivative when the observable parameter is the phase shift inducedby the cross-shaped object. Black points represent values of the complex correlation coefficient for the plane waves usedto construct the reflection matrix, and the red point represent the value of the complex correlation coefficient for the optimalincident state. b , Analogous to a when the observable parameter is the lateral displacement of the circular object.[2] Helstrom, C. W. Quantum detection and estimation theory. J. Stat. Phys. , 231–252 (1969).[3] Braunstein, S. L., Caves, C. M. & Milburn, G. J. Generalized Uncertainty Relations: Theory, Examples, and LorentzInvariance. Ann. Phys. (N. Y.) , 135–173 (1996).[4] Demkowicz-Dobrzański, R., Jarzyna, M. & Kołodyński, J. Chapter Four - Quantum Limits in Optical Interferometry. In
Progress in Optics , vol. 60, 345–435 (Elsevier, 2015).[5] Mandel, L. & Wolf, E.
Optical Coherence and Quantum Optics (Cambridge University Press, 1995).[6] Ambichl, P. et al.
Focusing inside Disordered Media with the Generalized Wigner-Smith Operator.
Phys. Rev. Lett. ,033903 (2017).[7] Lerosey, G., Rosny, J. d., Tourin, A. & Fink, M. Focusing Beyond the Diffraction Limit with Far-Field Time Reversal.
Science , 1120–1122 (2007).[8] Xu, X., Liu, H. & Wang, L. V. Time-reversed ultrasonically encoded optical focusing into scattering media.
Nat. Photonics , 154–157 (2011).[9] Judkewitz, B., Wang, Y. M., Horstmeyer, R., Mathy, A. & Yang, C. Speckle-scale focusing in the diffusive regime withtime reversal of variance-encoded light (TROVE). Nat. Photonics , 300–305 (2013).[10] Bosch, J., Goorden, S. A. & Mosk, A. P. Frequency width of open channels in multiple scattering media. Opt. Express , 26472–26478 (2016).[11] Zhou, E. H., Ruan, H., Yang, C. & Judkewitz, B. Focusing on moving targets through scattering samples. Optica ,227–232 (2014).[12] Ma, C., Xu, X., Liu, Y. & Wang, L. V. Time-reversed adapted-perturbation (TRAP) optical focusing onto dynamic objectsinside scattering media. Nat. Photonics , 931–936 (2014).[13] Ruan, H. et al. Focusing light inside scattering media with magnetic-particle-guided wavefront shaping.
Optica , 1337–1343 (2017).[14] Fan, S. & Kahn, J. M. Principal modes in multimode waveguides. Opt. Lett. , 135–137 (2005).[15] Carpenter, J., Eggleton, B. J. & Schröder, J. Observation of Eisenbud–Wigner–Smith states as principal modes in multi-mode fibre. Nat. Photonics , 751–757 (2015).[16] Xiong, W. et al. Spatiotemporal Control of Light Transmission through a Multimode Fiber with Strong Mode Coupling.
Phys. Rev. Lett.117