Spatiotemporal focusing through a multimode fiber via time-domain wavefront shaping

We shape fs optical pulses and deliver them in a single spatial mode to the input of a multimode fiber. The pulse is shaped in time such that at the output of the multimode fiber an ultrashort pulse appears at a predefined focus. We shape a set of pulses corresponding to a spatial grid of such nonlinear foci. Our result shows how to raster scan an ultrashort pulse at the output of a stiff piece of square-core step-index multimode fiber and in this way the potential for making a nonlinear fluorescent image of the scene behind the fiber, while the connection to the multimode fiber can be established via a thin and flexible single-mode fiber. The experimental results match our model well.


I. INTRODUCTION
All-optical imaging via multimode fibers (MMF) has the potential to become the method of choice for imaging in confined spaces, combining the smallest access diameter with the highest NA [1,2]. The most important application is minimally invasive endoscopy, but other use cases such as product inspection in an industrial setting are notable as well [3].
Multimode fibers support a large number of optical modes and hence transmit patterns from their input to their output facet. However, complex multimode interference makes it challenging to reconstruct the original input; bending of the MMF scrambles the multimode interference. Different methods have already been investigated overcoming this, such as spatial wavefront shaping [4][5][6], machine learning [7], and compressive sensing [8]. However, so far most MMF imaging methods are based on linear scattering or absorption.
In free-space microscopy, a plethora of special imaging modalities have been devised exploiting nonlinear imaging with ultrashort pulses. Nonlinear microscopy has multiple advantages: it reduces out-of-focus background and phototoxicity, allows to initiate highly localized photochemistry in thick samples, and provides optical sectioning that results in higher sensitivity and 3D imaging ability [9]. Considerable efforts have been put into the development of nonlinear endo-microscopy methods [10,11]. Unfortunately, combining ultrashort pulses with MMF imaging is non-trivial, as the modal interference and modal dispersion in a MMF results in a complex spatiotemporal output field [12,13]. Recently, several nonlinear optical imaging techniques through a single MMF probe have been demonstrated including two-photon excitation microscopy [14,15], 3D microfabrication based on two-photon polymerization [16], and coherent anti-Stokes Raman scattering (CARS) microscopy [17]. All these methods of nonlinear imaging require spatial-domain wavefront shaping and consequently control of many spatial modes on the MMF input.
Here we propose a new approach for imaging through a single MMF probe. We 'focus' light at any point on the distal fiber facet by using a single input mode utilizing light scrambling in an MMF, pulse shaping in time and a nonlinear optical detection. Our system allows control over the position of a nonlinearly focused beam in space on the MMF output facet by shaping an input pulse in a single spatial mode in time. We experimentally demonstrate grid scanning an ultrashort pulse over the output facet of a stiff piece of the MMF by temporally re-shaping the pulse. This new way of light control at the MMF output can help to avoid the perturbation sensitivity of MMF-based imaging probes.

II. THEORETICAL DESCRIPTION
With continuous-wave (CW) light, it is possible to spatially shape the input field of a MMF in such a way that a focus appears at the output facet of the fiber. However, timedomain shaping is necessary in order to allow the input to travel in a single spatial mode.
The output field of an MMF for a broadband pulsed input is also time-dependent, which can be exploited to do time-domain wavefront shaping. The principle is illustrated in Fig.   1. Two spots at the output, A and B, of the MMF are assumed to have an independent temporal response to a transform-limited input pulse. Inverting one of the responses and using that as the input pulse shape results in a transform-limited pulse in either spot A or B, depending on which response was inverted. This enables the output of a short pulse in many possible output spots, even though the input pulse is still in a single spatial mode.
We will now elaborate on why the output is time-dependent, how time-domain wavefront shaping is defined and how the system can be modelled.

Time-dependent output patterns
A multi-mode fiber supports many modes, each with its own propagation constant. For a focussed CW input, many of these modes are excited, with their amplitude given by the overlap integral of the mode field and the input field. Due to the different propagation constants, however, the modal amplitudes do not stay in phase after travelling through the fiber. As a result, the output field is a superposition of mode fields with approximately random phases, which leads to a speckled output field. The propagation constants and mode fields are in general also frequency-dependent, which leads to the desired time-dependent output field. Without taking into account polarization or mode mixing, the time-dependent output field of a MMF of length L with N modes is described by where ω are (discrete) optical frequencies [PP: at this moment not clear how these are defined/picked and why this is not an integral over the spectrum], β n (ω) the propagation constants, A n (ω) the initial mode amplitudes, and Ψ n (x, y) the (orthonormal) mode fields.
For simplicity we have assume that the mode profiles are not frequency-dependent. We take the input field to have a constant amplitude C(x, y), but with a phase shift θ(ω), so that E in (x, y, ω) = exp(iθ(ω))C(x, y). We can therefore approximate the mode amplitudes as Time-domain wavefront shaping By altering θ(ω), we can change the output field in time and target a specific output location to produce an ultrashort pulse there. In the center of the fiber for example, the output field at t = 0 is given by Eqs. (1) and (2): By setting the phase shifts to we have which is a strong peak due to the coherent sum. In general, the argument (and the amplitude) of the inner sum in Eq. (1) varies rapidly with x and y. Since θ(ω) is fixed and independent of x and y, the sum over all frequencies for positions away from (x, y) = (0, 0) is incoherent and the output is therefore not peaked in time there. To produce a peaked pulse in time at an arbitrary position (x, y) = (X, Y ), we can simply set In an experimental setting, however, the exact propagation constants are not known, and mode mixing further complicates the propagation through the fiber, so that the required phase shifts cannot be calculated a priori. Instead, the phase shifts can be optimized using an iterative algorithm.

Enhancement
An important parameter in wavefront shaping is the enhancement, defined as the ratio of (nonlinear) intensity in the wavefront shaping region after shaping and the (nonlinear) intensity for the same input wavefront but for a different, random sample. Similar to (linear) spatial wavefront shaping, we expect that the enhancement scales linearly with the number of controls in the pulse shaper and that the enhancement should go to zero for zero controls.
This expectation is perhaps counterintuitive, as the enhancement in our case can only be measured nonlinearly. But, even though the intensity peak in time grows quadratic with the number of frequencies and thus with the number of controls, the width of this peak shrinks linearly as well. This should result in a linear enhancement increase in a time-averaged nonlinear measurement after time-domain wavefront shaping. Mathetical support for the linear scaling of the enhancement can be found in Supplement 1, section S1.

Numerical model
We use a numerical model of the square core fiber for testing time-domain wavefront shaping algorithms and to simulate time-domain behaviour that cannot be measured in the lab. The full details of the model can be found in Supplement 1, section S2.

III. EXPERIMENTAL DETAILS
To experimentally verify the principle of time-domain wavefront shaping, we use the setup as illustrated in Fig. 2. The output of a mode-locked Ti:Sa laser with 13 nm bandwidth (≈100 fs) pulses, centred at 800 nm (Spectra Physics Tsunami, 80 MHz) is shaped in time with a 4f pulse shaper [18]. The pulse shaper can shift the phase of 640 pixels with a spectral resolution of 0.064 nm/pixel, and is calibrated using a spectrometer [19]. The output of the pulse shaper is focussed into a multi-mode fiber, after which the average output power is 50 mW.
The speckle-like output pattern cannot be directly measured in the time-domain because of the short timescales. To still detect temporal behaviour indirectly, the output pattern is imaged with a nonlinear imaging method. To this end, the output facet of the MMF is imaged into a 50 µm thin cuvette filled with a two-photon fluorescent medium (Rhodamine 6G in ethylene glycol). This medium does not have linear absorption for 800 nm pump light, but can absorb two 800 nm and emit a green fluorescence photon [20]. This two-photon process is sensitive to the square of the instantaneous optical power, so temporal compression can be made visible. The pump light is removed with a short-pass filter and the weak fluorescence is imaged with an EMCCD camera with high gain (Andor iXon DV885). Swapping the short-pass filter for an ND filter allows for linear imaging of the output intensity.

Optimization algorithm
As explained in the theoretical section, in the experiment the phase shifts that correspond to an ultrashort pulse at a specific output position are difficult to determine a priori. Instead, the phase shifts are found with an optimization procedure. On the camera, a circle with a 7-pixel radius (≈ 2.5 µm) is placed around the desired output position. In every step of the algorithm, 160 of the central 320 pulse shaper pixels are randomly selected and shifted relative to their current phase from 0 to 2π in increments of π/4. At each shift, the average nonlinear intensity in the circle is recorded. This intensity will show a sinusoidal variation with phase shift. A sine curve is therefore fitted and the phase shift that gives maximum nonlinear intensity is now set to each of the 160 pixels. This optimization of the phases of a random sample of 160 pixels is repeated for 5000 steps. This algorithm is inspired by spatial wavefront shaping, where it is known that this type of algorithm gives a good signal-to-noise ratio in determining the optimal phases [21]. The phases are initially set to random values, so that the algorithm is more likely to find a global optimum.
A complication is that the resolution of the nonlinear imaging method is poor, which means that multiple spots can be enhanced simultaneously because light and fluorescence from spots slightly outside the circle still leaks into the circle and thus contributes to the average nonlinear intensity in the circle. To improve the contrast of the optimization, we instead optimize on the ratio between the nonlinear intensity inside and outside the circle.
Only in the beginning we optimize the intensity directly and slowly change into contrast enhancement, by gradually changing the optimization parameter from intensity to ratio.
After 3000 steps, the optimization is fully based on the ratio.

Square-core multimode fiber and numerical model
The multimode fiber in the experiment is a 70 × 70 µm 2 square-core fiber (Ceramoptec, 0.22 NA). We found that a square core fiber has a flatter intensity profile at the output and shows less correlation between the input location and the output pattern in comparison to round MMFs, which indicates more mode mixing and better excitation of higher-order modes.
To characterize the frequency dependence of the output patterns, a tunable CW Ti:Sa laser (Coherent MBR-110) is used as input. The wavelength is scanned with a step size of 0.1 cm −1 (7 pm), limited by the resolution of our wavelength meter (Burleigh WA-10L).
The speckled (linear) intensity images are correlated with the output pattern at 799.50 nm with the Pearson correlation coefficient. The experimental intensity pattern at 799.50 nm is shown in Fig. 3(a), and the correlation coefficients are plotted in Fig. 3(b), together with the prediction from the model. The patterns fully decorrelate after a wavelength shift of approximately 40 pm, which matches remarkably well with our zero-free-parameter model.
The fact that the measured decorrelation width matches the model so well gives confidence that we understand the modal dynamics of this square-core fiber good enough to proceed using it for our imaging method. all. This is likely due to the difficulty of optimizing the pulse form with very few controls, especially in the presence of noise. This small offset is less relevant at a high number of controls. However, if we account for this offset by adding its absolute value to the measured enhancement data, the enhancement does scale linearly with the number of controls. The slope of this linear relationship is much lower in the experimental results, which we explain by the large amount of noise in the nonlinear imaging method.

Temporal compression
As stated before, due to the very short timescale of the output intensity pattern, it is difficult to experimentally measure temporal compression at an output spot directly.
However, we can use our numerical model and simulate time-domain wavefront shaping with it. Fig. 6(a) shows a magnification of the modeled nonlinear intensity both before and after time-domain wavefront shaping. The linear intensity over a time period of 4 ps in the highlighted regions for both cases is shown in Fig. 6(b). Before optimization, both regions show a random and broad distribution of light in time, which is expected due to the random phase-shifts of the frequencies that are present at the spots. After optimizing in region A, however, region B still shows a similar trace but region A now shows a high and narrow pulse of light. This confirms the idea we sketched in Fig. 1, namely that both spots have an independent temporal response, which can be selectively compressed by finding the corresponding optimal input pulse shapes.

Raster scanning
In order for time-domain wavefront shaping to have applications in nonlinear endoscopic imaging, a single optimization position is not sufficient. The simplest way to scan an ultrashort pulse over the entire output facet of the MMF is to define an optimization grid with many points and optimize the nonlinear intensity for each point individually. The variations in intensity enhancement are likely due to noise in the nonlinear imaging method, which can make it difficult for the algorithm to precisely determine the optimal phase. This effect of reduced enhancement due to noise is also known in spatial wavefront shaping [21]. For the model, we simulate shot noise with a similar amplitude as in the experiment, but if we artificially increase the noise further it can happen that no enhancement is ever found, which further supports this reasoning.The positional variations are likely due to the large optimization region we use in the algorithm. A spot can start to get enhanced anywhere in this region, which causes random variations in the final focus position. The contrast between an optimized spot and the background is much lower in the experiment, but this is likely caused by out-of-focus nonlinear fluorescence in the cuvette.
The grid optimization shows that it is possible to grid scan an ultrashort pulse over the output facet of a multi-mode fiber. Being single-spatial mode, the input pulse shapes for each position are not sensitive to spatial perturbations of its delivery channel. Therefore, if the MMF is stiff, the different shapes can be stored and reused many times. Moreover, the fluorescence from a single grid position can in principle be collected back through the same MMF. Time-domain wavefront shaping therefore enables deterministic and robust grid scanning of an ultrashort pulse over the output facet, and has applications in nonlinear endoscopic imaging.

V. CONCLUSION AND OUTLOOK
We have demonstrated spatial grid scanning of an ultrashort pulse at the output facet of a square-core multimode fiber by only changing the temporal shape of an ultrashort pulse at the input facet. The results match well with our numerical model, which can be used to directly show temporal behaviour at the ultrashort timescales. Delivery of the input pulse is done in a single spatial mode, which is insensitive to spatial perturbations. Our method has applications in nonlinear endoscopic imaging, since the pulse delivery and corresponding response readout can be done via the same fiber.

I. ENHANCEMENT DERIVATION
Ordinary scattering media are typically modeled with N spatial input modes and a complex-valued transmission matrix t mn that connects the field of the n th input mode to the m th output mode [1,2]. With the n th input field written as E in n = A n e iφn , we have Our model for a time-domain medium is analogous to this. We take a single mode as input, and M spatial output modes. Furthermore, we assume a discrete set frequencies Ω of size N and spacing δ.
Again, we use a complex matrix t mn to connect the input field of the n th frequency mode to the m th spatial output mode. With the input field of the n th frequency mode as E in n (t) = A n e i[(ω 0 +nδ)t+φn] , the output field in the m th spatial mode is given by Both the input field and output fields are time dependent and 2π/δ-periodic. For each ω ∈ Ω, we assume the transmission to be independent and random in phase for each output mode.
For simplicity, we consider only random phase, but fixed amplitude for each spatial output mode. Furthermore, the total transmission is taken to be unity. Under these assumptions, where f Θmn (θ mn ) is the probability density function for θ mn .
It is easy to see from equation (1) that when the input is phase shaped such that φ n = − arg(t mn ), the amplitude of E out m is maximized. This is called wavefront shaping [2]. Similarly, we can maximize E out m (t = 0) by setting φ n = − arg(t mn ) = −θ mn , as then all the frequency components are in phase (see equation (3)). Since the input field is in a single spatial mode and only depends on time, we call this time-domain wavefront shaping. For a non-shaped input (φ n = 0), the t mn coefficients are essentially the Fourier coefficients of the output signal. Letting φ n = − arg(t mn ) is the same as taking the unshaped output signal and conjugating its Fourier coefficients. Due to the properties of the Fourier transform, this is equal to the (complex conjugated) time reversed signal. This effect has an intuitive explanation, since taking the exact reverse of a time-scrambled signal and putting it through the scrambler again will unscramble it [3].
In wavefront shaping, the most important figure of merit is the enhancement η, which is defined as the ratio of the intensity in the shaping region after optimization, I N , and the intensity in the same region with the same optimized input, but ensemble-averaged over all possible samples, I 0 . So, where ... denotes the ensemble-averaged expected value. In spatial wavefront shaping, assuming circular complex Gaussian random t mn [2], The enhancement thus scales linearly with the number of controlled modes N .
We now derive the enhancement for time-domain wavefront shaping using the model defined above. Using equation (4), we find three useful ensemble-averaged expected values: 2π 0 e iθmn dθ mn = 0, Without any phase conjugation and with unity amplitude input E in n (t) = e i(ω 0 +nδ)t , the expected value for the output intensity in the m th output mode is given by If we phase conjugate for output mode j, then E in n (t) = e i[(ω 0 +nδ)t−arg(t jn )] = √ M t * jn e i(ω 0 +nδ)t . For any output mode k = j, where we have used the fact that t kn and t jn are independent for k = j. It is logical that equations (8) and (9) are the same if all t mn are uncorrelated, since averaging with an unshaped input wavefront is then the same as averaging with a random input wavefront.
For the output mode j, we find Based on the definition of the enhancement and using equations (9) and (10), we write Since I N (t) = N 2 /M for t → 0, the maximum enhancement in time is N . However, the temporal features in the input and output fields are of the order ∆t ∼ 1/N δ. For a physical system, the total bandwidth N δ can be several THz, giving temporal features in the femtosecond regime. This makes a direct, time-resolved measurement of I N (t) virtually impossible. Any physical detector will thus likely measure a time-averaged signal.
To emulate a time-averaged measurement, we only need to time-average over a single period 2π/δ, because the input and output fields are periodic in time. A linear detector will detect signals proportional to Both signals are the same, hence it is impossible to perform time-domain wavefront shaping with the feedback of a time-averaged linear detector. Because we can only shape an input pulse in the time domain, we cannot increase the average output energy in a spatial output mode. It is possible to use a linear detector and an interferometric measurement to reconstruct the output fields, but this method is slow and therefore not suitable for direct feedback for wavefront shaping.
A possibility to get a feedback signal to base our temporal wavefront shaping on is the use of a non-linear detector. Let's assume such a detector is sensitive to I 2 (t), allowing to detect signals proportional to With this detector, the enhancement becomes The enhancement expression is very similar to the result for spatial wavefront shaping (equation (6)). For large N , η ≈ 2N/3, so the enhancement scales linearly with N . Since I N (0) 2 /N 2 ∝ N 2 , one might expect the enhancement to grow quadratic and not linear.
However, the width of the central peak must decrease linearly with N because of energy conservation, so S 2 N /S 2 0 is linear in N . This effect is well-known in non-linear detection, where for constant average power the signal scales inversely proportional to the pulse width [4].
Our experimental implementation of a non-linear detector is discussed in the main paper.
So far we have assumed the time-domain medium to be loss free. In case it is is lossy, the modeling assumption that all |t mn | = 1/ √ M will no longer be valid. Making the amplitude of t mn also randomly distributed will change the average values from equation (7). This may alter the theoretical enhancement for small N , but for large N the enhancement should still be linear in N , which is the most important result of this section.
A complication in the experiment is that we cannot naturally vary the true number of independent frequency channels, as the resolution of our pulse shaper is fixed. Instead, we emulate smaller N by binning together pixels on the pulse shaper SLM. This does not result in a lower number of frequency channels in the sample, but does result in fewer controllable frequency channels, effectively reducing N .

II. SQUARE CORE FIBER MODEL
In the following we detail our model for the square core fiber. We will first describe the transverse fiber modes, then the propagation through the fiber and finally discuss the parameter choices to mimic the actual fiber used in the experiment.

A. Mode profiles
We only consider a single polarization (horizontal) in the square core fiber, for which one component of the electric field is given by [5] Here, for mode numbers p and q (both 1, 2, ...). Note that the coordinate system origin is centered on the square core fiber. The wave numbers can be found with the transcendental equations and Here, a is half the width of the fiber (i.e. the fiber is 2a-by-2a), and n co and n cl are the core and cladding refractive index, respectively, and k the wave number of the light.
Equations (17) and (18) are easily solvable using a few iterations of Newton's method, and the maximum values for p and q are those that still give a real solution for all wave numbers.
Finally, the propagation constant β inside the fiber can be found with If we neglect γ x and γ y , all wavelengths would have the same mode profile, but their propagation constants would be different. It is precisely this difference in propagation constants that gives rise to phase shifts and independent speckle patterns for a broadband input.

B. Mode propagation and mode mixing
Having found all possible values for p and q, there are a total of max(p) * max(q) modes, which can be enumerated with a single index n. Without any coupling between the modes, the evolution of the mode amplitudes c n is described by the matrix equation The initial amplitudes are found by calculating the overlap integral of the mode electric field profiles and the input electric field.
In a real multimode fiber there is mode coupling between the modes due to bending and refractive index variations [6]. The coupling due to a single bend with radius r may be described by dc n dz = iβ n c n − n co kξ r m A nm c m , where ξ is a correction factor (0.77 for silica), and A nm = E x,n | cos(θ)x + sin(θ)y|E x,m are the overlap integrals between modes n and m for a bend with projected angle θ with respect to the x-axis. Equation (21) can be solved as a matrix differential equation by using eigenvalue decomposition.

C. Implementation
Our fiber is modelled as a 1 meter long, 70-by-70 micron square core fiber, with a numerical aperture of 0.22 and a cladding refractive index of 1.4533 (assuming pure silica at 800 nm light). No dispersion of the numerical aperture or refractive index is taken into account, only dispersion due to equation (19). In the lab, the fiber is strongly bent by winding it around a cylinder a few times and tie-wrapped for stability. Because the exact bending is unknown, we randomly bend the modelled fiber 10 times, in a random direction and with a random radius of curvature between 1 and 4 centimeter. Table I shows the bend parameters in the final model. All bends are 10 centimeters long. For each bend, the propagation from beginning to end is calculated by solving equation (21). Figure 1 shows the difference for the linear output intensity for an ultrashort input pulse for the fiber model with and without mode coupling due to bending. Without bends, the output intensity shows intensity peaks at the input pulse location, because the modes are not coupled while propagating through the fiber. In contrast, the output intensity with bends is much more evenly distributed, without clearly visible patterns. A similar improvement is visible in the nonlinear output intensity.  (b) Linear output intensity for the same input pulse, but now for the fiber model with bends. The speckle pattern intensity is much more even due to mode coupling.
bins can be changed during time-domain wavefront shaping. The frequency amplitude array with square bins is then convoluted with a Gaussian with a FWHM of √ 0.61 bins in order to simulate the finite width of a single frequency in the Fourier plane of our experimental pulse shaper (61 µm compared to 100 µm wide pixels). The algorithm for wavefront shaping in the simulation is the same as the algorithm that is used in the experiment.