Correlating light fields through disordered media across multiple degrees of freedom

Speckle patterns are inherent features of coherent light propagation through complex media. As a result of interference, they are sensitive to multiple experimental parameters such as the configuration of disorder or the propagating wavelength. Recent developments in wavefront shaping have made it possible to control speckle pattern statistics and correlations, for example using the concept of the transmission matrix. In this article, we address the problem of correlating scattered fields across multiple degrees of freedom. We highlight the common points between the specific techniques already demonstrated, and we propose a general framework based on the singular value decomposition of a linear combination of multiple transmission matrices. Following analytical predictions, we experimentally illustrate the technique on spectral and temporal correlations, and we show that both the amplitude and the phase of the field correlations can be tuned. Our work opens up new perspectives in speckle correlation manipulation, with potential applications in coherent control.


I. INTRODUCTION
The propagation of coherent light through a scattering medium results in a complex interference pattern generally called "speckle" [1].Speckle patterns are very difficult to predict and are affected by many experimental parameters: the medium itself, the wavelength or angle of illumination, and, in the case of pulsed light, the propagation time [2,3].Because of their complexity and apparent randomness, speckle patterns are often perceived as an inconvenience.Scattering is a major limitation for imaging both through weakly and strongly diffusive media.Examples of the former are ubiquitous in astronomy [4], and those of the latter in bio-medical imaging [5].
To overcome difficulties due to scattering, present-day techniques frequently take advantage of correlations such as in the so-called "optical memory effect" (ME), which facilitate image reconstruction [6][7][8][9] such as in speckle scanning microscopy [10,11].Resulting from a small perturbation of the speckle pattern [12,13], the ME is responsible for maintaining the specific shape of a given speckle pattern when tuning an experimental parameter within a small range (memory effect range).In a disordered system, however, a whole variety of different parameters can be tuned, resulting in many variants of the ME.The best known ME for scattering media is the tilt-tilt ME [7,14], generalized by [15] to tilt-tilt and shift-shift MEs.But the ME also involves other degrees of freedom as illustrated by the recent observation of spectro-axial MEs [16,17].The ME's limited range, however, drastically restricts its applicability.This range depends on the medium and is inversely proportional to its thickness so that it is often too small to improve imaging through thick scattering media.Several previous studies have characterized the ME [7], aiming to extend its range either by input or output light selection [18][19][20] or by engineering specific metasurfaces [21,22].However, little effort has been made so far to engineer it.
A versatile approach to the study of complex media is to use wavefront shaping tools and techniques.Using spatial light modulators (SLM), the incoming beam can be shaped to control the propagation of light.These techniques allowed introducing non-normally occurring properties [23] into speckle patterns as tailoring their statistics through iterative algorithms [24,25] or focusing light behind them [26].Another approach is to measure the transmission matrix (TM) [27] of the system.The TM makes it possible to determine the input state that focuses light behind the medium [28] or to modulate the energy delivery distribution using its singular modes [29].Non-local correlations of the output field are also possible by appropriate filtering of the Fourier components of the TM [30].Furthermore, transmission control with a TM, initially developed in monochromatic light, can even be extended to pulses [31,32].
Recently, wavefront shaping has been used to create adaptable and reconfigurable MEs.To our knowledge, three techniques have been presented for this purpose.Related to the conventional tilt-tilt ME, the authors of [33] construct an operator to customize the angular memory effect.In [34][35][36][37], eigenmodes of the Wigner-Smith operator are used to create speckle patterns resistant to frequency variations [38] -a manner of frequency ME.Finally in [39], the authors manage to relate the output fields when light is propagating through air or a scattering medium.One thing all of the techniques have in common is that they rely on a specific TM-basedoperator which combines the information for two propagation schemes (i.e., angle, wavelength, and disorder, respectively).They then use the eigenvectors of this operator as wavefront inputs.However, so far this approach is limited to a single type of correlation and two output fields.For the imaging applications based on ME mentioned above, an extension of the current control over the scattered field would improve the speed and accuracy of image reconstruction.In this article, we address the lack of a robust and unified method for inducing speckle correlations behind a complex medium.We present a general method for on-demand control/modulation of correlations between any number of speckle patterns.To do so, we take advantage of the singular value decomposition (SVD) of a linear combination of different TMs.We first present and apply this technique to the hitherto unexplored temporal domain by correlating the fields for two distinct delays in a pulse (Sec.II).We present an analytical model that allows us to predict the magnitude of the correlation.We highlight its correspondence with experimental data in a monochromatic scheme to illustrate the generality of the method (Sec.III).Moreover, taking advantage of the Fourier relation between time and frequency, we investigate periodic correlation in the temporal (spectral) domain that arises when generating the input state with information from spectral (temporal) measurements (Sec.IV).Finally, we generalize this effect by presenting the correlations obtained when mixing spectral and temporal information by calculating the input state from the singular value decomposition of the sum of a monochromatic TM and a time-gated TM.

A. Concept
Let us start by considering the following problem: one measures two transmission matrices T 1 and T 2 associated with two different experimental conditions.How should one shape the input field such that in both cases it leads to the same output?In a suitable mode basis this output field is described by a vector that results from a corresponding input field vector X: We consider Eq. ( 1) as the problem of finding a vector in the kernel of T 1 − T 2 .When the matrices T 1 and T 2 are unitary, this is equivalent to finding an eigenvector of T 1 T † 2 , as studied in [39].However, in the general case (non-unitary), our method assures that we will find the solution if it exists.We compute the kernel of T 1 − T 2 through its singular value decomposition, which even allows us to deal with rectangular matrices, and yields two output fields that are perfectly correlated (corresponding to a correlation value of 1, see Eq. ( 2)).

B. Experimental implementation
We first test this hypothesis to generate temporal speckle correlations.To this end, we send short pulses (100 fs) to a scattering medium (TiO 2 layer suspended on a glass coverslip) and measure the intensity on the sample output plane with a CCD camera (see Fig. 1(a)).To extract the scattered and elongated transmitted pulse, part of the incoming pulse is diverted to a delay-line that we interfere with the transmitted light following the approach of [40,41].By varying the length of the delay-line, we measure the interference for specific delays on the camera and apply filtering to access the temporal field information.Using this time-gated experimental setup (the details of which are presented in SM Sec.S1) together with a SLM to control the light impinging on the scattering medium, we can measure time-gated TMs and obtain temporal control over the transmission of an incoming pulse of light [32].The TM measurement is performed following the approach of [27], where the camera images are binned such that one speckle grain (defined as the width of the autocorrelation of the amplitude speckle) corresponds to one output mode.The fields then used for the analysis are binned in the same manner, with no significant difference with respect to non-binned fields, see SM Sec.S6.

C. Results
Correspondingly, we measure two time-gated TMs, T 1 and T 2 , at delays of t 1 =0.9 ps and t 2 =1.7 ps, respectively.Once the TMs have been measured and normalized with respect to their intensity, we calculate their difference (the normalization allows the two TMs to have the same weight whatever their measured delay in the pulse) and calculate the singular value decomposition of T 1 − T 2 .We then select the singular vector associated with the smallest non-zero singular value (corresponding to a normalized value s = s/ ⟨s 2 ⟩ of 0.1), denoted by v min and display its phase on the SLM to modulate the incoming field.We extract the temporal evolution of the pulse for this modulation scheme.The spatially averaged amplitude of the recorded field is shown in Fig. 1(b).To study the relationship between the recorded electric fields E(t), we correlate the fields measured for all delays with the field measured at t 1 and denote the correlation obtained C t1 .We use the correlation function where E i ≡ E(t i ) are complex fields (to get C t1 , t i = t 1 and t j varies).We plot in Fig. 1(c, top) the evolution of the absolute value of C t1 .For the plane wave input (black dashed line), the value reaches 1 for t = t 1 and decreases to reach 0 for a delay longer than the temporal speckle grain width.However, for the shaped input v min , another peak is observed at t 2 indicating increased correlation.To highlight the characteristic of v min , Fig. visualization of the correlations between all delays t, t ′ we plot in Fig. 1(d) a 2D excess correlation graph (simply named 2D correlation graph in the following) with axis (⟨t⟩ = (t + t ′ )/2, ∆ = |t − t ′ |).The normalized excess correlation value is given by a color scale and its position in the 2D graph indicates the relative and absolute values between t and t ′ .Only one correlation peak is visible at the expected position (i.e., t 0 = (t 1 +t 2 )/2 = 1.3 ps and δt = |t 2 − t 1 | = 0.8 ps) confirming that only the correlation between these two times is enhanced.Note that, the correlation between t 1 and t 2 is lower than 1, despite what one could expect from Eq. (1).This is due to the experimental limitation of phase-only control of the SLM (see Fig. 3(a)).

III. THEORETICAL EXPECTATIONS AND EXPERIMENTAL VERIFICATION
The above experimental results are consistent with the search for the kernel of T 1 − T 2 .However, output fields of singular vectors associated with small singular values have reduced transmission and are more subjected to noise.For our approach to be relevant for applications (e.g., imaging), it needs to combine high transmission and high correlation.We thus consider the field correlation obtained by modulating the incoming wavefront with singular vectors associated with singular values larger than zero and especially the largest one, v max .Significant correlation values are also observed in this case, as visible in Fig. 2. Thus, to understand the observations, more accurate modeling is needed.In the following, we present a general study of the eigenvalue problem for T 1 + e iα T 2 , where α is an arbitrarily fixed phase (recovering the previous case for α = π).It is important to keep in mind that this mathematical approach corresponds, from a more physical point of view, to determining the wavefront of the input field that favors specific interferences.In this case, the interference leads, behind the complex media characterized by T 1 and T 2 , to related fields.This is quite similar to the usual wavefront shaping experiments while here the research is complicated by taking into account not only the propagation information contained in a single TM but in two of them.The singular value decomposition of T 1 + e iα T 2 is equivalent to the eigenproblem of the following matrix product: Describing TMs as random matrices [27], we derive in SM Sec.S7 A based on Eq. ( 3) the expected correlation value as a function of the singular values s (X s is the singular vector associated to the singular value s), where γ ≡ N CCD /N SLM is defined as the ratio between the degree of freedom and the number of controlled modes.We begin by comparing the analytical model with simulations (including the experimental phase-only constraint) and experiments in Fig. 2 for no dephasing between the TMs.We set α = 0 and plot the absolute value of the correlation.The experimental data (blue dots) correspond well to the analytical predictions (black dotted line) as well as to the simulations (orange solid line).(For a rapid comparison, it is easy to verify that the correlation cancels when s = √ γ as predicted by Eq. ( 4).)A perfect correspondence is not expected because Eq. ( 4) is a prediction for full modulation of the wavefront.The combination of experimental errors and the use of only the singular vector's phase to shape the incoming light leads to a lower measured correlation.
In addition to giving access to the correlation value, Eq. ( 3) illustrates an important aspect of the method: the existence of two sets of terms.The first set (T † 1 T 1 +T † 2 T 2 ) corresponds to a control at the dedicated time-delays t 1 and t 2 .The second set, consisting of the cross information of the TMs, (T † 1 T 2 + T † 2 T 1 ) contains the coupling between T 1 and T 2 , which provides the key to control the correlations (see SM Sec.S4).As the induced correlations come from the cross terms only, it is possible to construct the correlation operator directly from this expression.Note that this is a symmetrized form of the approach of [39].
So far, we have focused on the absolute value of the correlations, independently of their phase.However, Eq. ( 4) highlights the possibility of controlling the phase α between the two fields.To explore this aspect, we performed field measurements and simulations by varying the relative phase α between two combined TMs. Figure 3(a) presents simulations in agreement with Eq. ( 4).The real part of the correlation as a function of the normalized singular values is plotted for different values of the parameter α.For α = 0 or π the correlation is real (with opposite signs) and therefore the real part is extreme, for intermediate values an imaginary part exists.The possibility of reaching a correlated or anti-correlated field is experimentally illustrated in Fig. 3(b) where two monochromatic TMs are measured at λ 1 = 806 nm and λ 2 = 810 nm.They are combined with α varying between 0 and 2π.For all values of α, the SVD is calculated and the phase of v max is displayed on the SLM to shape the input beam.The wavelength is smoothly tuned from 800 nm to 816 nm while recording the field.The real part of the correlation with the field measured at 806 nm is plotted as a function of λ.We observe a correlation when α = 0 (brown solid line) and an anticorrelation for α = π (blue solid line).All intermediate values are plotted in the complex plane of Fig. 3(c) and compared to an input plane wave (black dot).We have thus demonstrated that the SVD approach allows us to control the amplitude of the correlation by selecting the singular vector according to the associated singular value (radius of the circle in Fig. 3(c)) and the relative phase α of the fields by fixing the relative phase of the TMs (position on the circle in Fig. 3(c)) for a given parameter of interest (time-delay, wavelength etc.).
A practical aspect of the above SVD approach is the symmetric roles that the two TMs naturally play.Based on this symmetrized operator form, an extension to more than two field correlations is straightforward.To illustrate this point, we measured three time-gated TMs (measured for delays of 1.3, 2 and 2.7 ps), sum them and calculate the SVD of the new matrix T 1 + T 2 + T 3 .The 2D correlation plot corresponding to the sending of the singular vector associated with the highest singular value (v max ) is presented Fig. 4. One can observe three correlation peaks (corresponding to the number of combinations of two TMs in a set of three) located at (t 0 = (t 1 + t 2 )/2 = 1.65 ps, δt = t 2 − t 1 = 0.7 ps), (t 0 = (t 1 + t 3 )/2 = 2 ps, δt = t 3 − t 1 = 1.4 ps) and (t 0 = (t 2 +t 3 )/2 = 2.3 ps, δt = t 3 −t 2 = 0.7 ps).Slight increases are also visible for near-zero |δt| values, resulting from an intensity-induced artefact detailed in SM Sec.S2.Note that increasing the number of correlated fields unsurprisingly decreases the degree of correlations achieved.The derivation of the scaling of the correlation with the number of TMs is presented in SM Sec.S7 B.
It should be noted that the experiments are either performed with pulsed light correlating different delays (Fig. 1, Fig. 2 and Fig. 4) or monochromatic light taking the wavelength as a variable parameter (Figs.3(b) and (c)) to illustrate the claimed versatility of the approach.Furthermore, additional experiments are presented in SM Sec.S6, showing that fields can also be correlated after propagation through different disorders, generalizing [39] for any pairs of scattering media.In the following, we show how to further customize the correlation of the fields.

IV. CROSS-EFFECTS
We now take advantage of the particular link between time and frequency, which are conjugate quantities.The idea is to create spectral (temporal) speckle correlations merely using temporal (spectral) information.We expect such an approach to be useful from the imaging perspective when one domain is easier to access or faster to characterize than the conjugate one.The information is obtained here by measuring the monochromatic or time-gated TMs and by calculating the input state of interest as described above.However, we no longer measure the field directly by varying the wavelength or the delay, but we change the laser settings so that we measure the evolution of the field for the conjugate quantity.In the case of operating from time (characterization of the TMs) to frequency (measurement of the correlations) the procedure is as follows: (i) measure the two time-gated TMs at delays t 1 and t 2 ; (ii) compute the SVD of their sum; (iii) modulate the field according to the singular vector of interest (v max ); (iv) change the laser settings and measure the field and correlations varying the illumination wavelength.The experimental results for both cases (i.e., frequency to time and time to frequency) are presented in Fig. 5.The application of beam shaping results in a periodic modulation of the field correlations, well observable in terms of the checkerboard appearing in the 2D correlation plot presented in Figs.5(a) and (c).The period of the modulation is related to the spacing of the frequencies (delays) of the TMs and is given by δλ = λ 2 0 /(cδt) (δt = λ 2 0 /(cδλ)).These correlations of the fields are accompanied by amplitude peaks of the same period.The application of beam shaping results in a periodic modulation of the field amplitude, as can be seen in Figs.5(b) and (d).Note that the correlation checkerboard pattern is the generalization to many correlated delays of the peaks observed in the three delays experiment presented in Fig. 4 but with the advantage of retaining a high correlation value.A very interesting aspect of the generation of multiple temporal correlations using frequency measurements is therefore both the drastic reduction of the measurement time and the fact that we retain a high correlation value where a sum of several TMs would have led to its non-negligible reduction.The latter advantage is also valid for the reciprocal experiment.The effect of the relative phase α between the TMs for the cross experiment is discussed in SM Sec.S5.To further extend the possibilities of this technique, we present in the following experiments where we relax the constraint of investigating invariant modes of only one type of ME at a time.The combined TMs can now be measured in different experimental configurations (not only varying a parameter in one given configuration).We choose here to combine a monochromatic TM and a time-gated TM.The monochromatic TM is measured at λ = 808 nm and the time-gated TM is measured at a 1.3 ps delay in the pulse for the same central wavelength.Their sum is performed and the phase of the leading singular vector is displayed on the SLM to modulate the incident beam.Field correlations are extracted and presented in Figs.5(e) and (f) from a pulse scan (e) and a monochromatic scan (f).In both cases, in the 2D correlation plot, correlation lines are observed.We observe a "V-shaped" correlation whose tip lies at t 1 (resp.λ 1 ) for δt = 0 (resp.δλ = 0).This particular correlation indicates that all speckles at different delays (resp.frequencies) are all partially correlated with the speckle at t 1 .Indeed, in this experiment, as for the cross effect, the correlations are spectrally and temporally delocalized.Here, the sum of the two TMs relates the field information for wavelength λ 1 and delay t 1 .The interpretation is therefore as follows: for all delays, a part of the speckle, with wavelength λ 1 , is correlated to the speckle at t 1 .The measured slope of 2 is expected because the horizontal axis shows t 0 instead of the delay itself.This result is particularly interesting because it means that a spectrally selective correlation is feasible in the pulse.
Extending the field correlation with the SVD to any type of TM makes our approach very practical and promising for future applications.For example, one could consider combining a TM with a reflection matrix.Such a choice of matrices provides us with a new experimental tool to investigate how the electric field gets distributed and eventually stored inside a complex medium [42] by constraining the reflected and transmitted fields to be the same.

V. DISCUSSION AND CONCLUSION
We present here a very general method to correlate two or more fields behind a scattering medium with the knowledge of TMs.Our approach relies on the natural mixing of information that occurs in the SVD of a linear combination of TMs and on the sorting of the singular vectors by the corresponding singular values.We show that the resulting correlation control is broadly applicable, grants full control of the real and imaginary parts of the correlation, and is determined by the distribution of singular values.With the modulation of phase and amplitude, even a perfect correlation (or anti-correlation) is possible.Due to its generality, our approach applies to different illumination configurations such as pulsed light or different wavelengths in monochromatic illumination; also extensions to other experimental platforms such as multi-mode optical fibers [43] can be considered.Moreover, our approach is not limited to classical sources since all the techniques used here are also applicable to single photons [44].For example, maintaining and engineering correlations between entangled photons, usually scrambled in scattering media, would be extremely valuable for understanding fundamental questions or for developing quantum technologies [45].Finally, because of the many analogies and techniques shared with acoustics or microwaves, our approach should not be restricted to optics [46][47][48].This flexible and wide extension of the range of field correlation paves the way for engineering the memory effect at will, with all the practical implications that this entails.Most promisingly, the correlations are achievable independently of the natural ME range of the mediumwhich is typically very small for thick scattering samples, like the one used here.Taking advantage of the phase of the correlations may therefore prove to be an asset for imaging biological tissues, for example in the context of structured illumination techniques [49].A central point in the field of scattering through complex media is that wavefront-shaping techniques can turn these media into arbitrary optical elements [26,[50][51][52].
Here the joint control of the intensity and correlations of the field at the output of the medium goes further by opening the way to coherent control.Our technique makes it possible to create a source with adjustable scattering properties.When used with pulsed light, the structured pulse and repetition of the same pattern at fixed intervals are ideal for pump-probe experiments, where the ultrafast switching of a shaped wavefront is highly desirable [53].For these latter experiments, it is very interesting to note that the presence of the scattering medium allows the illumination to be modulated at a higher rate than the initial pulse repetition rate.measuring a single time-gated TM displaying one vector of its SVD on the SLM and looking at the field correlation for monochromatic light when varying the wavelength.

S2. IMPACT OF INTENSITY ON CORRELATIONS
As visible in the correlation formula given in Eq. ( 2) of the main text, the correlated fields are normalized to avoid intensity-induced artefacts.However, this does not fully prevent intensity changes from influencing the correlation measure.For example, imagine a system with a TM that continuously changes in time.For a random input, the output state of this TM will show the same temporal correlations as the matrix itself.Yet, if we calculate the output generated by the first singular vector of the TM at a given point in time t 0 we see that its temporal correlations are wider (Fig. S2(a)).Does this mean that the output of the singular vector is more stable in time?It does, but the effect is only due to the intensity enhancement associated with the first singular vector (Fig. S2(b)).As the output speckle at this point in time has a much higher intensity it also has a higher contribution at neighbouring times.This is something the normalization in Eq. ( 2) does not account for.To understand that this is solely caused by intensity variations, let us consider an even simpler example: A TM that is composed of only two components with linearly changing weights T ∝ T 1 (1 − t) + T 2 t.Also here, a random input decorrelates faster than the first singular vector of T 1 (Fig. S2(c)).However, if we construct an artificial random output with the temporal evolution r = 2r 1 (1 − t) + r 2 t, given that r 1,2 are random vectors, we observe the same temporal correlations as for the singular vector output.This shows that simply by matching the intensity evolution of the singular vector output state we could reproduce the same correlation feature (Fig. S2(d)).Therefore, here the enhanced correlations are only due to intensity differences and not inherent to the SVD.
In the measurements, this effect is visible for small δt or δλ, for example, in Fig. 4 or Fig. S4.However, it only acts on scales of the temporal decorrelation length and does not affect the long-range correlations investigated in this work.This can be seen in Fig. 1 of the main text.There, field correlations are observed between two different delay times in the pulse that show a slightly reduced intensity.This leads to a decreased width of the speckle auto-correlations around these times, i.e., the inverse of the effect explained in Fig. S2.In the 2D representation of the correlations presented in Fig. 1(c), this manifests as two dark spots around t 0 = t 1,2 but does not affect the SVD-induced correlations at larger δt.
Note that, the intensity artefacts described here occur in complex media where the speckle pattern and the transmitted intensity vary.In the case of multimode fibres, where the transmission is generally constant, we do not expect these effects to play a significant role [37].

S3. EVOLUTION OF THE EFFECT WITH THE TEMPORAL DELAY BETWEEN TMS
To illustrate the variation of the position of the peak in the correlation plot, a set of experiments is performed keeping the same central time (t 0 = 2.6 ps) but varying the delay δt between the two time-gated TMs.The results are presented Fig. S3.The correlation peak moves at a fixed t 0 along |δt|.

S4. DIFFERENCE WITH MASKS SUM AND CORRELATION TECHNIQUE
To emphasise the importance of the cross terms in Eq. (3) in the main text, we compared the correlations and pulse shapes of the following two situations.For Fig. S4(a), the usual procedure of this work is applied, i.e., we measure two time-gated TMs at two different times in the pulse, add them together and perform the SVD of the sum.We measure both the correlations and the pulse shape when the first calculated singular vector is displayed on the SLM.One observes intensity peaks at the position two times of the time-gated TMs visible both on the pulse shape and on the 2D correlation plot (see Sec. S2) for small |δt|.There is also a correlation peak corresponding to the correlation between the field at the two time-gated TMs times, thus appearing at t 0 = (t 1 + t 2 )/2 = 1.3 ps and |δt| = t 2 − t 1 = 0.8 ps.On Fig. S4(b) the SVD of the two time-gated TMs are performed separately and their 1st singular vectors are summed to create the vector whose phase is displayed on the SLM, as for the multi-time control [31,32].In the latter case, the intensity enhancement is still present (both in the pulse shape and the 2D correlation plot) but no further field correlation is observable.
In Fig. 1(a) of the main text, due to both the weak control given by the ratio between the degree of freedom and the number of controlled modes γ ≡ N CCD /N SLM ≃ 0.36 and the use of the singular vector associated to the smallest singular value, no clear intensity variation is visible for this specific input state (green curve) compared to the field obtained for a plane wave input (black curve).

S5. PHASE IMPACT FOR THE CROSS EFFECT
The effect of the relative phase α between the two TMs has already been studied and the main text (Sec.III).It is presented in the case of two monochromatic TMs in Fig. 3.Here we present its impact on the cross effect.In Fig. S5 as in Figs.5(c) and (d) we measure two monochromatic TMs and extract the singular vectors of their sum for different values of α.A pulse scan is performed to extract the field.Changing α causes the complete correlation plot to be translated along t 0 .To visualise this, one can look at the shape of the extracted pulse when performing the frequency-to-time cross effect.On the top part of Fig. S5 we track the positions of the maxima for all values of α.For α from 0 to 2π, all the maxima are shifted by one period.

FIG. S5
. Impact of α on the frequency to time cross effect.The bottom plot represents the shape of the pulse for α = 0.The top part represents the tracking of the maxima of the pulse while varying α: the maxima are shifted.The data come from the same experiment as in Fig. 5, thus averaged over 4 realizations of the medium.

S6. SPATIAL EFFECT
As stated in the main text, the method we present is very general and does not depend on the parameter used to vary the TM.We have worked mainly with time and frequency but other parameters can be used.To show this explicitly, we present here the same experiment performed for different spatial positions.We measure two TMs linking the plane of the SLM to two different regions of interest (ROIs) R 1 and R 2 on the camera.TMs are, as usual, summed and the phase of the first singular vector obtained is displayed on the SLM.We extract the field extract for a large ROI of the camera including the two regions where the TMs were measured.We display in Figs.S6(a) and (b) the amplitude of the fields with or without binning of the camera pixels.We calculate the correlations between R 1 and all other possible ROIs using a moving window, as shown in Fig. S6(c).An increase in correlation can be observed when the moving window coincides with R 2 .The correlation values obtained for the binned and unbinned images are similar.In the case of the unbinned image, the correlation value is slightly lower, this behaviour is expected as the image contains information about the high k vectors that is not contained in the TM itself.In Fig. S6 the geometry presented is simple because the scan is only performed on one dimension.However, the result generalises to a two-dimensional scan.

A. Analytical formula for the correlation
In this section, we derive a formula for the correlation between the output fields of two random transmission matrices, when sending a singular vector of a linear combination.Formally, let T 1 and T 2 be n-by-m matrices where all coefficients are drawn independently from the same (complex) Gaussian distribution of mean 0 and standard deviation σ.Fix two parameters α 1 and α 2 and define An important property is that our Gaussian distribution is invariant to rotations, such as the one transforming (T 1 , T 2 ) into (M, ∆).More precisely, coefficients of (T 1 , T 2 ) are independent Gaussian variables of mean 0 and standard deviation σ, thus coefficients of (M, ∆) are also independent Gaussian variables of mean 0 and standard deviation σ.Therefore, we draw M which is now a fixed matrix, but ∆ is still random.We write the singular value decomposition M = U ΣV † , where U ∈ C n×n and V ∈ C m×m are unitary matrices and Σ ∈ R n×m + is diagonal.Let X ∈ C m and Y ∈ C n be the right and left singular vectors associated with a singular value µ.More precisely, µ is the j-th coefficient of Σ for some j, and X and Y are respectively the j-th column of V and U .In particular, we have Rewritting T 1 X and T 2 X using M and ∆, we obtain Recall that ∆ is random Gaussian and that X and Y are fixed unit vectors.Therefore z is also Gaussian, and a short calculation shows it has mean 0 and standard deviation σ.
We are finally ready to write the correlation between T 1 X and T 2 X.
Each inner sum is a Gaussian random variable of mean 0 and standard deviation σ, thus the expected value of ||∆X|| 2  2 is exactly equal to nσ 2 .Notice that up to this point, all the formulae are exact.We will now approximate random variables by their expected value (which can be formalized with concentration inequalities) and say that ||∆X|| 2  2 ≈ nσ 2 .Thus, terms involving z are negligible, and we write To conclude, we define the aspect ratio γ = n/m.We normalize μ2 = µ 2 /⟨µ 2 ⟩, computing the average over n singular values (including zeros if γ > 1).More precisely Observe that in the formula, both X and μ are invariant when multiplying M by a constant (here √ 2).Hence, we can instead compute the singular value decomposition of e iα1 T 1 + e iα2 T 2 , which exactly corresponds to Eq. ( 4) from the main text when α 1 = 0, α 2 = α, n = N CCD and m = N SLM .Finally, Marchenko-Pastur's law says that μ2 varies between (1 − √ γ) 2 and (1 + √ γ) 2 , so correlations vary between e i(α1−α2) In particular, computing the singular value decomposition of ℓ j=1 e iαj T j corresponds to coefficients r 1,j = e iαj / √ ℓ in the first row of R, in which case the correlation formula can be simplified into Using random matrix theory, Marchenko-Pastur law shows that μ2 q ranges from (1 − √ γ) 2 to (1 + √ γ) 2 , and the correlation ranges from e i(αj

C. Random unitary matrices
Let us now consider two unitary random matrices T 1 and T 2 , drawn uniformly at random (which corresponds to the distribution induced by the Haar measure over the group of unitary matrices).As in the previous section, we compute the singular value decomposition decomposition of M = e iα1 T 1 + e iα2 T 2 = U ΣV † , where α 1 and α 2 are two parameters.Let X and Y be the right and left singular vectors associated with a singular value µ.More precisely, µ is the j-th coefficient of Σ for some j, and X and Y are respectively the j-th column of V and U .We are interested in the correlation between T 1 X and T 2 X, that is As a product of unitary matrices, C is unitary.We are now going to show that C is also a diagonal matrix, which will imply that its diagonal coefficients are complex numbers of modulus 1.Using both the definition of M and its singular value decomposition, one can write

Combining both equations gives e
This shows that the real part of Z = e i(α2−α1) C = Σ 2 /2 − I is a diagonal matrix.Because T 1 and T 2 are random unitary matrices, then all coefficients from Z are distinct with probability 1.We are going to show that the imaginary part of Z is also diagonal.First, observe that the matrices (Z + Z † ) and (Z − Z † ) commute: More precisely, this means that ℜ(Z) and ℑ(Z) commute, and that both matrices stabilize the eigen-spaces of the other matrix.If ℜ(Z) is a diagonal matrix with distinct coefficients, then each of its eigen-spaces has dimension 1 and is spanned by one of the vectors of the canonical basis.Thus, each vector of the canonical basis is an eigenvector of ℑ(Z), which in turn is diagonal.Hence, Z is diagonal, and because it is also unitary then diagonal coefficients must be complex of modulus 1.
Going back to our correlations, we showed that the correlation between T 1 X and T 2 X is a complex of modulus 1, such that ℜ(e i(α2−α1) C(T 1 X, T 2 X)) = µ 2 /2 − 1.Thus, we have To perform the simulations presented in the article, we used a simple approach modelling the TM by a random matrix [27,30].When quantitatively comparing simulations to experimental results, the remaining grain size after binning is added to the random matrix by convolving each output dimension with the adequate Gaussian.

S8. COMPARISON OF THE DIFFERENT OPERATORS
We present in this section simulation results that compare the correlation achievements obtained for different operators: the SVD of T 1 + T 2 (eigenvalue of (T 1 + T 2 ) † (T 1 + T 2 )), the eigen decomposition of T † 1 T 2 , of T † 1 T 2 + T † 2 T 1 (its symmetric version) and T † 1 (T 2 − T 1 ) (equivalent for the Wigner-Smith approach where (T 2 − T 1 ) corresponds to the derivative part [38]).We only consider the case of square non-unitary matrices.A more complete comparison goes beyond the frame of this work.For all four operators, the input vectors X s are sorted by decreasing the absolute value of the singular value/ eigenvalue (it is noteworthy that the eigenvalues are not necessarily real).The resulting absolute values of the correlation (between T 1 X s and T 2 X s ) are plotted as a function of the number of the singular vector (No. in Fig. S7(a)).
Only the SVD approach and the Wigner-Smith approach allow us to fully modulate the correlations (|C| = 1 in case of a square TM) while the other methods access a more restricted range with a lower maximum correlation observed at 0.8 with the definition of Eq. ( 2).However, it should be noted that perfect correlation is expected for the lowest singular value, which is the most sensitive to noise or phase-only modulation.In the current experimental capabilities, the latter point diminishes the specific advantage of the SVD method, but continued progress in beam shaping control will overcome this purely experimental limitation.Using the Wigner-Smith operator grants similar results to the SVD approach while losing the benefit of the sorting due to real singular values, which is of practical interest.It also lacks easy and independent control of the absolute value and the phase of the correlation.
The methods of [39] and its symmetrized version give similar correlation values.The main difference between the latter two methods is visible in the transmitted intensities presented in Fig. S7(b).The ratio of the total transmitted intensities is displayed and one can observe that all symmetrized versions lead to balanced intensities (i.e., the symmetrized version of [39] and SVD method).This observation is not surprising as the two TMs play non-symmetric roles.

S9. EXTENSION TO MULTIMODE FIBRES
We also performed simulations to predict the results of the technique for unitary systems.A good example of a complex but lossless platform is multimode fibers (MMF).In addition, the finite number of propagation modes makes it possible to measure complete TMs.We performed the simulations using the code of [54] which solves the transverse scalar propagation equation and allows to calculate the modes for fibers with arbitrary index profiles and their TMs.We have applied it to measure the TMs of a 10 cm long step-index MMF with 0.22 numerical aperture and 25 µm radius illuminated with a light of either 700 nm or 800 nm.Due to the variation of the wavelength, the number of propagation modes is different.However, the code allows them to be expressed on the pixel basis, as in the experiments, so that their size is the same, allowing them to be added together.Thus, once the zero singular values corresponding to the rank difference between the two TMs are removed, the correlation results presented in Fig. S8 are obtained.S3) is presented by the black dashed line.The horizontal black dashed line represents the mean correlation for a random input vector and the grey-shaded area is its standard deviation.The data are not averaged over several realizations.
Two main points are worth noting: the absolute value of the correlation, equal to unity, does not depend on the singular vector and the real part of the correlation well follows the law predicted in Eq. (S3).Thus even if the resulting fields are perfectly correlated, the unitarity of the transformation prevents any variations.For this reason, losses and non-unitary transformations are often sought [55].

FIG. 1 .
FIG. 1. Experimental illustration of temporal speckle correlations engineering.(a) Schematic of the experimental setup.The light delivered by a laser is split into two paths with a half wave plate (λ/2) and a polarized beam splitter (PBS): (i) signal where the wavefront is modulated by a reflective phase-only SLM and passes through a layer of TiO2 (focused with microscope objectives of 0.4 numerical aperture (NA)); (ii) delay line introducing a delay τ .Both paths are recombined with a beam splitter (BS) and their interference, cleaned from the unscattered light with a polarizer (P) is imaged on a CCD camera.(b) Averaged temporal amplitude of the field measured behind a thick layer of TiO2 when illuminated by a laser pulse of 100 fs.The black curve corresponds to the field when a blank pattern is displayed on the SLM whereas the green curve results from the display of vmin (singular vector associated with the smallest singular value of T1 − T2).The time-gated TMs are measured with NSLM ≈ 620 modes on the SLM and NCCD = 225 pixels on the CCD-camera at t1 and t2.The camera pixels are always binned to have one speckle grain per pixel, except in the images of the field amplitude presented in the insets (for visualization purposes).(c) Top: Temporal evolution of the absolute value of the field correlation with respect to t1 for both the blank SLM pattern (black dashed curve) and for vmin.Bottom: subtraction of the reference correlation to highlight only the anomalous correlation increases.(d) 2D representation of the correlations to observe simultaneously all the correlations (with subtraction of the reference as in the bottom of Fig. 1(c)).For symmetry reasons, only ∆ = |t − t ′ | i.e., spacing between the two correlated delays, is represented as a function of ⟨t⟩ = (t + t ′ )/2, which is the mean value of the two correlated delays.A correlation increase is observed for ∆ = |δt| = |t2 − t1| and ⟨t⟩ = t0 = (t2 + t1)/2.Data are averaged over 4 realizations of the medium.

FIG. 3 .
FIG. 3. Impact of the relative phase α between the TMs on the field correlations.(a) Simulation using 1024 × 1024 random matrices.Plot of the real part of the correlation for many values of α for phase and amplitude control (solid line) and phase-only control (dotted line).The data are averaged over 4 matrix realizations.(b) Experimental result of two-wavelength correlations for α = 0 (brown), α = π (purple), and blank input pattern (dotted black).TMs are measured for NCCD = 225 and NSLM ≈ 550, the input vector is vmax and the data are averaged over 4 realizations of the medium.(c) Imaginary part of the correlations as a function of the real part for a set of values of α for the same experiment as in (b).The central black dot corresponds to the blank input pattern.

FIG. 5 .
FIG. 5. Experimental observation of cross-effects.(a,b) Time to frequency: measurement of two time-gated TMs at t1 = 0.9 ps and t2 = 1.7 ps for observing the spectral variations when displaying vmax and tuning the wavelength.The expected period of the modulation is δλ = λ 2 0 /(cδt) = 2.7 nm. (a) 2D correlation plot where (b) represents the amplitude modulation A1 of the singular vector divided by the amplitude A ref when sending a blank SLM pattern.The grey vertical dashed lines illustrate the expected spacing of the peaks.(c,d) Frequency to time: measurement of two monochromatic TMs at λ1 = 805.5 nm and λ2 = 810.5 nm and for temporal field control when the sample is illuminated by pulses.The expected period of the modulation is δt = λ 2 0 /(cδλ) = 0.4 ps (see vertical dashed lines).(c) 2D correlation plot where (d) represents the amplitude modulation (the black dotted line corresponds to the reference amplitude).In both cases (frequency to time and time to frequency) the TMs are measured with NSLM ≈ 2180 and NCCD = 225 with the data being averaged over 4 realizations of the medium.(e,f) Extension of the correlations by varying different parameters.Two TMs with NSLM ≈ 2220 and NCCD = 225 are measured, one is time-gated and one is acquired in monochromatic illumination.2D correlation plot for (e) a temporal scan and (f) a monochromatic scan when displaying vmax.The effect being weaker than for the observations in (a) and (b), the data presented are averaged over 8 realizations of the medium.

FIG. S2 .
FIG. S2.Correlations induced by intensity variations.The TMs are of size 256 × 256 and sampled every 0.1-time unit (a,b) or 0.03 (c,d).(a) Field correlations between t = 0 and other times for a random input (blue) and the singular vector associated with the largest singular value of the TM measured at t = 0 (1st SVD, red).(b) Total intensity for the same data as in (a).(c) Field correlations for a random input (blue), the singular vector associated with the largest singular value (red) and an artificial random input having the same time evolution as the TMs (dashed yellow).(d) Intensity for the same data as in (c).The data are averaged over 20 realisations of the disorder (a,b) and 10 realisations (c,d).
FIG. S3.Evolution of the appearance of the correlation plot when the delay time between the two measured time-gated TMs (NCCD = 225 and NSLM ≈ 680) is modified from 0.4 ps (a) to 3.2 ps (d) while keeping the same central position.The correlation peak appears at a fixed t0 but with a |δt| that varies with the relative position of the time-gated TMs.The full correlation is only plotted in (b) whereas on the other graphs, only the area of the dashed red rectangle is shown.The data are not averaged.
FIG. S4.Importance of the SVD terms coupling the two TMs.(a) SVD of (T1 + T2) and (b) coherent sum of the phase masks obtained from the SVD of (T1) and the SVD of (T2) performed separately.The top parts represent the pulse shape in both cases.The intensity-induced correlations are visible on the 2D correlation plot at the bottom.No further correlations are present in the coherent sum case (b) whereas one appears for |δt| = 0.8 ps in (a).TMs are measured NCCD = 225 and NSLM ≈ 570.Data averaged over 4 realizations of the medium.

FIG. S6 .
FIG. S6.Spatial field correlation (a) Unbinned image of the speckle.The TMs are measured on the left (R1) and central (R2) parts for λ = 808 nm, NCCD = 225 and NSLM ≈ 610.(b) Same speckle image but binned so that one pixel maps one speckle grain.(c) Correlation of the left ROI (R1) to all the other positions with a moving window.An increase in correlation is observed when the moving window reaches the position where the other TM was measured.Data averaged over 4 realizations of the medium.

FIG. S7 .
FIG. S7.Comparison of different methods for correlating fields (a)The absolute value of the correlation is plotted as a function of the input vector number (ordered by decreasing the absolute value of the singular/eigenvalue).The simulation is performed with TMs of size 1000 × 1000.(b) Ratio of the total transmitted intensities for the same data as in (a).

FIG. S8 .
FIG.S8.Simulation of field correlation for multimode fibers.Real part (blue dots) and absolute value (orange dots) of the correlations when summing two unitary TMs.The analytical prediction of Eq. (S3) is presented by the black dashed line.The horizontal black dashed line represents the mean correlation for a random input vector and the grey-shaded area is its standard deviation.The data are not averaged over several realizations.