This paper develops a 3D vector solution for the scattering of partially coherent laser-beam illumination from statistically rough surfaces. Such a solution enables a rigorous comparison to the well-known Priest and Meier polarimetric bidirectional reflectance distribution function (pBRDF) [Opt. Eng. 41(5), 988 (2002) [CrossRef] ]. Overall, the comparison shows excellent agreement for the normalized spectral density and the degree of polarization. Based on this agreement, the 3D vector solution also enables an extension to the Priest and Meier pBRDF that accounts for the effects of active illumination. In particular, the 3D vector solution enables the development of a closed-form expression for the spectral degree of coherence. This expression provides a gauge for the average speckle size based on the spatial-coherence properties of the laser source. Such an extension is of broad interest to long-range applications that deal with speckle phenomena.
1. Introduction
The Priest and Meier polarimetric bidirectional reflectance distribution function (pBRDF) is an oft-used model for rough-surface scattering [1]. At its core, this model assumes the use of Gaussian distributed and Gaussian correlated micro facets, in addition to the geometrical-optics approximation (i.e., ray-optics propagation in the limit that the wavelength goes to zero). Such simplifications enable a straightforward model to the problem at hand, but lack the rigor associated with Maxwell’s equations.
In this paper, we develop a partially coherent 3D vector solution. This solution assumes the use of Gaussian distributed and Gaussian correlated surface heights, in addition to the physical-optics (PO) approximation (i.e., wave-optics propagation with Kirchhoff boundary conditions). Such a solution then enables a rigorous comparison to the Priest and Meier pBRDF in terms of the normalized spectral density (SD) and the degree of polarization (DoP). It also enables an extension that properly accounts for the effects of active illumination.
When using active illumination, more often than not, a partially coherent laser beam illuminates an object of interest. This topic plays a key role in long-range applications that use the light scattered from distant objects for tactical purposes. Interestingly enough, limited literature exists pertaining to the scattering of partially coherent laser-beam illumination from rough surfaces.
In an effort to bridge this gap, recent publications derived a 2D scalar-equivalent solution for the scattering of partially coherent laser-beam illumination from statistically rough surfaces using the PO approximation [2,3]. Specifically, the analysis made use of a Gaussian Schell-model (GSM) form in creating the incident field cross-spectral density function (CSDF). This choice allows one to vary the size and spatial coherence properties of the incident radiation. In so doing, the analysis formulated closed-form expressions for the scattered field CSDF to observe the size and spatial coherence properties of the scattered radiation in the far zone. The analysis also validated these analytical expressions through computational simulations and showed good agreement between the theoretical predictions and the numerical results.
While the 2D scalar-equivalent solution is a convenient tool for gaining insight into rough surface scattering, a complete understanding of the problem requires a 3D vector solution. With that said, this paper makes use of the PO approximation to determine a 3D vector solution for the far-field scattering of laser-beam illumination with partial spatial coherence from statistically rough surfaces. By formulating the 3D vector solution in a manner consistent with Wolf’s unified theory of coherence and polarization [4,5], all physical implications inherent in Wolf’s work apply here.
The 3D vector solution formulated in this paper considers three-different material substrates: dielectrics, conductors, and a perfect electrical conductor (PEC). In turn, it goes well beyond the results of a recent conference paper [6], which only studied the ideal case of a PEC. Additionally, the 3D vector solution uses a GSM form in creating the incident field cross-spectral density matrix (CSDM). This choice allows for the development of a closed-form expression for the scattered field CSDM that is applicable to very rough surfaces.
As formulated, the 3D vector solution contains complicated functions of the source parameters (size and coherence properties) and the surface parameters (surface height standard deviation and correlation length). Under certain circumstances, it maintains a GSM form. Based on these circumstances, we develop a closed-form expression for the spectral degree of coherence (SDoC). This expression also contains complicated functions of both the source and the surface parameters. Like with the 2D scalar-equivalent solution and a paraxial solution [7,8], we demonstrate that for long-range applications, we can approximate/simplify the SDoC radius as a function of just the source parameters.
In what follows, Section 2 provides the methodology used to obtain the 3D vector solution. Here, we state all simplifying assumptions and explain their physical implications. Section 3 then provides a comparison to the Priest and Meier pBRDF, which shows excellent agreement for the normalized SD and the DoP. Based on this agreement, Section 4 formulates a closed-form expression for the SDoC and provides a comparison to the 2D scalar-equivalent solution and a paraxial solution. Throughout Sections 2–4, we regularly refer to Supplement 1 for more details. The conclusion for this paper follows thereafter in Section 5.
2. Methodology for the 3D vector solution
Figure 1 describes the geometry used to obtain a 3D vector solution. As shown, a zero mean 2D sample function $h = h({x,y} )$ describes the surface height at the rough interface with standard deviation ${\sigma _h}$ and correlation length ${\ell _h}$. This function gives rise to a statistically rough surface S. Spatially partially coherent laser-beam illumination (parameters given below) emanates from the source-plane specified by the coordinates $({x,u,v} )$, which are different from the surface-plane coordinates $({x,y,z} )$. As such, the vector $\boldsymbol{\mathrm{\rho}}=x\hat{{\textbf x}} + u\hat{{\textbf u}}$ points from the source-plane origin to a transverse beam location since $v = 0$ in the source-plane; the vector ${{\textbf r}_s} = 0\hat{{\textbf x}} - {y_s}\hat{{\textbf y}} + {z_s}\hat{{\textbf z}}$ points from the source-plane origin to the surface-plane origin; and the vector ${\textbf r} = x\hat{{\textbf x}} + y\hat{{\textbf y}} + z\hat{{\textbf z}}$ points from the surface-plane origin to an observation point. Here, we assume that the rough surface is statistically isotropic [9,10] with impedance $\eta $. Therefore, and without loss of generality, we align the x axes of the source and surface planes. Above, the medium is free space with impedance ${\eta _0}$.
2.1 Incident field cross-spectral density matrix
With Fig. 1 in mind, spatially partially coherent laser-beam illumination emanates from the source plane. As such, we use a Gaussian Schell-model (GSM) form for the incident field cross-spectral density matrix (CSDM) $\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{W}}{^{i}}({{\boldsymbol{\mathrm{\rho}}_1},{\boldsymbol{\mathrm{\rho}}_2}})$ [5], such that
(1)$$\begin{aligned}\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{W}}{^{i}}\left(\boldsymbol{\mathrm{\rho}}_1, \boldsymbol{\mathrm{\rho}}_2\right) & \equiv\left\langle\mathbf{E}^i\left(\boldsymbol{\mathrm{\rho}}_1\right) \mathbf{E}^{i \dagger}\left(\boldsymbol{\mathrm{\rho}}_2\right)\right\rangle \\ & =\left[\begin{array}{ll} \left\langle E_x^i\left(\boldsymbol{\mathrm{\rho}}_1\right) E_x^{i^*}\left(\boldsymbol{\mathrm{\rho}}_2\right)\right\rangle & \left\langle E_x^i\left(\boldsymbol{\mathrm{\rho}}_1\right) E_u^{i^*}\left(\boldsymbol{\mathrm{\rho}}_2\right)\right\rangle \\ \left\langle E_u^i\left(\boldsymbol{\mathrm{\rho}}_1\right) E_x^{i^*}\left(\boldsymbol{\mathrm{\rho}}_2\right)\right\rangle & \left\langle E_u^i\left(\boldsymbol{\mathrm{\rho}}_1\right) E_u^{i^*}\left(\boldsymbol{\mathrm{\rho}}_2\right)\right\rangle \end{array}\right] \\ & =W_{\mathrm{mn}}^i\left(\boldsymbol{\mathrm{\rho}}_1, \boldsymbol{\mathrm{\rho}}_2\right) \quad(\mathrm{m}=x, u ; \mathrm{n}=x, u) \\ & =A_{\mathrm{m}} \exp \left(-\frac{\left|\boldsymbol{\mathrm{\rho}}_1\right|^2}{4 w_s^2}\right) A_{\mathrm{n}} \exp \left(-\frac{\left|\boldsymbol{\mathrm{\rho}}_2\right|^2}{4 w_s^2}\right) B_{\mathrm{mn}} \exp \left(-\frac{\left|\boldsymbol{\mathrm{\rho}}_2-\boldsymbol{\mathrm{\rho}}_1\right|^2}{2 \ell_{\mathrm{mn}}^2}\right),\end{aligned}$$
where $\left\langle \cdots \right\rangle $ denotes ensemble average, $\dagger $ denotes Hermitian conjugate, and ${\ast} $ denotes complex conjugate. In Eq. (1), the element-based parameters ${A_\textrm{m}}$ and ${A_\textrm{n}}$ are the beam amplitudes in the x and u directions, respectively, ${w_s}$ is the source width, and ${B_{\textrm{mn}}}$ and ${\ell _{\textrm{mn}}} = {\ell _{\textrm{nm}}}$ are the correlation amplitude and correlation length, respectively. Note that ${B_{\textrm{mn}}}$ follows additional constraints [5,11]; namely, (2)$$\begin{array}{{ccc}} {\begin{array}{*{20}{c}} {{B_{\textrm{mn}}} = 1}&{\textrm{when m} = \textrm{n}} \end{array}}\\ {\begin{array}{*{20}{c}} {|{{B_{\textrm{mn}}}} |\le 1}&{\textrm{when m} \ne \textrm{n}} \end{array}}\\ {{B_{\textrm{mn}}} = B_{\textrm{nm}}^\ast } \end{array}. $$
Also note that, in general, the incident field ${{\textbf E}^i}$ and the parameters ${A_\textrm{m}}$, ${w_s}$, ${\ell _{\textrm{mn}}}$, and ${B_{\textrm{mn}}}$ are radian frequency $\omega $ dependent [5]; however, we omit this dependence for brevity in the notation.2.2 Scattered field
In what follows, we use the PO approximation to develop a far-field expression for the scattered field ${{\textbf E}^s}({\textbf r} )$. For this purpose, we can write the incident field ${{\textbf E}^i}({\textbf r} )$ in terms of its spectrum ${{\textbf T}^i}({k_x^i,k_u^i} )$ using the plane-wave spectrum representation [12]. Using the macro-scale scattering geometry given in Fig. 1, the following expressions result:
(3)$${{\textbf E}^i}({\textbf r} )= \frac{1}{{{{({2\pi } )}^2}}}\int\limits_{ - \infty }^\infty {\int\limits_{ - \infty }^\infty {{{\textbf T}^i}({k_x^i,k_u^i} ){\textrm{e}^{ - j{{\textbf k}^i} \cdot {{\textbf r}_s}}}{e ^{ - j{{\textbf k}^i} \cdot {\textbf r}}}dk_x^idk_u^i} }$$
and (4)$${{\textbf T}^i}({k_x^i,k_u^i} )= \int\limits_{ - \infty }^\infty {\int\limits_{ - \infty }^\infty {{{\textbf E}^i}(\boldsymbol{\mathrm{\rho}} ){\textrm{e}^{jk_x^ix}}{\textrm{e}^{jk_u^iu}}dxdu} }, $$
which are valid in the source-free half space where $v \ge 0$. In Eqs. (3) and (4), ${{\textbf k}^i} = {k_0}{\hat{{\textbf k}}^i} = k_x^i\hat{{\textbf x}} + k_u^i\hat{{\textbf u}} + k_v^i\hat{{\textbf v}}$ is the incident propagation vector, ${k_0} = {{2\pi } / {{\lambda _0}}}$ is the free-space wavenumber, and ${\lambda _0}$ is the free-space wavelength.Consistent with long-range applications, we assume that the observation points of interest are in the far zone. As such, the scattered electric field ${{\textbf E}^s}({\textbf r} )$ depends on the far-field vector potentials, ${\textbf L}({\textbf r} )$ and ${\textbf N}({\textbf r} )$, using the following relationships [13]:
(5)$${{\textbf E}^s}\left( {\textbf r} \right) \approx j{k_0}\frac{{{e ^{ - j{k_0}r}}}}{{4\pi r}}\left[ {\left(\boldsymbol{\hat{\mathrm{\phi}}}{\boldsymbol{\hat{\mathrm{\theta}}} - {\boldsymbol{\hat{\mathrm{\theta}}}\boldsymbol{\hat{\mathrm{\phi}}}}} \right) \cdot {\textbf L}\left( {\textbf r} \right) - {\eta _0}}\left( \boldsymbol{\hat{\mathrm{\theta}}}{\boldsymbol{\hat{\mathrm{\theta}}}} + {\boldsymbol{\hat{\mathrm{\phi}}}\boldsymbol{\hat{\mathrm{\phi}}}} \right) \cdot {\textbf N}\left( {\textbf r} \right) \right],$$
(6)$${\textbf L}({\textbf r} )= \int\!\!\!\int\limits_S {{\textbf M}({{\textbf r^{\prime}}} ){e ^{j{k_0}\hat{{\textbf r}} \cdot {\textbf r^{\prime}}}}ds^{\prime}}, $$
and (7)$${\textbf N}({\textbf r} )= \int\!\!\!\int\limits_S {{\textbf J}({{\textbf r^{\prime}}} ){e ^{j{k_0}\hat{{\textbf r}} \cdot {\textbf r^{\prime}}}}ds^{\prime}}. $$
In Eqs. (5)–(7), $\boldsymbol{\hat{\mathrm{\theta}}}$ and $\boldsymbol{\hat{\mathrm{\phi}}}$ are unit vectors in the polar (vertical polarization) and azimuth (horizontal polarization) directions, respectively, ${\textbf M}({{\textbf r^{\prime}}} )$ and ${\textbf J}({{\textbf r^{\prime}}} )$ are the equivalent surface-current densities, respectively, and the vector ${\textbf r^{\prime}} = x^{\prime}\hat{{\textbf x}} + y^{\prime}\hat{{\textbf y}} + h^{\prime}\hat{{\textbf z}}$ points from the surface-plane origin to a point on the statistically rough surface S. Note that in the far zone, $r > {{2{D^2}} / {{\lambda _0}}}$, where $D \approx 2L$ and $r \gg L$. Consequently, we neglect all contributions to the scattered field ${{\textbf E}^s}$ that are in the radial $\hat{{\textbf r}}$ direction because their contributions scale as ${1 / {{r^n}}}$, where $n = 2,3, \ldots $, and are negligible [13]. Also note that in the case of rough surfaces, one can relate the far zone to the correlation length of the surface, not its actual size. This outcome explains why lasers reflected from rough surfaces almost immediately result in speckle.Using the micro-scale scattering geometry given in Fig. S1 of Supplement 1 and the PO approximation [13,14],
(8)$$\mathbf{M}\left(\mathbf{r}^{\prime}\right) \approx-\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{M}}{'} \cdot \mathbf{E}^i\left(\mathbf{r}^{\prime}\right)$$
and (9)$$\mathbf{J}\left(\mathbf{r}^{\prime}\right) \approx \frac{1}{\eta_0} \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{J}}{'} \cdot \mathbf{E}^i\left(\mathbf{r}^{\prime}\right),$$
where $\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{J}}{'}$ and $\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{M}}{'}$ are dyadics, such that (10)$$\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{J}}{'} = {\hat{\textbf{n}}^{\prime}} \times [{{\hat{\textbf{s}}^{\prime}}({1 - {{r^{\prime}}_\parallel }} ){\hat{\textbf{p}}^{\prime}} - {\hat{\textbf{p}}^{\prime}}({1 - {r_ \bot }^\prime } ){\hat{\textbf{s}}^{\prime}}} ]$$
and (11)$$\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{M}}{'} = {\hat{\textbf{n}}^{\prime}} \times [{{\hat{\textbf{s}}^{\prime}}({1 + {{r^{\prime}}_ \bot }} ){\hat{\textbf{s}}^{\prime}} + {\hat{\textbf{p}}^{\prime}}({1 + {{r^{\prime}}_\parallel }} ){\hat{\textbf{p}}^{\prime}}} ]. $$
In Eqs. (10) and (11), ${\hat{\textbf{n}}^{\prime}} = \hat{{\textbf n}}({x^{\prime},y^{\prime}} )$ is the 2D unit outward normal vector given by (12)$${\hat{\textbf{n}}^{\prime}} = \frac{{{\textbf n^{\prime}}}}{{|{{\textbf n^{\prime}}} |}} = \frac{{ - {h_{x^{\prime}}}\hat{{\textbf x}} - {h_{y^{\prime}}}\hat{{\textbf y}} + \hat{{\textbf z}}}}{{\sqrt {1 + h_{x^{\prime}}^2 + h_{y^{\prime}}^2} }}, $$
where (13)$$h_{x^{\prime}}=\frac{\partial h^{\prime}}{\partial x^{\prime}}=\frac{\partial h\left(x^{\prime}, y^{\prime}\right)}{\partial x^{\prime}} \quad h_{y^{\prime}}=\frac{\partial h^{\prime}}{\partial y^{\prime}}=\frac{\partial h\left(x^{\prime}, y^{\prime}\right)}{\partial y^{\prime}}.$$
Furthermore, ${\hat{\textbf{s}}^{\prime}}$ and ${\hat{\textbf{p}}^{\prime}}$ are the unit perpendicular and parallel vectors, whereas ${r_ \bot }^\prime $ and ${r_\parallel }^\prime $ are the corresponding Fresnel reflection coefficients, respectively. Referencing the micro-scale scattering geometry in Fig. S1 of Supplement 1, the following relationships result: (14)$$\hat{\mathbf{s}}=\frac{\hat{\mathbf{k}}^i \times \hat{\mathbf{n}}}{\left|\hat{\mathbf{k}}^i \times \hat{\mathbf{n}}\right|} \quad \hat{\mathbf{p}}=\hat{\mathbf{s}} \times \hat{\mathbf{k}}^i \quad \hat{\mathbf{p}}^r=-\hat{\mathbf{s}} \times \hat{\mathbf{k}}^r.$$
Thus, in arriving at the expressions found in Eqs. (10) and (11), we use the geometrical-optics (GO) approximation [14]; specifically, the law of reflection, such that $-\hat{\mathbf{n}} \cdot \hat{\mathbf{k}}^i=\hat{\mathbf{n}} \cdot \hat{\mathbf{k}}^r$. Please see Section S1 of Supplement 1 for more details on the GO and PO approximations.Based on Eq. (12), it is important to note that the integration in Eqs. (6) and (7) is over the parameterized rough surface (i.e., $ds^{\prime} = |{{\textbf n^{\prime}}} |dx^{\prime}dy^{\prime}$). Consequently, using Eq. (3) and substituting Eqs. (8)–(12) into Eqs. (6) and (7), we obtain the following expressions:
(15)$$\mathbf{L}(\mathbf{r})=-\frac{1}{(2 \pi)^2} \int_{-L}^L \int_{-L}^L \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{L}}{'} \cdot \mathbf{T}^i\left(k_x^i, k_u^i\right) \mathrm{e}^{-j \mathbf{k}^i \cdot \mathbf{r}_s} \mathrm{e}^{j \mathbf{q} \cdot \mathbf{r}^{\prime}} d k_x^i d k_u^i d x^{\prime} d y^{\prime}$$
and (16)$$\mathbf{N}(\mathbf{r})=\frac{1}{(2 \pi)^2 \eta_0} \int_{-L}^L \int_{-L}^L \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{N}}{'} \cdot \mathbf{T}^i\left(k_x^i, k_u^i\right) \mathrm{e}^{-j \mathbf{k}^i \cdot \mathbf{r}_s} \mathrm{e}^{j \mathbf{q} \cdot \mathbf{r}^{\prime}} d k_x^i d k_u^i d x^{\prime} d y^{\prime},$$
where $\mathbf{q}=k_0\left(\hat{\mathbf{r}}-\hat{\mathbf{k}}^i\right)=q_x \hat{\mathbf{x}}+q_y \hat{\mathbf{y}}+q_z \hat{\mathbf{z}}, \quad \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{L}}{'}=\left|\mathbf{n}^{\prime}\right| \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{M}}{'}, \quad {\textrm {and}} \quad \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{N}}{'}=\left|\mathbf{n}^{\prime}\right| \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{J}}{'}$. Without further simplifications, no closed-form expression exists for the far-field vector potentials, ${\textbf L}({\textbf r} )$ and ${\textbf N}({\textbf r} )$. This outcome is because the integrands in Eqs. (15) and (16) are complicated functions of surface height and surface slopes; namely, $h^{\prime}$, ${h_{x^{\prime}}}$, and ${h_{y^{\prime}}}$ with respect to the integrals over the parameterized rough surface. We simplify these integrals using the stationary-phase (SP) approximation [14,15], where (17)$$\begin{array}{{ccc}} {\frac{\partial }{{\partial x^{\prime}}}({{\textbf q} \cdot {\textbf r^{\prime}}} )\approx 0}&{\frac{\partial }{{\partial y^{\prime}}}({{\textbf q} \cdot {\textbf r^{\prime}}} )\approx 0} \end{array}. $$
As a result, the relationships found in Eq. (13) simplify, such that (18)$$\begin{array}{{ccc}} {{h_{x^{\prime}}} \approx{-} \frac{{{q_x}}}{{{q_z}}}}&{{h_{y^{\prime}}} \approx{-} \frac{{{q_y}}}{{{q_z}}}} \end{array}, $$
and in turn, $\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{L}}{'} \approx \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{L}}$ and $\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{N}}{'} \approx \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{N}}$ in Eqs. (15) and (16). Similar to the PO approximation, the SP approximation physically dictates that reflection from the rough surface is locally specular and excludes all local diffraction effects [14,15]. Note that in using Eqs. (6)–(18), we assume that all observation points are in the far zone, the PO approximation holds, and the effects of shadowing/masking and multiple scattering are negligible [9,16].2.3 Scattered field cross-spectral density matrix
Here, we develop closed-form expressions for the scattered field CSDM $\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{W}}{^{s}}({{{\textbf r}_1},{{\textbf r}_2}})$. In general, $\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{W}}{^{s}}({{{\textbf r}_1},{{\textbf r}_2}})$ depends on the scattered field ${{\textbf E}^s}({\textbf r} )$, such that in the far zone
(19)$$\begin{aligned} \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{W}}{^{s}}\left(\mathbf{r}_1, \mathbf{r}_2\right) & \equiv\left\langle\mathbf{E}^s\left(\mathbf{r}_1\right) \mathbf{E}^{s \dagger}\left(\mathbf{r}_2\right)\right\rangle \\ & =\left[\begin{array}{ll} \left\langle E_\theta^s\left(\mathbf{r}_1\right) E_\theta^{s^*}\left(\mathbf{r}_2\right)\right\rangle & \left\langle E_\theta^s\left(\mathbf{r}_1\right) E_\phi^{s^*}\left(\mathbf{r}_2\right)\right\rangle \\ \left\langle E_\phi^s\left(\mathbf{r}_1\right) E_\theta^{s^*}\left(\mathbf{r}_2\right)\right\rangle & \left\langle E_\phi^s\left(\mathbf{r}_1\right) E_\phi^{s^*}\left(\mathbf{r}_2\right)\right\rangle \end{array}\right] . \\ & =W_{\mathrm{pq}}^s\left(\mathbf{r}_1, \mathbf{r}_2\right) \quad(\mathrm{p}=\theta, \phi ; \mathrm{q}=\theta, \phi) \end{aligned}$$
Using Eq. (5), we determine the matrix elements found in Eq. (19) as (20)$$\begin{aligned} W_{\theta \theta }^s({{{\textbf r}_1},{{\textbf r}_2}} )&= {\Omega _{12}}\sum\limits_{\textrm{i} = x,y,z} {\sum\limits_{\textrm{j} = x,y,z} {\left[ {({{{\boldsymbol{\hat{\mathrm{\phi}}}}_1} \cdot \hat{{\textbf i}}} )({{{\boldsymbol{\hat{\mathrm{\phi}}}}_2} \cdot \hat{{\textbf j}}} )\left\langle {{L_\textrm{i}}({{{\textbf r}_1}} )L_\textrm{j}^\ast ({{{\textbf r}_2}} )} \right\rangle } \right.} } \\ &\quad + {\eta _0}({{{\boldsymbol{\hat{\mathrm{\phi}}}}_1} \cdot \hat{{\textbf i}}} )({{{\boldsymbol{\hat{\mathrm{\theta}}}}_2} \cdot \hat{{\textbf j}}} )\left\langle {{L_\textrm{i}}({{{\textbf r}_1}} )N_\textrm{j}^\ast ({{{\textbf r}_2}} )} \right\rangle \\ &\quad + {\eta _0}({{{\boldsymbol{\hat{\mathrm{\theta}}}}_1} \cdot \hat{{\textbf i}}} )({{{\boldsymbol{\hat{\mathrm{\phi}}}}_2} \cdot \hat{{\textbf j}}} )\left\langle {{N_\textrm{i}}({{{\textbf r}_1}} )L_\textrm{j}^\ast ({{{\textbf r}_2}} )} \right\rangle \\ &\quad\left. {\textrm{ } + \eta_0^2({{{\boldsymbol{\hat{\mathrm{\theta}}}}_1} \cdot \hat{{\textbf i}}} )({{{\boldsymbol{\hat{\mathrm{\theta}}}}_2} \cdot \hat{{\textbf j}}} )\left\langle {{N_\textrm{i}}({{{\textbf r}_1}} )N_\textrm{j}^\ast ({{{\textbf r}_2}} )} \right\rangle } \right], \end{aligned}$$
(21)$$\begin{aligned} W_{\theta \phi }^s({{{\textbf r}_1},{{\textbf r}_2}} )&= {\Omega _{12}}\sum\limits_{\textrm{i} = x,y,z} {\sum\limits_{\textrm{j} = x,y,z} {\left[ { - ({{{\boldsymbol{\hat{\mathrm{\phi}}}}_1} \cdot \hat{{\textbf i}}} )({{{\boldsymbol{\hat{\mathrm{\theta}}}}_2} \cdot \hat{{\textbf j}}} )\left\langle {{L_\textrm{i}}({{{\textbf r}_1}} )L_\textrm{j}^\ast ({{{\textbf r}_2}} )} \right\rangle } \right.} } \\ &\quad + {\eta _0}({{{\boldsymbol{\hat{\mathrm{\phi}}}}_1} \cdot \hat{{\textbf i}}} )({{{\boldsymbol{\hat{\mathrm{\phi}}}}_2} \cdot \hat{{\textbf j}}} )\left\langle {{L_\textrm{i}}({{{\textbf r}_1}} )N_\textrm{j}^\ast ({{{\textbf r}_2}} )} \right\rangle \\ &\quad - {\eta _0}({{{\boldsymbol{\hat{\mathrm{\theta}}}}_1} \cdot \hat{{\textbf i}}} )({{{\boldsymbol{\hat{\mathrm{\theta}}}}_2} \cdot \hat{{\textbf j}}} )\left\langle {{N_\textrm{i}}({{{\textbf r}_1}} )L_\textrm{j}^\ast ({{{\textbf r}_2}} )} \right\rangle \\ &\quad\left. {\textrm{ } + \eta_0^2({{{\boldsymbol{\hat{\mathrm{\theta}}}}_1} \cdot \hat{{\textbf i}}} )({{{\boldsymbol{\hat{\mathrm{\phi}}}}_2} \cdot \hat{{\textbf j}}} )\left\langle {{N_\textrm{i}}({{{\textbf r}_1}} )N_\textrm{j}^\ast ({{{\textbf r}_2}} )} \right\rangle } \right], \end{aligned}$$
(22)$$\begin{aligned} W_{\phi \theta }^s({{{\textbf r}_1},{{\textbf r}_2}} )&= {\Omega _{12}}\sum\limits_{\textrm{i} = x,y,z} {\sum\limits_{\textrm{j} = x,y,z} {\left[ { - ({{{\boldsymbol{\hat{\mathrm{\theta}}}}_1} \cdot \hat{{\textbf i}}} )({{{\boldsymbol{\hat{\mathrm{\phi}}}}_2} \cdot \hat{{\textbf j}}} )\left\langle {{L_\textrm{i}}({{{\textbf r}_1}} )L_\textrm{j}^\ast ({{{\textbf r}_2}} )} \right\rangle } \right.} } \\ &\quad - {\eta _0}({{{\boldsymbol{\hat{\mathrm{\theta}}}}_1} \cdot \hat{{\textbf i}}} )({{{\boldsymbol{\hat{\mathrm{\theta}}}}_2} \cdot \hat{{\textbf j}}} )\left\langle {{L_\textrm{i}}({{{\textbf r}_1}} )N_\textrm{j}^\ast ({{{\textbf r}_2}} )} \right\rangle \\ &\quad + {\eta _0}({{{\boldsymbol{\hat{\mathrm{\phi}}}}_1} \cdot \hat{{\textbf i}}} )({{{\boldsymbol{\hat{\mathrm{\phi}}}}_2} \cdot \hat{{\textbf j}}} )\left\langle {{N_\textrm{i}}({{{\textbf r}_1}} )L_\textrm{j}^\ast ({{{\textbf r}_2}} )} \right\rangle \\ &\quad\left. {\textrm{ } + \eta_0^2({{{\boldsymbol{\hat{\mathrm{\phi}}}}_1} \cdot \hat{{\textbf i}}} )({{{\boldsymbol{\hat{\mathrm{\theta}}}}_2} \cdot \hat{{\textbf j}}} )\left\langle {{N_\textrm{i}}({{{\textbf r}_1}} )N_\textrm{j}^\ast ({{{\textbf r}_2}} )} \right\rangle } \right], \end{aligned}$$
and (23)$$\begin{aligned} W_{\phi \phi }^s({{{\textbf r}_1},{{\textbf r}_2}} )&= {\Omega _{12}}\sum\limits_{\textrm{i} = x,y,z} {\sum\limits_{\textrm{j} = x,y,z} {\left[ {({{{\boldsymbol{\hat{\mathrm{\theta}}}}_1} \cdot \hat{{\textbf i}}} )({{{\boldsymbol{\hat{\mathrm{\theta}}}}_2} \cdot \hat{{\textbf j}}} )\left\langle {{L_\textrm{i}}({{{\textbf r}_1}} )L_\textrm{j}^\ast ({{{\textbf r}_2}} )} \right\rangle } \right.} } \\&\quad - {\eta _0}({{{\boldsymbol{\hat{\mathrm{\theta}}}}_1} \cdot \hat{{\textbf i}}} )({{{\boldsymbol{\hat{\mathrm{\phi}}}}_2} \cdot \hat{{\textbf j}}} )\left\langle {{L_\textrm{i}}({{{\textbf r}_1}} )N_\textrm{j}^\ast ({{{\textbf r}_2}} )} \right\rangle \\&\quad - {\eta _0}({{{\boldsymbol{\hat{\mathrm{\phi}}}}_1} \cdot \hat{{\textbf i}}} )({{{\boldsymbol{\hat{\mathrm{\theta}}}}_2} \cdot \hat{{\textbf j}}} )\left\langle {{N_\textrm{i}}({{{\textbf r}_1}} )L_\textrm{j}^\ast ({{{\textbf r}_2}} )} \right\rangle \\&\quad\left. {\textrm{ } + \eta_0^2({{{\boldsymbol{\hat{\mathrm{\phi}}}}_1} \cdot \hat{{\textbf i}}} )({{{\boldsymbol{\hat{\mathrm{\phi}}}}_2} \cdot \hat{{\textbf j}}} )\left\langle {{N_\textrm{i}}({{{\textbf r}_1}} )N_\textrm{j}^\ast ({{{\textbf r}_2}} )} \right\rangle } \right], \end{aligned}$$
where (24)$${\Omega _{12}} = k_0^2\frac{{{e ^{ - j{k_0}{r_1}}}{e ^{j{k_0}{r_2}}}}}{{{{({4\pi } )}^2}{r_1}{r_2}}}. $$
In addition, using Eqs. (6)–(18), we determine the element-based correlations found in Eqs. (19)–(23) from the following relationships:
(25)$$\begin{aligned} &\left\langle\mathbf{L}\left(\mathbf{r}_1\right) \mathbf{L}^{\dagger}\left(\mathbf{r}_2\right)\right\rangle= \frac{1}{(2 \pi)^4} \\ & \int_{-L}^L \int_{-L}^L \int_{-L}^L \int_{-L}^L \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{L}}_1 \cdot\left\langle\mathbf{T}^i\left(k_{x 1}^i, k_{u 1}^i\right) \mathbf{T}^{i \dagger}\left(k_{x 2}^i, k_{u 2}^i\right)\right\rangle \cdot \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{L}}{}_{2}^{\dagger} \\ & \mathrm{e}^{-j k_{v 1}^i r_s} \mathrm{e}^{j k_{v 2}^i r_s} \mathrm{e}^{j q_{x 1} x_1^{\prime}} \mathrm{e}^{-j q_{x 2} x_2^{\prime}} \mathrm{e}^{j q_{y 1} y_1^{\prime}} \mathrm{e}^{-j q_{y 2} y_2^{\prime}}\left\langle\mathrm{e}^{j q_{z 1} h_1^{\prime}} \mathrm{e}^{-j q_{z 2} h_2^{\prime}}\right\rangle \\ & d k_{x 1}^i d k_{x 2}^i d k_{u 1}^i d k_{u 2}^i d x_1^{\prime} d x_2^{\prime} d y_1^{\prime} d y_2^{\prime}, \end{aligned}$$
(26)$$\begin{aligned} &\left\langle\mathbf{L}\left(\mathbf{r}_1\right) \mathbf{N}^{\dagger}\left(\mathbf{r}_2\right)\right\rangle=-\frac{1}{(2 \pi)^4 \eta_0} \\ & \int_{-L}^L \int_{-L}^L \int_{-L}^L \int_{-L}^L \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{L}}_1 \cdot\left\langle\mathbf{T}^i\left(k_{x 1}^i, k_{u 1}^i\right) \mathbf{T}^{i \dagger}\left(k_{x 2}^i, k_{u 2}^i\right)\right\rangle \cdot \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{N}}{}_{2}^{\dagger} \\ & \mathrm{e}^{-j k_{v 1}^i r_s} \mathrm{e}^{j k_{v 2}^i r_s} \mathrm{e}^{j q_{x 1} x_1^{\prime}} \mathrm{e}^{-j q_{x 2} x_2^{\prime}} \mathrm{e}^{j q_{y 1} y_1^{\prime}} \mathrm{e}^{-j q_{y_2} y_2^{\prime}}\left\langle\mathrm{e}^{j q_{z 1} h_1^{\prime}} \mathrm{e}^{-j q_{z 2} h_2^{\prime}}\right\rangle \\ & d k_{x 1}^i d k_{x 2}^i d k_{u 1}^i d k_{u 2}^i d x_1^{\prime} d x_2^{\prime} d y_1^{\prime} d y_2^{\prime}, \end{aligned}$$
(27)$$\begin{aligned} &\left\langle\mathbf{N}\left(\mathbf{r}_1\right) \mathbf{L}^{\dagger}\left(\mathbf{r}_2\right)\right\rangle=-\frac{1}{(2 \pi)^4 \eta_0} \\ & \int_{-L}^L \int_{-L}^L \int_{-L}^L \int_{-L}^L \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{N}}_1 \cdot\left\langle\mathbf{T}^i\left(k_{x 1}^i, k_{u 1}^i\right) \mathbf{T}^{i \dagger}\left(k_{x 2}^i, k_{u 2}^i\right)\right\rangle \cdot \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{L}}{}_{2}^{\dagger} \\ & \mathrm{e}^{-j k_{v 1}^i r_s} \mathrm{e}^{j k_{v 2}^i r_s} \mathrm{e}^{j q_{x 1} x_1^{\prime}} \mathrm{e}^{-j q_{x 2} x_2^{\prime}} \mathrm{e}^{j q_{y 1} y_1^{\prime}} \mathrm{e}^{-j q_{y 2} y_2^{\prime}}\left\langle\mathrm{e}^{j q_{z 1} h_1^{\prime}} \mathrm{e}^{-j q_{z 2} h_2^{\prime}}\right\rangle \\ & d k_{x 1}^i d k_{x 2}^i d k_{u 1}^i d k_{u 2}^i d x_1^{\prime} d x_2^{\prime} d y_1^{\prime} d y_2^{\prime}, \end{aligned}$$
and (28)$$\begin{aligned} &\left\langle\mathbf{N}\left(\mathbf{r}_1\right) \mathbf{N}^{\dagger}\left(\mathbf{r}_2\right)\right\rangle=\frac{1}{(2 \pi)^4 \eta_0^2} \\ & \int_{-L}^L \int_{-L}^L \int_{-L}^L \int_{-L}^L \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{N}}_1 \cdot\left\langle\mathbf{T}^i\left(k_{x 1}^i, k_{u 1}^i\right) \mathbf{T}^{i \dagger}\left(k_{x 2}^i, k_{u 2}^i\right)\right\rangle \cdot \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{N}}{}_{2}^{\dagger} \\ & \mathrm{e}^{-j k_{v 1}^i r_s} \mathrm{e}^{j k_{v 2}^i r_s} \mathrm{e}^{j q_{x 1} x_1^{\prime}} \mathrm{e}^{-j q_{x 2} x_2^{\prime}} \mathrm{e}^{j q_{y 1} y_1^{\prime}} \mathrm{e}^{-j q_{y 2} y_2^{\prime}}\left\langle\mathrm{e}^{j q_{z 1} h_1^{\prime}} \mathrm{e}^{-j q_{z 2} h_2^{\prime}}\right\rangle \\ & d k_{x 1}^i d k_{x 2}^i d k_{u 1}^i d k_{u 2}^i d x_1^{\prime} d x_2^{\prime} d y_1^{\prime} d y_2^{\prime}, \end{aligned}$$
where ${r_s} = |{{{\textbf r}_s}} |$. Inherent in Eqs. (25)–(28) is the assumption that the incident field plane-wave spectrum is statistically independent of the rough surface. This assumption is physically intuitive; thus, Eqs. (25)–(28) contain two separate correlations.The first correlation is with respect to the incident field plane-wave spectrum. This correlation is equivalent to a dyadic [cf. Equation (4)]; namely,
(29)$$\left\langle\mathbf{T}^i\left(k_{x 1}^i, k_{u 1}^i\right) \mathbf{T}^{i \dagger}\left(k_{x 2}^i, k_{u 2}^i\right)\right\rangle=\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{\Phi}}{^{i}},$$
where $\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{\Phi}}{^{i}}=\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{\Phi}}{^{i}}\left(k_{x 1}^i, k_{x 2}^i, k_{u 1}^i, k_{u 2}^i\right)$. The second correlation is with respect to the parameterized rough surface. This correlation is a joint-characteristic function $\chi ^{\prime}$ of the random variables $h_1^{\prime}=h\left(x_1^{\prime}, y_1^{\prime}\right)\; {\textrm{and}}\; h_2^{\prime}=h\left(x_2^{\prime}, y_2^{\prime}\right)$, such that (30)$$\left\langle\mathrm{e}^{j q_{z 1} h_1^{\prime}} \mathrm{e}^{-j q_{z 2} h_2^{\prime}}\right\rangle=\chi^{\prime}.$$
In practice, we must choose a form for this joint-characteristic function.A very common choice for the statistical distribution of the rough surface is to assume that the surface heights are Gaussian distributed and Gaussian correlated. In so doing, the joint probability density function $p^{\prime}=p\left(h_1^{\prime}, h_2^{\prime}\right)$ of the random variables ${h^{\prime}_1}$ and ${h^{\prime}_2}$ takes the following form [17]:
(31)$$p^{\prime}=\frac{1}{2 \pi \sigma_h^2 \sqrt{1-\Gamma^{\prime 2}}} \exp \left[-\frac{h_1^{\prime 2}+h_2^{\prime 2}-2 \Gamma^{\prime} h_1^{\prime} h_2^{\prime}}{2 \sigma_h^2\left(1-\Gamma^{\prime 2}\right)}\right],$$
where $\Gamma^{\prime}=\Gamma\left(x_1^{\prime}-x_2^{\prime}, y_1^{\prime}-y_2^{\prime}\right)$ is the surface autocorrelation function, such that (32)$$\Gamma^{\prime}=\exp \left[-\frac{\left(x_1^{\prime}-x_2^{\prime}\right)^2}{\ell_h^2}\right] \exp \left[-\frac{\left(y_1^{\prime}-y_2^{\prime}\right)^2}{\ell_h^2}\right].$$
History shows that one typically chooses Gaussian-Gaussian (G-G) models for analytical convenience [18]; however, other models exist in practice. For example, the stretched exponential-stretched exponential (SE-SE) model better characterizes surfaces roughened by random industrial processes [19]. Basu et al. highlighted this point with profilometer measurements of sandblasted metallic surfaces [20,21]. Unfortunately no general analytical form exists for the SE joint-characteristic function; nevertheless, the analysis of Basu et al. also showed that G-G models were still fairly good approximations for sandblasted metallic surfaces [20,21]. Thus, Fourier transforming the joint probability density function $p^{\prime}$ in Eq. (31) yields the desired joint-characteristic function $\chi ^{\prime}$ [17], such that (33)$$\begin{aligned} \chi^{\prime} & =\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} p^{\prime} \mathrm{e}^{j q_{z 1} h_1^{\prime}} \mathrm{e}^{-j q_{z 2} h_2^{\prime}} d h_1^{\prime} d h_2^{\prime} \\ & =\exp \left[-\frac{\sigma_h^2}{2}\left(q_{z 1}^2+q_{z 2}^2\right)\right] \exp \left(\sigma_h^2 q_{z 1} q_{z 2} \Gamma^{\prime}\right),\end{aligned}$$
where $\chi^{\prime}=\chi\left(k_{x 1}^i, k_{x 2}^i, k_{u 1}^i, k_{u 2}^i ; x_1^{\prime}-x_2^{\prime}, y_1^{\prime}-y_2^{\prime}\right)$. Note that throughout the literature, numerous other surface models exist in addition to G-G and SE-SE models; for example, a recent publication explored the use of non-Gaussian surface autocorrelation functions [22].Using the relationships found in Eqs. (29)–(33), the integrands in Eqs. (25)–(28) still contain complicated functions with respect to the source and surface parameters. To simplify things, we can separate these complicated functions into amplitude and phase terms, viz.
(34)$$\begin{aligned} \left\langle\mathbf{L}\left(\mathbf{r}_1\right) \mathbf{L}^{\dagger}\left(\mathbf{r}_2\right)\right\rangle & =\frac{1}{(2 \pi)^4} \\ & \int_{-L}^L \int_{-L}^L \int_{-L}^L \int_{-L}^L\left(\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{f}}{^\mathcal{L}} \mathrm{e}^{j k_0 g} d k_{x 1}^i d k_{x 2}^i d k_{u 1}^i d k_{u 2}^i\right) d x_1^{\prime} d x_2^{\prime} d y_1^{\prime} d y_2^{\prime},\end{aligned}$$
(35)$$\begin{aligned} \left\langle\mathbf{L}\left(\mathbf{r}_1\right) \mathbf{N}^{\dagger}\left(\mathbf{r}_2\right)\right\rangle= & -\frac{1}{(2 \pi)^4 \eta_0} \\ & \int_{-L}^L \int_{-L}^L \int_{-L}^L \int_{-L}^L\left(\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{f}}{^\mathcal{K}} \mathrm{e}^{j k_0 g} d k_{x 1}^i d k_{x 2}^i d k_{u 1}^i d k_{u 2}^i\right) d x_1^{\prime} d x_2^{\prime} d y_1^{\prime} d y_2^{\prime},\end{aligned}$$
(36)$$\begin{aligned} \left\langle\mathbf{N}\left(\mathbf{r}_1\right) \mathbf{L}^{\dagger}\left(\mathbf{r}_2\right)\right\rangle= & -\frac{1}{(2 \pi)^4 \eta_0} \\ & \int_{-L}^L \int_{-L}^L \int_{-L}^L \int_{-L}^L\left(\int_{-\infty}^L \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{f}}{^\mathcal{M}} \mathrm{e}^{j k_0 g} d k_{x 1}^i d k_{x 2}^i d k_{u 1}^i d k_{u 2}^i\right) d x_1^{\prime} d x_2^{\prime} d y_1^{\prime} d y_2^{\prime},\end{aligned}$$
and (37)$$\begin{aligned} \left\langle\mathbf{N}\left(\mathbf{r}_1\right) \mathbf{N}^{\dagger}\left(\mathbf{r}_2\right)\right\rangle & =\frac{1}{(2 \pi)^4 \eta_0^2} \\ & \int_{-L}^L \int_{-L}^L \int_{-L}^L \int_{-L}^L\left(\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{f}}{^\mathcal{N}} \mathrm{e}^{j k_0 g} d k_{x 1}^i d k_{x 2}^i d k_{u 1}^i d k_{u 2}^i\right) d x_1^{\prime} d x_2^{\prime} d y_1^{\prime} d y_2^{\prime}.\end{aligned}$$
Here, (38)$$\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{f}}{^{\mathcal{L}}}=\chi^{\prime}\left(\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{L}}{_{1}} \cdot \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{\Phi}}{^{i}} \cdot \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{L}}{}_{2}^{\dagger}\right), $$
(39)$$\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{f}}{^{\mathcal{K}}}=\chi^{\prime}\left(\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{L}}{_{1}} \cdot \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{\Phi}}{^{i}} \cdot \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{N}}{}_{2}^{\dagger}\right), $$
(40)$$\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{f}}{^{\mathcal{M}}}=\chi^{\prime}\left(\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{N}}{_{1}} \cdot \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{\Phi}}{^{i}} \cdot \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{L}}{}_{2}^{\dagger}\right), $$
and (41)$$\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{f}}{^{\mathcal{N}}}=\chi^{\prime}\left(\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{N}}{_{1}} \cdot \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{\Phi}}{^{i}} \cdot \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{N}}{}_{2}^{\dagger}\right), $$
are amplitude dyadics that contains all of the amplitude terms, and (42)$$\begin{aligned} g & =\left[r_s+(\hat{\mathbf{y}} \cdot \hat{\mathbf{v}}) y_2^{\prime}\right] k_{v 2}^i-\left[r_s+(\hat{\mathbf{y}} \cdot \hat{\mathbf{v}}) y_1^{\prime}\right] k_{v 1}^i \\ & +\frac{1}{k_0}\left(x_2^{\prime} k_{x 2}^i-x_1^{\prime} k_{x 1}^i\right)+\frac{(\hat{\mathbf{y}} \cdot \hat{\mathbf{u}})}{k_0}\left(y_2^{\prime} k_{u 2}^i-y_1^{\prime} k_{u 1}^i\right) \\ & +\left(\hat{\mathbf{x}} \cdot \hat{\mathbf{r}}_1 x_1^{\prime}-\hat{\mathbf{x}} \cdot \hat{\mathbf{r}}_2 x_2^{\prime}\right)+\left(\hat{\mathbf{y}} \cdot \hat{\mathbf{r}}_1 y_1^{\prime}-\hat{\mathbf{y}} \cdot \hat{\mathbf{r}}_2 y_2^{\prime}\right) \end{aligned}$$
is a common phase function that contains the phase terms. Without further simplifications, no closed-form expressions exist for the integral relationships given in Eqs. (34)–(37).To simplify the integrals found in parenthesis in Eqs. (34)–(37), we use an asymptotic mathematical technique known as the method of stationary phase (MoSP) [23,24]. In so doing, we assume that the amplitude terms are slowly varying in the interval $({ - \infty ,\infty } )$. We also assume that the phase terms are rapidly oscillating in the interval $({ - \infty ,\infty } )$ except near special points where the rate of change is zero or “stationary.” These special points are critical points of the first kind [24]. Away from these points, the phase terms are rapidly oscillating and the positive and negative contributions of the integrand in Eq. (34)–(37) effectively cancel out.
Using the MoSP to simplify Eqs. (34)–(37) has two implications with regards to the macro-scale scattering geometry given in Fig. 1. The first implication is with respect to the v component of the incident propagation vector ${{\textbf k}^i}$; namely,
(43)$$k_v^i = \sqrt {k_0^2 - {{({k_x^i} )}^2} - {{({k_u^i} )}^2}} . $$
In particular, we assume that $k_v^i \gg k_x^i$ and $k_v^i \gg k_u^i$. As a result, (44)$$k_{v 1,2}^i \approx\left\{\begin{array}{cc} k_0 & \text { in } \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{f}}{}_{0}^{\mathcal{L}}, \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{f}}{}_{0}^{\mathcal{K}}, \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{f}}{}_{0}^{\mathcal{K}}, \text { and } \mathbf{\overset{\scriptscriptstyle\leftrightarrow}{f}}{}_{0}^{\mathcal{N}} \\ k_0-\frac{\left(k_{x 1,2}^i\right)^2}{2 k_0}-\frac{\left(k_{u 1,2}^i\right)^2}{2 k_0} & \text { in } g \end{array}\right..$$
This physically implies that the incident electromagnetic fields are highly directional being predominately directed along the v direction in Fig. 1. The second implication is that the distance from the source-plane origin to the surface-plane origin must be much greater than half the surface length (i.e., ${r_s} \gg L$), which is typically the case for long-range applications. Applying the MoSP to evaluate these integrals physically means that the surface plane is in the far zone of the source plane. As such, using the MoSP to simplify Eqs. (34)–(37) results in the following relationships: (45)$$\left\langle {{\textbf L}({{{\textbf r}_1}} ){{\textbf L}^\dagger }({{{\textbf r}_2}} )} \right\rangle \approx \frac{{k_0^2}}{{{{({2\pi } )}^2}r_s^2}}{\overset{\scriptscriptstyle\leftrightarrow}{\mathcal{L}}}, $$
(46)$$\left\langle {{\textbf L}({{{\textbf r}_1}}){{\textbf N}^\dagger}({{{\textbf r}_2}})} \right\rangle \approx{-} \frac{{k_0^2}}{{{{({2\pi })}^2}{{\eta}_0}r_s^2}}{\overset{\scriptscriptstyle\leftrightarrow}{\mathcal{K}}},$$
(47)$$\left\langle {{\textbf N}({{{\textbf r}_1}}){{\textbf L}^\dagger}({{{\textbf r}_2}})} \right\rangle \approx{-} \frac{{k_0^2}}{{{{({2\pi })}^2}{{\eta}_0}r_s^2}}{\overset{\scriptscriptstyle\leftrightarrow}{\mathcal{M}}},$$
and (48)$$\left\langle {{\textbf N}({{{\textbf r}_1}} ){{\textbf N}^\dagger }({{{\textbf r}_2}} )} \right\rangle \approx \frac{{k_0^2}}{{{{({2\pi } )}^2}\eta _0^2r_s^2}}{\overset{\scriptscriptstyle\leftrightarrow}{\mathcal{N}}}, $$
where $\overset{\scriptscriptstyle\leftrightarrow}{\mathcal{L}}$, $\overset{\scriptscriptstyle\leftrightarrow}{\mathcal{K}}$, $\overset{\scriptscriptstyle\leftrightarrow}{\mathcal{M}}$, and $\overset{\scriptscriptstyle\leftrightarrow}{\mathcal{N}}$ are dyadics that contain all of the amplitude and phase terms evaluated at the critical points of the first kind, which we determine as (49)$$k_{x 1,2}^i \approx \frac{k_0}{r_s} x_{1,2}^{\prime} \quad k_{u 1,2}^i \approx \frac{k_0(\hat{\mathbf{y}} \cdot \hat{\mathbf{u}})}{r_s} y_{1,2}^{\prime}.$$
Please see Section S2 of Supplement 1 for more details. Please also see Section S3, where we explicitly formulate $\overset{\scriptscriptstyle\leftrightarrow}{\mathcal{L}}$, $\overset{\scriptscriptstyle\leftrightarrow}{\mathcal{K}}$, $\overset{\scriptscriptstyle\leftrightarrow}{\mathcal{M}}$, and $\overset{\scriptscriptstyle\leftrightarrow}{\mathcal{N}}$ for different material substrates (i.e., dielectrics, conductors, and a PEC).Provided Sections S2 and S3, we are still left with integrals with respect to the parameterized rough surface. These integrals take the following element-based form:
(50)$$\begin{aligned} \Psi_{\mathrm{mn}} & =\int_{-L}^L \int_{-L}^L \int_{-L}^L \int_{-L}^L \Phi_{\mathrm{mn}}^i\left(\frac{k_0}{r_s} x_1^{\prime}, \frac{k_0}{r_s} x_2^{\prime}, \frac{k_0(\hat{\mathbf{y}} \cdot \hat{\mathbf{u}})}{r_s} y_1^{\prime}, \frac{k_0(\hat{\mathbf{y}} \cdot \hat{\mathbf{u}})}{r_s} y_2^{\prime}\right) \\ & \chi^{\prime}\left(\frac{k_0}{r_s} x_1^{\prime}, \frac{k_0}{r_s} x_2^{\prime}, \frac{k_0(\hat{\mathbf{y}} \cdot \hat{\mathbf{u}})}{r_s} y_1^{\prime}, \frac{k_0(\hat{\mathbf{y}} \cdot \hat{\mathbf{u}})}{r_s} y_2^{\prime} ; x_1^{\prime}-x_2^{\prime}, y_1^{\prime}-y_2^{\prime}\right) \\ & \exp \left[j k_0\left(\hat{\mathbf{x}} \cdot \hat{\mathbf{r}}_1 x_1^{\prime}-\hat{\mathbf{x}} \cdot \hat{\mathbf{r}}_2 x_2^{\prime}\right)\right] \exp \left[j k_0\left(\hat{\mathbf{y}} \cdot \hat{\mathbf{r}_1} y_1^{\prime}-\hat{\mathbf{y}} \cdot \hat{\mathbf{r}}_2 y_2^{\prime}\right)\right] \\ & \exp \left[-j \frac{k_0}{2 r_s}\left(x_1^{\prime 2}-x_2^{\prime 2}\right)\right] \exp \left[-j \frac{k_0(\hat{\mathbf{y}} \cdot \hat{\mathbf{u}})^2}{2 r_s}\left(y_1^{\prime 2}-y_2^{\prime 2}\right)\right] \\ & \exp \left[-j k_0(\hat{\mathbf{y}} \cdot \hat{\mathbf{v}})\left(y_1^{\prime}-y_2^{\prime}\right)\right] d x_1^{\prime} d x_2^{\prime} d y_1^{\prime} d y_2^{\prime}, \end{aligned}$$
where $({\textrm{m} = x,u;\textrm{ n} = x,u} )$. In Eq. (50), $\Phi _{\textrm{mn}}^i$ is equivalent to the Fourier transform of the incident field CSDM elements found in Eq. (1), that is (51)$$\begin{aligned} \Phi _{\textrm{mn}}^i &= \int\limits_{ - \infty }^\infty {\int\limits_{ - \infty }^\infty {\int\limits_{ - \infty }^\infty {\int\limits_{ - \infty }^\infty {W_{\textrm{mn}}^i({{\boldsymbol{\mathrm{\rho}}_1},{\boldsymbol{\mathrm{\rho}}_2}} )} } } } {e ^{jk_{x1}^i{x_1}}}{e ^{ - jk_{x2}^i{x_2}}}{e ^{jk_{u1}^i{u_1}}}{e ^{ - jk_{u2}^i{u_2}}}d{x_1}d{x_2}d{u_1}d{u_2}\\ &= \frac{{{\pi ^2}{A_\textrm{m}}{A_\textrm{n}}{B_{\textrm{mn}}}}}{{({a_{\textrm{mn}}^2 - b_{\textrm{mn}}^2} )}}\textrm{exp} [{ - {{\tilde{a}}_{\textrm{mn}}}{{({k_{x1}^i} )}^2} - {{\tilde{a}}_{\textrm{mn}}}{{({k_{x2}^i} )}^2} + 2{{\tilde{b}}_{\textrm{mn}}}k_{x1}^ik_{x2}^i} ]\\ &\textrm{ }\textrm{exp} [{ - {{\tilde{a}}_{\textrm{mn}}}{{({k_{u1}^i} )}^2} - {{\tilde{a}}_{\textrm{mn}}}{{({k_{u2}^i} )}^2} + 2{{\tilde{b}}_{\textrm{mn}}}k_{u1}^ik_{u2}^i} ]\end{aligned}, $$
and (52)$$\begin{array}{{ccc}} {{a_{\textrm{mn}}} = \frac{1}{{4w_s^2}} + {b_{\textrm{mn}}}}&{{b_{\textrm{mn}}} = \frac{1}{{2\ell _{\textrm{mn}}^2}}}\\ {{{\tilde{a}}_{\textrm{mn}}} = \frac{{{a_{\textrm{mn}}}}}{{4({a_{\textrm{mn}}^2 - b_{\textrm{mn}}^2} )}}}&{{{\tilde{b}}_{mn}} = \frac{{{b_{\textrm{mn}}}}}{{4({a_{\textrm{mn}}^2 - b_{\textrm{mn}}^2} )}}} \end{array}. $$
We can reduce the integrals found in Eq. (50) into a closed-form expression. For this purpose, we first perform the following variable transformations: (53)$$\begin{array}{ll} x_d=x_1^{\prime}-x_2^{\prime} & x_a=x_1^{\prime}+x_2^{\prime} \\ y_d=y_1^{\prime}-y_2^{\prime} & y_a=y_1^{\prime}+y_2^{\prime}, \end{array}$$
so that Eq. (50) simplifies into the following expression: (54)$$\begin{aligned} {\Psi _{mn}} &= \frac{{{\pi ^2}{A_\textrm{m}}{A_\textrm{n}}{B_{\textrm{mn}}}}}{{4({a_{\textrm{mn}}^2 - b_{\textrm{mn}}^2} )}}\textrm{exp} \left[ { - \frac{{k_0^2\sigma_h^2}}{2}({\vartheta_{z1}^2 + \vartheta_{z2}^2} )} \right]\\ &\textrm{ }\int\limits_{ - 2L}^{2L} {\int\limits_{|{{y_a}} |- 2L}^{2L - |{{y_a}} |} {\int\limits_{ - 2L}^{2L} {\int\limits_{|{{x_a}} |- 2L}^{2L - |{{x_a}} |} {\textrm{exp} \left[ {k_0^2\sigma_h^2{\vartheta_{z1}}{\vartheta_{z2}}\textrm{exp} \left( { - \frac{{x_d^2}}{{\ell_h^2}}} \right)\textrm{exp} \left( { - \frac{{y_d^2}}{{\ell_h^2}}} \right)} \right]} } } } \\ &\textrm{ }\textrm{exp} \left[ { - \frac{{k_0^2}}{{2r_s^2}}({{{\tilde{a}}_{\textrm{mn}}} + {{\tilde{b}}_{\textrm{mn}}}} )x_d^2} \right]\textrm{exp} \left[ { - \frac{{k_0^2}}{{2r_s^2}}({{{\tilde{a}}_{\textrm{mn}}} - {{\tilde{b}}_{\textrm{mn}}}} )x_a^2} \right]\\ &\textrm{ }\textrm{exp} \left[ { - \frac{{k_0^2{{({\hat{{\textbf y}} \cdot \hat{{\textbf u}}} )}^2}}}{{2r_s^2}}({{{\tilde{a}}_{\textrm{mn}}} + {{\tilde{b}}_{\textrm{mn}}}} )y_d^2} \right]\textrm{exp} \left[ { - \frac{{k_0^2{{({\hat{{\textbf y}} \cdot \hat{{\textbf u}}} )}^2}}}{{2r_s^2}}({{{\tilde{a}}_{\textrm{mn}}} - {{\tilde{b}}_{\textrm{mn}}}} )y_a^2} \right]\\ &\textrm{ }\textrm{exp} \left[ {j\frac{{{k_0}}}{2}({{\vartheta_{x1}} + {\vartheta_{x2}}} ){x_d}} \right]\textrm{exp} \left[ {j\frac{{{k_0}}}{2}({{\vartheta_{x1}} - {\vartheta_{x2}}} ){x_a}} \right]\\ &\textrm{ }\textrm{exp} \left[ {j\frac{{{k_0}}}{2}({{\vartheta_{y1}} + {\vartheta_{y2}}} ){y_d}} \right]\textrm{exp} \left[ {j\frac{{{k_0}}}{2}({{\vartheta_{y1}} - {\vartheta_{y2}}} ){y_a}} \right]\\ &\textrm{ }\textrm{exp} \left[ { - j\frac{{{k_0}}}{{2{r_s}}}{x_a}{x_d}} \right]\textrm{exp} \left[ { - j\frac{{{k_0}{{({\hat{{\textbf y}} \cdot \hat{{\textbf u}}} )}^2}}}{{2{r_s}}}{y_a}{y_d}} \right]d{x_d}d{x_a}d{y_d}d{y_a} \end{aligned}, $$
and (55)$$\begin{array}{ccc} {{\vartheta _{x1,2}} = \hat{{\textbf x}} \cdot {{\hat{{\textbf r}}}_{1,2}}}&{{\vartheta _{y1,2}} = \hat{{\textbf y}} \cdot {{\hat{{\textbf r}}}_{1,2}} - \hat{{\textbf y}} \cdot \hat{{\textbf v}}}&{{\vartheta _{z1,2}} = \hat{{\textbf z}} \cdot {{\hat{{\textbf r}}}_{1,2}} - \hat{{\textbf z}} \cdot \hat{{\textbf v}}} \end{array}. $$
From here, we must handle the exponential term containing the surface-autocorrelation function (i.e., the first exponential term inside the integrals above).In the literature, this term has been handled in two ways. The first is to expand it in a Taylor series and proceed with the evaluation of the integrals [17]. Mathematically, this approach is applicable to all surfaces; however, because the series is slowly convergent, we limit this approach to smooth-to-moderately rough surfaces. The other approach involves expanding the surface autocorrelation function [cf. Equation (33)] in a Taylor series and retaining only the first and second order terms [17]. This treatment is applicable to very rough surfaces and is the approach taken in this paper. In practice, the criterion ${\sigma _h} \ge 0.5{\lambda _0}$ helps in discerning the transition point from the smooth-to-moderately rough surface regime to the very rough surface regime and is an empirically determined relationship. This criterion is quite physical and says that very rough surfaces are surfaces where the average optical-path difference is greater than one wave. For such surfaces, we should expect complete constructive/destructive interference of the scattered field. This outcome stands in contrast to smooth-to-moderately rough surfaces, where we get partial interference.
2.4 Very rough surfaces
When considering very rough surfaces, we expand the surface autocorrelation function found inside the first exponential term in Eq. (54). We first rewrite the joint-characteristic function found in Eq. (33) in an alternative form, such that
(56)$$\begin{array}{l} \chi ^{\prime}\left( {\frac{{{k_0}}}{{{r_s}}}{{x^{\prime}}_1},\frac{{{k_0}}}{{{r_s}}}{{x^{\prime}}_2},\frac{{{k_0}({\hat{{\textbf y}} \cdot \hat{{\textbf u}}} )}}{{{r_s}}}{{y^{\prime}}_1},\frac{{{k_0}({\hat{{\textbf y}} \cdot \hat{{\textbf u}}} )}}{{{r_s}}}{{y^{\prime}}_2};{x_d},{y_d}} \right) = \\ \textrm{ }\textrm{exp} \left( { - k_0^2\sigma_h^2{\vartheta_{z1}}{\vartheta_{z2}}\left[ {\frac{{{\vartheta_{z1}}}}{{2{\vartheta_{z2}}}} + \frac{{{\vartheta_{z2}}}}{{2{\vartheta_{z1}}}} - \textrm{exp} \left( { - \frac{{x_d^2}}{{\ell_h^2}}} \right)\textrm{exp} \left( { - \frac{{y_d^2}}{{\ell_h^2}}} \right)} \right]} \right) \end{array}. $$
For very rough surfaces, (57)$$k_0^2\sigma _h^2{\vartheta _{z1}}{\vartheta _{z2}} \gg 1, $$
and the alternative form found in Eq. (56) maintains significant value when (58)$$\frac{{{\vartheta _{z1}}}}{{2{\vartheta _{z2}}}} + \frac{{{\vartheta _{z2}}}}{{2{\vartheta _{z1}}}} - \textrm{exp} \left( { - \frac{{x_d^2}}{{\ell_h^2}}} \right)\textrm{exp} \left( { - \frac{{y_d^2}}{{\ell_h^2}}} \right) \approx 0. $$
Since all of the observation points of interest are in the far zone, (59)$$\frac{{{\vartheta _{z1}}}}{{2{\vartheta _{z2}}}} \approx \frac{{{\vartheta _{z2}}}}{{2{\vartheta _{z1}}}}, $$
and the relationship found in Eq. (58) is only possible for small ${x_d}$ and ${y_d}$.With Eqs. (56)–(59) in mind, it makes sense to expand the exponential functions found in Eq. (58) and retain only the first and second order terms, so that
(60)$$\textrm{exp} \left( { - \frac{{x_d^2}}{{\ell_h^2}}} \right)\textrm{exp} \left( { - \frac{{y_d^2}}{{\ell_h^2}}} \right) \approx 1 - \frac{{x_d^2}}{{\ell _h^2}} - \frac{{y_d^2}}{{\ell _h^2}}. $$
Substituting Eq. (60) into Eq. (54) allows us to again separate the integrals over the parameterized rough surface. As mentioned before, this allows for the development of a closed-form expression for very rough surface conditions without having to convert to polar coordinates. Carrying out the subsequent integrations, a closed-form expression results. Specifically, we must complete the square in the exponential terms and use the following integral relationship [24, p. 266]: (61)$$\int_{ - \infty }^\infty {\textrm{exp} ({ - a{t^2}} )} \textrm{exp} ({ - jbt} )dt = \sqrt {{\pi / a}} \textrm{exp} [{ - {{{b^2}} / {({4a} )}}} ]$$
where $a > 0$. In so doing, (62)$$\begin{aligned} {\Psi _{\textrm{mn}}} &= \frac{{4{\pi ^4}r_s^4\ell _h^2{A_\textrm{m}}{A_\textrm{n}}{B_{\textrm{mn}}}}}{{k_0^2|{\hat{{\textbf y}} \cdot \hat{{\textbf u}}} |({a_{\textrm{mn}}^2 - b_{\textrm{mn}}^2} )}}\textrm{exp} \left[ { - \frac{{k_0^2\sigma_h^2}}{2}{{({{\vartheta_{z1}} - {\vartheta_{z2}}} )}^2}} \right]\frac{1}{{\sqrt {{\cal A}_{\textrm{mn}}^x{\cal A}_{\textrm{mn}}^y} }}\\ &\textrm{exp} \left[ { - \frac{{k_0^2r_s^2\ell_h^2}}{{4{\cal D}_{\textrm{mn}}^x{\cal A}_{\textrm{mn}}^x}}({{\cal A}_{\textrm{mn}}^x - r_s^2\ell_h^2} )({\vartheta_{x1}^2 + \vartheta_{x2}^2} )} \right]\\ &\textrm{exp} \left[ { - \frac{{k_0^2r_s^2\ell_h^2}}{{4{\cal D}_{\textrm{mn}}^y{\cal A}_{\textrm{mn}}^y}}({{\cal A}_{\textrm{mn}}^y - r_s^2\ell_h^2{{({\hat{{\textbf y}} \cdot \hat{{\textbf u}}} )}^2}} )({\vartheta_{y1}^2 + \vartheta_{y2}^2} )} \right]\\ &\textrm{exp} \left[ { - \frac{{k_0^2r_s^2}}{{{\cal A}_{\textrm{mn}}^x}}({\ell_h^2{{\tilde{b}}_{\textrm{mn}}} + r_s^2\sigma_h^2{\vartheta_{z1}}{\vartheta_{z2}}} ){{({{\vartheta_{x1}} - {\vartheta_{x2}}} )}^2}} \right]\\ &\textrm{exp} \left[ { - \frac{{k_0^2r_s^2}}{{{{({\hat{{\textbf y}} \cdot \hat{{\textbf u}}} )}^2}{\cal A}_{\textrm{mn}}^y}}({\ell_h^2{{\tilde{b}}_{\textrm{mn}}}{{({\hat{{\textbf y}} \cdot \hat{{\textbf u}}} )}^2} + r_s^2\sigma_h^2{\vartheta_{z1}}{\vartheta_{z2}}} ){{({{\vartheta_{y1}} - {\vartheta_{y2}}} )}^2}} \right]\\ &\textrm{exp} \left( {j\frac{{{k_0}r_s^3\ell_h^2}}{{2{\cal A}_{\textrm{mn}}^x{\cal A}_{\textrm{mn}}^y}}[{{\cal A}_{\textrm{mn}}^y({\vartheta_{x1}^2 - \vartheta_{x2}^2} )+ {\cal A}_{\textrm{mn}}^x({\vartheta_{y1}^2 - \vartheta_{y2}^2} )} ]} \right) \end{aligned}, $$
and (63)$$\begin{array}{l} {\cal D}_{\textrm{mn}}^x = k_0^2\ell _h^2({{{\tilde{a}}_{\textrm{mn}}} + {{\tilde{b}}_{\textrm{mn}}}} )+ 2k_0^2r_s^2\sigma _h^2{\vartheta _{z1}}{\vartheta _{z2}}\\ {\cal D}_{\textrm{mn}}^y = k_0^2\ell _h^2({{{\tilde{a}}_{\textrm{mn}}} + {{\tilde{b}}_{\textrm{mn}}}} ){({\hat{{\textbf y}} \cdot \hat{{\textbf u}}} )^2} + 2k_0^2r_s^2\sigma _h^2{\vartheta _{z1}}{\vartheta _{z2}}\\ {\cal A}_{\textrm{mn}}^x = r_s^2\ell _h^2 + 4({{{\tilde{a}}_{\textrm{mn}}} - {{\tilde{b}}_{\textrm{mn}}}} ){\cal D}_{\textrm{mn}}^x\\ {\cal A}_{\textrm{mn}}^y = r_s^2\ell _h^2{({\hat{{\textbf y}} \cdot \hat{{\textbf u}}} )^2} + 4({{{\tilde{a}}_{\textrm{mn}}} - {{\tilde{b}}_{\textrm{mn}}}} ){\cal D}_{\textrm{mn}}^y \end{array}. $$
Provided the closed-form expression given in Eq. (62), we must satisfy the following conditions: (64)$$\begin{array}{lll} \begin{array}{{ccc}} {L > \frac{{{r_s}}}{{{k_0}{w_s}}}\sqrt {\frac{{ - \ln ({{\delta_x}} )}}{2}\left( {1 + \frac{4}{{\alpha_{\textrm{mn}}^2}}} \right)} }&{L > \frac{{{r_s}}}{{{k_0}{w_s}|{\hat{{\textbf y}} \cdot \hat{{\textbf u}}} |}}\sqrt {\frac{{ - \ln ({{\delta_y}} )}}{2}\left( {1 + \frac{4}{{\alpha_{\textrm{mn}}^2}}} \right)} } \end{array}\\ \end{array}, $$
where (65)$${\alpha _{mn}} = \frac{{{\ell _{mn}}}}{{{w_s}}}$$
is the element-based source ratio, and ${\delta _x}$ and ${\delta _y}$ are again user-defined parameters, respectively. They denote the points at which the exponential functions with respect to ${x_a}$ and ${y_a}$ no longer maintain significant value in Eq. (54). The conditions given in Eq. (64) physically mean that the projected partially coherent incident beam must “fit” on the rough surface. In satisfying these conditions, the ${x_a}$ and ${y_a}$ limits of integration extend to $({ - \infty ,\infty } )$ in Eq. (54). Additionally, if ${\delta _x} = {\delta _y} = \delta $, the second condition in Eq. (64) becomes the most stringent because of the $|{\hat{{\textbf y}} \cdot \hat{{\textbf u}}} |= \cos ({{\theta^i}} )$ in the denominator.The closed-form expression obtained in Eq. (62) is remarkably physical. For instance, the exponential terms on the second and third lines of Eq. (62) are predominately responsible for the angular extent of the scattered SD. These exponential terms are functions of the sum of the squares of the observation projections: ${\vartheta _{x1,2}}$ and ${\vartheta _{y1,2}}$. On the other hand, the exponential terms on the fourth and fifth lines of Eq. (62) determine the angular extent of the scattered SDoC. Note that these terms depend on the differences of the observation projections (i.e., $|{{\vartheta_{x1}} - {\vartheta_{x2}}} |$ and $|{{\vartheta_{y1}} - {\vartheta_{y2}}} |$). In turn, we can state that the closed-form expression obtained in Eq. (62) allows the scattered field CSDM to maintain a GSM form with respect to ${\vartheta _{x1,2}}$ and ${\vartheta _{y1,2}}$. We examine these points more closely in the sections that follow.
3. Comparison to the Priest and Meier pBRDF
In this section, we compare the 3D vector solution to a pBRDF developed by Priest and Meier [1]. In general, the Priest and Meier pBRDF assumes fully incoherent illumination; thus, there is no coherence information contained within the solution. Instead, the Priest and Meier pBRDF provides the Mueller matrices for statistically rough surfaces that are characterized by G-G PDFs [cf. Equation (31)]. As a result, if we assume unpolarized illumination, then the first column of the Mueller matrix (given by the pBRDF solution) becomes the scattered Stokes vector [25]. Provided this scattered Stokes vector, we can then compare the normalized SD and the DoP between the pBRDF solution and the 3D vector solution.
By formulating the 3D vector solution in a manner consistent with Wolf’s unified theory of coherence and polarization [4,5], all physical implications inherent in Wolf’s work apply here. Accordingly, we can readily formulate the scattered normalized SD $S_N^s({\textbf r} )$ and the scattered DoP ${P^s}({\textbf r} )$ from the closed-form expression developed above for the scattered field CSDM $\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{W}}{^{s}}({{{\textbf r}_1},{{\textbf r}_2}})$. For this purpose, we use the following relationships [5]:
(66)$$S_N^s(\mathbf{r})=\frac{\operatorname{Tr}\left\{\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{W}}{^{s}}(\mathbf{r}, \mathbf{r})\right\}}{\max \left[\operatorname{Tr}\left\{\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{W}}{^{s}}(\mathbf{r}, \mathbf{r})\right\}\right]},$$
and (67)$$P^s(\mathbf{r})=\sqrt{1-\frac{4 \operatorname{Det}\left\{\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{W}}{^{s}}(\mathbf{r}, \mathbf{r})\right\}}{\left(\operatorname{Tr}\left\{\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{W}}{^{s}}(\mathbf{r}, \mathbf{r})\right\}\right)^2}},$$
where $\textrm{Tr}\{\cdots \}$ denotes the trace operation, $\textrm{Det}\{\cdots \}$ denotes the determinant operation, and ${\textbf r} = {{\textbf r}_{1,2}}$ corresponds to a single observation point. These relationships contain measurable quantities and serve as metrics with which to compare the 3D vector solution to the Priest and Meier pBRDF.For this comparison, we used a $5.08\textrm{ cm } \times \textrm{ 5}\textrm{.08 cm}$ Labsphere Infragold coupon [26]. We also used a realizable laboratory setup [27], where ${\lambda _0} = 10.6$ µm, ${r_s} = {r_1}_{,2} = 185$ cm, and ${w_s} = 1.9$ mm. As such, the Labsphere Infragold coupon maintains the following complex index of refraction: $n = \textrm{13}.\textrm{45} - j\textrm{63}.\textrm{62}$ [28]. Note that a KLA Tencor Alpha-Step IQ Surface Profiler [29] determined the surface statistics of the Labsphere Infragold coupon as ${\sigma _h} = 11.09$ µm and ${\ell _h} = 116.9$ µm using four 1 cm scans (step size $0.2\mathrm{\ \mu m}$). These surface statistics relate to very rough surface conditions [cf. Equation (57)].
3.1 Normalized SD comparison
Assuming unpolarized illumination (i.e., ${A_x} = {A_u}$, and ${B_{xu}} = {B_{ux}} = 0$), Fig. 2 shows the comparison for the scattered normalized SD at ${\theta ^i} = {20^ \circ }$. Note the excellent agreement between the Priest and Meier pBRDF and 3D vector solution for various scattering geometries. In practice, we can obtain the scattered normalized SD from the pBRDF solution by cosine correcting the first term of the scattered Stokes vector (i.e., multiplying by $\cos ({{\theta^s}} )$ [30]) and thereafter dividing by the max value.
3.2 DoP comparison
Again, assuming unpolarized illumination, Fig. 3 shows the comparison for the scattered DoP at ${\theta ^i} = {20^ \circ }$. Note the excellent agreement between the Priest and Meier pBRDF and the 3D vector solution for various scattering geometries. In general, the scattered DoP exists where light exists, which makes sense considering that, by definition, the scattered DoP depends on the scattered SD [cf. Equations (66) and (67)].
4. Closed-form expression for the SDoC
Once again, by formulating the 3D vector solution in a manner consistent with Wolf’s unified theory of coherence and polarization [4,5], all physical implications inherent in Wolf’s work apply here. Accordingly, we can readily formulate the scattered SDoC ${\mu ^s}({{{\textbf r}_1},{{\textbf r}_2}} )$ from the closed-form expression developed above for the scattered field CSDM $\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{W}}{^{s}}\left(\mathbf{r}_1, \mathbf{r}_2\right)$. For this purpose, we use the following relationship [5]:
(68)$$\mu^s\left(\mathbf{r}_1, \mathbf{r}_2\right)=\frac{\operatorname{Tr}\left\{\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{W}}{^{s}}\left(\mathbf{r}_1, \mathbf{r}_2\right)\right\}}{\sqrt{\operatorname{Tr}\left\{\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{W}}{^{s}}\left(\mathbf{r}_1, \mathbf{r}_1\right)\right\}} \sqrt{\operatorname{Tr}\left\{\mathbf{\overset{\scriptscriptstyle\leftrightarrow}{W}}{^{s}}\left(\mathbf{r}_2, \mathbf{r}_2\right)\right\}}},$$
where again, $\textrm{Tr}\{\cdots \}$ denotes the trace operation and ${\textbf r} = {{\textbf r}_{1,2}}$ corresponds to a single observation point. This relationship contains a measurable quantity and serves as a metric with which to compare the 3D vector solution to the 2D scalar-equivalent solution [2,3] and a paraxial solution [7,9].4.1. Comparison to the 2D scalar-equivalent solution
In order to compare the 3D vector solution to the previously validated 2D scalar-equivalent solution [2,3], we assume horizontally polarized (s-pol) illumination and an in-plane scattering geometry (i.e., ${A_u} = {B_{ux}} = {B_{xu}} = 0$ and ${\phi ^s} = \phi _{1,2}^s = {90^{^\circ} }$). This provides the setup needed to make a fair comparison between the two solutions.
For this comparison, we can formulate a closed-form expression for the angular SDoC radius, which describes the angular extent over which the scattered field is correlated. In general, the angular SDoC radius provides a gauge for the average speckle size observed in the far zone and is a quantity of importance for long-range applications that deal with speckle phenomena. Note that the formulation presented here is highly analogous to that performed for the 2D scalar-equivalent solution [2,3].
Assuming that Eq. (64) holds, so that the incident illumination “fits” on the rough surface, the exponential terms on the fourth and fifth lines of Eq. (62), in general, determine the angular extent of the scattered SDoC. Provided s-pol illumination and an in-plane scattering geometry, only ${\Psi _{xx}}$ exists, and the difference of the observation projections simplify, such that
(69)$$|{{\vartheta_{x1}} - {\vartheta_{x2}}} |= |{\sin ({\theta_1^2} )\cos ({\phi_1^s} )- \sin ({\theta_2^2} )\cos ({\phi_2^s} )} |= 0$$
and (70)$$|{{\vartheta_{y1}} - {\vartheta_{y2}}} |= |{\sin ({\theta_1^s} )\sin ({\phi_1^s} )- \sin ({\theta_2^s} )\sin ({\phi_2^s} )} |= |{\sin ({\theta_1^s} )- \sin ({\theta_2^s} )} |. $$
Consequently, the following “correlation” exponential $\gamma$ results from Eq. (62) for ${\Psi _{xx}}$: (71)$$\begin{aligned} \gamma &= \textrm{exp} \left[ { - \frac{{k_0^2r_s^2}}{{{{({\hat{{\textbf y}} \cdot \hat{{\textbf u}}} )}^2}{\cal A}_{xx}^y}}({\ell_h^2{{\tilde{b}}_{xx}}{{({\hat{{\textbf y}} \cdot \hat{{\textbf u}}} )}^2} + r_s^2\sigma_h^2{\vartheta_{z1}}{\vartheta_{z2}}} ){{[{\sin ({\theta_1^s} )- \sin ({\theta_2^s} )} ]}^2}} \right]\\ &\approx \textrm{exp} \left[ { - \frac{{r_s^2{\cal D}_{xx}^y}}{{2{{({\hat{{\textbf y}} \cdot \hat{{\textbf u}}} )}^2}{\cal A}_{xx}^y}}{{[{\sin ({\theta_1^s} )- \sin ({\theta_2^s} )} ]}^2}} \right] \end{aligned}$$
Upon setting $\gamma$ equal to ${1 / \textrm{e}}$, the following expression results: (72)$$\begin{aligned} {|{\sin ({\theta_1^s} )- \sin ({\theta_2^s} )} |_{{1 / \textrm{e}}}} &\approx \frac{{|{\hat{{\textbf y}} \cdot \hat{{\textbf u}}} |}}{{{r_s}}}\sqrt {\frac{{2A_{xx}^y}}{{{\cal D}_{xx}^y}}} \\ &= \frac{{|{\hat{{\textbf y}} \cdot \hat{{\textbf u}}} |}}{{{r_s}}}\sqrt {\frac{{8w_s^2\ell _{xx}^2}}{{\ell _{xx}^2 + 4w_s^2}} + \frac{{2r_s^2\ell _h^2{{({\hat{{\textbf y}} \cdot \hat{{\textbf u}}} )}^2}}}{{k_0^2[{\ell_h^2w_s^2{{({\hat{{\textbf y}} \cdot \hat{{\textbf u}}} )}^2} + 2r_s^2\sigma_h^2{\vartheta_{z1}}{\vartheta_{z2}}} ]}}} \end{aligned}$$
Because the magnitude of the argument of $\gamma$ is large (specifically the $k_0^2r_s^4\sigma _h^2{\vartheta _{z1}}{\vartheta _{z2}}$ term), $\sin ({\theta_1^s} )\approx \sin ({\theta_2^s} )$ for Eq. (71) to have significant value. This implies that $\theta _1^s \approx \theta _2^s$ and that $\gamma$ is approximately a function of $\Delta {\theta ^s} = \theta _1^s - \theta _2^s$. Using this insight, Eq. (72) simplifies, since (73)$$\begin{aligned} \sin ({\theta_1^s} )&= \sin ({\theta_2^s + \Delta {\theta^s}} )\\ &= \sin ({\theta_2^s} )\cos ({\Delta {\theta^s}} )+ \cos ({\theta_2^s} )\sin ({\Delta {\theta^s}} )\\ &\approx \sin ({\theta_2^s} )+ \cos ({\theta_2^s} )\Delta {\theta ^s} \end{aligned}, $$
(74)$$\begin{aligned} \cos ({\theta_1^s} )&= \cos ({\theta_2^s + \Delta {\theta^s}} )\\ &= \cos ({\theta_2^s} )\cos ({\Delta {\theta^s}} )- \sin ({\theta_2^s} )\sin ({\Delta {\theta^s}} )\\ &\approx \cos ({\theta_2^s} )\end{aligned}, $$
and (75)$$\frac{{{\vartheta _{z1}}{\vartheta _{z2}}}}{{{{({\hat{{\textbf y}} \cdot \hat{{\textbf u}}} )}^2}}} = \frac{{[{\cos ({\theta_1^s} )+ \cos ({{\theta^i}} )} ][{\cos ({\theta_2^s} )+ \cos ({{\theta^i}} )} ]}}{{{{\cos }^2}({{\theta^i}} )}} \approx {\left[ {1 + \frac{{\cos ({\theta_2^s} )}}{{\cos ({{\theta^i}} )}}} \right]^2}. $$
After some simple algebra, the expression for the angular SDoC radius becomes (76)$${|{\Delta {\theta^s}} |_{{1 / \textrm{e}}}} \approx \frac{1}{{\varpi {r_s}}}\sqrt {\frac{{8w_s^2}}{{1 + {{({{2 / {{\alpha_{xx}}}}} )}^2}}} + \frac{2}{{k_0^2\sigma _{h^{\prime}}^2{{({1 + \varpi } )}^2}\left( {1 + {{\left[ {\frac{{{\Omega _s}}}{{\sigma_{h^{\prime}}^2({1 + \varpi } )}}} \right]}^2}} \right)}}}, $$
where ${\alpha _{xx}} = {{{\ell _{xx}}} / {{w_s}}}$ is a source ratio [cf. Equation (65)], $\varpi = {{\cos ({\theta_2^s} )} / {\cos ({{\theta^i}} )}}$ is a projection ratio, ${\Omega _s} = {{{w_s}} / {{r_s}}}$ is the source half angle (viewed from the rough surface), and (77)$${\sigma _{h^{\prime}}} = \sqrt 2 \frac{{{\sigma _h}}}{{{\ell _h}}}$$
is the surface slope standard deviation [17]. Based on the assumptions used, valid surface slope standard deviations must satisfy the condition ${\sigma _{h^{\prime}}} \le 0.25\textrm{ rad}$ [15,31].For all intents and purposes, we can neglect the term involving ${{{\Omega _s}} / {{\sigma _{h^{\prime}}}}}$ in Eq. (76). This ratio results in values on the order of ${10^{ - 4}}$ for long-range applications. With that said, we can also claim that the source term contained in the radical above is much greater than the surface term. Thus, factoring out the source term and using the binomial approximation yields
(78)$${|{\Delta {\theta^s}} |_{{1 / \textrm{e}}}} \approx \frac{{{\Omega _s}}}{\varpi }\left( {2\sqrt {\frac{2}{{1 + {{({{2 / {{\alpha_{xx}}}}} )}^2}}}} + \frac{{\sqrt {2[{1 + {{({{2 / {{\alpha_{xx}}}}} )}^2}} ]} }}{{4k_0^2w_s^2\sigma_{h^{\prime}}^2{{({1 + \varpi } )}^2}}}} \right) \approx \frac{{2{\Omega _s}}}{\varpi }\sqrt {\frac{2}{{1 + {{({{2 / {{\alpha_{xx}}}}} )}^2}}}}$$
For most cases of interest, we can neglect the second term contained within the parenthesis in Eq. (78). It only provides a “small” correction to the angular SDoC radius due to the surface parameters. As a result, the angular SDoC radius becomes a function of only the source parameters. This is highly analogous to the result obtained by the 2D scalar-equivalent solution [2,3]. It is also consistent with the classic, narrow-band, fully coherent illumination result derived by Goodman [32].In the comparison that follows, we used a nominal laboratory setup, where ${\lambda _0} = 1.064$ µm, ${r_s} = {r_1}_{,2} = 2$ m, and ${w_s} = 2$ mm. We also used partially coherent illumination of NKB7 glass ($n = \textrm{1}\textrm{.507}$ [28]) at ${\theta ^i} = {0^ \circ }$ and ${\theta ^i} = {56.4^ \circ }$ (i.e., Brewster’s angle). In support, we varied the coherence of the incident illumination, so that ${\alpha _{xx}} = 2\textrm{ and 0}\textrm{.5}$ [cf. Equation (65)]. These values relate to a relatively coherent source and to a relatively incoherent source, respectively. In addition, we varied the surface roughness of the NKB7 glass, so that ${\sigma _{h^{\prime}}} = \textrm{0}\textrm{.1}$ rad and 0.2 rad [cf. Equation (77)], where ${\ell _h} = 100{\lambda _0}$ and ${\sigma _h} = {{{\sigma _{h^{\prime}}}{\ell _h}} / {\sqrt 2 }}$. These values both correspond to very rough surface conditions, respectively.
Assuming partially coherent s-pol illumination (i.e., ${A_x} = 1$ and ${A_u} = {B_{xu}} = {B_{ux}} = 0$), Fig. 4 shows the comparison for the magnitude of the scattered SDoC. Note the excellent agreement between the 2D scalar-equivalent solution and the 3D vector solution. Also note the excellent agreement for the angular SDoC radius [cf. Equation (78) and vertical lines in Fig. 4]. This agreement supports the finding that the angular SDoC radius approximately depends on only the source parameters.
With Fig. 4 in mind, Fig. S2 of Supplement 1 shows the comparison for the scattered normalized SD. In general, the results deviate between the 2D and 3D solutions as ${\sigma _{h^{\prime}}}$ increases. This deviation is deterministic in practice. For s-pol illumination and an in-plane scattering geometry, the 2D and 3D solutions have the same functional dependence in the exponential term which predominately drives the angular extent of the scattered SD. In practice, we can use the angular SD radius $\theta _{{1 / \textrm{e}}}^s$ to identify the ${1 / \textrm{e}}$ location for the scattered normalized SD. Assuming near normal incidence,
(79)$$\theta _{{1 / \textrm{e}}}^s \approx {\cos ^{ - 1}}\left( {\frac{2}{{1 + 2\sigma_{h^{\prime}}^2}} - 1} \right), $$
which states that the angular SD radius is approximately a function of only the surface parameters. Please see Section S4 of Supplement 1 for more details.Provided Section S4, there is a small difference contained in the amplitude terms when comparing the 2D and 3D solutions. As formulated, the 3D vector solution contains a ${1 / {\sqrt {{\cal A}_{xx}^x{\cal A}_{xx}^y} }}$ amplitude factor [cf. Equation (63)], whereas the 2D scalar-equivalent solution contains only a ${1 / {\sqrt {{\cal A}_{xx}^y} }}$ amplitude factor. These amplitude factors (in addition to polarization terms) appear in front of the exponential term which predominantly drives the angular extent of the scattered normalized SD. Functionally, these amplitude factors tend to “push” the scattered normalized SD to the right, whereas the polarization terms tend to “pull” the scattered normalized SD to the left (for positive incident angles—the opposite is true for negative incident angles).
As shown in Fig. S2(c) and (d), as ${\sigma _{h^{\prime}}}$ increases, the SD peak for the 2D scalar-equivalent solution shifts more and more to the left when compared to the 3D vector solution. The full-wave 2D method-of-moments (MoM) solution of Basu et al. shows similar behavior [20,21]. This outcome is due to the lack of the aforementioned amplitude factor in the 2D solutions. Note that comparisons with a full-wave 3D MoM solution are unrealizable at optical wavelengths with the setup used in Figs. 4 and S2—the computational sampling and memory requirements are too great for modern-day desktop computers.
4.2 Comparison to a paraxial solution
In practice, we can reformulate the 3D vector solution within a Cartesian coordinate system. With that said, we must use the following relationships:
(80)$$\begin{array}{ccc} \begin{array}{{ccc}} {{\vartheta _{x1,2}} = \frac{{{x_{1,2}}}}{{{r_{1,2}}}}}&{{\vartheta _{y1,2}} = \frac{{{y_{1,2}}}}{{{r_{1,2}}}} - \sin ({{\theta^i}} )}&{{\vartheta _{z1,2}} = \frac{{{z_{1,2}}}}{{{r_{1,2}}}} + \cos ({{\theta^i}} )} \end{array}\\ \begin{array}{*{20}{c}} {{{\boldsymbol{\hat{\mathrm{\theta}}}}_{1,2}} \cdot \hat{{\textbf x}} = \frac{{{z_{1,2}}{x_{1,2}}}}{{{r_{1,2}}{\rho _{1,2}}}}}&{{{\boldsymbol{\hat{\mathrm{\theta}}}}_{1,2}} \cdot \hat{{\textbf y}} = \frac{{{z_{1,2}}{y_{1,2}}}}{{{r_{1,2}}{\rho _{1,2}}}}}&{{{\boldsymbol{\hat{\mathrm{\theta}}}}_{1,2}} \cdot \hat{{\textbf z}} ={-} \frac{{{\rho _{1,2}}}}{{{r_{1,2}}}}} \end{array}\\ \begin{array}{*{20}{c}} {{{\boldsymbol{\hat{\mathrm{\phi}}}}_{1,2}} \cdot \hat{{\textbf x}} ={-} \frac{{{y_{1,2}}}}{{{\rho _{1,2}}}}}&{{{\boldsymbol{\hat{\mathrm{\phi}}}}_{1,2}} \cdot \hat{{\textbf y}} = \frac{{{x_{1,2}}}}{{{\rho _{1,2}}}}} \end{array}\\ \begin{array}{*{20}{c}} {{\rho _{1,2}} = \sqrt {x_{1,2}^2 + y_{1,2}^2} }&{{r_{1,2}}} \end{array} = \sqrt {x_{1,2}^2 + y_{1,2}^2 + z_{1,2}^2} \end{array}. $$
Provided these relationships, we can compare the 3D vector solution to an ABCD-matrix approach outlined by Korotkova et al. [7,8]. In general, the ABCD-matrix approach describes paraxial-wave propagation through any complex optical system. When modeling rough surface scattering using the ABCD-matrix approach, a phase-screen transmittance function characterizes the scattering process for very rough surface conditions. The inclusion of a soft-Gaussian aperture in the model also accounts for the size of the scattering surface and accompanying diffraction effects.In the comparison that follows, we used horizontally polarized (p-pol) illumination, so that ${A_u} = 1$, ${A_x} = {B_{xu}} = {B_{ux}} = 0$, ${\lambda _0} = 1.064$ µm, ${r_s} = {r_1}_{,2} = 10$ km, ${w_s} = {{2.54} / 2}$ cm, ${\ell _{xx}} = {\ell _{uu}} = 0.25{w_s}$, and ${\theta ^i} = {0^ \circ }$. This nominal long-range setup corresponds to partially coherent illumination at normal incidence. For a fair comparison, we removed the soft-Gaussian aperture in the ABCD solution to ensure that the incident illumination “fits” on the rough surface. In addition, we assumed an idealized perfectly reflecting rough surface, where ${\sigma _h} = 10{\lambda _0}$ and ${\ell _h} = 100{\lambda _0}$. Such surface statistics correspond to very rough surface conditions [cf. Equation (57)]. Provided this setup, Fig. 5 shows a comparison between the ABCD paraxial solution and the 3D vector solution.
The magnitude of the scattered SDoC [Fig. 5(a)] shows excellent agreement. Here, the vertical line depicts the angular SDoC radius [cf. Equation (78)] multiplied by the propagation distance to the observation plane. Note that it correctly identifies the ${1 / \textrm{e}}$ location for the scattered SDoC. In addition, the normalized scattered SD [Fig. 5(b)] shows relatively good agreement between the two solutions. The wings of the power distribution associated with the 3D vector solution tend to be wider than those associated with the ABCD solution; nonetheless, the results show good agreement out to the vertical lines, which depict the angular SD radius [cf. Equation (79)] multiplied by the propagation distance to the observation plane. Also note that it correctly identifies the ${1 / \textrm{e}}$ location for the scattered normalized SD. Furthermore, the scattered DoP [Fig. 5(c)] shows excellent agreement between the two solutions. Based on the assumptions used in this paper, no de-polarization occurs upon scattering from the perfectly reflecting rough surface [14], at least for isotropic beam parameters (i.e., when ${\ell _{xx}} = {\ell _{uu}}$). This outcome is also the case for s-pol illumination at normal incidence and is consistent with the 2D scalar-equivalent solution [2,3]. Please see Section S5 of Supplement 1 for an example where this outcome is not the case. Also see Code 1 (Ref. [33]) for MATLAB code to plot all figures in the main document (Figs. 2–5) and Supplement 1 (Figs. S2–S4).
5. Conclusion
This paper formulated a 3D vector solution for the far-field scattering of spatially partially coherent laser-beam illumination from statistically rough surfaces. Compared to previous efforts, the 3D vector solution formulated here significantly extends the rough surface scattering literature, since most of the theory developed to date only considered the effects of fully coherent or fully incoherent illumination in the formation of solutions. Moreover, the 3D vector solution considered three different material substrates: dielectrics, conductors, and a PEC. By incorporating the effects of the material substrates, the 3D vector solution transcends previous efforts, which included the effects of partially coherent beam illumination but not the effects of the material substrates.
To develop the 3D vector solution, we used the physical optics approximation (Kirchhoff boundary conditions). We also used a GSM form for the incident field CSDM. This choice allowed for the development of a closed-form expression for the scattered field CSDM that is applicable to very rough surfaces. The closed-form expression is extremely physical, and under certain circumstances, maintained a GSM form. This outcome agrees with published results valid only in the paraxial regime.
By using the CSDM notation within the 3D vector solution, all aspects inherent in Wolf’s unified theory of coherence and polarization apply here. Accordingly, the 3D vector solution enabled a rigorous comparison to the Priest and Meier pBRDF in terms of the normalized SD and the DoP. Overall, the comparison showed excellent agreement for the scattered normalized SD density and the scattered DoP. Based on this agreement, the 3D vector solution also enabled an extension to the Priest and Meier pBRDF to account for the effects of active illumination. In particular, the 3D vector solution enabled the development of a closed-form expression for the scattered SDoC. This expression provides a gauge for the average speckle size based on the spatial-coherence properties of the laser source. As such, this extension is of broad interest to long-range applications that deal with speckle phenomena.
Acknowledgments
The Authors would like to thank O. Korotkova for her helpful hints regarding the ABCD paraxial solution. Note that the results of this paper go well beyond the initial results contained in a conference paper by the authors [6]. With that said, some of the verbiage used in both papers are appropriately similar. The views expressed are those of the authors and do not necessarily reflect the official policy or position of the Department of the Air Force, the Department of Defense, or the U.S. government.
Disclosures
The authors declare no conflicts of interest.
Data availability
No data were generated or analyzed in the presented research.
Supplemental document
See Supplement 1 for supporting content.
References
1. R. G. Priest and S. R. Meier, “Polarimetric microfacet scattering theory with applications to absorptive and reflective surfaces,” Opt. Eng. 41(5), 988–993 (2002). [CrossRef]
2. M. W. Hyde, S. Basu, S. J. Cusumano, and M. F. Spencer, “Scalar wave solution for the scattering of a partially coherent beam from a statistically rough metallic surface,” Proc. SPIE 8550, 85503A (2012). [CrossRef]
3. M.W. Hyde IV, S. Basu, M. F. Spencer, S. J. Cusumano, and S. T. Fiorino, “Physical optics solution for the scattering of a partially-coherent wave from a statistically rough material surface,” Opt. Express 21(6), 6807–6825 (2013). [CrossRef]
4. E. Wolf, “Unified theory of coherence and polarization of random electromagnetic beams,” Phys. Lett. A 312(5-6), 263–267 (2003). [CrossRef]
5. E. Wolf, Introduction to the Theory of Coherence and Polarization of Light (Cambridge, 2007).
6. M. F. Spencer, M. W. Hyde, S. Basu, and M. A. Marciniak, “The scattering of partially coherent electromagnetic beam illumination from a statistically rough surface modeled as a perfect electrical conductor,” Proc. SPIE 9205, 92050J (2014). [CrossRef]
7. O. Korotkova, Y. Cai, and E. Watson, “Stochastic electromagnetic beams for LIDAR systems operating in turbulent atmosphere,” Appl. Phys. B 94(4), 681–690 (2009). [CrossRef]
8. O. Korotkova, Random Light Beams: Theory and Applications (CRC Press, 2013).
9. J. A. Ogilvy, Theory of Wave Scattering From Random Rough Surfaces (IOP, 1991).
10. A. A. Maradudin, ed., Light Scattering and Nanoscale Surface Roughness (Springer, 2007).
11. H. Roychowdhury and O. Korotkova, “Realizability conditions for electromagnetic Gaussian Schell-model sources,” Opt. Commun. 249(4-6), 379–385 (2005). [CrossRef]
12. T. B. Hansen and A. D. Yaghjian, Plane-Wave Theory of Time-Domain Fields Near Field Scanning Applications (IEEE Press, 1999).
13. C. A. Balanis, Advanced Engineering Electromagnetics, 2nd ed. (Wiley, 2012).
14. P. Beckmann, The Depolarization of Electromagnetic Waves (Golem Press, 1968).
15. F. T. Ulaby, R. K. Moore, and A. K. Fung, Microwave Remote Sensing: Active and Passive (Artech, 1986).
16. M. Nieto-Vesperinas, Scattering and Diffraction in Physical Optics, 2nd ed. (World Scientific, 2006).
17. P. Beckmann and A. Spizzichino, The Scattering of Electromagnetic Waves form Rough Surfaces (Pergamon Press, 1963).
18. J. W. Goodman, Speckle Phenomena in Optics: Theory and Applications (Roberts and Company, 2007).
19. C. A. Mack, “Analytic form for the power spectral density in one, two, and three dimensions,” J. Micro/Nanolithogr. MEMS MOEMS 10(4), 040501 (2011). [CrossRef]
20. S. Basu, M. W. Hyde IV, S. J. Cusumano, M. A. Marciniak, and S. T. Fiorino, “Validity of using a Gaussian Schell model for extended beacon studies,” Proc. SPIE 8380, 83800E–83800E-13 (2012). [CrossRef]
21. S. Basu, M. W. Hyde IV, S. J. Cusumano, M. A. Marciniak, and S. T. Fiorino, “Examining the validity of using Gaussian Schell-model source to model the scattering of a fully coherent Gaussian beam from a rough impedance surface,” Opt. Eng. 52(2), 038001 (2013). [CrossRef]
22. S. Chatterjee, V. C. Vani, and R. K. Banyal, “Intensity profile of light scattered from a rough surface,” Appl. Opt. 52(24), 6000–6010 (2013). [CrossRef]
23. A. Ishimaru, Electromagnetic Wave Propagation, Radiation, and Scattering (Prentice Hall, 1990).
24. L. Mandel and E. Wolf, Optical Coherence and Quantum Optics (Cambridge, 1995).
25. E. Collett, Polarized Light in Fiber Optics (PolaWave Group, 2003).
26. LabSphere, Inc., “Reflectance Coatings and Materials,” (LabSphere, Inc., 2023), https://www.labsphere.com/wp-content/uploads/2023/06/ReflectanceCoatingsMaterialsTechGuide.pdf.
27. J. C. Vap, S. E. Nauyoks, and M. A. Marciniak, “Optimization of a dual-rotating-retarder polarimeter as applied to a tunable infrared Mueller-matrix scatterometer,” Meas. Sci. Technol. 24(5), 055901 (2013). [CrossRef]
28. M. N. Polyanskiy, Refractive Index Database (RefractiveIndex.Info, 2014), http://refractiveindex.info.
29. KLA-Tencor Corporation, “Alpha-step IQ Surface Profiler” (KLA-Tencor Corporation, 2012), https://www.kla.com/products/instruments/innovation-history.
30. J. C. Stover, Optical Scattering Measurements and Analysis, 3rd ed. (SPIE Press, 2012).
31. A. Ishimaru, C. Le, Y. Kuga, L. A. Sengers, and T. K. Chan, “Polarimetric scattering theory for high slope rough surfaces,” Prog. Electromagn. Res. 14, 1–36 (1996). [CrossRef]
32. J. W. Goodman, “Statistical properties of laser speckle patterns,” in Laser Speckle and Related Phenomena, Vol. 9 of Topics in Applied Physics, J. C. Dainty, ed. (Springer-Verlag, 1975, pp. 9–75).
33. M. F. Spencer, “Active-illumination extension to the Priest and Meier pBRDF: MATLAB code,” figshare (2023), https://doi.org/10.6084/m9.figshare.23915856.
Open Access
Add to BibSonomy
