The Small Perturbation Method for 2D scattering from Gaussian random rough surface

Contents

The scattering amplitude

For an incident plane wave impinging from air a $y$-invariant surface $z=\eta(x)$ with an angle $\theta_i$ from the $z$-axis, the field in the air above the surface region (that is $z>\max\eta)$ writes as the plane wave sum : $$ \psi(x,z)=e^{i(k_0x-q_0z)}+\int\frac{s(k,k_0)}{q}e^{i(kx+qz)}dk $$ for an implicit $e^{-i\omega t}$ time dependency. With wavenumber $K_0=\frac{\omega}{c}=\frac{2\pi}{\lambda}$, $k_0=K_0\sin\theta_i$ and $q_0=K_0\cos\theta_i$ define the incident wavevector. For the scattering one, $q$ is chosen so that $k^2+q^2=K_0^2$ with $0\le\arg q\le\pi/2$.

lambda=1;       % the EM wavelength is chosen as length unit
K0=2*pi/lambda;
qf=@(k,ep) sqrt(K0^2*ep-k.^2); % ep=1 in air

The perturbative expansion

Assuming $\epsilon=\eta/\lambda\ll 1$ and following the small perturbation method (SPM), the scattering amplitude $s(k,k_0)$ is expanded into a series : $$ s(k,k_0)=B_0(k_0)\delta(k-k_0)+B_1(k,k_0)\tilde\eta(k-k_0)+\int B_2(k,x,k_0) \tilde\eta(k-x)\tilde\eta(x-k_0)dx+\int B_3(k,x,y,k_0) \tilde\eta(k-x)\tilde\eta(x-y)\tilde\eta(y-k_0)dxdy+O(\epsilon^4) $$ where $\tilde\eta$ is the Fourier transform of the height function $\eta$ $$ \tilde\eta(k)=\frac{1}{2\pi}\int\eta(x) e^{ikx}dx $$ and the coefficients $B_n$ are roughness-independent.

The Dirichlet boundary condition

Enforcing $\psi(x,\eta(x))=0$ on the surface leads to :

BD0=@(k0) -qf(k0,1);
BD1=@(k,k0) 2*1i*qf(k,1).*qf(k0,1);
BD2=@(k,x,k0) 2*qf(k,1).*qf(x,1).*qf(k0,1);

The Neumann boundary condition

Enforcing $(\mathbf{\hat n\cdot grad\,}\psi(x,\eta(x))=0$ on the surface leads to :

BN0=@(k0) +qf(k0,1);
BN1=@(k,k0) -2*1i*(K0^2-k.*k0);
BN2=@(k,x,k0) -2*(K0^2-k.*x).*(K0^2-k0.*x)./qf(x,1);

The transmission problem for S-polarized electromagnetic waves

Assuming continuity of both the field $\psi$ and its normal derivative $\mathbf{\hat n\cdot grad\,}\psi$ when passing from air (wavenumber K_0) to the scattering medium (wavenumber $K_0\sqrt{\varepsilon}$), we first define reflection $r_s$ coefficient for the plane : $$ r_s(\theta)=\frac{\cos\theta-\sqrt{\varepsilon-\sin^2\theta}}{\cos\theta-\sqrt{\varepsilon-\sin^2\theta}} $$ and associated transmission coefficient $t_s(\theta)=1+r_s(\theta)$. Then,

This corresponds to the transmission problem of an electromagnetic wave which electric field is parallel to the invariance direction $y$, between air and a nonmagnetic medium of relative permittivity $\varepsilon$.

rs=@(k,ep) (qf(k,1)-qf(k,ep))./(qf(k,1)+qf(k,ep));
ts=@(k,ep) 1+rs(k,ep);
BS0=@(k0,ep) +qf(k0,1).*rs(k0,ep);
BS1=@(k,k0,ep) 0.5*1i*K0^2*(ep-1).*ts(k,ep).*ts(k0,ep);

The transmission problem for P-polarized electromagnetic waves

Now considering an electromagnetic wave wich electric field is in the incidence plane, that is perpendicular to the invariance direction $y$, it is the quantity $\frac{1}{\varepsilon}\mathbf{\hat n\cdot grad\,}\psi$ that is continuous besides filed $\psi$. When passing from air (wavenumber K_0) to the scattering medium (wavenumber $K_0\sqrt{\varepsilon}$), we first define reflection $r_p$ coefficient for the plane : $$ r_p(\theta)=\frac{\varepsilon\cos\theta-\sqrt{\varepsilon-\sin^2\theta}}{\varepsilon\cos\theta-\sqrt{\varepsilon-\sin^2\theta}} $$ and associated transmission coefficient $t_p(\theta)=1+r_p(\theta)$. Then,

This corresponds to the transmission problem of an electromagnetic wave which electric field is parallel to the invariance direction $y$, between air and a nonmagnetic medium of relative permittivity $\varepsilon$.

rp=@(k,ep) (ep*qf(k,1)-qf(k,ep))./(ep*qf(k,1)+qf(k,ep));
tp=@(k,ep) 1+rp(k,ep);
BP0=@(k0,ep) +qf(k0,1).*rp(k0,ep);
BP1=@(k,k0,ep) 0.5*1i*((ep-1)*(ep*k.*k0-qf(k,ep).*qf(k0,ep))/ep^2).*tp(k,ep).*tp(k0,ep);

Stationary random rough surface with normal height distribution and Gaussian roughness spectrum

The surface $z=\eta(x)$ is a stationary centered random process with normally-distributed height. The roughess spectrum $\tilde\rho$, Fourier transform of the two-points correlation function $\langle \eta(x_1)\eta(x_2)\rangle=\rho(x_2-x_1)$, with brackets denoting statistical averaging, is a Gaussian function with two parameters $h$ and $\ell$ the height root mean quare and correlation length, respectively. $$ \tilde\rho(k)=\frac{h^2\ell}{2\sqrt{\pi}}e^{-\frac{k^2\ell^2}{4}} $$

h=lambda/50;
l=lambda/2;
spec=@(k) (0.5*h^2*l/sqrt(pi))*exp(-0.25*k.^2*l^2);

*** THE COHERENT REFLEXION COEFFICIENT ***

Stationarity implies that the mean field writes in air above the surface region : $$ \langle\psi(x,z)\rangle=e^{i(k_0x-q_0z)}+r_c(\theta_i)e^{i(k_0x+q_0z)} $$ as the sum of the incident field and a specularly-reflected plane wave. The amplitude of that reflection defines the coherent reflexion coefficient. It is related to the scattering amplitude by : $$ \frac{\langle s(k,k_0) \rangle}{q_0}=r_c(\theta_i)\delta(k-k_0) $$ in the general case, and applying the third-order SPM expansion to Gaussian process, $$ r_c(\theta_i)=\frac{B_0(k_0)}{q_0}+\int \frac{B_2(k_0,x+k_0,k_0)}{q_0}\tilde\rho(x)dx+O(\epsilon^4) $$ since, due to the Gaussian moment theorem, first- and third-orders do not contribute.

xmax=6*(sqrt(2)/l); % integration bound for Gaussian spectrum
rc=@(k0,B0,B2) B0(k0)./qf(k0,1)+integral(@(x)B2(k0,k0+x,k0).*spec(x)./qf(k0,1),-xmax,xmax);

The Dirichlet boundary condition

tid=[0:35 37:89];
for ii=1:length(tid)
    spm2rcd(ii)=rc(K0*sind(tid(ii)),BD0,BD2);
    spm2rcn(ii)=rc(K0*sind(tid(ii)),BN0,BN2);
end
figure(1)
subplot(121)
[AX,H1,H2]=plotyy(tid,abs(spm2rcd),tid,angle(spm2rcd)*180/pi);grid on
set(AX(1),'FontSize',24,'XLim',[0 90],'XTick',[0:15:90])
set(AX(2),'FontSize',24,'XLim',[0 90],'XTick',[0:15:90])
set(H1,'LineWidth',2)
set(H2,'LineWidth',2)
xlabel('Incidence angle (deg)')
title('Coherent reflexion coefficient - Dirichlet')
subplot(122)
[AX,H1,H2]=plotyy(tid,abs(spm2rcn),tid,angle(spm2rcn)*180/pi);grid on
set(AX(1),'FontSize',24,'XLim',[0 90],'XTick',[0:15:90])
set(AX(2),'FontSize',24,'XLim',[0 90],'XTick',[0:15:90])
set(H1,'LineWidth',2)
set(H2,'LineWidth',2)
xlabel('Incidence angle (deg)')
title('Coherent reflexion coefficient - Neumann')
legend('Modulus','Argument (deg)','Location','Best')

The incoherent bistatic coefficient

With $\Delta s(k,k_0)=s(k,k_0)-\langle s(k,k_0) \rangle$ the incoherent scattered amplitude, the incoherent bistatic coefficent $\sigma(k,k_0)$ is defined by $$ \langle \Delta s^\ast(k,k_0)\Delta s(k',k_0) \rangle=\sigma(k,k_0)\delta(k-k') $$ For a Gaussian process, the SPM expansion leads to : $$ \sigma(k,k_0)=|B_1(k,k_0)|^2\tilde\rho(k-k_0)+O(\epsilon^4) $$ since, due to the Gaussian moment theorem, the third-order contribution vanishes.

*** THE DIFFERENTIAL REFLECTIVITY ***

In air, the time-averaged power density through an horizontal plane is proportional to $$ w(x,z)=\Im m\left(\psi^\ast(x,z)\frac{d\psi}{dz}(x,z)\right) $$ It writes $w_i=-q_0$ for the incident plane wave and $$ \langle w_d \rangle=q_0|r_c|^2+\Re e\int\frac{\sigma(k,k_0)}{q}dk=q_0|r_c|^2+\int_{-\pi/2}^{+\pi/2}\sigma(k,k_0)d\theta_d $$ for the statistically averaged scattered field. Note that only the propagative part of the scattered field contributes. The reflectivity $R=-\langle w_d \rangle/w_i$ thus writes : $$ R(\theta_i)=|r_c(\theta_i)|^2+\int_{-\pi/2}^{+\pi/2}r_d(\theta_d,\theta_i)d\theta_d $$ introducing $r_d(\theta_d,\theta_i)=\sigma(k,k_0)/q_0$ the (incoherent) differential reflectivity : $$ \frac{\langle \Delta s^\ast(k,k_0)\Delta s(k',k_0) \rangle}{q_0}=r_d(\theta_d,\theta_i)\delta(k-k') $$

rd=@(k,k0,B1) abs(B1(k,k0)).^2.*spec(k-k0)./qf(k0,1);

Monostatic diagrams

tid=[0:89];
k0=K0*sind(tid);
k=-k0; % backscattering configuration
spm1rdd=rd(k,k0,BD1); % Dirichlet
spm1rdn=rd(k,k0,BN1); % Neumann
ep=2.25; % dielectric constant of glass at optical frequencies
spm1rds=rd(k,k0,@(k,k0)BS1(k,k0,ep)); % S-pol.
spm1rdp=rd(k,k0,@(k,k0)BP1(k,k0,ep)); % P-pol.
dB=@(x) 10*log10(x);
figure(2)
H1=plot(tid,dB(spm1rdd),tid,dB(spm1rdn),tid,dB(spm1rds),tid,dB(spm1rdp));grid on
set(gca,'FontSize',24,'XLim',[0 90],'XTick',[0:15:90])
set(H1,'LineWidth',2)
xlabel('Monostatic angle (deg)')
title('Differential reflectivity - dB')
legend('Dirichlet','Neumann','S-pol.','P-pol.','Location','Best')

Bistatic diagrams

ep=2.25; % dielectric constant of glass at optical frequencies
tid=atand(sqrt(ep)); % The BREWSTER angle
fprintf('The Brewster angle is %gdeg for index %g\n',tid,sqrt(ep))
k0=K0*sind(tid);
tsd=[-89:89];
k =K0*sind(tsd);
spm1rdd=rd(k,k0,BD1); % Dirichlet
spm1rdn=rd(k,k0,BN1); % Neumann
spm1rds=rd(k,k0,@(k,k0)BS1(k,k0,ep)); % S-pol.
spm1rdp=rd(k,k0,@(k,k0)BP1(k,k0,ep)); % P-pol.
dB=@(x) 10*log10(x);
figure(3)
H1=plot(tsd,dB(spm1rdd),tsd,dB(spm1rdn),tsd,dB(spm1rds),tsd,dB(spm1rdp));grid on
set(gca,'FontSize',24,'XLim',[-90 90],'XTick',[-90:15:90])
set(H1,'LineWidth',2)
xlabel('Scattering angle (deg)')
title('Differential reflectivity - dB')
legend('Dirichlet','Neumann','S-pol.','P-pol.','Location','Best')

The Brewster angle is 56.3099deg for index 1.5

Summary - 2D scattering

The rough surface reflectivity decomposes as $$ R(\theta_i)=|r_c(\theta_i)|^2+\int_{-\pi/2}^{+\pi/2}r_d(\theta_d,\theta_i)d\theta_d $$ with $r_c(\theta_i)$ the coherent reflexion coefficient $$ \frac{\langle s(k,k_0) \rangle}{q_0}=r_c(\theta_i)\delta(k-k_0) $$ and $r_d(\theta_d,\theta_i)$ the (incoherent) differential reflectivity $$ \frac{\langle \Delta s^\ast(k,k_0)\Delta s(k',k_0) \rangle}{q_0}=r_d(\theta_d,\theta_i)\delta(k-k') $$ both obtained from the scattering amplitude $s(k,k_0)$, itself defined through $$ \psi(x,z>\max\eta)=e^{i(k_0x-q_0z)}+\int\frac{s(k,k_0)}{q}e^{i(kx+qz)}dk $$


Gabriel Soriano