Distortions of an optical wave propagating through the atmosphere, for example, manifested as the twinkle of stars [1], have been observed long before the invention of lasers. These distortions affect both the intensity and phase of the optical wave. The origins of wavefront deteriorations are inhomogeneities in the temperature and pressure, which cause variations of the refractive index across and along the transmission path [2]. Although laser beam propagation through free space has been extensively studied, it is a subset of the more general study of optical wave propagation through spatiotemporally random media [1,3]. Multimode optical fibers (MMFs) are another important example of such random media. In MMFs, inherent imperfections and external perturbations introduce pulse broadening and distortion [4]. The most important free-space optics (FSO) applications include communication and LiDAR. Although the functionalities of these systems are quite different, their performances all depend on the energy throughput from the source to the receiver [5,6]. MMFs, on the other hand, are potentially the best candidate for medical endoscopy since only passive components enter the body and do not suffer from the pixelated images obtained by multi-core fibers (MCFs) [7]. However, wavefront distortion presents a major obstacle toward imaging through MMFs that is yet to be overcome completely. The growing applications of FSO and MMFs necessitate a solution for wavefront distortion more than ever before.

For applications such as FSO communications where the operator has access to both the transmitter and the receiver/target, one of the existing approaches to combat turbulence and wavefront distortion is the combination of wavefront sensors and adaptive optics (AO) in conjunction with feedback controls [8]. This is not necessarily the case for LiDAR and imaging through MMF because the targets are either inaccessible or noncooperative. Therefore, the appropriate solution for the measurement of wavefront distortion should use channel estimation techniques rather than conventional optics methods.

Channel estimation has been a major task in communications and there have been many different approaches developed through the years [9]. The key to utilizing the communication channel estimation methods to resolve the distortions is to synthesize wavefront correction using modal superposition rather than directly deploying AO. In this vein, light contained in the space subtended from the transmitter to the receiver/target can be decomposed as a superposition of spatial modes, such as Hermite–Gaussian (HG) or Laguerre–Gaussian (LG) modes in free space and linearly polarized (LP) modes in MMFs [10]. To estimate the coherent transfer matrix (CTM) of an optical channel without access to its distal end, the reflections of the beam should be investigated. Rayleigh scattering is the main reflection source in free space (depending on the weather) and optical fibers. In both cases, Rayleigh scattering is caused by small-scale inhomogeneities of the local electric susceptibility, which act as induced dipole oscillators. Rayleigh scattering has been studied for several decades, and its theoretical foundations are well-established [11–13]. It was originally utilized in an optical time-domain reflectometer (OTDR) for evaluating the attenuation of optical fibers, and further utilized in a polarization OTDR (POTDR), coherent OTDR (COTDR), and phase OTDR (Φ-OTDR) for vibration sensing [14–19].

Here, for the first time (to the best of our knowledge), we characterize the CTM of an optical channel using a technique similar to COTDR. The theoretical foundation of coherent channel estimation based on distributed reflections is demonstrated in Section 2. Since the reflections originate from the naturally occurring Rayleigh scattering, the theoretical discussion can be deployed for free space, MMFs, single-mode fibers (SMFs), or any other optical channel where Rayleigh scattering exists. In Section 3, the derivation of the CTM of an SMF for two polarization modes is explained, and its corresponding experimental design and results are presented in Sections 4 and 5. Furthermore, a polarization imaging system is implemented with a very promising error of less than ∼5%.

We describe the proposed CTM characterization method using an endoscopy system, shown in Fig. 1(a), as an example. The input laser beam is split into two arms, one is used as the local oscillator (LO), and the other one is sent toward the object through an MMF. The reflected light from the object will be returned through the same MMF and combined with LO to measure its phase and amplitude deploying two optical photoreceivers [e.g., a balanced photodetector (BPD)]. Therefore, the laser beam propagates through the MMF twice and inevitably interacts with the structural fluctuations of the fiber. Such interaction with the inhomogeneities inside the fiber results in the reflection/backscattering of some photons toward the light source, which can then be deployed to extract information about the fiber properties.

The input light is the polarized output of a narrow linewidth laser; therefore, the light passing through the MMF is either polarized parallel to the scattering plane or unpolarized after passing through enough length of the fiber and experiencing polarization mode coupling [20]. So, Rayleigh scattering is the main source of reflections from the fiber [21]. Exploiting the distributed Rayleigh reflections shown in Fig. 1(b), we propose a method to estimate the CTM of the optical transmission media without access to its distal end.

Specifically, the transfer matrix of a fiber (H) is the multiplication of transfer matrices of its concatenated segments, H=∏n=1NHn,(1)where Hn denotes the transfer matrix of the nth segment as shown in Fig. 1(b). Also, S represents the coherent scattering matrix due to Rayleigh scattering. It is worth noting that any variations and couplings between different modes in forward propagation are modeled in the transfer matrices of different segments (H).

To find the transfer matrix of each segment, S is considered to be known and constant on average for all different segments of the fiber. Following the same logic as Ref. [12], suppose that electric field Ein is launched into the fiber propagating in the z direction, Ein(x,y,z,t)=E0e^in(x,y)cos[β(z−vpt)]e−12αz,(2)where e^in(x,y) is the unit amplitude field distribution of the input field, α is the power attenuation constant, β is the phase constant, and vp is the phase velocity. The scattering caused by small-scale inhomogeneity Δχ(x,y,z) of the local electric susceptibility acts as induced dipole oscillators. The inhomogeneous part of the dipole moment per unit volume is P*(x,y,z,t)=ε0Δχ(x,y,z)Ein(x,y,z,t).(3)

Consider the scattered light from a length element dzs located at Zs from the input facet of the fiber with a volume of dVs=dxdydzs. Assuming dzs≪λ so that all the dipoles are excited in phase, the dipole moment of each element will be dp=P*dVs=[ε0E0e^in(x,y)Δχ(x,y,zs)e−12αzsdVs]⏟dp^cosβ(z−vpt).(4)

Therefore, the amplitude of the scattered electric field (Es) with the unit amplitude field distribution of e^s(x,y) is |Es|∝e^s(x,y)e^in(x,y)Δχ(x,y,zs).(5)

Since Δχ can be regarded as white spatial noise independent of the input and scattering spatial mode distribution, the average scattered electric field from the cross section of the fiber is ⟨Es⟩=0.(6)

Although the white noise is zero on average, it has a finite power (proportional to its variance); therefore, Eq. (5) mathematically supports, and it is also experimentally reported before in Refs. [12,17,22] that the Rayleigh backscattering power from an optical fiber is not zero. As a result, to find the transfer matrix of each segment, we need to look for second-order equations with nonzero elements of ⟨S2⟩.

Coupling between the HG modes in free space and spatial modes in MMFs can be investigated using a generalization of principal states for polarization in SMFs [23]. The generalized principal states are the best means to describe optical beam distortions [4,24]. In the absence of turbulence, an input fundamental mode will maintain its profile as it propagates toward the receiver/target. However, turbulence is inevitable, and the input fundamental spatial mode will evolve into a superposition of many spatial modes. Beamforming techniques developed for wireless communication offer a recipe for constructing a desired wavefront at the input, which compensates for the turbulence effects through the channel and provides the desired wavefront at the target/receiver [25,26]. In what follows, we illustrate the concept of distributed coherent channel estimation, which applies to any transmission medium with an arbitrary number of vector modes.

The CTM of a transmission channel, such as an optical fiber, can be found relying on the naturally occurring Rayleigh backscattering as explained in Section 2. The overall scheme of the method is shown in Fig. 2. For simplicity, it illustrates the concept of distributed coherent channel estimation for an SMF with two polarization modes, which can be extended to an MMF supporting many vector modes. The polarization of an optical beam E=[|Ex|∡φx|Ey|∡φy] can be represented by a point on the Poincare sphere with ρ=|Ex|/|Ey| and Ψ=φx−φy. A pulse with a width of a few nanoseconds (Δt) and a known polarization state (Ein) is used as the input into the fiber-under-test. Reflections from the fiber can be divided into pulses with the width of Δt corresponding to a fiber segment with the length of Δz=cΔt/2n in which c is the speed of light and n is the effective index of the fiber. Therefore, each pulse in the reflected signal is corresponding to a segment of the fiber depending on its arrival time. The CTM of the optical fiber is the product of the CTMs of the concatenated segments shown in Fig. 2 and can be found from Eq. (1). Although a smaller Δt provides better length resolution, a larger Δt will result in higher peak power in the backreflected signal.

For the first segment, the equation should be written as ⟨Er1Er1†⟩=⟨H1TRH1Ein1Ein1†H1†R†H1T†⟩=⟨H1RH1Ein1Ein1†H1†R†H1†⟩,(7)where ⟨⟩ denotes averaging over the cross section of the fiber-under-test (performed already through photodetection), and the superscripts T and † denote transpose and conjugate transpose operators, respectively. Er1 and Ein1 are the reflected and the input electric fields of the first segment, respectively, and R is the Rayleigh reflection coefficient matrix of one segment (a specific case of the scattering matrix S mentioned in Section 2). Since this experiment is based on the Rayleigh reflections for which the scattering angle is π, R is a coefficient of the 2×2 identity matrix I [27], R=r×I=[r00r].(8)

The coupling between orthogonal polarizations in the Rayleigh backscattered light has also been studied for silica independently, concluding that the off-diagonal elements of R can be considered as insignificant and approximately zero [28,29]. In addition, Δt and Δz are selected to be small enough so that the CTM of each segment Hn is near identity, Hn=[1+δ11δ12δ211+δ22]n,|δij|≪1,δijδkh≈0.(9)

Further practical assumptions on Hn can help to solve Eq. (7). The CTM of each fiber segment is unitary (Hn†Hn=I) and symmetric since the loss is the same and negligible for both polarizations in a short piece of SMF. Therefore, it can be shown that Re{δij}=0,δ11=−δ22, and δ12=δ21. Applying the defined R and H matrices of Eqs. (8) and (9), after straightforward algebraic manipulations, Eq. (7) can be cast into solving for two unknowns δ11, δ12 besides the Rayleigh reflection coefficient ⟨|r|2⟩pulse, which depends on both the fabrication process of the fiber and the input pulse width or equivalently the length of each fiber segment, $$\begin{gathered} ⟨|r|2⟩pulse=xr=2Im{Exr1Eyr1*}Im{Exin1Eyin1*}−(|Eyin1|2−|Exin1|2)|Exr1|2+2Re{Exin1Eyin1*}Re{Exr1Eyr1*}−(|Eyin1|2−|Exin1|2)|Exin1|2+2Re{Exin1Eyin1*}Re{Exin1Eyin1*}+2Im{Exin1Eyin1*}Im{Exin1Eyin1*}, \\ Im{δ11}=Re{Exr1Eyr1*}−xrRe{Exin1Eyin1*}−4xrIm{Exin1Eyin1*},Im{δ12}=|Exr1|2−xr|Exin1|2−4xrIm{Eyin1Exin1*}. \end{gathered}$$(10)

We measured ⟨|r|2⟩pulse for the first few segments and used their average as a constant coefficient for the rest of the fiber to avoid error propagation through the entire channel estimation process. Once the CTM of the first segment is found, a similar procedure can be used to find the CTM of the next segment, and so on. Accordingly, for the nth segment, Eq. (7) can be modified as ⟨ErnErn†⟩=⟨(∏i=1nHi)R(∏i=n1Hi)⁢Ein1Ein1†(∏i=1nHi)†R†(∏i=n1Hi)†⟩,(11)where H1,2,…,n−1 are already obtained based on the previous equations and Hn is the unknown. The variable substitutions will result in similar solutions as shown in Eq. (10).

To find the CTM of an SMF with two polarization modes, the experimental setup shown in Fig. 3(a) is designed. The fiber-under-test is 96 m of AcoustiSens Wideband Vibration Sensor Fiber with 10–15 dB higher Rayleigh backscattering. Since the polarization beat length of an SMF (the length scale on which the orientation of the fiber’s axes of birefringence changes randomly) is in the order of 10 m [30], a 96 m optical fiber has a completely random polarization state at the output. As the source with a large coherence length, the NKT BASIK E15 CW laser with a linewidth of less than 100 Hz is deployed to minimize the phase noises from the source and improve the stability of the interference measurements. The laser output splits into two parts: one serves as LO and the other one passes through an electro-optic modulator (EOM), followed by an erbium-doped fiber amplifier (EDFA). Light passing through the EOM is modulated as 3 ns pulses repeating every 1.1 μs with an extinction ratio of 27.5 dB. The modulated beam passes through a high-power circulator with 56 dB directivity to isolate it from the reflections. The polarization controller (PC) on the LO path was used to maintain equal powers in the x and y polarizations. The reflected signals from the test fiber and the LO interfere in a dual-polarization 90 deg optical hybrid, and each pair of outputs is connected to a BPD to capture the electric fields in the two quadratures and two polarizations simultaneously using a four-channel real-time oscilloscope (RTO). Since the power received from the Rayleigh reflections is about 60 dB lower than the launched input power, the received signal is averaged over 3000 periods for accuracy.

Figure 3(b) depicts the obtained signal on RTO where different time slots are labeled as follows. (1)Reflection from the input facet: A 3 ns pulse with a relatively higher peak compared with the rest of the signal waveform.(2)Reflection from the output facet: A 3 ns pulse with a time difference of about 0.96 μs from the input facet reflection indicating twice the length of the fiber-under-test.(3)Rayleigh reflections from the optical fiber: Between the two mentioned peaks, a comparatively lower-level signal is observed corresponding to Rayleigh reflections from the inhomogeneities through the fiber-under-test.(4)No signal (noise): Although noise is distributed everywhere all over the signal, more specifically an almost flat level of noise can be seen between one output facet reflection and the next input facet reflections. Accordingly, there is a clear difference between the Rayleigh reflection signals and the noise-originated variations of the signal that are suppressed by averaging over multiple periods of the detected light.

Processing of the captured data results in the characterization of the SMF empowering the important application of polarization imaging through SMFs. The research on the applications of polarization imaging has a long history [31]. Polarization-based imaging and polarimetry can provide many diagnostic capabilities in clinical and preclinical studies [32] in addition to many other active sensing applications such as the determination of a target shape and orientation, or simply visualizing it through random media [33,34]. Due to the partially polarized nature of the optical and microwave fields used to form images, one can also utilize polarization information to characterize an image more completely [35]. Here, we demonstrated the experimental setup shown in Figs. 3(c) and 3(d) to validate the functionality of the proposed method of coherent channel estimation.

At first, the CTM of the SMF is estimated using the proposed distributed channel estimation method described in Sections 2 and 3. Afterward, a polarization object with known characteristics is set at the end of the fiber-under-test leading to an additional reflected pulse in the detected signal on the RTO, compared with the captured data in Fig. 3(b), corresponding to the reflection from the object. The polarization state of the light at the receiver is affected by both the polarization object and the polarization evolution through the imaging SMF. The estimated CTM of the SMF will provide enough information to acquire the polarization state of light before entering the fiber. Two separate polarization objects are designed for this experiment including a quarter-wave plate (QWP) shown in Fig. 3(c) and a vortex half-wave plate (V-HWP) shown in Fig. 3(d). In what follows, the experiments are described in more detail, and the results will be discussed.

In the experimental setup shown in Fig. 3(a), suppose that Ein=Exin∡φxinx^+Eyin∡φyiny^, Eout=Exout∡φxoutx^+Eyout∡φyouty^, and ErN=ExrN∡φxrNx^+EyrN∡φyrNy^ are the input, output, and Nth reflection electric fields in Fig. 2, respectively. Deploying Eq. (10) for all different segments of the fiber, H can be found as a symmetric matrix, and Eout is consequently estimated with two independent methods, Eout_1=H×EinandEout_2=H−1×ErN.(12)

The first one (Eout_1) utilizes the calculated H and the measured Ein while the second one finds Eout_2 based on the calculated H and the measured ErN. The two conditions in Eq. (12) must be independently satisfied, which can be used to determine the precision of the estimated CTM of the fiber. Ideally, the polarization state is expected to be the same for both Eout_1 and Eout_2, and the error from the ideal case determines the precision of the proposed channel estimation method. The two sets of values for ρ and Ψ are compared and the errors are defined as ρerr=ρ1−ρ2ρ2×100%andΨerr=Ψ1−Ψ2.(13)

We recorded 100 sets of data, and Figs. 4(a) and 4(b) show distributions of ρerr and Ψerr, respectively, for all sets of measurements. The CTM is, thus, estimated with a precision of ±0.3dB in amplitude and ±8.5° in phase based on the standard deviations of the distributions.

In the configuration of Fig. 3(c), the reflected light from the mirror passes through a linear polarizer (LP) and then transmits through a QWP before coupling into the SMF-under-test. The rotation of either QWP or LP will provide us with a controllable deterministic polarization state. With straightforward Jones matrix calculations, the degree of linear polarization of light (POL) POL=I∥−I⊥I∥+I⊥(14)is shown to be |cos(2α)| for two consecutive LP and QWP, where α is the angle between the axes of these two wave plates. If [ExEy∡φ] is the reflected signal from the object after passing through the SMF with an estimated CTM of H, the expected reflected light from the object is H−1×[ExEy∡φ],(15)where H is found using the proposed distributed channel estimation method. The expected (blue line) and measured (black dots) degree of linear polarization for the object in Fig. 3(c) is shown in Fig. 4(c) illustrating an error of less than 5.2% (normalized standard deviation of the difference) through the whole range of angle α.

In the configuration of Fig. 3(d), the reflected light from the mirror passes through an LP and then transmits through a zero-order V-HWP. Thorlabs’ liquid crystal polymer (LCP) vortex retarders are half-wave retarders designed to affect the radial and azimuthal polarization of optical fields. A vortex retarder has a constant retardance across the clear aperture while its fast axis rotates continuously over the area of the optic. The ideal fast axis angle distribution of the vortex HWP is depicted in Fig. 4(d), where the fast axis orientation rotates 180 deg over the surface of the wave plate.

Straightforward calculations will show that the Jones vector of the beam after passing through the LP and HWP is [cos2βsin2β], where β is the angle between the fast axis of the HWP and the LP axis. Therefore, by measuring the polarization state of the light and utilizing Eq. (15), we can image the distribution of the HWP fast axis orientation. The results of our measurements are shown in Fig. 4(e), where we have reported the results of imaging the POL passing through eight slices of the circular HWP vortex. Compared with the reported graph in Thorlabs specs (WPV10L-1550) and shown in Fig. 4(d), there was supposed to be a 22.5 deg rotation in the direction of the fast axis between two adjacent slices while the average rotation angle is found to be 22.1 deg in our experiment.

Reporting all these promising results for imaging and characterization, it is also worthwhile to mention that the distributed channel estimation is a single-ended method and runs with access to only the input of the optical channel independent of the distal end. This unique feature makes our approach suitable for all biomedical and military applications, and we envision having a major impact on future studies in the field.

Distributed channel estimation allows the characterization of the CTM of a coupled multimode optical channel such as optical fibers or free space. Here we deployed a pulsed signal to investigate different segments of the fiber and studied the Rayleigh reflections from each segment to obtain the CTM. The promising results for the case of two polarization modes in an SMF suggest that it is possible to generalize the same approach for the characterization of the CTM of an MMF or free-space link.

The main difference among SMF, MMF, and free space is the number of modes and, consequently, the number of equations. Solving those equations needs measurements of all or at least several of the spatial and polarization modes. However, since the origin of the reflections is Rayleigh scattering, the same as in an SMF, the relations and mathematical methods are almost the same, albeit with more tedious numerical calculations. The experimental setup will also become more complex. There are two approaches that we would propose to run such experiments. (1)Using an optical switch to map the information in different modes to be received at different time slots. The research on similar concepts of wavelength-to-time and space-to-time mapping is mature and can help the demonstration of a mode-to-time mapping scheme [36]. Since the optical channels-under-test are changing slower than a few milliseconds (kHz), the additional time consumption of the mode-to-time mapping is tolerable specifically since the speed of optical switches is getting faster every day and the commercially available ones are currently on the order of MHz.(2)Multimode and few-mode devices designed for integrated photonics are currently hot topics of research and have shown promising results [37]. Additionally, dual-polarization high-speed optical coherent transceivers have been commercialized and are available now. The combination of spatial mode demultiplexers and dual-polarization coherent receivers can provide a suitable integrated photonics platform for measurements of the distributed coherent channel estimation for optical channels with more than two modes.

The implementation of coherent channel estimation for MMFs and free space with only access to the proximal end of the channel will open up a vast area of applications including the following. (1)Imaging through multimode fibers: Modes tend to couple to each other in multimode fibers, rendering imaging through a multimode fiber with distortion. The coherent OTDR technique to be established here will enable the complete characterization of mode coupling in multimode fibers, and as a result, image distortion can be compensated with the complete knowledge of mode coupling along the multimode fibers.(2)LiDARs: An adaptive optical beam control system optimized for laser focusing through a turbulent atmosphere is one of the main challenges for this application. There have been only a limited number of solutions that address this critical issue. The proposed method enables single-ended characterization and pre-compensations of the transmitted signal and presents a paradigm shift that will yield solutions that complement and enhance the capability of currently available methods.