Spatial, spectral, temporal and polarisation resolved state tomography of light

In standard Stokes polarimetry illustrated in Fig. 1(a),

the unknown polarisation state of light, mathematically described by a density matrix $\hat{\rho}_{\text{in}}$, is interrogated by a sequence of polarisation analyser states that are defined by the orientation of the polarisation optics. In the depicted example, the power coupled into the SMF is equal to the expectation value of the unknown input state ($\hat{\rho}_{\text{in}}$) in the horizontal polarisation state. This horizontal polarisation state would exist at the input plane if the light was hypothetically launched from the SMF. Detecting the power with an SMF for multiple analyser states enables the recovery of the Stokes vector ($\mathbf{S}$), with its elements $S_{n}$ signifying weights of Pauli matrices $\hat{\sigma}_{n}$ in a sum that reconstructs the density matrix ($\hat{\rho}_{\text{measured}}$) of the unknown state. The density matrix encapsulates complete information about the unknown input state, including whether it represents a mixture of up to $N=2$ states or a pure state. The situation is mathematically equivalent for spatial states (Fig. 1(b)), however, to perform the projective measurements experimentally in this case, we have to substitute the polarisation optics for a spatial light modulator, which acts as a reconfigurable spatial analyser state device. The power coupled into the SMF in this scenario is equal to the expectation value ($|\text{field overlap}|^{2}$) of the unknown input ($\hat{\rho}_{\text{in}}$) in the state defined by LP11 mode with lobes oriented horizontally, as illustrated in Fig. 1(b). We can again imagine that we hypothetically launch the light from the SMF towards the SLM which transforms the fundamental LP01 mode of the SMF into LP11 projective state at the source plane. We can mathematically overlap the source and the hypothetical LP11 field at the source plane but since no detector exists there the overlap value cannot be determined experimentally at this plane. However, the reciprocity of light ensures conservation of field overlap value in the system, which in turn allows its experimental recovery at the SMF plane (more details in Supplementary Note 4).

Fig. 2(a) shows a simplified schematic of the experimental setup.

The unknown, pure or mixed state of light that generally depends on wavelength and time ($\hat{\rho}_{\text{in}}\left(\lambda,t\right)$) is separated by a polarisation beam splitter (PBS) into two paths to allow independent interrogation of vertical (V) and horizontal (H) polarisation by a single, large-area SLM that handles both polarisations. The orthogonal polarisation paths are recombined with another PBS and coupled into an SMF. The SMF subsequently routes the light to an optical spectrum analyser (OSA) and oscilloscope to analyse the input state with respect to wavelength and temporal degrees of freedom. More detailed experimental setup configuration can be found in Supplementary Note 1, along with the calibration and alignment procedures of the system in Supplementary Note 2.

The unknown state, represented by the density matrix ($\hat{\rho}_{\text{in}}\left(\lambda,t\right)$), can have a spatial profile with fine features due to the existence of multiple coherent and incoherent spatial fields at each wavelength and time. In order to extend the technique to more complicated spatial fields with more features, necessitates the extension of Stokes formalism from the limited two-dimensional case (Fig. 1(b)), where the analyser states are only sensitive to two modal components, into much higher dimension, until the Stokes states span all the modal components present in the beam. In our experimental case, a Stokes space dimension of $N=21$ satisfies this condition, and actually exceeds it. When no a priori information about the beam is known, the Stokes dimensionality $N$ should be chosen with some redundancy, at the price of carrying more projective measurements.

The natural generalisation of Pauli matrices, used for the reconstruction of the density matrix in two-dimensional scenario, into higher dimensions is facilitated by the generalised Gell-Mann matrices⁴¹ that span the space of observables of an $N$-dimensional Hilbert space. A vector in such an $N$-dimensional Hilbert space, specifically its vector elements, can represent complex superposition coefficients of an optical field in any orthogonal spatial basis of interest. The choice of this basis can be arbitrary. To demonstrate that the tomographic measurement can be performed using a random orthogonal spatial basis, we generate the spatial analyser states in the following fashion. We take as a starting point a Laguerre-Gaussian (LG) basis (Fig. 2(b)), with the waist of $w_{0}=3\,\rm{\mu m}$ reflecting the approximate extent of the investigated beam, and generate a single new spatial state composed of a random complex superposition of $N$ states in the LG basis (Fig. 2(c)). We subsequently apply the Gram-Schmidt orthogonalisation process to create the remaining $N-1$ orthogonal random spatial modes (Fig. 2(c)) to form the full random Laguerre-Gauss (RLG) basis with the dimension $N$.

The number of modes $N$ in the initial selected basis should reflect the expected level of spatial details in the probed system. For instance, the spatial analyser states can sense high spatial frequency components only if the initial Stokes basis and the related random basis contains such high spatial frequency components. The choice of LG basis is convenient in cases with expected rotational symmetry, but the method is not limited to LG basis and works with Hermite-Gauss, position, Bessel, Hadamard or Fourier basis as per application needs.

To generate the spatial analyser states that perform the tomographic measurements, we first find a set of unique eigenvectors of the Gell-Mann matrices. The vector elements of each eigenvector represent the complex superposition coefficients in the RLG basis and make up a given spatial analyser state (Fig. 2(d)). By sequentially applying the spatial analyser states using the SLM and measuring the intensity as a function of wavelength and time, using the OSA and oscilloscope respectively, we can reconstruct the generalised Stokes vector as both a function of wavelength and time. A blueprint for the calculation of $2N^{2}-N=861$ analyser masks for Stokes dimension of $N=21$ that generate the projective states, along with the mask alignment requirements is in Supplementary Note 2. Similar to the formalism presented in Fig. 1(b), by summing the Gell-Mann matrices with weights given by the Stokes vector, we reconstruct the spectrally and temporally resolved density matrices of the unknown beam. The eigenvectors and eigenvalues of the reconstructed density matrices (expressed in the RLG basis), are the spatial eigenstates and their probabilities in the investigated wavelength or temporal bin. Since both the spectrum and the time trace are collected in a single sweep for one specific spatial analyser state, only $2N^{2}-N$ measurements in total are required for the full reconstruction of the spatial, spectral and temporal information of the beam. A detailed theory of the high-dimensional Stokes analysis is presented in Supplementary Note 3.

We chose the vertical-cavity surface-emitting laser (VCSEL) diode as the unknown state $\hat{\rho}_{\text{in}}\left(\lambda,t\right)$ to demonstrate the method. VCSEL is a spatially, spectrally and temporally interesting source that represents a class of optical beams difficult to analyse using existing techniques. VCSEL cavity supports several spatial modes, with distributed spectral peaks (Fig. 3(a)).

The position of the spectral peaks and the number of spatial modes varies based on the VCSEL driving conditions like the bias current and the temperature (see Supplementary Note 7 for modally resolved light-current curves and details of driving conditions). The number of spatial modes and their relative power impacts the overall spatial profile of the VCSEL, which influences its performance in applications such as face recognition and 3D scene scanning for augmented reality ⁴². The ability to measure the spatio-spectral profile can provide useful feedback during VCSEL engineering and fabrication to improve the VCSEL performance in those applications.

In this section, we use the measured spectrally resolved density matrix to distinguish the amplitude and phase of each spatial mode for each spectral slice. Uniquely, the density matrix unlocks the ability to resolve multiple mutually incoherent modes within the same spectral bin, which typically occurs in situations when spectral peaks of modes are closely packed and thus beyond the resolution of the spectrometer. We confirm the validity of the results by comparing the probability weighted sum of reconstructed mode intensities with the intensity profiles obtained from a simple raster-scan, with results obtained very similar in both cases.

To perform the spatio-spectral analysis, we first collect the spectrum for each spatial analyser state in Fig. 2(d) and find the Stokes vector and the density matrix for each wavelength bin (resolution bandwidth of $70\,{\rm pm}$ given by the spectrometer) in the RLG basis. The measured density matrix at selected spectral peaks is depicted in Fig. 3(b). The eigenvalues and eigenvectors of the density matrix represent the probabilities and the spatial states at each wavelength, respectively. We note that the acquired density matrix should theoretically be positive-semidefinite to ensure positive probability for each supported state. However, due to measurement induced noise, this is not generally the case, and we have to apply Higham algorithm that finds the nearest positive semidefinite matrix in the Frobenius norm sense ⁴³.

The mutually incoherent spatial states with the three highest probabilities at each spectral peak of Fig. 3(a) are shown in Fig. 3(c-e). The spatial states were constructed from the eigenvectors of the density matrix, with eigenvectors representing the complex superposition coefficients in the RLG basis. Despite the randomness of the RLG spatial basis, the reconstructed modes have the well-known linearly polarised (LP) mode spatial profiles typical for a circular laser cavity. Some spectral peaks have a significant second dominant state (Fig. 3(d)) and another weaker third state (Fig. 3(e)). This is most pronounced for the spatial modes corresponding to LP02 (green column), LP21 (red column) and LP11 (purple column) but also faintly visible for LP31 (orange column). In the LP11 case, the second dominant state is present simply due to leakage of the powerful neighbouring LP11 mode with a perpendicular orientation of lobes (brown). In the LP21 case, the wavelength bin contains a mixture of mutually incoherent states that are spectrally very close, within the $70\,{\rm pm}$ resolution of our spectrometer. Remarkably, even though we cannot spectrally resolve the modes, we can utilise the density matrix to recover not only the individual spatial profiles in the mixture but also their relative power. The ability of the technique to resolve spectrally overlapping spatial modes is one of its hallmark features. On the basis of mutual incoherence, it enables spatial differentiation of features beyond the spectral resolution. We explore this compelling attribute of Stokes analysis in detail in Supplementary Note 8, where we intentionally decrease the resolution bandwidth of the OSA to $10\,{\rm nm}$, fully spectrally mixing the spatial modes in a single detected wavelength bin, yet the Stokes analysis yields the correct spatial states and their relative powers in the mixture. For the LP21 case, the recovered probabilities in the mixture are $63\%$ for the LP21a and $27\%$ for the LP21b state. Apart from the two dominant modes totalling $90\%$ of power in the investigated wavelength bin, there is also approximately $2\%$ of LP02 mode as visible on a partially distorted spatial profile in Fig. 3(e). The $8\%$ of unaccounted power distributes among the remaining 18 spatial eigenstates (not displayed) of the density matrix, all of them with a random, speckle-like spatial profile. There is no physical reason to expect that such random spatial states exist in the laser cavity. We instead associate the $8\%$ of the remaining power to noise that inevitably occurs during the hundreds of sequential projections, either due to power fluctuation of the VCSEL or the thermally induced spectral offset.

The noise not only affects the amount of light that cannot be attributed to physical spatial modes but also influences the reconstruction fidelity of the spatial profile of the relatively weak modes. This is noticeable for the weak LP02 (Fig. 3(e), red column) and also evident in the green column where both Figs. 3(d,e) are likely the weak, leaking LP21 modes. Nevertheless, despite these noise-induced effects, the mixed state reconstruction in Fig. 3(f), obtained as a probability weighted intensity sum of all 21 spatial eigenstates, is in almost perfect match with the intensity profiles (Fig. 3(g)) obtained by raster scan, which implies accurate determination of the spatial modes and their relative powers. We perform the raster scan by deflecting the beam using linear phase ramps (tilt) at the SLM, which scans the beam around the fixed core of the single-mode fibre. For each tilt, we collect the corresponding spectrum to generate the $(x,y,\lambda)$ datacube, with relevant wavelength slices plotted in Fig. 3(g). The influence of the noise level on the Stokes analysis is explored in-depth in Supplementary Note 9 for different levels of noise, RLG basis dimensionality, for pure and mixed states and also exploring the effect of orthogonality of input spatial states.

The number of spatial modes, their spectral peak position and relative power depends on the VCSEL temperature and bias. The applied bias does not have to be constant. VCSELs are often modulated at speeds of tens of GHz in short-reach interconnects that form the backbone of modern data centers and local area networks. The quality and speed of modulation can be influenced by mode competition within the cavity and the dynamic modal content of the beam impacts the data transmission through multimode fibres due to modal and chromatic dispersion. Having access to the temporal dynamics of VCSELs on an individual spatial mode basis is therefore of paramount importance in the quest to design faster optical interconnects.

The measured temporally resolved density matrix distinguishes individual mutually incoherent modes, all evolving simultaneously at high speed within the beam. In the spatio-spectral case, the method distinguished between mutually incoherent modes residing in the same spectral bin. Here, the method resolves between all mutually incoherent fields of the beam in each temporal bin, which allows tracking their power (probability) as a function of time. The ability of the method to differentiate mutually incoherent fields despite no way of separating them in time or spectrum is the added layer of flexibility that empowers the method. Whenever there is a lack of temporal or spectral resolution, the method can still spatially resolve the fields, which can be seen as a form of spontaneous spatial filtering.

We modulate the VCSEL via a pseudorandom binary sequence (PRBS) generator that is synchronised with an oscilloscope (more details of optical setup in the Supplementary Note 1). The PRBS enables repeatable temporal modulation of the VCSEL between the high ($7\,{\rm mA}$) and the low ($5\,{\rm mA}$) bias states. PRBS modulation in tandem with a synchronised oscilloscope facilitates the sequential application of several hundreds of spatial analyser states on the optical output with identical temporal modulation. We emphasize that the SLM refresh rate does not limit our temporal resolution as we collect the full temporal trace at $22\,{\rm GHz}$ sampling rate – given by the bandwidth of the photodiode used – for each spatial analyser state displayed on the SLM.

The modulation of VCSEL disrupts the equilibrium in the cavity, which triggers dynamic changes to modes, with their relative power evolving and redistributing. Our technique allows observation of these changes for all spatial mode fields, at all delays. Fig. 4 shows two temporal snapshots of the VCSEL at times $t=0.32\,{\rm ns}$ (yellow box, high bias) and $t=1.08\,{\rm ns}$ (purple box, low bias). The selected times correspond to the opposite extremes of bias applied to the VCSEL with values of $7\,{\rm mA}$ and $5\,{\rm mA}$ respectively. The overall VCSEL response (dark-grey curve, obtained from the SLM raster scan and integrated spatially for each temporal bin) closely follows the PRBS drive signal (light-grey curve) modulating the VCSEL at $10\,{\rm GHz}$. The remaining coloured curves in Fig. 4(a) are the probabilities of the spatial states that are given by the eigenvectors of the density matrices (Fig. 4(h,q)) in a given temporal bin.

The reconstructed spatial profiles of the first six dominant states are in Figs. 4(b-g, k-p). Interestingly, during each transition from high to low bias state and vice-versa, the probability distribution of the spatial eigenstates change. The observed change is most evident on the LP21 (Fig. 4(c)) and LP01 (Fig. 4(d)) eigenstates. For high bias state, LP21 (Fig. 4(c)) dominates over LP01 (Fig. 4(d)) by about $5\%$. The situation is completely reversed for the low bias case, with LP01 (Fig. 4(l)) having $10\%$ higher probability compared to LP21 (Fig. 4(m)). This is an expected behaviour for laser diodes, with low bias favouring the low order modes, in this case the fundamental mode LP01, at the expense of higher order modes. Due to this effect, the overall spatial intensity profile of VCSEL is not constant during modulation as can be seen by comparing Figs. 4(i) and Figs. 4(r) obtained by SLM raster scan. The high-bias intensity profile in Fig. 4(i) is brighter and has more pronounced side lobes due to the higher probability of LP11 and LP21 modes compared to the low bias case in Fig. 4(r) which is not only dimmer but also more influenced by the LP01 mode. The observed eigenstate probability changes as a function of bias and follows similar trends to measurements of modally resolved light-current curves presented in Supplementary Note 7. We also note that the SLM raster scan profiles (Figs. 4(i,r)) are in agreement with the probability-weighted sum of intensity profiles of all reconstructed spatial states (Figs. 4(j,s)), obtained using the high-dimensional Stokes method. This manifests the viability of the method for temporal analysis of the beam.

Compared to the spatio-spectral analysis case, the reconstructed profiles of spatial modes appear to have distorted wavefronts. To determine the cause of this effect, we measured the variation of the VCSEL power for a fixed SLM mask (around $5\%$) and performed the numerical simulations of the Stokes analysis for a mixture of modes with a relative power given by the relative strength of the spectral peaks in Fig. 3(a). A total of $6$ modes detected in the spectral analysis (the LP31a and LP31b are not counted as they are weak) is received by the temporal detector at all times, which means that the Stokes analysis has to distinguish a mixture of $6$ modes for each temporal bin (compared to 2 or 3 in spectral analysis). The numerical simulation yields similar wavefront distortion of the reconstructed modes, with increasing level of noise negatively affecting the field overlap value between the input and the reconstructed modes (see Supplementary Note 9 for detailed analysis of the noise effect on Stokes analysis). An additional source of noise can be attributed to the spatial modes of the real-world cavity being not completely orthogonal due to current injection patterns and manufacturing defects. This is evidenced by LP11 odd and even modes not having perfectly orthogonal orientation of lobes in SLM raster scans of Fig. 3(g). In spectral analysis case, the effect of this imperfect orthogonality on wavefront reconstruction is not so evident because of limited number of modes in each spectral bin compared to the temporal case. We numerically simulate the effect of imperfect orthogonality on reconstructed modes in Supplementary Note 9.

The reconstructed modes and their temporally resolved probabilities/relative powers provide unique insight into the temporal dynamics within the cavity on a mode-by-mode basis. This result can be useful to explore VCSELs of varying aperture scale and shape to develop higher yield of optimum mode dynamics during fabrication. The typically used raster-scanning approaches only obtain the intensity variation of the VCSEL during modulation similar to Figs. 4(i,r). While such information is useful, it does not elucidate the role of individual modes. Our tomographic approach provides not only the overall intensity profile (Fig. 4(j,s)), but crucially, the amplitude, phase and the relative powers of all spatial modes at each time.