Superluminal Motion-Assisted Four-Dimensional Light-in-Flight Imaging

Advances in high-speed imaging techniques have opened new possibilities for capturing ultrafast phenomena such as light propagation in air or through media. Capturing light in flight in three-dimensional xyt space has been reported based on various types of imaging systems, whereas reconstruction of light-in-flight information in the fourth dimension z has been a challenge. We demonstrate the four-dimensional light-in-flight imaging based on the observation of a superluminal motion captured by a new time-gated megapixel single-photon avalanche diode camera. A high-resolution light-in-flight video is generated without laser scanning, camera translation, interpolation, or dark noise subtraction. An unsupervised machine-learning technique is applied to analyze the measured spatiotemporal data set. A theoretical formula is introduced to perform least-square regression for numerically solving a nonlinear inverse problem, and extra-dimensional information is recovered without prior knowledge. The algorithm relies on the mathematical formulation equivalent to the superluminal motion in astrophysics, which is scaled by a factor of a quadrillionth. The reconstructed light-in-flight trajectory shows good agreement with the actual geometry of the light path. Applicability of the reconstruction approach to more complex scenes with multiple overlapped light trajectories is verified based on a data set generated by Monte Carlo simulations. Our approach could potentially provide novel functionalities to high-speed imaging applications such as non-line-of-sight imaging and time-resolved optical tomography.


I. INTRODUCTION
Progress in high-speed imaging techniques has enabled observation and recording of light propagation dynamics in free space as well as in transparent, translucent, and scattering media.Various approaches for capturing the light in flight have been demonstrated: holographic techniques [1][2][3], photonic mixers [4], time-encoded amplified imaging [5], streak cameras [6,7], intensified charge-coupled devices [8], and silicon complementary metal-oxide-semiconductor image sensors [9].In recent years, one-and two-dimensional arrays of single-photon avalanche diodes (SPAD) have been adopted for the light-in-flight imaging systems [10][11][12][13][14][15].These sensors boost data acquisition speed by employing pixel-parallel detection of time stamping with picosecond time resolution and single-photon sensitivity.However, the sensors are capable of sampling only three-dimensional spatiotemporal information (x, y, t) of the target light, and hence the tracking of light-in-flight breaking outside an xy plane is a challenge.
A method to reconstruct extra-dimensional position information z from measured spatiotemporal information has recently been investigated by researchers [16,17].The authors remark that different propagation angles with respect to the xy plane result in the different apparent velocity of light.Thus, comparing the measured spatiotemporal data set with this theory could give an estimation of the z component in the light propagation vector.Yet, the analysis is limited to a simplified case with either single straight light paths within the xz plane passing a fixed point or folded light paths with the angles between the light propagation vector and focal plane no larger than 10 degrees; however, a complete four-dimensional reconstruction of light-in-flight in arbitrary paths remains to be verified.In addition, spatial resolution of the detector array is critical for accurate reconstruction of the four-dimensional trace; the reported SPAD sensor resolutions-e.g., 64 × 32 [15], 32 × 32 [10], and 256 × 1 [13] pixels-are not sufficient, whereas significant improvement of the estimation error for (x, y, z, t) is expected with orders-of-magnitude-larger arrays, typically in the megapixels range.Complementary metal-oxidesemiconductor sensors routinely achieve megapixel resolution; e.g., 0.3-megapixel resolution is reported for highspeed imaging [9]; however, the necessary timing resolution is not reached, and it is limited to 10 nanoseconds.In this paper, we demonstrate the four-dimensional light-in-flight imaging based on the time-gated megapixel SPAD camera [18].In contrast to conventional time-correlated singlephoton counting (TCSPC) approaches, where power-and area-consuming time-to-digital converter (TDC) circuits restrict the scaling of the pixel array [19][20][21], our timegating approach [18] achieves a much more compact pixel circuit, suitable for large-scale arrays.Our camera achieves a megapixel format, while at the same time ensuring a time resolution comparable to that of a TDC.Owing to an unprecedented high-spatiotemporal resolution, our camera successfully captured an exotic behavior of laser pulse propagation.An object moving at the speed of light can appear to be moving faster than light in a camera image if one neglects to account for the travel time of light from the object to the camera.This "apparent" phenomenon is referred to as superluminal motion.We then introduce an equation to fit the measured data set in three-dimensional space for solving a nonlinear inverse problem.A four-dimensional point cloud is computationally reproduced without prior knowledge, exhibiting good agreement with the actual light propagation vectors.The reconstruction method is further extended to more complex scenes, where trajectories of multiple beams appear to cross each other in the field of view.

II. PRINCIPLE AND EXPERIMENTAL SETUP
To illustrate our method, we introduce a Minkowski space that describes the spatiotemporal propagation of a laser pulse and scattered photons towards fixed observation points.Figure 1(a) shows a schematic illustration of the light-in-flight observation in Minkowski space at y ¼ 0. A straight light propagation along the x axis is depicted as ðx; y; zÞ ¼ ðct; 0; 0Þ, where c is the speed of light.In Minkowski space, this equation corresponds to a single tilted line crossing the origin [purple arrow in Fig. 1(a)].A camera is located at position C i (i ¼ 1,2,5).The observation of the light in flight at the camera is mediated by scattered light at each point on the path (colored dots on the purple arrow).The propagation of scattered light from the points is described as light cones.Here, we define t as the time at which the light propagating from a laser source reaches a point (x, y, z), called the scattering point on the light path, and t 0 as a time at which the light scattered by the medium at (x, y, z) reaches the camera.Because of the finite speed of light, a time difference from t to t 0 is proportional to the spatial distance between the scattering point and the camera location.In Minkowski space, the observation time t 0 can be visualized as an intersection point of the light cone and a vertical line L passing C i (dotted line).The relation between the light position x and the corresponding observation time t 0 is shown in Fig. 1(b).The behavior is highly dependent on the camera position with respect to the scattering point; i.e., at certain points it appears to be faster than in others.For example, at C 1 where the light is coming towards the camera, the scattering light from all the points on the light path reaches the camera at the same time.The corresponding apparent velocity, defined as v ap ¼ dx=dt 0 , is infinite, or t 0 appears not to change irrespective of the distance of the scattering point (x, y, z) from the camera.At C 5 , in contrast, the light propagating away from the camera is observed with the apparent velocity of c=2.For other scattering points corresponding to relative camera positions C 2 , C 3 , and C 4 , one observes intermediate apparent velocities, depending on the relative position and angle of the line of sight between the scattering point and the camera.Hence, analyzing the apparent velocity provides an estimated light propagation vector with extra-dimensional information.The discussion can be generalized for higher dimensions; the detailed analysis of the measured (x, y, t) data set and its local apparent velocity would enable a four-dimensional reconstruction of the light-in-flight history.A 2D-projected intuitive visualization for the propagation of scattered light is provided in Video 1 in Supplemental Material [22].Note that angle-dependent apparent velocity of light can be seen in nature, for instance, in an astronomical phenomenon called superluminal motion.In this case, relativistic jets of matter are emitted from radio galaxies, quasars, and other celestial objects, appearing to travel faster than light to the observer on a timescale of months to years.Light echoes are also known to produce superluminal motion [23,24].In addition, observation of superluminal motion in the scattering of a wavefront from a diffusive surface and diffusive media was recently reported [7,25].In Ref. [7], compensation of the finite time offset from a scattering event to the detection point in a camera was performed based on separately acquired 3D scanning data of the scattering surface, whereas our approach does not require any additional measurement and relies only on the assumption that light propagates at constant speed.
Figure 1(c) is the experimental setup to verify the fourdimensional light-in-flight imaging.A 510-nm laser (average power, 2 mW; frequency, 40 MHz; optical pulse width, 130 ps; PicoQuant GmbH, Berlin, Germany) and the megapixel SPAD camera are synchronized by a pulse generator.The emitted laser pulses are reflected by four mirrors to construct a three-dimensional trajectory.The megapixel SPAD camera can be operated in the intensityimaging mode to capture the background scene and in the time-gating mode to capture laser propagation [18].

III. THREE-DIMENSIONAL LIGHT-IN-FLIGHT IMAGING
Figure 2 describes a sequence of images showing threedimensional light in flight in an actual setup.The principle of time-gating operation in the megapixel SPAD camera is depicted in Ref. [18].The data acquisition time for full gate scanning (250 slices) is 15 to 30 seconds, much shorter than in previous works (e.g., hours [6,7] to 10 minutes [10]).This time can be reduced further by introducing a higher power laser.
Figure 2(a) shows the procedure to generate a threedimensional light-in-flight video.The 250 captured slices of time-gating images form a photon intensity distribution for each pixel as a function of gate position.A rising edge position τ 0 and high level of photon intensity I in the intensity profile can be extracted for each pixel.Time-ofarrival information t 0 over the array is obtained by subtracting independently measured pixel-positiondependent timing skews from τ 0 .The intensity profile is given as a convolution of the arriving photon distribution and the gating profile, and the laser pulse width appears to be elongated by 3.8 ns in the raw intensity profile.The laser pulse width can be narrowed to the picosecond scale by replacing the intensity profile with a Gaussian distribution having a mean value of t 0 , a standard deviation of 150 ps (corresponding to the combined jitter of the laser and the detector), and integrated photon counts of I. Figure 2(b) shows the selected frames of a reconstructed threedimensional light-in-flight video, where each frame is superimposed with a two-dimensional intensity image independently captured by the same SPAD camera.We start to observe a laser pulse around t 0 ¼ 0 ns, which is reflected by multiple mirrors and goes out of sight around t 0 ¼ 5.832 ns.The laser beam at t 0 ¼ 5.112 ns appears to be stretched compared to the beam at t 0 ¼ 2.232 ns, which implies an enhanced apparent velocity of the beam coming towards the camera at t 0 ¼ 5.112 ns.Note that our system requires no mechanical laser scanning, spatial interpolation, or dark noise image subtraction owing to its high spatial resolution array and low dark count rate [18].
A full three-dimensional light-in-flight video with a frame interval of 36 ps is available in Video 2 in Supplemental Material [22].

IV. FOUR-DIMENSIONAL RECONSTRUCTION OF LIGHT IN FLIGHT
Figure 3(a) shows a color-coded plot of the observation time t 0 for the propagating light.To reconstruct the four-dimensional light path, the data set needs to be subdivided into straight paths.We adopt the 2D Gaussian mixture model (GMM) fitting for data clustering, which is a common unsupervised machine-learning technique.Figure 3(b) shows the data points separated into five clusters.Note that the data clustering is performed based on two-dimensional spatial information (x and y), whereas the time information t 0 is not taken into account.
Figure 3(c) depicts the coordinate system for the light-inflight analysis, based on the assumption that all photons are scattered exactly once.Here, we derive a theoretical formula to perform the least-square regression for estimation of the four-dimensional light trajectory: ð⃗ r 0 ; θ; ϕÞ ¼ arg min ⃗ r 0 ;θ;ϕ where N is the total number of measurement data points, (x i p , y i p ) the corresponding pixel position on the focal plane for the ith data point, −z p the focal length (12.5 mm), t 0i the measured observation time for the ith data point, and ⃗ r 0 ¼ ðx 0 ; y 0 ; z 0 Þ the position of the beam at t ¼ 0; the normalized light propagation vector is represented as ⃗ n ¼ ðsin θ cos ϕ; sin θ sin ϕ; cos θÞ using a polar coordinate system.A single straight light trajectory can be expressed as ⃗ r ¼ ⃗r 0 þ ct⃗ n, defined by five fitting parameters, x 0 , y 0 , z 0 , θ, and ϕ.These parameters are estimated by solving the above optimization problem for each data cluster.
Figures 3(d) and 3(e) show the measured observation time t 0 as a function of x p and y p , where (x p , y p ) is the position of the pixel on the focal plane.Different colors of the data points correspond to the different data clusters.The result of fitting for each data cluster is shown as a dashed line, exhibiting good agreement with the measured points for both x p and y p .Note that the fitting involves five independent parameters and is sensitive to the variation of the data.Hence, the spatiotemporal resolution of the camera is a critical factor determining the accuracy of the four-dimensional reconstruction.
As shown in Fig. 4(a), the four-dimensional point cloud is reconstructed without prior knowledge.The time evolution t of the beam is shown as colors of points, with a certain offset from the observation time t 0 due to the finite traveling time from the scattering point to the camera.The fully recovered four-dimensional information for each point enables the visualization of the point cloud from an arbitrary viewpoint, apart from the actual location of the SPAD camera.In contrast to the three-dimensional light-inflight video where the apparent velocity of light changes as a function of propagation angle, the four-dimensional point cloud indicates the uniform speed of light irrespective of the beam position and angle.Relatively larger variation of the data points is observed for clusters 4 and 5.This variation originates from the enhanced apparent velocity of light coming towards the camera, leading to the reduced fitting accuracy and increased deviation.
The accuracy of the reconstruction is evaluated in Figs. 4  A detailed analysis of the reconstruction performance for various parameter combinations is investigated based on Monte Carlo simulation in Supplemental Material [22].The simulation shows that with the horizontal resolution of 1024 pixels, the estimation errors of the fitting parameters improve by a factor of 6 to 9 with respect to previously reported 32 × 32 pixel SPAD cameras, thereby justifying the necessity of megapixel resolution for high-precision four-dimensional light-in-flight imaging.The results also suggest that a further increase of the pixel array size could have a non-negligible impact on the stability and accuracy of light-in-flight measurements, along with the potential for further expansion of the measurement in distance and field of view.
To quantify the superluminal motion in this experiment, we introduce apparent transverse velocity v T ap , which is often used to characterize the deviation of the observed velocity from the actual velocity in astrophysics.The apparent transverse velocity is defined as the apparent velocity projected onto the focal plane and can be expressed as where ζ is the angle between a light propagation vector and a line-of-sight vector.Figure 5 shows the reconstructed apparent transverse velocity distribution in xyz space.As expected, v T ap higher than c is observed for light paths coming towards the SPAD camera.A gradual change of v T ap in a single straight path is caused by the fact that the laser position is dynamically shifted during the propagation, thereby changing the angle ζ.The range of v T ap for each data cluster is summarized in Table I.
Note that precise observation of superluminal motion in astrophysics often takes months, whereas the observed laboratory-scale superluminal motion takes place in a nanosecond, corresponding to a 10 15 -times reduction with respect to the astronomical observation.This case indicates that the high-speed SPAD camera enables us to reproduce the astrophysical phenomenon in a spatiotemporal dimension scaled by a factor of a quadrillionth.

V. GENERALIZATION
In the previous section, we experimentally demonstrated the four-dimensional reconstruction for the trajectory of a single laser beam based on 2D GMM fitting.The measured data set consists of multiple straight paths, which are not spatially overlapped with each other.However, the assumption of constant speed of light over all trajectories allows the reconstruction approach to be extended to more complex situations, such as multiple beams propagating independently, possibly appearing to cross each other in the field of view.
Experimentally, both time-gating-based and TDC-based approaches are capable of recording multiple time stamps for each pixel with a single measurement, which ensures the acquisition of a full spatiotemporal data set for the apparently crossing paths.A promising approach to perform data clustering for such a complex data set is to utilize temporal information in addition to two-dimensional xyspace information.A straightforward extension of the proposed method is to perform 3D GMM fitting.When the data set can be expressed as U ¼ fu 1 ; u 2 ; …; u i ; …g with u i ¼ ðx i p ; y i p ; t 0i Þ, the 3D GMM can be formulated as where K is the number of Gaussian components, π k the kth mixing coefficient, N the multivariate Gaussian distribution, μ k the three-dimensional mean vector for the kth Gaussian component, and Figures 6(c) and 6(d) show the Monte Carlo-generated virtual data sets.Gaussian-distributed noise is added to the measured time stamps.In this setup, some of the pixels around the center of the field of view detect two distinct time stamps originating from the scattering events in two different paths.Figure 6(c) shows the distribution of the time stamps for primary events.The black region indicates that there is no measured time stamp.Time stamps for the secondary events detected later than the primary events are plotted in Fig. 6(d).The full data set for data clustering is formed by merging all the time stamps for the primary and secondary events.Although the analysis in this paper is conducted for propagation of pointlike light pulses, the proposed reconstruction approach can be partially extended to cases where the propagating wavefront forms a continuous line or curve at a fixed time, provided that time evolution of the wavefront can be analytically formulated.For instance, the trajectory of an arclike wavefront generated from a laser pulse through a one-dimensional diffuser could be reconstructed based on the extended model.In contrast, the reconstruction approach cannot be directly extended to cases where propagating light forms a continuous surface, i.e., a laser pulse scattered by a two-dimensional diffuser.
The reconstruction approach implicitly assumes that temporal distribution of incident photons for each pixel can be approximated by a superposition of finite numbers of Dirac delta functions; this ensures that the measured data set can be expressed as the discrete point cloud U ¼ fu 1 ; u 2 ; …; u i ; …g, with u i ¼ ðx i p ; y i p ; t 0i Þ.However, when capturing the propagation of a wavefront with a continuous surface (or some exceptional cases for a wavefront with a continuous curve), the temporal distribution of incident photons can be continuous, and the resulting data set is no longer discrete.This result gives additional complexity to the model because the continuous distribution function of the incident photons is dependent on local light intensity at a scattering point and its scattering angle towards the camera.Introducing some assumptions on light intensity distribution at the wavefront surface and the scattering model, e.g., Rayleigh scattering or Mie scattering, could be useful for building more sophisticated models to potentially address this issue, which lies beyond the scope of this paper.

VI. DISCUSSION
Our approach of reconstructing extra-dimensional lightin-flight information with a high spatiotemporal resolution single-photon camera can be applied to a wide range of high-speed imaging techniques.One of the promising applications is non-line-of-sight imaging [26][27][28][29][30][31]; instead of the direct detection of photons reflected by multiple surfaces and arriving back at the camera, introducing scattering media in the NLOS setup enables us to trace the propagation of emitted and reflected laser beams or wavefronts.Our reconstruction approach enables the estimation of the light propagation vector, in addition to the scattered time and position, even for cases where the photons reflected by the target object do not directly return to the camera and only the secondary photons scattered in air or media are detected by the camera.This technique brings extended functionalities for imaging around and behind obstructing objects, e.g., in a foggy environment or under water.Another potential application is the combination with optical tomography techniques [32][33][34]; generalization of our theory can be useful for noninvasive three-dimensional monitoring inside a target object structure as well as nondestructive measurement of the threedimensional distribution of physical parameters such as the refractive index and transmittance.

5 C 1 C 2 C 3 C 4 FIG. 1 .
FIG. 1. Principle and experimental setup for four-dimensional light-in-flight imaging.(a) Schematic illustration of light propagation and scattering in Minkowski space at y ¼ 0. Propagating light is depicted as a purple dashed arrow, and scattered light is shown as a light cone.(b) Schematic plot of t 0 dependence of x.The slope of the dashed line indicates the apparent velocity.(c) Experimental setup for light-in-flight imaging.The picosecond laser and SPAD camera are both controlled by a pulse generator.When performing the light-in-flight measurements, all the mirrors are confined by a large transparent acrylic box (not shown), and a fog generator forms a small amount of mist in the box to enhance the scattering of the laser in air.The origin of the Cartesian coordinate is set at the optical center of the lens of the SPAD camera.

FIG. 2 .
FIG. 2. Video generation of three-dimensional light-in-flight imaging.(a) Data processing flow for generation of light-in-flight video.(b) Laser pulse propagating over time, shown along the dashed line in the top-left image.Superluminal motion of the laser pulse is observed in the latter frames.The light-in-flight video is taken under dark conditions, and each frame is superimposed with a 12-bit two-dimensional intensity image independently taken with the same SPAD camera under room light.The image resolution is 1024 × 500 pixels.

FIG. 3 .
FIG. 3. Analysis of extracted spatiotemporal data set.(a) Color-coded plot of observation time for the light path.Red segments are observed earlier, and blue segments are observed later.Dark regions have no timing data.(b) Light path data subdivided into five data clusters of straight paths using 2D GMM fitting.(c) Schematic view of the coordinate system for light-in-flight reconstruction.The laser beam image is projected onto the focal plane by an objective lens in a flipped geometry.(d) Observation time t 0 as a function of x p , the x position on the focal plane.Each data cluster is individually fitted using a theoretical formula.Dashed lines are the fitting curves.(e) t 0 as a function of y p .
(b) and 4(c).The figures show the fitted light propagation angles θ and ϕ as a function of the actual

FIG. 4 .
FIG. 4. Reconstructed four-dimensional light-in-flight observation.(a) Four-dimensional point cloud reproduced from the measured three-dimensional spatiotemporal data.The origin of the Cartesian coordinate system is defined at the optical center of the lens for the SPAD camera.The (x, y, z, t) coordinates of each data point are reproduced using the obtained fitting parameters for the corresponding data cluster.(b,c) Comparison of fitted light propagation angles θ and ϕ, respectively, with the actual angles for five data clusters.The fitted results show good agreement with the actual angles.

1 FIG. 5 .
FIG. 5. Reconstructed apparent transverse velocity distribution in xyz space.The color of each dot corresponds to the locally estimated apparent transverse velocity, normalized by the speed of light in vacuum.

Figures 7 (
a) and 7(b) show the result of 3D GMM-based data clustering, where the full data set is successfully classified into cluster 1 (red) and cluster 2 (blue).For each data cluster, the data are fitted by Eq. (1) to reproduce the four-dimensional point cloud.Figures7(c) and 7(d) show the reconstructed four-dimensional point cloud projected onto the xz plane and yz plane, respectively.The reconstructed light trajectories are consistent with the virtual configuration shown in Figs.6(a) and 6(b).The result confirms that the reconstruction method can be naturally extended to more complex situations where multiple light trajectories are apparently crossing in the field of view.

FIG. 6 .FIG. 7 .
FIG. 6. Virtual setup and data set generated by Monte Carlo simulation of four-dimensional light-in-flight imaging for two separate light trajectories.(a) Schematic illustration for top view of the virtual setup.(b) Schematic illustration for side view of the virtual setup.(c) Time stamps for primary events detected at each pixel, generated by Monte Carlo simulation.The black region indicates "no data."(d) Generated time stamps for secondary events.

TABLE I .
Summary of apparent transverse velocity for each data cluster.