Parallel random LiDAR with spatial multiplexing of a many-mode laser

We propose and experimentally demonstrate parallel LiDAR using random intensity fluctuations from a highly multimode laser. We optimize a degenerate cavity to have many spatial modes lasing simultaneously with different frequencies. Their spatio-temporal beating creates ultrafast random intensity fluctuations, which are spatially demultiplexed to generate hundreds of uncorrelated time traces for parallel ranging. The bandwidth of each channel exceeds 10 GHz, leading to a ranging resolution better than 1 cm. Our parallel random LiDAR is robust to cross-channel interference, and will facilitate high-speed 3D sensing and imaging.


I. INTRODUCTION
LiDAR (Light Detection and Ranging) is a key technology in autonomous driving, robotics, augmented and virtual reality. High-speed acquisition of distance information over a wide field of view is crucial for such applications [1]. The most common LiDAR scheme is based on time-of-flight measurement: the echo of an optical pulse tells the distance of a target from the time lag [2]. Raster-scanning of a probe beam is commonly employed for three-dimensional (3D) ranging [3,4]. However, the mechanical scanning rate limits the acquisition speed. Parallel LiDAR based on multiple-input multipleoutput (MIMO) scheme [5] facilitates high-speed 3D mapping. To generate multiple probe beams in parallel, spectral and/or temporal multiplexing of broadband light sources [6][7][8][9][10][11] are employed. With total bandwidth fixed, a larger number of spectral channels is accompanied by a narrower bandwidth of each channel, which means a worse resolution. Similarly, temporal multiplexing relies on interleaved transmission/reception, which limits the detection speed. Spatial multiplexing can avoid such problems while providing many more channels by using a 2D array of lasers [12,13]. However, all spatial channels have similar temporal/spectral waveforms, which will cause channel interference and ranging ambiguity.
Here we propose and demonstrate parallel random Li-DAR by using a single multimode laser for the simultaneous creation of many uncorrelated random time traces. Instead of using chaotic lasing dynamics, we resort to a different physical process-the spatiotemporal beating of many lasing modes in a single cavity-to generate ultrafast intensity fluctuations that vary spatially. In particular, a near-degenerate cavity can support lasing in many transverse and longitudinal modes. A slight detuning of the degenerate cavity lifts the frequency degeneracy of transverse modes, leading to a dense RF spectrum. A large number of transverse lasing modes offers several hundreds of probe beams with uncorrelated intensity fluctuations for parallel ranging. The resolution of an individual channel is determined by the lasing emission bandwidth, which is about 100 GHz for the solidstate laser used here. With sufficiently fast photodetection, the ideal resolution can be as fine as 1 mm, higher than that of chaotic semiconductor lasers. We further show that our parallel LiDAR scheme is robust to crossinterference between channels.

II. MULTIMODE LASING
Instead of building a 2D array of separate lasers, we combine many lasers in a single cavity. This is realized with a degenerate cavity [37], as shown schematically in Fig. 1(a). It consists of two lenses in a 2f − 2f telescope arrangement between two flat mirrors [38]. The focal length f of both lenses is 100 mm, which yields a cavity length L = 4f = 400 mm. A thin disk of Nd:YVO 4 crystal, located at one end of the cavity, provides optical gain when optically pumped by a high-power laser diode at 808 nm (Coherent, FAP800-40W). The other end of the cavity is an output coupler with 99.6% reflectivity at 1064 nm.
In principle, the degenerate cavity has the self-imaging property, where any point on one end mirror is imaged back onto itself after light propagates one round trip inside the cavity [ Fig. 1(a)]. This property, combined with optical pumping, allows simultaneous lasing  in many transverse modes with degenerate frequency and quality factor [39,40]. Each transverse mode creates a diffraction-limited area on the gain disk, and the total number of lasing modes is given by the ratio of the pump area over the diffraction area. The latter is determined by the numerical aperture of the cavity as well as optical misalignment. Small misalignments, aberrations, and thermal lensing induce a small degree of inherent detuning [41], therefore in practice the frequency degeneracy is slightly lifted with limited spectral spread.
Continuous-wave lasing is achieved in the degenerate cavity with optical pumping at room temperature. The pump beam has a diameter of 2 mm on the gain disk. The optical pump power is 22 W, well above the lasing threshold of 4.1 W. The output emission power from the degenerate cavity is 1.1 W. Since the optical gain bandwidth (> 100 GHz) is much wider than the frequency spacing of longitudinal modes (free spectral range ∆ν = 0.375 GHz), lasing occurs for hundreds of longitudinal modal groups, each containing many near-degenerate transverse modes [42].
The beating of transverse and longitudinal lasing modes creates intensity variations in space and time. If the temporal fluctuation is random at each location and uncorrelated with other locations, these intensity traces can be used for parallel LiDAR via spatial multiplexing. We characterize the laser emission intensity fluctuation using an amplified InGaAs photodetector (Newport, 818-BB-30A, 1.5 GHz bandwidth) and a real-time high-speed oscilloscope (Keysight UXR0204A, 20 GHz bandwidth). We note that 1.5 GHz (20 GHz) is the minimum 3-dB bandwidth tabulated by the manufacturer, beyond which the signal is attenuated electrically (digitally). Fourier transform F of the time trace of emission intensity I t (t) gives the radio-frequency (RF) spectrum. Figure 1(b) shows the experimental RF spectrum of the nearly degenerate cavity laser, obtained by taking the Fourier magnitude squared |F[I t (t)]| 2 . It features periodic modulation, with dominant RF components clustered around the integer multiples of ∆ν. We also compute the temporal correlation of emission intensity, C t (∆t) = I t (t) I t (t+∆t) t , where ... t denotes averaging over time t. In Fig. 1(c), C t (∆t) exhibits a series of correlation peaks spaced by ∆t = 2.7 ns, which corresponds to the round-trip time in the cavity. Such long-range correlation reflects periodic oscillation of emission intensity, which will significantly degrade the performance of our random LiDAR.
To suppress the long-range correlation in time, we need to generate a flatter RF spectrum, which is also a key feature in chaotic LiDAR [28]. This is achieved by further breaking the frequency degeneracy between different transverse lasing modes with cavity detuning [43]. As shown in Fig. 1(a), we slightly offset the axial location of the output coupler by about 1 cm. The transverse modes, represented by diffraction-limited spots at the original plane of the output coupler, will diffract and partially overlap at the new position. The spatial overlap leads to mode coupling, which lifts the frequency degeneracy. Temporal interference between these transverse modes will create additional beat notes that fill the gaps between integer multiples of ∆ν in the RF spectrum [ Fig. 1(d)]. Consequently, long-range temporal correlations are significantly reduced, as shown in Fig. 1(e). Hence, detuning the degenerate cavity geometry randomizes the temporal variation of laser intensity, which can be used for LiDAR.

III. AXIAL RESOLUTION
The axial resolution of random LiDAR is determined by the time scale of intensity fluctuation, or the bandwidth of the RF spectrum. To investigate how broad the RF spectrum is for the detuned cavity, we replace the aforementioned amplified photodetector with one having an order-of-magnitude higher bandwidth (Newport, 818-BB-36, 22 GHz bandwidth) combined with a low-noise amplifier (AT Microwave, AT-LNA-0018-3825H, 38-dB gain). The far-field emission out of the laser cavity is collected and sampled by a real-time high-speed oscilloscope (Keysight UXR0204A, 20 GHz bandwidth) at a sampling rate of 128 GS/s. As shown in Fig 2(a), the measured emission intensity exhibits rapid random fluctuation in time. The RF spectrum in Fig. 2(b) features densely packed peaks caused by many-mode beating. Each peak has a narrow linewidth (∼5 kHz measured from 1 ms-long trace, frequency resolution of 1 kHz). The numerous peaks are distributed more or less uniformly over a broad range in the RF domain. These peaks have irregular spacings and varying amplitudes, constituting a broad, dense RF spectrum. The sharp drop at 20 GHz is due to the limited bandwidth of the oscilloscope.
The wide RF spectrum corresponds to a short time scale of emission intensity fluctuation. Figure 2(c) shows the temporal autocorrelation of intensity trace from the detuned degenerate cavity laser. The correlation width ∆t c , defined by full-width-at-half-maximum (FWHM) of the temporal correlation function C t (∆t), is 45 ps. It determines the axial resolution of the random Li-DAR, namely, two targets separated axially by a distance ∆D = c∆t c /2 = 7 mm, can be resolved with ∆t c = 45 ps.
The broad RF spectrum and fast temporal decorrelation of emission intensity from the detuned degenerate cavity laser will significantly improve the axial resolution of random LiDAR. The experimentally measured RF bandwidth and decorrelation time are still limited by the electrical bandwidth of our photodetection setup. The intrinsic RF spectrum width is dictated by the laser emission spectral width at optical frequency, which sets the largest beat notes (frequency difference) of lasing modes. We measure the optical spectrum of laser emission using a spectrometer with a wavelength resolution of 0.04 nm (Horiba TRIAX550 equipped with CCD3000, 1200 gr/mm). The wavelength width (FWHM) of the multimode laser emission is about 0.4 nm, corresponding to a frequency bandwidth of 100 GHz [46]. With sufficiently fast photodetection, the ultimate axial resolution of our LiDAR can be about 1 mm, which exceeds that of chaotic LiDAR.

IV. SPATIAL MULTIPLEXING
The number of spatial channels available for parallel random LiDAR is determined by the number of transverse lasing modes. To ensure a large number of trans- verse modes lasing, the detuning of the degenerate cavity is kept small, so that the reduction of quality factors for high-order transverse modes is less significant. The beating of many transverse and longitudinal lasing modes creates complex intensity fluctuations in space and time.
The number of distinct random waveforms generated at different spatial locations can be used as separate channels for parallel LiDAR.
To find the number of uncorrelated time traces generated by the detuned degenerate cavity laser, we characterize the spatially-resolved temporal fluctuations of emission intensity. The laser beam profile at the output coupler is enlarged by a factor of 6.7 with a set of imaging optics. Figure 3(a) shows the emission intensity profile measured by a CCD camera (Allied Vision, Mako G-131B). The time-integrated output pattern is a superposition of many transverse lasing modes. The effective beam diameter, defined by FWHM of azimuthallyaveraged intensity, is 9.8 mm. The central region with strong emission within this effective beam diameter will be investigated for parallel LiDAR.
To measure intensity fluctuations at different spatial locations simultaneously, we divide the output laser beam using a 50/50 beamsplitter and place two amplified In-GaAs photodetectors (Newport, 818-BB-30A, 1.5 GHz bandwidth) at each arm. The two arms have an identical length so that the two photodetectors simultaneously measure the intensity fluctuations at different locations of the output beam. By varying the transverse positions of the two photodetectors, we obtain two time traces with different separations across the beam. Figure 3(b) shows temporal fluctuations of emission intensity at two locations with a transverse distance of 1 mm, denoted as α and β. The distinct time traces can be used as independent probes for parallel LiDAR.
To find how close two spatial channels can be, we characterize the spatial correlation of intensity fluctuations. For two transverse positions (x, y) and (x + ∆x, y + ∆y), the cross-correlation of intensity fluctuations is given by C s (∆x, ∆y; ∆t) = I(x, y, t) I(x + ∆x, y + ∆y, t + ∆t) t . Experimentally, one photodetector is fixed in space, as the other one is scanned laterally across the output beam. For every offset (∆x, ∆y), the maximal magnitude of |C s (∆x, ∆y; ∆t)| among all time delay ∆t is plotted in Fig. 3(c). The maximal cross-correlation has the largest value when the lateral positions of two photodetectors coincide (∆x = ∆y = 0). With increasing offset in the lateral position, the time traces become different, and the maximal correlation drops rapidly. Eventually, at a large offset, the maximal cross-correlation levels off to a constant value. We attribute this residual correlation to the finite number of beat notes in the RF spectrum.
The width of maximal cross-correlation sets the lateral width of a spatial channel of parallel LiDAR. Figure 3(d) shows the azimuthally averaged maximal crosscorrelation as a function of radial distance ∆r = [(∆x) 2 + (∆y) 2 ] 1/2 . Each curve is fit by a Gaussian function with a constant background. The three curves represent the cross-correlations measured at different locations across the beam: at the center (x = 0 mm, y = 0), at the effective beam radius (x = 5 mm, y = 0), and near the beam edge (x = 7 mm, y = 0). The peak value of correlation at ∆r = 0 decreases from the beam center to the edge, as the emission intensity gets weaker and the signal-tonoise ratio (SNR) becomes lower (10 dB, 4.5 dB, 2.3 dB for x = 0, 5, 7 mm, respectively). To find the correlation width, we normalize the correlation functions in the inset of Fig. 3(d). The half-width-at-half-maximum (HWHM) increases from 0.15 mm at the beam center to 0.18 mm near the edge. The channel width ∆r c , defined as twice the HWHM, varies from 0.3 mm to 0.36 mm across the beam. To have spatial channels with uncorrelated intensity fluctuations, we set the channel spacing to 0.6 mm. Given the effective beam diameter of 9.8 mm [ Fig. 3(a)], the output beam of our single laser will provide approximately 300 independent spatial channels for parallel LiDAR.

V. ISOLATING SPATIAL CHANNELS
One issue for long-distance ranging is the mixing of spatial channels outside the laser cavity. Upon axial propagation, individual channels diffract and overlap with other channels in space. Their interference will make temporal intensity fluctuations vary with propagation distance. Consider a spatial channel at transverse location (x, y), its intensity trace at different longitudinal positions (along z-axis) is not related simply by a temporal offset, i.e., I(x, y, z, t) = I(x, y, z + d, t + d/c), where c is the speed of light. Instead, the temporal profile completely changes along z.
Typically for ranging applications, the reference beam of one channel propagates a much shorter distance to a detector than the probe beam reflected from an object. Their intensity traces would be uncorrelated, as confirmed experimentally below. A multimode laser beam is split to reference and probe beams by a beamsplitter [see Fig. 4(a)]. The reference beam is directed to a photodetector. An iris is placed in front of the lens to select a single spatial channel. The probe beam is reflected from a mirror and goes to a second detector. Another iris in front of the second lens selects the same spatial channel as the reference. The path length difference between the reference and probe beams is 1.35 m. Their intensity traces are completely different [ Fig. 4(c)], and exhibit no cross-correlation peak at any time lag [black dashed curve in Fig. 4(e)]. This result indicates the temporal profile of intensity fluctuation in a single spatial channel becomes uncorrelated with free-space propagation, and thus cannot be used for ranging applications.
To resolve this issue, the spatial channels must be isolated to stop their beating outside the laser cavity. There are several methods of isolation. One is coupling the multimode laser emission into a fiber bundle by a lens array. The fields in individual single-mode fibers no longer couple, and the time traces of intensity remain invariant with propagation (except for a time lag). Another way is using a two-dimensional array of pinholes to isolate spatial channels. As a proof-of-concept demonstration, we filter a single spatial channel with a pinhole. Below we show its temporal profile remains invariant with axial propagation, and will use it for ranging. Figure 4(b) is a schematic of our experimental LiDAR setup with an isolated spatial channel. The output beam from the detuned degenerate cavity laser is attenuated to 0.3 W. Then, to select a single spatial channel, in the beam we place a pinhole with 0.5 mm diameter, which is slightly smaller than the channel spacing of 0.6 mm. The power of this separated single channel is 0.36 mW. The filtered beam passes through a 50/50 beamsplitter to form the reference and the probe. The reference beam directly goes to a photodetector R. As an object for LiDAR demonstration, we use a flat aluminium mirror with a reflectivity of 85% at 1064 nm. The probe beam is reflected from the object and subsequently measured by a second photodetector S. To increase the collection efficiency, we use a plano-convex lens (focal length = 40 mm) to focus light onto the active area of the amplified photodetector (Newport, 818-BB-30A, 1.5 GHz bandwidth, 100 µm diameter). The probe beam power (after the 50/50 beamsplitter) is 0.18 mW, well below the safety regulation limit (class-1 lasers, 2 mW for emission at 1064 nm) for commercial LiDAR applications [44,45]. Figure 4(d) shows the time traces measured by two photodetectors S and R. The signal trace S exhibits nearly identical intensity fluctuation as the reference R, except for an offset in time. The offset is given by 2D/c, where D is the distance from the beamsplitter to the object, as the distance from the beamsplitter to S is equal to that to R. The object distance D can be retrieved from the cross-correlation between the signal and reference traces. Figure 4(e) shows a dominant peak of the cross-correlation (solid red curve) at the time lag corresponding to 2D/c. Experimentally we move the object so that D varies from 7 cm to 66 cm, corresponding to an interval from 0.46 ns to 4.4 ns time delay. For all distances, a cross-correlation peak stands out of the fluctuating background [ Fig. 4(f)]. Its time lag gives the object distance D, which agrees well with the actual distance (vertical lines).

VI. RANDOM LIDAR PERFORMANCE
We characterize the axial resolution of the ranging measurement. The resolution is a measure to distinguish two adjacent peaks, and thus determined by the peak width of the cross-correlation. From Fig. 4(f), the FWHM of the correlation peak ∆t c = 0.35 ns yields the axial resolution c ∆t c /2 = 5.3 cm. We note that the axial resolution is limited by the bandwidth of our photodetectors (1.5 GHz), and it can be improved with highbandwidth photodetectors.
We further measure the precision and detection capability of our LiDAR. The target is a high-reflectivity Aluminium mirror (reflectivity of 85% at 1064 nm) placed at a distance D = 30 cm. The sampling rate is 128 GS/s, and the time interval to compute cross-correlation is 200 ns. To characterize the LiDAR performance as a function of the SNR, we vary the signal strength by inserting neutral density filters into the beam path. To find the precision, we perform 100 consecutive ranging measurements and obtain statistics of the time lag for the maximal cross-correlation. Its standard deviation remains less than δt = 10 ps, which implies that the precision is c δt/2 = 1.5 mm, for SNR higher than 3 dB [ Fig. 5(a)]. However, the precision quickly becomes worse, once the SNR falls below 3 dB.
The detection capability of LiDAR is characterized by the peak-sidelobe ratio, which indicates how strong the peak is compared to the background of cross-correlation. It is defined as the ratio between the maximal peak height and three times the standard deviation of the background. We measure the peak-sidelobe ratio for 100 consecutive measurements, and the averaged results are plotted as a function of SNR in Fig. 5(b). With SNR higher than 3 dB, the peak-sidelobe ratio remains larger than 5 dB, above the 3 dB limit for successful detection [34].
Lastly, we investigate the detection probability [ Fig. 5(c)]. We repeat the measurement one hundred times and examine whether the target can be detected. The detection probability is defined by the fraction of successful detection over the total number of measurements. The detection probability remains nearly 100% for the SNR higher than 4 dB. However, it rapidly drops at lower SNR.

VII. PARALLEL LIDAR
Next, we conduct simultaneous ranging of different targets using multiple spatial channels. As a proof-ofprinciple demonstration, we use two spatial channels. In order to generate two separate beams whose time traces are invariant with axial propagation, we use a doublepinhole mask P to isolate two spatial channels [ Fig. 6(a)]. Each pinhole diameter is 0.5 mm, and the edge-to-edge distance between the two pinholes is 2.0 mm. Since the double-pinhole mask is placed at the same plane where the spatial correlation function of laser emission in Fig. 3(d) is measured, the spatial correlation width of 0.3 mm is much smaller than the separation of two pinholes, and the intensity fluctuations in two filtered channels are uncorrelated. Both channels are split to reference and probe beams. Two targets (aluminum mirrors with a reflectivity of 85% at 1064 nm) are axially separated by 37 cm. Channel 1 (2) probes the first (second) target O 1 (O 2 ), and the difference in propagation distance between

Max cross-correlation
FIG. 6. Parallel LiDAR with spatial multiplexing. (a) Schematic of ranging two objects O1 and O2 using two spatial channels. The laser emission passes through a dual-pinhole mask P, which generates two beams with time traces invariant to axial propagation and uncorrelated to each other. A portion of beam serves as reference (R1,2), while another portion probes the objects O1,2 and return (S1,2). (b) The crosscorrelation of time traces between every combination of S1,2 and R1,2. A peak is present only between the same channel. The distance between peaks 2∆D/c reveals the distance between two objects ∆D = 37 cm. (c) The measured distance with respect to the actual distance ∆D, showing an excellent agreement. (d) For a wide range of ∆D, the maximal crosscorrelation between S and R remains high (low) for the same (different) channels.
signal and reference beams is D 1 = 16 cm (D 2 = 53 cm). Both signal and reference of two spatial channels are measured simultaneously by four photodetectors (Newport, 818-BB-30A, bandwidth 1.5 GHz). Figure 6(b) shows the cross-correlations between the same or different channels. Between the signal and the reference of the same channel (S 1 − R 1 and S 2 − R 2 ), their cross-correlation clearly exhibits a peak. The peaks for the two channels are located at different time lags, as the distances to O 1 and O 2 are different. The spacing between two peaks is given by 2∆D/c, where ∆D = |D 2 − D 1 | is the axial distance between two objects. Figure 6(b) reveals that ∆D is 37 cm, which matches the actual distance between O 1 and O 2 .
To show that our method can detect the objects with a wide range of distances, we vary ∆D by shifting one object O 2 axially, while keeping the other object O 1 fixed. Figure 6(c) shows that the measured distance ∆D agrees well with the actual distance.
Lastly, we verify that two spatial channels are uncorrelated. Figure 6(b) shows negligible cross-correlations between the signal and reference from different channels (S 1 − R 2 and S 2 − R 1 ). This is further confirmed at a varying distance between two objects ∆D. As shown in Fig. 6(d), over a wide range of ∆D, the maximal crosscorrelation for different channels remains low, while the maximal cross-correlation between signal and reference of the same channel is significantly higher. The absence of correlation between different spatial channels is a key to parallel random LiDAR.

VIII. ROBUSTNESS TO INTERFERENCE
Finally, we investigate the anti-interference capability of our parallel random LiDAR. Interference may occur when multiple beams from different spatial channels overlap in space due to a long free-space propagation distance, which can degrade the ranging performance. We start with the simplest case of overlapping fields E 1 and E 2 from two spatial channels, where I 1 (t) = |E 1 (t)| 2 and I 2 (t) = |E 2 (t)| 2 are field intensity of two channels, and I 1,2 (t) = 2Re[E * 1 (t) E 2 (t)] represents their interference. For ranging with channel 1, the signal in Eq. (1) is cross-correlated with channel 1 reference [proportional to I 1 (t)]. Here I 2 (t) and I 1,2 (t) cause additional fluctuations in channel 1 signal. Since the interference term I 1,2 (t) contains E 1 (t), it may lead to unwanted correlation with channel 1 reference I 1 (t).
We first experimentally investigate the interference of two spatial channels. Their respective probe beams, after being reflected from separate targets, become partially overlapped at the detector. In order to enhance their spatial overlap, we use two pinholes with 0.5 mm edge-to-edge distance to select two spatial channels, each pinhole diameter is still 0.5 mm. The experimental setup is modified from that of two-channel LiDAR in Fig. 6(a). The axial distance from the beamsplitter to object O 1 is 17 cm, and to object O 2 is 26 cm. Thus the two objects are separated axially by 9 cm. The two signal beams hit O 1 and O 2 separately. The reflected beams S 1 and S 2 diffract and partially overlap at a photodetector. The measured signal S 1+2 is a mixture of S 1 and S 2 .
As shown in Fig. 7(a), the time trace of S 1+2 features a larger fluctuation amplitude than S 1 , which is the echo of O 1 with O 2 removed. The SNR of two mixed channels Cross-correlation Interference of two spatial channels. (a) Experimentally measured time trace of mixed echoes S1+2 of two objects O1 and O2, probed by separate spatial channels 1 and 2 (upper panel). The magnitude of intensity fluctuation is larger than the echo S1 of channel 1 from O1 without O2 (lower panel). (b) Cross-correlation between the mixed signal S1+2 and the reference of individual channel R1, R2. Despite notable channel interference, strong correlation peaks reveal the axial distances of O1 and O2. The time interval for computing the correlation is T = 200 ns. (c) Maximal correlations of numerically simulated time traces, showing auto-correlation of intensity fluctuation in a single channel Iα Iα = 1, and cross-correlation between two spatial channels Iα I β 0. The correlation between a single-channel intensity Iα and its interference with another channel I α,β is close to 0. All correlations are computed within 200 ns interval and averaged over 15 spatial channels. Inset: Simulated effect of detection noise on correlation as a function of SNR. Iα Iα (solid black) decreases with SNR, but Iα I β (dotted blue) and Iα I α,β (dashed red) do not change significantly. (d) Simulated maximal correlations Iα I β 0 and Iα I α,β , decaying with increasing accumulation time T of correlation in the absence of detection noise. The straight lines represent a linear fit in the log-log scale. S 1+2 is 12 dB, while that for one channel S 1 is 9 dB, which are both far below the 33 dB SNR at the saturation limit of our photodetector. The difference between S 1+2 and S 1 is 3 dB, implying that the fluctuation power for S 1+2 is twice of that for S 1 . It indicates that the echoes S 1 and S 2 from two objects O 1 and O 2 are mixed with nearly identical amplitude of fluctuation. Figure 7(b) shows cross-correlations between measured time traces. First of all, the two references R 1 and R 2 have negligible correlation, confirming the two spatial channels are uncorrelated. The mixed signal S 1+2 has correlation peaks with both references R 1 and R 2 , but at different time lags, which correspond to the axial distances of objects O 1 and O 2 . Hence, both objects are clearly distinguished despite the spatial overlap of their echoes.
To understand the experimental results, we numerically simulate cross-channel interference. A multimode laser beam is constructed by linear superposition of many transverse and longitudinal modes (see Appendix for details). The transverse mode frequencies are randomly distributed over the free spectral range. We choose 15 spatial channels at different positions α = 1, · · · , 15, and compute the time trace of the optical field in each channel E α (t). From it, we calculate the intensity trace I α (t) = |E α (t)| 2 as well as its interference with another channel I α, (1)]. The mean intensity I α (t) t is set to be identical for all channels. For every interference term I α,β (t), two channels α and β are mixed with equal amplitude.
We calculate the cross-correlation for every pair of spatial channels, and average the maximal cross-correlation over all possible pairs of 15 channels. The accumulation time for correlation is set to 200 ns, to be consistent with the experiment. Figure 7(c) shows that cross-correlation between different channels I α I β is close to zero, while autocorrelation of every channel I α I α is unity. Moreover, the correlation between the intensity trace of a single channel and its interference with another channel, I α I α,β , also vanishes. Hence, the interference between spatial channels does not create an additional peak in the correlation between signal (I α + I β + I α,β ) and reference (I α ) of the individual channel.
The channel interference effects diminish as a result of averaging correlation over a sufficient time interval T . Even when the interference term I α,β (t) has a comparable magnitude to the signal I α (t), its correlation with the reference I α (t) decreases monotonically with increasing accumulation time T . In Fig. 7(d), maximal correlations of both I α I β and I α I α,β decay with a similar slope in the double logarithmic scale.
The numerical results in Figs. 7(c,d) are obtained without noise, yet in real measurements the detection noise plays a significant role. Based on the measured noise of our photodetector, we simulate the detection noise with uncorrelated Gaussian intensity distribution. The inset of Fig. 7(c) shows that the noise has a notable impact on the autocorrelation (I α I α ), which decreases with reducing SNR. The cross-correlation (I α I β ) and the interference term (I α I α,β ) are less affected by the noise, as they are already uncorrelated before adding the noise. The entire curves of maximal cross-correlation in Fig. 7(d) shift downward in the presence of noise but its slopes are unchanged (not shown).
When the number of overlapping channels is more than two, the mixed echo with equal contribution from M channels is expressed as: I inc (t) denotes an intensity sum of all channels except α, and I coh (t) includes field interference terms among all channels. Correlating the mixed signal with channel α reference produces three terms. The correlation peak of the first term gives the distance of an object probed by channel α, and the second and third terms add fluctuating backgrounds to cross-correlation. Our numerical results reveal that the maximum of the second term I α I inc scales as √ M for M 1, and the maximum of the third term I α I coh scales as approximately M . Both drop with accumulation time T , making background correlations negligible for sufficiently long T . Hence, our parallel LiDAR is robust to cross-channel interference.

IX. DISCUSSION AND CONCLUSION
In summary, we experimentally demonstrate parallel random LiDAR by utilizing spatiotemporal intensity fluctuations from a highly multimode laser. By detuning a degenerate cavity, we lift the frequency degeneracy of transverse modes without significant reduction of quality factors, so that lasing occurs in several hundreds of transverse modes simultaneously. Thanks to their frequency diversity, the spatio-temporal interference of all lasing modes creates numerous beat frequencies densely packed in the RF domain. The complex intensity fluctuations in space and time produce several hundreds of probe beams for parallel ranging. Our method does not need any optical modulator or external feedback. Instead, a single free-running laser can generate a large set of random waveforms in parallel. For the proof-of-concept demonstration, we have used an existing degenerate cavity laser of length 40 cm. For applications such as autonomous vehicles, robots, drones, compact size is essential, and our laser cavity length may be halved by placing a mirror at the cavity center, and further reduced to a few centimeters by using a lens of shorter focal length, while keeping the same number of transverse lasing modes. The output power of our degenerate cavity laser can reach 10 W with sufficient cooling [38]. To further increase the power, a single multimode amplifier may be used to amplify the multimode laser emission.
To keep the temporal profiles of intensity fluctuations in individual spatial channels invariant with free-space propagation, we stop the spatio-temporal mode beating outside the laser cavity by isolating spatial channels with pinholes. The number of parallel spatial channels with uncorrelated intensity fluctuations is limited by the number of transverse lasing modes [46]. We expect 300 spatial channels would be available from a single laser. In this proof-of-principle demonstration, we use only two spatial channels, but it will be straightforward to further increase the number of spatial channels by using a 2D array of pinholes or coupling the emission into a fiber bundle. M collimated beams may be steered simultaneously by a single mechanical element for block scanning. Compared to the temporal or spectral multiplexing scheme which uses a single photodetector in the receiver, our spatial multiplexing scheme requires an array of M photodetectors to measure the time traces of individual spatial channels simultaneously. Alternatively, a single photodetector may record the signals from multiple spatial channels, thanks to uncorrelated intensity fluctuations in different channels [ Figs. 7(a,b)].
The broad and dense RF bandwidth of highly multimode laser emission enables a superior axial resolution of LiDAR. In our two-channel demonstration, the photodetectors have a bandwidth of 1.5 GHz, which limits the axial resolution to 5.3 cm. By using the 22-GHz bandwidth photodetector, the axial resolution can be further improved to 0.7 cm. The axial resolution is ultimately limited by the optical spectral bandwidth of the laser emission. Given the measured multimode emission bandwidth of 100 GHz, the axial resolution of our LiDAR can be as high as 1 mm. In practical situations of limited detection bandwidth, a narrowband spectral filter may be inserted into the degenerate cavity to reduce the laser emission linewidth in order to increase the detection efficiency.
We also demonstrate the robustness of our parallel ranging scheme to cross-channel interference. Although the mixing of signals from different spatial channels significantly changes the temporal profile of intensity fluctuations, cross-correlation with the reference of the individual channel can distinguish different targets. This is because the cross-channel interference terms have a negligible correlation with the intensity fluctuation of individual channels. Such anti-interference characteristics will be particularly useful in the presence of many-channel interference, making our method suitable for massively parallel LiDAR.
In addition to LiDAR, our scheme is applicable to parallel RADAR. Direct measurement of laser emission from a detuned degenerate cavity by a 2D array of photodiodes will generate distinct random waveforms to feed an array of antennas for spatially demultiplexed RADAR. Our scheme has advantages over other demultiplexing schemes of chaotic signals in radio-frequency [47], optical frequency [8,48,49], polarization [50,51], and time domain [9,52]. Time-division demultiplexing inevitably re-quires a longer acquisition time for more channels, while for spatial demultiplexing the number of channels is independent of acquisition time. Compared to spectral demultiplexing, spatial demultiplexing is able to provide many channels without sacrificing the bandwidth and resolution of the individual channel. Moreover, our source has an intrinsic bandwidth of 100 GHz and can thus generate simultaneously several independent channels in various RF bands of interest for RADAR, such as the K u band (12-18 GHz), the K a band (27)(28)(29)(30)(31)(32)(33)(34)(35)(36)(37)(38)(39)(40) and the W band (75-110 GHz), where the last is not accessible with temporal chaos from semiconductor lasers under optical feedback [53]. Without those trade-offs, our parallel ranging scheme based on a high-power many-mode laser will facilitate high-speed 3D sensing and imaging.

Multimode lasing
To simulate the laser beam in a detuned degenerate cavity, we calculate a linear superposition of many spatial modes with randomly distributed frequencies. Considering the cylindrical symmetry of the laser cavity, we choose Laguerre-Gaussian modes as the modal basis. The spatial profile of a Laguerre-Gaussian mode in a polar coordinate r = (ρ, θ) is given by [54], where p and l are the radial and azimuthal mode numbers, L |l| p (·) is the generalized Laguerre polynomial, and w 0 is the beam radius for the fundamental transverse mode. Every modal profile is normalized bý |E pl (r)| 2 dr = 1. We sort the modes E pl (r) from loworder to high-order transverse modes in ascending order of 2p + |l|, and label them by E m (r) with the transverse mode index m. For instance, m = 0, 1, 2, 3, 4, 5 correspond to (p, l) of (0,0), (0,1), (0,-1), (1,0), (0,2), (0,-2), respectively.
A linear superposition of many transverse and longitudinal modes is expressed as: E(r, t) = 0≤m<M qmin≤q<qmax E m (r)e −i(2πνm,qt+φm,q) , (4) where q is the longitudinal mode index, ν m,q is the resonant mode frequency, and φ m,q is the random phase between 0 and 2π. We sum over M = 300 Laguerre-Gaussian transverse modes, and q max −q min = 40 longitudinal modal groups. The lowest frequency ν 0,qmin = c/λ 0 is set by the emission wavelength λ 0 = 1064 nm. The free spectral range ∆ν q = ν 0,q+1 − ν 0,q is given by 375 MHz.
To simulate the detuning of a degenerate cavity, we randomize the distribution of transverse mode frequencies. The spacing between adjacent transverse modes ∆ν m,q = ν m+1,q − ν m,q follows a uniform probability distribution between 0 and 2 ∆ν m,q . The mean transverse mode spacing ∆ν m,q is set to ∆ν q /M , which makes the transverse modes almost uniformly distributed within one free spectral range.
After calculating the near-field profile E(r, t), we extract 15 time traces at a normalized distance r/w 0 = 2.9 from the cavity axis. These locations have the largest number of overlapping transverse modes, which leads to minimal long-range temporal correlations of intensity fluctuations. The time step of sampling is set to 20 ps. For Figs. 7(c-d), we average over 10 random distributions of mode frequencies.

Cross-correlation
We first normalize the field E α (t) at every channel α such that the variance of the intensity fluctuation δI α (t) = I α (t)− I α (t) t becomes unity, i.e. [δI α (t)] 2 t = 1, where I α (t) = |E α (t)| 2 is the intensity time trace. The maximal cross-correlation between two intensity time traces I α (t) and I β (t) are then calculated as, where ... t is the average within the time window of length T . The cross-correlation with any other intensity trace I x (t) [I α,β , I inc , and I coh of the main text] is also calculated by Eq. (5) but with δI β (t) replaced by δI x (t) = I x (t) − I x (t) t . The intensity time trace of noise is simulated by the independent Gaussian distribution, with the mean µ = 0 and the standard deviation σ = 10 −(s/20) , where s is the SNR in decibels. It is added to the intensity trace of signal I x (t), and the variance of the intensity fluctuation is normalized before plugging into Eq. 5.