Single-pulse, Kerr-effect Mueller matrix LiDAR polarimeter

Christian K. Keyser; Richard K. Martin; Helena Lopez-Aviles; Khanh Nguyen; Arielle M. Adams; Demetrios Christodoulides

doi:10.1364/OE.388565

1. Introduction

Stokes vector and Mueller matrix polarimetry has been applied broadly in both laboratory and remote sensing research. Polarimeters can be used to measure the state of polarization of light (Stokes vector polarimetry) or to characterize the polarization dependent reflection, transmission, and scattering of materials (Mueller matrix polarimetry). For example, laboratory polarimetry has found utility in a diverse range of applications including biological tissue analysis [1,2], nanostructure and metamaterial characterization [3,4] and materials characterization [5,6]. In the remote sensing domain, polarimetry has been employed for active imaging contrast enhancement of scattering targets [7–9], characterization of airborne particles [10–12], and vegetation classification [13,14].

Traditional Mueller matrix polarimeter designs, such as illustrated in Fig. 1(a), employ linear polarizers and rotating quarter wave plates to measure the Mueller matrix of a sample [15,16]. In this design, a series of laser pulses are transmitted through the PSG with each pulse exiting the PSG with a different polarization state due to the rotating quarter wave plate. Likewise, in the PSA, different return polarization states are analyzed for each pulse due to the rotating retarder. A sample’s Mueller matrix is estimated through serial measurements of the signal transmitted/scattered from the sample, with each measurement employing a different illumination polarization state and receiver analyzer state. In general, to measure the full Mueller matrix, a minimum of 16 individual measurements is required and often results in a measurement duration of seconds to minutes, depending on the polarimeter design. More recent designs have employed electro-optic crystals and liquid crystals to produce the diversity of polarized transmit and analyzer states necessary to measure the Mueller matrix [17–20]. Although these alternative designs provide shorter measurement times compared to the rotating quarter wave design, they are ultimately limited by the need for multiple laser pulses for serial measurements. For measurements requiring high temporal resolution or high spatial resolution from a moving platform, multi-pulse measurement times may be unacceptable. Another drawback of conventional systems is the complex opto-mechanical designs and control electronics which drive up cost, size, weight, and power of sensors.

Fig. 1. (a) Conventional dual rotating wave plate polarimeter and (b) Proposed single pulse LiDAR polarimeter design. The PSG exploits the Kerr effect by passing the pump laser through a birefringent optical fiber. The detector processing chain uses a beam splitter (if necessary) followed by one or more sub-PSA, each followed by a detector (D1, D2, D3). Each sub-PSA may be composed of a rotated stationary quarter wave plate followed by a linear polarizer or just a linear polarizer.

Download Full Size | PDF

We propose a new architecture that simultaneously reduces hardware complexity and reduces measurement time to nanoseconds, or less, depending on laser pulse width and PSA detector bandwidth. The basic architecture is illustrated in Fig. 1(b); the PSG transmits a probe pulse with a well-defined distribution of polarization states in time and the PSA analyzes the signal scattered from the sample in time. If one knows the temporal polarization distribution of the probe pulse and the temporal polarization distribution of the scattered signal measured through the PSA, the intervening scattering sample Mueller matrix can be estimated. The probe pulse is produced by coupling a linearly polarized pump laser to a birefringent optical fiber; through the nonlinear Kerr effect, the time varying instantaneous pump power induces a time varying retardation in the fiber, which produces a temporal distribution of polarization states. Recently, we analyzed a more conventional way of producing this probe pulse through the use of an electro-optic retarder [21]. In this approach a 20 ns high-voltage ramp is applied to a Pockels cell as a laser pulse passes through. The voltage ramp induces a time varying retardation, which produces a probe pulse with a time varying polarization state. In contrast to the Pockels cell approach, instead of applying a time varying voltage, we apply a time varying pump laser power to produce a time varying laser polarization. This is particularly interesting for a PSG as it is a very simple opto-mechanical device with no electrical control and can be fusion spliced to a fiber laser to create a monolithic PSG. In addition, as will be shown below, the fiber length, pump power, and pump polarization angle can be tailored to optimize the probe pulse polarization distribution relative to measurement accuracy.

The general PSA system is of the division of amplitude type and is shown in Fig. 1(b); it employs three sub-PSAs, which are generally sufficient to characterize the Stokes vector [22]. For example, one may employ a linear polarizer at zero degrees, a linear polarizer at 45 degrees, and a circular analyzer. In Fig. 1(b), each detector block measures both the transmitted and rejected signals from each polarizer so that information about the two orthogonal states is available to construct each Stokes vector element. The detector bandwidth is chosen to allow high fidelity measurements of the polarization modulation at each sub-PSA. Due to the static birefringence axes of the fiber and Pockels cell, both of these systems are limited to measuring 12 of 16 sample Mueller matrix elements [21]. As will be shown, reducing the PSA system to a single elliptical analyzer, measuring right and left elliptically polarized light, enables measurement of six Mueller matrix elements. An elliptical analyzer is a quarter wave plate at angle theta followed by a linear polarizer. Although using just a single elliptical analyzer limits measurement of the full Mueller matrix in general, using this concept in a monostatic configuration is consistent with measuring the Mueller matrix of scattering surfaces that might be encountered in a remote sensing application as their matrices tend to be diagonal [23–28]. In what follows, we will begin with a general PSA design employing three sub-PSAs but then narrow the treatment to the case focused on diagonal Mueller matrices. This is relevant to LADAR as the world is generally composed of scattering materials having diagonal Mueller matrices and is also relevant to laboratory applications requiring characterization of surfaces.

In the remainder of this paper we will develop the underlying theory describing the polarimeter operation and experimentally analyze the concept. In Section 2 we will derive the Stokes vector generated by the Kerr effect in the birefringent optical fiber and compare it to the Stokes vector generated with a Pockels cell. In Section 3 we will develop a polarimeter system model and consider the constraints linear algebra place on Mueller matrix measurement for a single PSA receiver. In Section 4, the Cramer-Rao Lower Bound is employed to estimate measurement accuracy and optimize system parameters. In Section 5. Monte-Carlo simulations are used to quantify measurement accuracy versus signal to noise. Lastly, Section 6 outlines measurement from an experimental demonstration.

2. Nonlinear optical polarization state generator

The Kerr effect is a well-known third order nonlinear optical process associated with processes such as self-focusing, phase modulation, and four-wave mixing. Although considerable research has focused on understanding the influence of self and cross-phase modulation on the polarization state of a laser transmitted through an optical fiber, it is not clear that this knowledge has transitioned to a useful application. As will be demonstrated below, phase modulation can be leveraged as a polarization state generator to produce an optical pulse with time varying polarization for an opto-mechanically simple, single pulse polarimeter PSG.

Consider the polarization maintaining fiber in Fig. 2; it has an indices of refraction $n_x$ and $n_y$ along the $x$ and $y$ axes, respectively. The inherent birefringence in the fiber, $\delta n = n_x - n_y$, results in a phase shift proportional to $\delta n$ and fiber length $L$; it is independent of the laser power. At sufficiently high pump laser power there occurs a nonlinear phase shift proportional to the instantaneous optical power. Given the temporally varying power distribution, one can expect the pump laser to induce a temporally varying relative phase shift, or polarization.

Fig. 2. Schematic of the experimental setup to obtain evolution of polarization along time using a polarization maintaining fiber.

Download Full Size | PDF

Now consider a field of the following form propagating in the fiber:

(1)$$\vec{E}(x,y,z,t) = \left[ { \hat{x} U(z,t) + \hat{y} V(z,t) e^{{-}ik_0 (n_x-n_y)z} } \right] e^{ik_x z-i\omega_0 t}$$

where $U$, $V$ represent slowly varying envelopes of $z$, $t$, $k_0=2\pi /\lambda _0$, $k_x=k_0 n_x$, and $k_y=k_0 n_y$. We will use the following set of coupled differential equations corresponding to the limit of strong birefringence:

(2)$$\begin{aligned} & i \left( {\frac{\partial U}{\partial z} + \frac{1}{v_{gx}}\frac{\partial U}{\partial t}} \right) - \frac{\beta_{0x}^{\prime \prime}}{2} \frac{\partial^{2} U}{\partial t^{2}} + k_0 n_2 \left[ {\left| {U} \right|^{2} + A \left| {V} \right|^{2}} \right] U = 0 \end{aligned}$$

(3)$$\begin{aligned} & i \left( {\frac{\partial V}{\partial z} + \frac{1}{v_{gy}}\frac{\partial V}{\partial t}} \right) - \frac{\beta_{0y}^{\prime \prime}}{2} \frac{\partial^{2} V}{\partial t^{2}} + k_0 n_2 \left[ {\left| {V} \right|^{2} + A \left| {U} \right|^{2}} \right] V = 0 \end{aligned}$$

where $v_{gx}=c/n_{gx}$, $v_{gy}=c/n_{gy}$, with $n_{gx}$ and $n_{gy}$ representing the group refractive indices along the x,y axes. $\beta _{0x,y}^{\prime \prime }=d^{2} k/d\omega ^{2}$ denotes the dispersion coefficient evaluated at the carrier frequency $\omega _0$, $n_2$ stands for the nonlinear Kerr coefficient associated with silica glass, and $A=2/3$ for glass fibers.[29] We assume a monochromatic laser with a hyperbolic secant pulse shape of the form

(4)$$P(t) = P_0 \ \textrm{sech}^{2}\left( {\frac{t}{\tau_0}} \right)$$

where $P_0 = A_0^{2}$, $A_0$ is the peak electric field strength, and $\tau _0$ is a temporal half-width parameter. The pulse’s Full Width at Half Maximum (FWHM) is related to $\tau _0$ by

(5)$$\textrm{FWHM} = 2 \ \tau_0 \ \textrm{asech}\left( {\sqrt{0.5}} \right) \approx 1.763 \ \tau_0 .$$

We will consider an optical pulse of approximately 10 ns in duration [FWHM], so the dispersion coefficient in Eqs. (2) and (3) can be neglected. Now, we will analyze two cases: (1) a system ignoring walk-off effects, and (2) a system including walk-off effects.

2.1 System without walk-off effects

In this case, the solution of Eqs. (2) and (3) is given by:

(6)$$\begin{aligned}& U(z,\tau) = A_0 \textrm{sech}(\tau) \cos \alpha \exp\left[ {i k_0 n_2 \left| {A_0 \textrm{sech}(\tau)} \right|^{2} \left( {1 - \frac{1}{3} \sin^{2} \alpha} \right) z} \right] \end{aligned}$$

(7)$$\begin{aligned}& V(z,\tau) = A_0 \textrm{sech}(\tau) \sin \alpha \exp\left[ {i k_0 n_2 \left| {A_0 \textrm{sech}(\tau)} \right|^{2} \left( {1 - \frac{1}{3} \cos^{2} \alpha} \right) z} \right] \end{aligned}$$

where $\tau = t/\tau _0 - z/(\overline {v}_g \tau _0)$, where $\tau _0$ is related to the FWHM width of the power envelope through FWHM $\sim 1.76 \tau _0$, $\overline {v}_g$ is the average group velocity, and $\alpha$ is the angle between the linearly polarized field and the fiber x-axis (herein designated the “slow axis”) at the fiber input facet.

2.2 System with walk-off effects

In this case, the solution of Eqs. (2) and (3) is given by:

(8)$$\begin{aligned}& U(z,\tau) = A_0 \textrm{sech}(\tau-\epsilon) \cos \alpha e^{i\psi} \end{aligned}$$

(9)$$\begin{aligned}& V(z,\tau) = A_0 \textrm{sech}(\tau+\epsilon) \sin \alpha e^{i\phi} \end{aligned}$$

where $\psi$ and $\phi$ are given by:

(10)$$\begin{aligned}&\psi = l_w k_0 n_2 \left| {A_0} \right|^{2} \\ & \ \ \times \left[ {(\tau+\epsilon) \textrm{sech}^{2}(\tau-\epsilon) \cos^{2} \alpha + A \tanh(\tau+\epsilon) \sin^{2} \alpha + F(\tau-\epsilon)} \right]\end{aligned}$$

(11)$$\begin{aligned} &\phi = l_w k_0 n_2 \left| {A_0} \right|^{2} \\ & \ \ \times \left[ {(\tau-\epsilon) \textrm{sech}^{2}(\tau+\epsilon) \sin^{2} \alpha + A \tanh(\tau-\epsilon) \cos^{2} \alpha + G(\tau+\epsilon)} \right] \end{aligned}$$

where $l_w=(\tau _0 c)/(n_x-n_y)$ represents the walk-off length, $\epsilon =z\delta n/\tau _0 c$, and $F(\tau -\epsilon )$, $G(\tau +\epsilon )$ are given by:

(12)$$\begin{aligned}& F(\tau-\epsilon) ={-}(\tau-\epsilon) \textrm{sech}^{2}(\tau-\epsilon) \cos^{2} \alpha - A \tanh(\tau-\epsilon) \sin^{2} \alpha, \end{aligned}$$

(13)$$\begin{aligned}& G(\tau-\epsilon) = (\tau+\epsilon) \textrm{sech}^{2}(\tau+\epsilon) \sin^{2} \alpha + A \tanh(\tau+\epsilon) \cos^{2} \alpha . \end{aligned}$$

Figures 3 and 4 show the evolution of the Stokes parameters, defined as:

(14)$$\begin{aligned}S_0 =& \ \left| {E_x} \right|^{2} + \left| {E_y} \right|^{2} \end{aligned}$$

(15)$$\begin{aligned}S_1 =& \ \left| {E_x} \right|^{2} - \left| {E_y} \right|^{2} \end{aligned}$$

(16)$$\begin{aligned}S_2 =& \ E_x E_y^{*} + E_x^{*} E_y \end{aligned}$$

(17)$$\begin{aligned}S_3 =& \ i\left( {E_x E_y^{*} - E_x^{*} E_y} \right) \end{aligned}$$

where $E_x$, $E_y$ are obtained from Eq. (1). In Fig. 3, we observe that the walk-off effect is a slight shift of the values of the Stokes parameters along time, and it does not represent a problem for such large pulse duration (10 ns). Figure 4 shows the pulse polarization content on the Poincaré sphere.

Fig. 3. Evolution of the Stokes parameters (a) $S_1/S_0$, (b) $S_2/S_0$, and (c) $S_3/S_0$ in time at the output of a 10 m PM fiber ($\delta n=3.0 \times 10^{-4}$), for two cases: ignoring walk-off effects and including walk-off effects. The initial conditions are: pulse duration 10 ns, peak power 500 W, angle 28.6$^{o}$, beam size 3.85 $\mu$m, and wavelength of 1.06 $\mu$m.

Download Full Size | PDF

Fig. 4. Polarization content generated through Kerr effect in a 10 m PM fiber ($\delta n=3.0 \times 10^{-4}$) with pulse duration 10 ns, peak power 500 W, angle 28.6$^{o}$, beam size 3.85 $\mu$m, and wavelength of 1.06 $\mu$m.

Download Full Size | PDF

In the limit of no walk off, the Stokes vector can be written as

(18)$$\begin{aligned} \textbf{S} =& \ P(\tau) \left[ {\begin{array}{c} 1 \\ \cos(2 \alpha) \\ \cos\psi \sin\left( {2 \alpha} \right) \\ -\sin\psi \sin\left( {2 \alpha} \right) \end{array}} \right] \end{aligned}$$

(19)$$\begin{aligned} \psi =& \ \frac{k_0 z}{3} \left( {3 \delta n + P(\tau) n_2 \cos{2 \alpha} } \right) \end{aligned}$$

where we have written the instantaneous optical power as $P(\tau ) = Ao^{2} \textrm {sech}^{2}\tau$ and all Stokes components are normalized to $S_0$. The transverse coordinate $z$ should be evaluated at the fiber length $L$ to obtain the Stokes vector at the fiber output. The term $S_1$ is dependent only on the optical power $P(\tau )$ and the input polarization orientation with respect to the fiber slow axis $\alpha$. $S_2$ and $S_3$ terms include phase shift due to the birefringence of the PM fiber ($\delta n$) and a time varying phase term due to instantaneous power induced phase modulation ($P(\tau ) n_2 \cos (2 \alpha$)). If a Pockels cell is used to generate the time varying polarization probe pulse, the resulting PSG Stokes vector will have a form similar to Eq. (18) but the phase term in $S_2$ and $S_3$ will be $\psi (t) = k_0 \delta n z + \pi V(t) / V_{\pi }$. Here $V(t)$ is the time varying voltage applied to the crystal and $V_{\pi }$ is the half-wave voltage. The Kerr effect uses the time varying pulse power $P(t)$ instead of voltage $V(t)$ to achieve time varying polarization [30].

Through the Kerr phase dependence on $P(\tau ) * z$ in Eq. (19), the temporal polarization distribution can be controlled through peak power, pulse shape, and fiber length. For Mueller matrix measurement, it is important to choose $\alpha$ to balance the energy partitioned between probe beam Stokes components in order to sample the target Mueller matrix components with sufficient sensitivity. The polarization angle, $\alpha$, yielding equal power in each of the Stokes polarization states of Eq. (18) can be found from the conservation of energy $\sqrt {S_1^{2} + S_2^{2} + S_3^{2}} = S_0$ and by setting $S_1=S_2=S_3$. From this consideration one finds $\alpha \approx 27^{o}$. Lastly, although Eq. (18) is written in the coordinate system of the fiber birefringence axes, it may be desirable to rotate the coordinates by $\alpha$ into the laboratory reference frame [24,31], leading to

(20)$$\textbf{S}_{\textrm{lab}} =P(\tau) \left[ {\begin{array}{c} 1 \\ \cos^{2}(2 \alpha) + \sin^{2}\left( {2 \alpha} \right) \cos\psi \\ -\frac{1}{2}\sin\left( {4 \alpha} \right) + \frac{1}{2}\sin\left( {4 \alpha} \right) \cos\psi \\ -\sin\psi \sin\left( {2 \alpha} \right) \end{array}} \right] .$$

2.3 Model of receiver hardware

The receiver diagram in Fig. 1 allows for the possibility of using a beam splitter followed by multiple PSAs. There is a wide variety of PSAs that could be used; we consider only the following options for brevity:

1. Horizontal/Vertical linear polarizers (denoted H/V).
2. Diagonal ($+45^{o}$ and $-45^{o}$) linear polarizers (denoted P/M).
3. A Quarter-Wave Plate (QWP) rotated by $\theta =45^{o}$ relative to the lab frame followed by an H/V linear polarizer (not rotated); this is effectively a left/right circular polarization analyzer (denoted L/R) [24,31].

In each pair of outputs, the first (e.g. “H”) is the polarization which is transmitted through the cube polarizer and the other (e.g. “V”) is rejected from the polarizer and is propagated on a separate optical path to another detector. These three options provide direct measurement of the Stokes vector elements. We will explore all three options when practical, though the nominal configuration will be just option 3 with no beam splitter, with the QWP rotation angle $\theta$ left unspecified, which allows for generalization to an elliptical L/R PSA rather than simply a circular L/R PSA.

Considering the atmospheric transmission $A(R)$, the transmitter efficiency $T_t(p)$, the receiver efficiency $T_r(p)$, and the detector response $D_{\textrm {det}}$, the terms affecting the received LiDAR signal can be grouped as

(21)$$k(R,p) = \frac{1}{R^{2}} A(R) \cdot T_t(p) \cdot T_r(p) \cdot D_{\textrm{det}}$$

where $R$ is the one-way range to target and $p$ denotes the detector polarization channel (H, V, P, M, R, L). Note that the detector response here is the optical response; the temporal detector response will be accounted for later. The terms $T_t(p)$, $T_r(p)$, and $D$ can be calibrated at the instrument level.

Let $D = 2R/c$ indicate the round-trip travel time to target, where $c$ is the speed of light. The signal exiting polarization channel $p$ is

(22)$$\textbf{S}_{\textrm{rx}}(t,p) = k(R,p) Q(R) \cdot \textbf{M}_p \cdot \textbf{M} \ \cdot \textbf{S}_{\textrm{lab}}\left( t - D \right)$$

where $Q(R)$ is the laser beam-target overlap integral, $\textbf {M}_p$ is the known Mueller matrix [24,31] of polarization analyzer $p$, and $\textbf {M}$ is the unknown Mueller matrix of the target. The scalar $Q(R)$ and the factor $A(R)$ within $k(R,p)$ are potentially unknown but can be canceled out later by normalization the estimated Mueller matrix. Unless otherwise noted, throughout this section $\textbf {S}_{\textrm {lab}}$ can represent the transmitted Stokes vector from either the Kerr-based or Electro-Optic (EO)-based transmitter.

Leaving $\theta$ unspecified for option 3, the Mueller matrices for the PSAs considered here are [24,31]

(23)$$\textbf{M}_{\textrm{h/v}} = \frac{1}{2} \begin{bmatrix} 1 & \pm 1 & 0 & 0 \\ \pm 1 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix}$$

(24)$$\textbf{M}_{\textrm{p/m}} = \frac{1}{2} \begin{bmatrix} 1 & 0 & \pm 1 & 0 \\ 0 & 0 & 0 & 0 \\ \pm 1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix}$$

(25)$$\textbf{M}_{\textrm{r/l}} = \frac{1}{2} \begin{bmatrix} 1 & \pm \cos^{2} 2 \theta & \pm \cos 2 \theta \sin 2 \theta & \mp \sin 2 \theta \\ \pm 1 & \cos^{2} 2 \theta & \cos 2 \theta \sin 2 \theta & -\sin 2 \theta \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix}$$

Now we can find the output intensities (i.e. the first element of the Stokes vector) from each possible PSA pair. Let $N \in \left \{{2,4,6}\right \}$ be the number of signals to be recorded by the fast detectors in Fig. 1; for example, $N=2$ for any single analyzer. We assume that the beam splitter evenly divides the input intensity into $N/2$ equal streams, with a fraction $2/N$ of the intensity going into each PSA pair. Using (23)–(25) the possible output intensities are

(26)$$\begin{aligned} & y_{\textrm{h/v}}(t) = \frac{k_{\textrm{h/v}}Q}{N} \sum_{j = 0}^{3} S_{\textrm{lab},j}(t-D) \left( M_{0j} \pm M_{1j} \right) \end{aligned}$$

(27)$$\begin{aligned} & y_{\textrm{p/m}}(t) = \frac{k_{\textrm{p/m}}Q}{N} \sum_{j = 0}^{3} S_{\textrm{lab},j}(t-D) \left( M_{0j} \pm M_{2j} \right) \end{aligned}$$

(28)$$\begin{aligned} & y_{\textrm{r/l}}(t) = \frac{k_{\textrm{r/l}}Q}{N} \sum_{j=0}^{3} S_{\textrm{lab},j}(t-D)\\ & \ \ \ \ \times \left( { M_{0j} \pm M_{1j} \cos^{2} 2 \theta \pm M_{2j} \cos 2 \theta \sin 2 \theta \mp M_{3j} \sin 2 \theta } \right) \end{aligned}$$

The estimation process must be carried out using a discrete time representation. We will model the conversion from continuous to discrete time by filtering the detector output though a filter with impulse response $h(t)$ and frequency response $H(f)$ and then instantaneously sampling the filter output with sampling period $T$. We assume that the filter’s frequency response has cut-off frequency $F_c$; ideally $H(f)$ will be a rectangle, but we need not assume that. Herein, we use the term “detector bandwidth” to refer to $F_c$.

We will use $k$ to index the sample time and $K$ samples are collected per polarization channel. The effects of the detector filter, sampling, and additive noise $v[k]$ can be represented by

(29)$$y_p[k] = {\left. {\left[ y_p(t) \star h(t) \right]} \right|_{t = kT}} + v[k]$$

where $\star$ denotes convolution. To account for the time dependence and the effects of this convolution, define

(30)$$s_j(t,D) = S_{\textrm{lab},j}(t-D) \star h\left( t \right), \ \ \ j \in \left\{{0,1,2,3}\right\}$$

using the transmitted Stokes vector from (18), where $j$ is the index of the four elements of the Stokes vector. Then define sampled versions, stacked into four $K \times 1$ vectors,

(31)$$\textbf{s}_j(D) = \left[ s_j\left( {0 \cdot T} \right), s_j\left( {1 \cdot T} \right),\ldots, s_j\left( {\left[ {K - 1} \right] \cdot T} \right) \right]^\textrm{T}, \ \ \ j \in \left\{{0 \ldots 3}\right\} .$$

The samples from the polarization channel outputs can be represented in vector form by

(32)$$\begin{aligned} \textbf{y}_p =& \ \left[ y_p(0),\ldots, y_p([K-1]T) \right]^\textrm{T} \end{aligned}$$

(33)$$\begin{aligned} \textbf{y} =& \ \left[ \textbf{y}_{p_1}^\textrm{T}, \ldots, \textbf{y}_{p_N}^\textrm{T} \right]^\textrm{T} \end{aligned}$$

Assume (temporarily) that all 6 polarization channels are used and that all 16 Mueller matrix elements are to be estimated. Then counting both outputs of each of (26)–(28), there are six output channels that are each sampled $K$ times by the detectors. Stacking these samples per (32)–(33), the overall system could be written as a linear model with respect to the unknown elements $M_{ij}$, in the form

(34)$$\begin{aligned} \textbf{y} =& \ \textbf{H}(D) \ \cdot \ \textbf{x} + \textbf{v} \end{aligned}$$

(35)$$\begin{aligned} \textbf{x} =& \ Q\left[ M_{00}, M_{01}, M_{02}, M_{03}, M_{10}, M_{11}, \ldots, M_{33} \right]^\textrm{T} , \end{aligned}$$

where $\textbf {v}$ is a vector of samples of the additive noise process. We have included the unknown scale factor $Q$ within the unknown vector $\textbf {x}$; the method for dealing with $Q$ will depend on the application for which the estimates of $\textbf {x}$ will be used. A detailed inspection of (26)–(28) is necessary to reveal the structure of $\textbf {H}$. It takes the form

(36)$$\begin{aligned} & \textbf{H}(D) = \textbf{H}_k \otimes \textbf{H}_s(D) \end{aligned}$$

(37)$$\begin{aligned} & \textbf{H}_k = \frac{1}{N} \begin{bmatrix} k_{\textrm{h}} & k_{\textrm{h}} & 0 & 0 \\ k_{\textrm{v}} & -k_{\textrm{v}} & 0 & 0 \\ k_{\textrm{p}} & 0 & k_{\textrm{p}} & 0 \\ k_{\textrm{m}} & 0 & -k_{\textrm{m}} & 0 \\ k_{\textrm{r}} & k_{\textrm{r}} \cos^{2} 2 \theta & k_{\textrm{r}} \cos 2 \theta \sin 2 \theta & -k_{\textrm{r}} \sin 2 \theta \\ k_{\textrm{l}} & -k_{\textrm{l}} \cos^{2} 2 \theta & -k_{\textrm{l}} \cos 2 \theta \sin 2 \theta & k_{\textrm{l}} \sin 2 \theta \\ \end{bmatrix} \end{aligned}$$

(38)$$\begin{aligned} & \textbf{H}_s(D) = \left[ \textbf{s}_0(D), \ \textbf{s}_1(D), \ \textbf{s}_2(D), \ \textbf{s}_3(D) \right] \end{aligned}$$

where $\otimes$ denotes a Kronecker product. In a simpler scenario where only $N$ of the 6 possible polarization channels are used and only $M$ of the 16 possible Mueller matrix elements are assumed to be non-zero, we can omit rows and columns of $\textbf {H}$ via

(39)$$\textbf{y} = \ \textbf{H}(D) \cdot \textbf{x}_{\textrm{reduced}} + \textbf{v}$$

(40)$$\begin{aligned} \textbf{H}(D) =& \ \left[ {\left( { \mathbb{1}_P \textbf{H}_k} \right) \otimes \textbf{H}_s(D)} \right] \mathbb{1}_X^{T} \end{aligned}$$

(41)$$\begin{aligned}\textbf{x}_{\textrm{reduced}} =& \ \mathbb{1}_X \textbf{x} \end{aligned}$$

where $\mathbb{1}$ denotes an “indicator” matrix. The $N \times 6$ polarization channel indicator $\mathbb{1}_P$ is obtained by deleting all but the $N$ relevant rows of a $6 \times 6$ identity matrix, and the $M \times 16$ Mueller matrix indicator $\mathbb{1}_X$ is obtained by deleting all but the $M$ relevant rows of a $16 \times 16$ identity matrix. Left-multiplying by $\mathbb{1}_P$ or $\mathbb{1}_X$ will select only the $P$ or $M$ rows of interest, and right-multiplying by $\mathbb{1}_X^{T}$ will select only the $M$ columns of interest.

Note that $\textbf {H}_k$ depends strongly on the rotation of the receiver’s QWP, $\theta$; and $\textbf {H}_s$ depends heavily on the transmitter parameters (i.e. for a Kerr system, the optical fiber parameters $\delta _n$, $n_2$, and $\alpha$; or for an EO system, the rotation of the transmitter’s EO retarder, $\alpha$, as well as the voltage ramp parameters $V_{\textrm {max}}$ and $\beta$, the ramp rate). Specifically, ignoring the detector’s temporal impulse response,

(42)$$\begin{aligned} \textbf{H}_s = & \ \begin{bmatrix} \textbf{P} & \textbf{P} \odot \cos {\boldsymbol{\psi}} & \textbf{P} \odot \sin {\boldsymbol{\psi}} \end{bmatrix}\\ & \ \ \ \ \times \begin{bmatrix} 1 & \cos^{2} 2 \alpha & \frac{1}{2} \sin 4 \alpha & 0 \\ 0 & \sin^{2} 2 \alpha & \frac{-1}{2} \sin 4 \alpha & 0 \\ 0 & 0 & &0 & \sin 2 \alpha \end{bmatrix} \end{aligned}$$

where $\odot$ is Hadamard (element-by-element) multiplication; and $\textbf {P}$ and ${\boldsymbol {\psi}} $ contain $K$ samples of the instantanteous optical power $P(\tau )$ and the phase argument from the $S_2$ and $S_3$ components of (18), respectively. The time samples are evaluated at times

(43)$$t_n = (n-1) T - D, \ \ \ n \in \left\{{1,\ldots,K}\right\} .$$

2.4 Noise model

Since the detector output voltage is proportional to incident power, we define the Signal to Noise Ratio (SNR) as a ratio of detector voltage values. However, additive Gaussian noise present within the detector can cause the detector voltage to go negative. As such, we compute a ratio of absolute values of average detector voltages, leading to

(44)$$\begin{aligned} \textrm{SNR} =& \ \frac{ \frac{1}{NKT} \sum_{p=1}^{N} \int y_p(t) dt}{{\textrm{E}\left\{ {\left| {v(t)} \right|} \right\}}}\\ \approx & \ \frac{\frac{1}{N} \sum_{p=1}^{N} \sum_{k=0}^{K-1} y_p[k] }{\sqrt{\frac{2}{\pi}} K \sigma } \end{aligned}$$

where $\sigma$ is the standard deviation of the detector noise $v[k]$. Since $\sigma$ is in the same units as the detector voltage, which is proportional to incident power, the denominator of (44) uses $\sigma$ rather than $\sigma ^{2}$, which ensures that this SNR definition is a power ratio. We have also made use of the fact that the average magnitude of a Gaussian random variable is $\sqrt {2/\pi }$ times the standard deviation of that Gaussian. Note that this SNR definition is with respect to the received signal, so factors such as atmospheric transmission and target reflectivity are included in the numerator.

In order to gain intuition as to how the various parameters affect the SNR, assume that during the observation window, the detector captures all but a negligible amount of energy of the received pulse, and that $M_{00} >> M_{ij}$ for $i,j \ne 0$ in (26). The first assumption is necessary anyway for accurate parameter estimation; the latter assumption is not strictly accurate, but $M_{00}$ is always the largest element and all of the other elements in the SNR computation will get multiplied by various sines and cosines which will further attenuate their effects, so this assumption is reasonable for the purpose of getting a sense of how various parameters affect the SNR. Using these assumptions,

(45)$$\textrm{SNR} \approx \frac{Q M_{00}}{N K T \sigma} \sqrt{\frac{\pi}{2}} \frac{\sum_p k_p}{N} \int \left| {S_\textrm{lab,0}(t)} \right| dt .$$

(There are two factors of $N$ in the denominator; one is from within each $y_p(t$) due to the beam splitter, and the other is because the SNR numerator is defined as an average across the $N$ polarization channels.) The zeroth element in the transmitted stokes vector is just the pulse power. For the sech$^{2}$ pulse model in (4),

(46)$$\int_{-\infty}^{\infty} \left| {S_\textrm{lab,0}(t)} \right| dt = P_0 \ \int_{-\infty}^{\infty} \textrm{sech}^{2}\left( {\frac{t}{\tau_0}} \right) dt = 2 P_0 \tau_0 .$$

Thus, the SNR can be approximated as

(47)$$\begin{aligned} & \textrm{SNR} \approx \frac{\sum_p k_p}{N} Q M_{00} \sqrt{\frac{\pi}{2}} \frac{2 \tau_0}{N K T} \frac{P_0}{\sigma}\\ & \ \ \approx \underbrace{\frac{1.42}{N} \cdot \frac{\sum_{p=1}^{N} k_p}{N} \cdot Q \cdot M_{00}}_{A} \cdot \underbrace{\frac{\textrm{FWHM}}{K T}}_{B} \cdot \underbrace{\frac{P_0}{\sigma}}_{C} . \end{aligned}$$

This holds for either the proposed Kerr system or for a similar EO-based system, so long as it uses the same pulse shape. Each of the terms $A$, $B$, and $C$ is unitless. Term $A$ is essentially the propagation loss through the optical path, term $B$ is a ratio of the laser pulse width over the duration of the range gate, and term $C$ is the ratio of the pulse peak power over the standard deviation of the detector noise. To maintain a unitless ratio, both $P_0$ and $\sigma$ should be represented in the same units, such as via the detector voltages they induce. The propagation loss terms $k(R,p)$ from (21) should also be unitless. Due to the assumption that $M_{00} >> M_{ij}$ for $i,j \ne 0$, the SNR expression in (47) is only correct to within perhaps a factor of two, but it can be used to gauge how key parameter values will affect the SNR. The simulations in Section 5 will use (44) to evaluate the SNR for a given discrete time waveform.

3. Linear algebra analysis

We now consider linear algebra concerns that potentially limit receiver operation. Throughout, we consider the special case where only the elliptical (R/L) polarization channels are used. In that case,

(48)$$\begin{aligned}& \mathbb{1}_P = \left[ { \textbf{0}_{2 \times 4} \ \textbf{I}_{2 \times 2} } \right] \end{aligned}$$

(49)$$\begin{aligned}& \mathbb{1}_P \textbf{H}_k = \end{aligned}$$

(50)$$\begin{aligned}& \ \ \ \ \frac{1}{2} \begin{bmatrix} k_{\textrm{r}} & k_{\textrm{r}} \cos^{2} 2 \theta & k_{\textrm{r}} \cos 2 \theta \sin 2 \theta & -k_{\textrm{r}} \sin 2 \theta \\ k_{\textrm{l}} & -k_{\textrm{l}} \cos^{2} 2 \theta & -k_{\textrm{l}} \cos 2 \theta \sin 2 \theta & k_{\textrm{l}} \sin 2 \theta \\ \end{bmatrix} \end{aligned}$$

Note that $\mathbb{1}_P \textbf {H}_k$ has a rank of 2; but regardless of the value of $\alpha$ (for either an EO or Kerr-based system), $\textbf {H}_s$ only has three linearly independent columns rather than four. This comes from (42), since the matrix has a factor with only 3 columns and thus has a maximum possible rank of 3. Thus, the overall rank of $\textbf {H}$ is $2 \times 3 = 6$, meaning that the maximum number of Mueller matrix elements we can recover is six. In order to transform this into a well-posed problem, redundant columns of $\textbf {H}$ must be removed by assuming some of the elements of the target Mueller matrix are zero. In the remainder of this section, we discuss the possible ways this can be done. Alternatively, the problem could also be resolved by adding additional polarization analyzers, at the cost of increased receiver complexity.

Inspection of the 16 columns of $\textbf {H}$ reveals that they are spanned by the six vectors listed in the top row of Table 1. The terms $\cos {\boldsymbol {\psi}}$ and $\sin {\boldsymbol {\psi}}$ are shorthand for the length-$K$ vectors of samples of $\cos \psi (t)$ and $\sin \psi (t)$, and the scale factors of $P(t)$ are omitted for notational brevity. Table 1 also shows the dependence of each column of $\textbf {H}$ on the six spanning vectors. A necessary condition for a well-posed problem is that the elements depend on as many spanning vectors as there are elements to be estimated. For example, we cannot estimate both $M_{13}$ and $M_{23}$ because they both depend only on the same spanning vector, $\textbf {b}_6$.

Table 1. Dependence of the columns of $ {\textbf H}$ on the six basis vectors $\left\{{ {\textbf b}_n}\right\}$. The $M_{ij}$ notation indicates the Mueller matrix elements corresponding to the indicated columns of $ {\textbf H}$. Green entries are the elements to be estimated when a diagonal Mueller matrix is assumed. Regarding the blue elements, refer to the text between (52) and (53).

View Table

An obvious choice for reducing the parameter space is to assume that $\textbf {M}$ is purely diagonal, i.e. that all elements of $\textbf {x}$ are zero except for #1, #6, #11, and #16. Columns of $\textbf {H}$ corresponding to null elements of $\textbf {x}$ can be deleted, resulting in

(51)$$\begin{aligned} & \textbf{H}_{\textrm{diag}} = \\ & \ \ \frac{1}{2} \begin{bmatrix} k_{\textrm{r}} \textbf{s}_0 & k_{\textrm{r}} \textbf{s}_1 \cos^{2} 2 \theta & k_{\textrm{r}} \textbf{s}_2 \cos 2 \theta \sin 2 \theta & -k_{\textrm{r}} \textbf{s}_3 \sin 2 \theta \\ k_{\textrm{l}} \textbf{s}_0 & -k_{\textrm{l}} \textbf{s}_1 \cos^{2} 2 \theta & -k_{\textrm{l}} \textbf{s}_2 \cos 2 \theta \sin 2 \theta & k_{\textrm{l}} \textbf{s}_3 \sin 2 \theta \end{bmatrix}\\ \end{aligned}$$

(52)$$\begin{aligned} & \textbf{x}_{\textrm{diag}} = Q\left[ M_{00}, M_{11}, M_{22}, M_{33} \right]^\textrm{T} . \end{aligned}$$

As mentioned in Section 1, it is generally accepted that the Mueller matrix for a scattering sample is diagonal, so this is a reasonable assumption for analysis of scattering samples. The resulting spanning vector dependence is highlighted in green in Table 1; limiting the problem to those four elements results in dependence on four unique spanning vectors, so the problem is well posed. This is why the special case considered in this section is our default configuration – it allows recovery of the full Mueller matrix diagonal with only a single PSA and a single pair of detectors for its two outputs.

Table 1 can also be used to show that it is (surprisingly) absolutely essential to retain both outputs of the linear polarizer, even under the assumption of a diagonal Mueller matrix. If only one polarization channel is retained, the basis vectors are each halved to only either the first $K$ or last $K$ elements, which would make $\left \{{ \textbf {b}_4, \textbf {b}_5, \textbf {b}_6}\right \}$ redundant with $\left \{{ \textbf {b}_1, \textbf {b}_2, \textbf {b}_3}\right \}$. As such, the rank of $\textbf {H}$ would be 3 yet there would be 4 unknowns, making the problem ill-posed.

We may wish to recover additional Mueller matrix elements, since up to six can be recovered. Inspection shows that the potential choices are extremely limited – we can include $M_{03}$, and we can include either $M_{01}$ or $M_{02}$; these elements are highlighted in blue in Table 1. Any other choices lead to an ill-conditioned problem. These two choices can be formulated as

(53)$$\begin{aligned} \textbf{H}_{\textrm{option2}} =& \ \textbf{H}\left( {:,[1,2,4,6,11,16]} \right) \end{aligned}$$

(54)$$\begin{aligned} \textbf{x}_{\textrm{option2}} =& \ Q\left[ M_{00}, M_{01}, M_{03}, M_{11}, M_{22}, M_{33} \right]^\textrm{T} \end{aligned}$$

(55)$$\begin{aligned} \textbf{H}_{\textrm{option3}} =& \ \textbf{H}\left( {:,[1,3,4,6,11,16]} \right) \end{aligned}$$

(56)$$\begin{aligned} \textbf{x}_{\textrm{option3}} =& \ Q\left[ M_{00}, M_{02}, M_{03}, M_{11}, M_{22}, M_{33} \right]^\textrm{T} \end{aligned}$$

For readability, our notation assumes that one computes all 16 columns of $\textbf {H}$ using the Kronecker product structure of (40) and then deletes the columns corresponding to the removed elements of $\textbf {x}$, though that is not necessarily the most computationally efficient approach.

In summary, for the default single PSA configuration, the only choices of Mueller matrix elements to estimate that lead to a well-conditioned problem are the three choices in (52), (54), or (56); or any subset thereof. All other Mueller matrix elements are implicitly assumed to be zero. In Section 5 we will evaluate the effects of making that assumption when the assumed-to-be-null values are small but non-zero.

4. Cramer-Rao bound analysis

We can use the Cramer-Rao Lower Bound (CRLB) to approximate the Mueller matrix estimation performance. Throughout this section, we assume that the Mueller matrix is diagonal, and that only the elliptical (R/L) channels are used.

As in prior work [21,32], we will assume the round-trip delay $\tau$ is known for the purposes of the CRLB, though it will be considered unknown later. For a model of the form $\textbf {y} = \textbf {H} \textbf {x} + \textbf {v}$ where $\textbf {v}$ is white Gaussian noise, the CRLB is the inverse of the Fisher Information Matrix (FIM) $\textbf {J}$, given by [33]

(57)$$\begin{aligned} Cov\left[\hat{ \textbf{x}}\right] \ge & \ \textbf{J}^{{-}1} , \end{aligned}$$

(58)$$\begin{aligned} \textbf{J} =& \ \sigma^{{-}2} \ \textbf{H}^{T} \textbf{H} , \end{aligned}$$

where $\sigma ^{2}$ is the variance of the noise from (44).

We will use the CRLB to gain some intuition regarding the effects of the angle $\alpha$ (either the EO retarder rotation or the Kerr fiber rotation) and the QWP rotation $\theta$. Assuming $k_{\textrm {h}} = k_{\textrm {v}}$ and ignoring constants such as $k_{\textrm {h}}$, $k_{\textrm {v}}$, $1/2$, and $\sigma ^{2}$, we have

(59)$$\textbf{H} = \begin{bmatrix} \textbf{s}_0 & \textbf{s}_1 \cos^{2} 2 \theta & \textbf{s}_2 \cos 2 \theta \sin 2 \theta & - \textbf{s}_3 \sin 2 \theta \\ \textbf{s}_0 & - \textbf{s}_1 \cos^{2} 2 \theta & - \textbf{s}_2 \cos 2 \theta \sin 2 \theta & \textbf{s}_3 \sin 2 \theta \end{bmatrix} .$$

For brevity of notation, define

(60)$$\begin{array}{cc} c_{ \theta} = \cos 2 \theta , & s_{ \theta} = \sin 2 \theta . \end{array}$$

Using these definitions,

(61)$$\textbf{J} = 2 \begin{bmatrix} \textbf{s}_0^{T} \textbf{s}_0 & \textbf{s}_0^{T} \textbf{s}_1 \ c_{\theta}^{4} & \textbf{s}_0^{T} \textbf{s}_2 \ \ c_{\theta} s_{\theta} & - \textbf{s}_0^{T} \textbf{s}_3 \ s_{\theta} \\ \textbf{s}_0^{T} \textbf{s}_1 \ c_{\theta}^{4} & \textbf{s}_1^{T} \textbf{s}_1 \ c_{\theta}^{4} & \textbf{s}_1^{T} \textbf{s}_2 \ \ c_{\theta}^{3} s_{\theta} & - \textbf{s}_1^{T} \textbf{s}_3 \ c_{\theta}^{2} s_{\theta} \\ \textbf{s}_0^{T} \textbf{s}_2 \ \ c_{\theta} s_{\theta} & \textbf{s}_2^{T} \textbf{s}_1 \ c_{\theta}^{3} s_{\theta} & \textbf{s}_2^{T} \textbf{s}_2 \ c_{\theta}^{2} s_{\theta}^{2} & - \textbf{s}_2^{T} \textbf{s}_3 \ c_{\theta} s_{\theta}^{2} \\ - \textbf{s}_0^{T} \textbf{s}_3 \ s_{\theta} & - \textbf{s}_3^{T} \textbf{s}_1 \ c_{\theta}^{2} s_{\theta} & - \textbf{s}_3^{T} \textbf{s}_2 \ c_{\theta} s_{\theta}^{2} & \textbf{s}_3^{T} \textbf{s}_3 \ s_{\theta}^{2} \end{bmatrix} .$$

Now we can gain some intuition by noting the dependence between the $\textbf {s}$ vectors via the linear dependence among the elements of (20). For example,

(62)$$\begin{aligned}\textbf{s}_0^{T} \textbf{s}_1 =& \ \sum_k \left( {P(kT)} \right) \left( {-P(kT) \sin \phi(kT) \sin 2 \alpha} \right) \end{aligned}$$

(63)$$\begin{aligned}=& \ - \sin 2 \alpha \sum_k P^{2}(kT) \sin \phi(kT) \approx 0 \end{aligned}$$

assuming the $\sin \phi (t)$ term oscillates faster than the rise and fall of the pulse shape. Using the generic symbol $\#$ to indicate a non-zero value,

(64)$$\textbf{J} \approx \begin{bmatrix} \# & 0 & 0 & 0 \\ 0 & \# & \# & 0 \\ 0 & \# & \# & 0 \\ 0 & 0 & 0 & \# \end{bmatrix} .$$

That is, elements $M_{11}$ and $M_{22}$ are coupled with each other, whereas $M_{00}$ and $M_{33}$ are uncoupled. This is because element $i,j$ of $\textbf {J}^{-1}$ is a (very tight) bound on the covariance between $M_{i,i}$ and $M_{j,j}$, and the inverse of a block diagonal matrix maintains the same block diagonal structure. This will be helpful in interpreting simulation results later.

Figures 5 and 6 show a numerical evaluation of the CRLB for the EO and Kerr systems, respectively. The target Mueller matrix is known to be diagonal. The Mueller matrix data used was for oak bark at $\lambda = 1.06 \ \mu$m, at normal incidence [6]. The SNR was set to 10 dB using (44); since the SNR is system parameter dependent, the waveforms used to compute the noise variance $\sigma ^{2}$ were generated using globally optimum values ($\alpha = 27.0^{o}$ and $\theta = 18.8^{o}$ for the EO system and $\alpha = 28.6^{o}$ and $\theta = 17.2^{o}$ for the Kerr system). Both systems used FWHM = 10 ns. The nominal values of the EO parameters were $V_{\textrm {max}} = 1344$ V, and $\beta = 0.1$ ns$^{-1}$, based on Quantum Technology Starfire driver and LiTaO3 crystal. The nominal values of the Kerr parameters were refractive index $n_0 = 1.45$, bifrefringence $\delta _n = 3.5 \cdot 10^{-4}$, nonlinear Kerr coefficient $n_2 = 1.2 \cdot 10^{-22}$ m$^{2}$/V$^{2}$, fiber length $L=1$ m, fiber radius $w_i=3.5 \ \mu$m, and power $P=5$ kW.

Fig. 5. CRLB vs. EO retarder rotation $\alpha$ and QWP rotation $\theta$, for SNR = 10 dB, $V_{\textrm {max}} = 1344$ V, and $\beta = 0.1$ ns$^{-1}$. This assumes a purely diagonal Mueller matrix. The middle and right plots are cross-sections of the contour data in the left plot. The plots are symmetric every $45^{o}$ in each angle, so only $0 \le \alpha ^{o} \le 45^{o}$, $0 \le \theta ^{o} \le 45^{o}$ is shown.

Download Full Size | PDF

Fig. 6. CRLB vs. Kerr fiber rotation $\alpha$ and QWP rotation $\theta$, for SNR = 10 dB. This assumes a purely diagonal Mueller matrix. The middle and right plots are cross-sections of the contour data in the left plot. The plots are symmetric every $45^{o}$ in each angle, so only $0 \le \alpha ^{o} \le 45^{o}$, $0 \le \theta ^{o} \le 45^{o}$ is shown.

Download Full Size | PDF

The first plot shows contours of the CRLB, averaged across the main diagonal of $\textbf {J}^{-1}$ (corresponding to the four Mueller matrix diagonal elements), as a function of $\alpha$ and $\theta$. The second and third plots show cross-sections of the first plot, taken across the global minimum. Note that all plot values will scale linearly with $\sigma$; i.e. if the SNR was increased by a factor of two (or 3 dB) because sigma was halved, then the mantissa would be halved but the curve shape would remain the same. The plots are symmetric every $45^{o}$ for each angle, so only ranges of $[0,45^{o}]$ are shown. Comparison of Figs. 5 and 6 indicates the EO (Pockels cell) and Kerr based PSG systems perform similarly if proper system parameters are chosen.

5. Performance simulations

5.1 Parameter estimation algorithm

The parameter estimation method that will be employed is in the family of least squares (LS) and Maximum Likelihood (ML) estimation; specifically, a separable least-squares style of estimation can be used [33]. The idea is to assume candidate values for any unknown parameters that have a nonlinear relationship with the observations (in this case, $D$), perform an LS or ML estimate of the parameters with a linear relationship with the observations (in this case, $\textbf {x}$), and iterate this approach across a range of values for $D$ to find the value leading to the best overall fit. For the case of additive white Gaussian noise, the ML estimate will always reduce to some variant of LS fit (e.g. standard LS or nonlinear LS), so we use the two terms synonymously here.

Assuming white Gaussian noise, the ML estimate is

(65)$$\left( {\hat {\textbf{x}},\hat D } \right) = \arg \mathop {\min }_{\textbf{{x}},D } \left\| {\textbf{{y}} - \textbf{{H}}\left( D \right)\textbf{{x}}} \right\|^{2} .$$

For a known (or tentatively assumed) value for $D$, Eq. (65) can be minimized over $\textbf {x}$ using standard multivariate calculus. Then Eq. (65) can be minimized over $D$ via a line search. This leads to

(66)$$\begin{aligned} \hat{D} =& \ \arg \mathop {\min }_{D} \left\| { \textbf{y} - \textbf{H}(D) \hat {\textbf{x}}(D)} \right\|^{2} \end{aligned}$$

(67)$$\begin{aligned} \hat {\textbf{x}}(D) =& \ {\left[ { \textbf{H}^\textrm{T}(D) \textbf{H}(D)} \right]^{ - 1}}{ \textbf{H}^\textrm{T}}(D) \textbf{y} . \end{aligned}$$

The choice for $\textbf {H}$ will be as in (51) if the target Mueller matrix is assumed to be diagonal, or (53) or (55) if the six elements in (54) or (56) are assumed to be nonzero, respectively.

5.2 Feature estimation error

In the remainder of this section, we will simulate the transmitted waveforms and estimate different subsets of Mueller matrix elements. The parameter settings are $\textrm {FWHM} = 10$ ns, a sampling rate and detector bandwidth of $1$ GHz, $V_{\pi } = 192$ V, $V_{\textrm {max}} = 1344$ V, $\beta = 0.1$ ns$^{-1}$, $\alpha _{EO} = 30^{o}$, $\theta _{EO} = 18.3^{o}$, wavelength $\lambda = 1.06~\mu$m, a fiber length of $L$=1 m, birefringence $\delta _n = 3.5 \cdot 10^{-4}$, nonlinear coefficient $n_2 = 1.2\cdot 10^{-22}$ m$^{2}$/V$^{2}$, $\alpha _{Kerr} = 26.9^{o}$, $\theta _{Kerr} = 17.6^{o}$, $P_0 = 5$ kW, a spot size of $w_i = 3.5 \cdot 10^{-6}$ m, and a core refractive index of $n_0 = 1.45$.

The Mueller matrix was from oak bark, measured at $1.0607\mu$m at normal incidence, and normalized relative to $M_{00}$:

(68)$$\textbf{M}_{\textrm{oakbark}} = \left[ { \begin{array}{rrrr} 1.0000 & 0.0079 & -0.0021 & 0.0024 \\ -0.0132 & 0.1657 & 0.0011 & -0.0008 \\ -0.0006 & 0.0000 & -0.1746 & -0.0025 \\ 0.0011 & 0.0010 & 0.0006 & -0.0230 \end{array} } \right]$$

For both systems, we assume $\textbf {M}$ is diagonal and only estimate the four elements along the diagonal. In order to remove the ambiguity caused by the unknown scale factor $Q$ in (52), the estimated elements of $\textbf {x}$ are divided by the estimate corresponding to $M_{0,0}$, resulting in three normalized estimates:

(69)$$\widehat{M}_{i,i} = \frac{\hat{x}_i}{\hat{x}_0} , \ \ \ i \in \left\{{1,2,3}\right\}.$$

The performance metric will be the Normalized Root Mean Squared Error (NRMSE). Given that the true Mueller matrix and the estimates have each already been normalized by $M_{0,0}$ and $\widehat {M}_{0,0}$, respectively, the NRMSE is

(70)$$\begin{aligned}\hbox{NRMSE of }M_{i,i} =& \ \sqrt{ \frac{1}{N_{tr}} \sum_{n_{tr}=1}^{N_{tr}} \left( {\widehat{M}_{i,i,n_{tr}} - M_{i,i}} \right)^{2} } , \ \ \ i \in \left\{{1,2,3}\right\} \end{aligned}$$

(71)$$\begin{aligned}\hbox{NRMSE (avg)} =& \ \sqrt{ \frac{1}{3} \sum_{i=1}^{3} \left( {\hbox{NRMSE of }M_{i,i}} \right)^{2} } , \end{aligned}$$

where $n_{tr} \in \left \{{1,\ldots ,N_{tr}}\right \}$ indexes the Monte Carlo trials that average over noise realizations.

Figure 7 shows the NRMSE averaged over $N_{tr} = 1000$ Monte Carlo trials, for an EO system [21] and the proposed Kerr-based system. There are two pairs of curves figures; in the two labeled “(diag)”, the off-diagonal elements of $\textbf {M}$ were set to zero when generating the waveforms; whereas in the other two curves, the full Mueller matrix of (68) was used, so the diagonal assumption is invalid. Either way, the estimator assumed the Mueller matrix was diagonal and only estimated the four diagonal elements. Figure 7 shows that, in this case, the error resulting from the assumption of a diagonal sample Mueller matrix has negligible effect on the measurement accuracy, and the estimation accuracy of a truly diagonal sample Mueller matrix has some marginal improvement as SNR increases. When the assumption of a diagonal Mueller matrix is met, performance improves monotonically as SNR increases; whereas the presence of residual off-diagonal elements causes modeling error that manifests as a “floor” on estimation performance and high SNR. The EOM and Kerr results are extremely similar, though the Kerr error in the $M_{1,1}$ estimate is about 46% higher than that of the EOM system (except in the presence of the error floor effect). However, this is offset by the fact that the EOM error is about 12% higher in the $M_{2,2}$ estimate, and even though that is a smaller percentage increase, the error in that estimate is higher so the error averaged across all three elements is very comparable between the EOM and Kerr systems. This comparison shows that the more compact and cost-effective Kerr-based system produces results comparable to those from the larger and more cumbersome EO system.

Fig. 7. Comparison of the NRMSE between the EO and Kerr systems, broken out by diagonal elements of the target Mueller matrix. In this case, the true Mueller matrix was not diagonal.

Download Full Size | PDF

Note that the error floor is higher in elements $M_{1,1}$ and $M_{2,2}$ than $M_{3,3}$. This is related to the analysis surrounding (64). The performance of the estimators for these two elements are coupled, so when one of them has increased error, so does the other.

6. Experimental demonstration

In this section, we demonstrate the Kerr PSG with an experimental setup similar to Fig. 1, but we implement the three PSAs serially, one pulse at a time. We employ a Keopsys fiber laser with wavelength = 1.064 $\mu$m, 10 ns FWHM pulse width and linear polarization. We used a Newport PM fiber with nominal parameters $\delta n = 3.0 \cdot 10^{-4}$ and mode field diameter 7.7 $\mu$m. The experimental demonstration is constrained by laser coherence length and the fiber damage threshold. The optical path difference for the orthogonal polarization states in the fiber is OPD$= L \cdot \delta n$ whereas the coherence length is $L_{\textrm {coh}} = (2 \ln 2 / \pi )^{1/2} \lambda ^{2} / (n \cdot \delta \lambda )$ for a Gaussian spectral distribution. For loss of coherence to be negligible, we need fiber OPD$\ll L_{\textrm {coh}}$, or $L \ll 18$ m for our laser with linewidth $\delta \lambda \sim 0.07$ nm. We chose $L = 5$ m and pump power 500 W to achieve a reasonable level of polarization modulation with limited loss of coherence and risk of damage in the fiber. From Section 2.2 we calculate $\epsilon \sim 6$ ps which is $\sim 0.1\%$ of the pulse width; this places the interaction in the “no walk-off” regime. Lastly, it should be pointed out that the experimental laser pulse shape was not $\textrm {sech}^{2}(t)$ but had a more complicated shape as shown by the grey line in the Fig. 8 $S_1$ plot. The solution to the coupled differential equations depends on the exact pulse shape employed. To get a first order estimate of the veracity of the model we replace the $\textrm {sech}^{2}(t)$ power distribution in Eq. (4) with the measured pump power envelope.

Fig. 8. Comparison of measured and modeled Stokes elements; the elements were normalized to the peak power of the incident pump pulse. In the $S_1$ plot, the gray line is pump envelope plotted on the negative $S_1$ axis for comparison to $S_1$. The red and blue lines denote data and model fits, respectively.

Download Full Size | PDF

A half-wave plate was used to orient the incident linearly polarized laser at $68^{o}$ to the fiber slow axis. Three PSAs were used to characterize the emission from the fiber; there were two linear PSAs oriented at $0^{o}$ and $45^{o}$ to measure the $S_1$ and $S_2$ Stokes elements, respectively, and a circular polarization analyzer was used to measure $S_3$. Two, 1 GHz bandwidth, Thorlabs pin detectors were used at each polarizer to record the orthogonal polarization signals and then construct the Stokes elements from the difference of measured signals. Figure 8 shows the recorded Stokes elements in comparison with the model predictions. The model was based on Eq. (19) and (20) and also allowed for some small angular misalignment of the fiber, wave plate, and polarizer. Reasonably good agreement between model and measurement is observed. The quantities $\sqrt {S_1^{2} + S_2^{2} + S_3^{2}}$ and $S_0$ should be identical if the pulse is fully coherent. Inspection of these quantities yields noticeable differences in pulse shape suggesting there is some loss of coherence due to the long fiber length. Coherence could be maintained by employing a fiber with end caps so that shorter fiber length and higher pump power could be employed. Future work will employ a pulse shape more consistent with the $\textrm {sech}^{2}(t)$ pump envelope and a short length of end-capped fiber to maintain coherence.

7. Conclusion

We have developed and analyzed a radically new, single pulse Mueller matrix Light Detection and Ranging (LiDAR) polarimeter based on a temporally multiplexed architecture and a PSG based on the Kerr effect in a birefringent optical fiber. We have derived, to our knowledge, the first solution describing the Kerr effect in a polarization maintaining fiber that includes walk-off effects. A temporally multiplexed LiDAR polarimeter performance model was constructed for a PSG based on the Kerr effect in PM fiber and a Pockels cell. The primary difference in the description of the two PSGs is the phase function; proper selection of fiber and pump parameters can enable both systems to produce similar temporal polarization distributions. Application of Cramer-Rao lower bound and Monte-Carlo simulations enabled identification of optimum measurement parameters and measurement accuracy dependence on signal to noise ratio. Simulations showed that samples with a diagonal Mueller matrix could be measured with high accuracy measurements in a single $\sim 10$ ns pulse at sufficient signal to noise ratio. Lastly, experimental demonstration of the Kerr PSG showed good agreement with model predictions. Future work will employ the single-pulse Kerr effect PSG in a polarimeter for Mueller matrix measurements.

Disclosures

The authors declare no conflicts of interest.

References

1. S. L. Jacques, J. C. Ramella-Roman, and K. Lee, “Imaging skin pathology with polarized light,” J. Biomed. Opt. 7(3), 329–340 (2002). [CrossRef]

2. J. T. Bruulsema, J. E. Hayward, T. Farrell, M. S. Patterson, L. Heinemann, M. Berger, T. Koschinsky, J. Sandahl-Christiansen, H. Orskov, M. Essenpreis, G. Schmelzeisen-Redeker, and D. B ocker, “Correlation between blood glucose concentration in diabetics and noninvasively measured tissue optical scattering coefficient,” Opt. Lett. 22(3), 190–192 (1997). [CrossRef]

3. T. W. H. Oates, H. Wormeester, and H. Arwin, “Characterization of plasmonic effects in thin films and metamaterials using spectroscopic ellipsometry,” Prog. Surf. Sci. 86(11-12), 328–376 (2011). [CrossRef]

4. I. S. Nerbø, S. LeRoy, M. Foldyna, M. Kildemo, and E. Sondergard, “Characterization of inclined GaSb nanopillars by Mueller matrix ellipsometry,” J. Appl. Phys. 108(1), 014307 (2010). [CrossRef]

5. D. H. Goldstein, L. L. Deibler, and B. B. Wang, “Measurement of small birefringence in sapphire and quartz plates,” Proc. SPIE 4819, 20 (2002). [CrossRef]

6. D. Jones, D. H. Goldstein, and J. C. Spaulding, “Reflective and polarimetric characteristics of urban materials,” in Proc. SPIE 6240, Polarization: Measurement, Analysis, and Remote Sensing VII, 62400A, (2006), pp. 1–9.

7. G. Anna, F. Goudail, and D. Dolfi, “Optimal discrimination of multiple regions with an active polarimetric imager,” Opt. Express 19(25), 25367–25378 (2011). [CrossRef]

8. F. Goudail and J. S. Tyo, “When is polarimetric imaging preferable to intensity imaging for target detection?” J. Opt. Soc. Am. A 28(1), 46–54 (2011). [CrossRef]

9. M. Alouini, F. Goudail, J. Bourderionnet, D. Dolfi, A. Bénière, I. Baarstad, T. Løke, P. Kaspersen, X. Normandin, and G. Berginc, “Near-infrared active polarimetric and multispectral laboratory demonstrator for target detection,” Appl. Opt. 48(8), 1610–1618 (2009). [CrossRef]

10. M. Del Guasta, E. Vallar, O. Riviere, F. Castagnoli, V. Venturi, and M. Morandi, “Use of polarimetric lidar for the study of oriented ice plates in clouds,” Appl. Opt. 45(20), 4878–4887 (2006). [CrossRef]

11. A. Shimizu, N. Sugimoto, I. Matsui, K. Arao, I. Uno, T. Murayama, N. Kagawa, K. Aoki, A. Uchiyama, and A. Yamazaki, “Continuous observations of Asian dust and other aerosols by polarization lidars in China and Japan during ACE-Asia,” J. Geophys. Res. 109(D19), D19S17 (2004). [CrossRef]

12. V. Noel, H. Chepfer, G. Ledanois, A. Delaval, and P. H. Flamant, “Classification of particle effective shape ratios in cirrus clouds based on the lidar depolarization ratio,” Appl. Opt. 41(21), 4245–4257 (2002). [CrossRef]

13. S. Tan and R. M. Narayanan, “A multiwavelength airborne polarimetric lidar for vegetation remote sensing: Instrumentation and preliminary test results,” in Proc. IEEE Int. Geosci. Remote Sens. Symp., vol. 5 (2002), pp. 2675–2677.

14. S. Tan and A. Haider, “A comparative study of polarimetric and nonpolarimetric lidar in deciduous-coniferous tree classification,” in Proc. IEEE Int. Geosci. Remote Sens. Symp., vol. 5 (2010), pp. 1178–1181.

15. R. M. A. Azzam, “Photopolarimetric measurement of the Muller matrix by Fourier analysis of a single detected signal,” Opt. Lett. 2(6), 148–150 (1978). [CrossRef]

16. D. H. Goldstein, “Mueller matrix dual-rotating retarder polarimeter,” Appl. Opt. 31(31), 6676–6684 (1992). [CrossRef]

17. R. M. A. Azzam, “Stokes-vector and Mueller-matrix polarimetry [Invited],” J. Opt. Soc. Am. A 33(7), 1396–1408 (2016). [CrossRef]

18. Y. Y. Markushin, N. P. Calvano, G. S. Pati, and R. Tripathi, “A polarimetric scanning LADAR: system development and performance analysis,” in Proc. SPIE, (2014), pp. 908004–1–908004–9.

19. G. Martínez-Ponce, C. Solano, and C. Pérez-Barrios, “Hybrid complete mueller polarimeter based on phase modulators,” Opt. Lasers Eng. 49(6), 723–728 (2011). [CrossRef]

20. A. De Martino, Y. K. Kim., E. Garcia-Caurel, B. Laude, and B. Drévillon, “Optimized mueller polarimeter with liquid crystals,” Opt. Lett. 28(8), 616–618 (2003). [CrossRef]

21. C. K. Keyser, R. K. Martin, P. K. Nguyen, and A. M. Adams, “Single-Pulse Mueller Matrix LiDAR Polarimeter: Modeling and Demonstration,” IEEE Trans. Geosci. Remote Sensing 57(6), 3296–3307 (2019). [CrossRef]

22. R. Azzam, “Division-of-amplitude Photopolarimeter (DOAP) for the Simultaneous Measurement of All Four Stokes Parameters of Light,” Opt. Acta 29(5), 685–689 (1982). [CrossRef]

23. R. C. Thompson, J. R. Bottiger, and E. S. Fry, “Measurement of polarized light interactions via the Mueller matrix,” Appl. Opt. 19(8), 1323–1332 (1980). [CrossRef]

24. D. H. Goldstein, Polarized Light (CRC Press, 2011), 3rd ed.

25. S. Breugnot and P. Clémenceau, “Modeling and performances of a polarization active imager at λ = 806nm,” Opt. Eng. 39(10), 2681–2688 (2000). [CrossRef]

26. H. C. Van de Hulst, Light Scattering by Small Particles (Dover, New York, 1981).

27. S. Clémenceau, P. Breugnot, and L. Collot, “Polarization diversity active imaging,” Proc. SPIE 3380, 284–291 (1998). [CrossRef]

28. J. Ellis, P. Caillard, and A. Dogariu, “Off-diagonal Mueller matrix elements in backscattering from highly diffusive media,” J. Opt. Soc. Am. A 19(1), 43–48 (2002). [CrossRef]

29. G. Agrawal, Nonlinear Fiber Optics (Academic Press, 2007), 5th ed.

30. A. Yariv and P. Yeh, Optical Waves in Crystals (John Wiley and Sons, 1984).

31. R. Chipman, “Chapter 22: Polarimetry,” in Handbook of Optics. Vol 2: Devices, Measurements, & Properties, 2ed, ed. M. Bass, McGraw Hill Education, 1995, URL: http://photonics.intec.ugent.be/education/IVPV/res_handbook/v2ch22.pdf.

32. R. K. Martin, C. Keyser, L. Ausley, and M. Steinke, “Pixel classification with a temporally multiplexed spectropolarimetric ladar system,” IEEE Trans. Geosci. Remote. Sens. 56(7), 3735–3746 (2018).

33. S. M. Kay, Fundamentals of Statistical Signal Processing: Estimation theory (Prentice Hall, 1993).

	$b_{1} = [\begin{matrix} 1 \\ 1 \end{matrix}]$	$b_{2} = [\begin{matrix} \cos ψ \\ \cos ψ \end{matrix}]$	$b_{3} = [\begin{matrix} \sin ψ \\ \sin ψ \end{matrix}]$	$b_{4} = [\begin{matrix} 1 \\ - 1 \end{matrix}]$	$b_{5} = [\begin{matrix} \cos ψ \\ - \cos ψ \end{matrix}]$	$b_{6} = [\begin{matrix} \sin ψ \\ - \sin ψ \end{matrix}]$
$H [:, 1] \leftrightarrow M_{00}$	$x$
$H [:, 2] \leftrightarrow M_{01}$	$x$	$x$
$H [:, 3] \leftrightarrow M_{02}$	$x$	$x$
$H [:, 4] \leftrightarrow M_{03}$			$x$
$H [:, 5] \leftrightarrow M_{10}$				x
$H [:, 6] \leftrightarrow M_{11}$				$x$	$x$
$H [:, 7] \leftrightarrow M_{12}$				x	x
$H [:, 8] \leftrightarrow M_{13}$						x
$H [:, 9] \leftrightarrow M_{20}$				x
$H [:, 10] \leftrightarrow M_{21}$				x	x
$H [:, 11] \leftrightarrow M_{22}$				$x$	$x$
$H [:, 12] \leftrightarrow M_{23}$						x
$H [:, 13] \leftrightarrow M_{30}$				x
$H [:, 14] \leftrightarrow M_{31}$				x	x
$H [:, 15] \leftrightarrow M_{32}$				x	x
$H [:, 16] \leftrightarrow M_{33}$						$x$

Single-pulse, Kerr-effect Mueller matrix LiDAR polarimeter

Abstract

1. Introduction

2. Nonlinear optical polarization state generator

2.1 System without walk-off effects

2.2 System with walk-off effects

2.3 Model of receiver hardware

2.4 Noise model

3. Linear algebra analysis

4. Cramer-Rao bound analysis

5. Performance simulations

5.1 Parameter estimation algorithm

5.2 Feature estimation error

6. Experimental demonstration

7. Conclusion

Disclosures

References

Cited By

Figures (8)

Tables (1)

Equations (71)

Optics Express