Terahertz Emission Modeling of Lunar Regolith

Abstract

We investigate the terahertz (THz) scattering and emission properties of lunar regolith by modeling it as a random medium with rough top and bottom boundaries and a host medium situated beneath. The total scattering and emission arise from three sources: the rough boundaries, the volume, and the interactions between the boundaries and the volume. To account for these sources, we model their respective phase matrices and apply the matrix doubling approach to couple these phase matrices to compute the total emission. The model is then used to explore insights into lunar regolith scattering and emission processes. The simulations reveal that surface roughness is the primary contributor to total scattering, while dielectric contrasts between the volume and the boundaries dominate total emission. The THz emissivity is highly sensitive to the regolith dielectric constant, particularly its imaginary part, making it a promising alternative for identifying previously undetected water ice in the lunar polar regions. The THz emissivity model developed in this study can be readily applied to invert the surface parameters of lunar regolith using THz observations.

1. Introduction

The Moon has been extensively explored over the decades, with missions including flybys, orbiters, landers, and impactors. Early laboratory analyses of lunar samples returned from the Apollo missions [1] have revealed the geology and formation of the Moon. Recent re-examinations of lunar samples returned from the Apollo missions have implied that water may be encased in glass beads of lunar regolith [2]. Comparative analyses of lunar samples returned by the latest Chang’E-5 mission have provided direct evidence of water content in lunar minerals [3,4]. In parallel, technological advancements have progressively revolutionized our understanding of the Moon through remote sensing techniques. Accordingly, X-ray and Gamma-ray spectrometers [5,6,7,8] have characterized the chemistry of the lunar surface and demonstrated global maps of the elemental composition of the lunar surface. Laser altimeters (Kaguya LALT, 2007 and NASA LOLA, 2009) [9,10] have created digital elevation models of the lunar topography, while terrain cameras (Kaguya TC, 2007 and NASA LROC, 2009) [11,12] have refined the geographic features of the lunar surface. Infrared spectrometers (NASA Cassini, 1997; Chandrayaan-1 M3, 2008; NASA Deep Impact, 2010; NASA/DLR SOFIA, 2021) [13,14,15,16,17] have discovered a widespread distribution of molecules on the lunar surface, including the detection of water molecules $\mathrm{O}\mathrm{H}/{\mathrm{H}}_2\mathrm{O}$ in the sunlit areas [13,14,15]. The impacting campaigns of the Lunar Reconnaissance Orbiter–Lyman Alpha Mapping Project (LRO-LAMP, 2009) and Lunar Crater Observation and Sensing Satellite (LCROSS, 2009) have provided the ground truth, including the detection of water frost in the crater Cabeus near the lunar south pole by impacting and simultaneously observing the ejected plume [18,19]. Radars, including the Arecibo radar telescope (1963), the Clementine bistatic radar (1994), the Chandrayaan-1 Mini-SAR (2008), and the LRO Mini-RF (2009), have detected subsurface structures on the Moon that suggest the possible presence of water ice deposits in the lunar regolith; however, these findings remain controversial [20,21,22,23]. Passive radiometers, such as Chang’E-1 (2007), NASA Diviner (2009), and Chang’E-2 (2010), have mapped the global brightness temperature on and under the lunar surface, identifying cold traps and potential water ice deposits in the lunar poles [24,25,26,27]. Despite the extensive knowledge gained from past missions about the Moon’s physical, chemical, mineralogical, and geographic characteristics, significant ambiguities remain concerning lunar regolith, water ice, and rare minerals. Consequently, the design of current and upcoming missions is focused on resolving these uncertainties and uncovering new insights.

A new phase of space exploration is underway, with a primary focus on the Moon and Cislunar space. Leading this era are prominent missions such as NASA’s Artemis program, CNSA’s Chang’E mission, ISRO’s Chandrayaan mission, and JAXA’s LUPEX mission, all aimed at searching for water ice and other resources on the Moon to support future human and robotic activities. To accurately locate and map the distribution of water ice on the Moon, comprehensive surface and subsurface observations are essential. In this pursuit, the passive radiometer [28,29,30,31] has been recognized as a powerful orbital instrument for investigating the geothermal and electrical properties of the Moon and other planetary bodies in the solar system. Notably, the Chang’E-1/-2 Microwave Radiometer (MRM) effectively monitored the lunar surface and subsurface brightness temperature using a multi-frequency (3 GHz∼37 GHz) microwave radiometer. Similarly, NASA’s Lunar Reconnaissance Orbiter, equipped with the Diviner Lunar Radiometer Experiment (DLRE), successfully measured the lunar surface brightness temperature using a multi-channel infrared (7.8 $\mathsf{\mu }$m∼400 $\mathsf{\mu }$m) radiometer. Despite these measurements, a significant probing gap remains in lunar exploration [32] within the terahertz (THz) band [33], which occupies a relatively untapped portion of the electromagnetic spectrum, lying in the transition region between far-infrared waves and microwaves. To explore the feasibility of this concept, a radiative transfer model is being developed specifically for THz wavelengths, aiming to provide comprehensive insights into the THz wave brightness temperature from the lunar surface and subsurface.

Utilizing THz waves for remote sensing holds great promise for detecting and measuring ever-finer details of the emitted and reflected waves that contain the “signatures” of the medium with which they interact, opening new windows of opportunity for planetary surface remote sensing [34]. THz waves can frequently and continuously monitor the lunar surface under all weather conditions. In addition, THz waves have the capability to sense deeper than infrared waves, allowing for high recognizability of the lunar subsurface. Moreover, THz waves have higher spatial resolution than microwaves, providing a more refined structuration of the lunar surface. These unique properties of THz waves are crucial in unveiling the mysteries of the physical-chemical composition and geometrical volume properties of the lunar regolith. This includes identifying new craters and regolith overturn [35], quantifying ilmenite abundance [36], and locating water ice in permanently shadowed regions—the primary objectives in lunar exploration. THz radiative transfer models are pivotal in understanding the physical mechanisms of interactions between electromagnetic waves and the medium, providing essential guidance for THz sensor design and aiding the development of algorithms to decipher the physical parameters of interest. For atmospheric remote sensing, numerous radiative transfer models have been developed, including the Atmospheric Radiative Transfer Simulator (ARTS) [37,38,39], Atmospheric TeraHertz Radiation Analysis and Simulation (AMATERASU) [40,41], and Earth Observing System Microwave Limb Sounder (EOSMLS) [42]. These models, created by various research groups, are used to observe Earth and planetary atmospheres in the THz frequency region. They serve as practical tools for simulating and analyzing electromagnetic radiation interactions with the atmosphere, studying atmospheric composition, cloud properties, and the planetary radiation budget. To the best of our knowledge, a THz radiative transfer model specifically for planetary surfaces, particularly the Moon, does not yet exist. The lunar regolith consists of granular particles with diameters ranging in the hundreds of microns [1], comparable to THz wavelengths, causing rough surface scattering due to boundary discontinuities. In addition, inhomogeneities [43] within the lunar regolith result in volume scattering, further complicating the THz signature and posing challenges in accurately estimating the physical parameters of interest. Modeling such a lunar regolith requires an appropriate combination of surface scattering and volume scattering methods. In this article, we present a theoretical model that accounts for surface scattering, volume scattering, and boundary interactions to investigate THz radiative properties of lunar regolith.

This article is structured as follows: Section 2 briefly introduces the THz remote sensing of the planetary surface, which is modeled, in the context of radiative transfer, as an inhomogeneous random volume bounded by randomly rough boundaries hosted by a background half-space, including the parameterization of physical properties such as surface roughness, dielectric properties, bulk density and porosity, and physical temperature. Section 3 presents the model development and the incorporation of the respective surface and volume phase matrices. Section 4 presents the model simulations to explore the THz radiative properties of lunar regolith. Section 5 summarizes the key findings and provides future outlooks.

2. Lunar Regolith as an Emission Layer

Structurally and dielectrically, we model the lunar regolith as an inhomogeneous layer of random particles, sitting on a host medium, as illustrated in Figure 1.

The inhomogeneous layer has a physical temperature and an effective dielectric constant, which is the average of the dielectric constants of the scatterers and the background medium. A lunar regolith layer will contain pebble and rock particles as scatterers and under-bed rock as the host medium.

In describing an inhomogeneous layer, we assume an optical depth of the layer, denoted as $\tau$. The extinction coefficient, $k_e$, is the sum of the absorption coefficient, $k_a$, and the scattering coefficient, $k_s$, due to linear processes, given by $k_e=k_a+k_s$.

The absorption coefficient is determined by the free-space wavenumber and dielectric property:

$$k_a=2k\Im \left(\sqrt{{ϵ}_l}\right),$$

(1)

with ℑ denoting the imaginary part and ${ϵ}_l$ being the complex dielectric constant of the lunar regolith. Here, k is the wavenumber in free-space.

The albedo is defined as:

$${\omega }_0={\displaystyle \frac{k_s}{k_a+k_s}}={\displaystyle \frac{k_s}{k_e}},$$

(2)

thus, $0\le {\omega }_0\le 1$.

The optical depth is given by:

$$\tau ={\int }_z^0\left(k_a+k_s\right)dz={\int }_z^0{k}_{e}dz.$$

(3)

2.1. Roughness Description

The lunar surface roughness characterizes its surface features and serves as one of the predominant variables determining radar scattering [44], thermal emission [45], optical reflection [46], and spectral radiance [47]. Therefore, it is essential to characterize the lunar surface as rough terrain, especially at THz frequencies.

For a given geologic unit of the lunar surface [48], a power-law roughness spectrum [49] can be generalized to the form:

$$W(K_x,K_y)={\displaystyle \sum _i}{\displaystyle \frac{2s_i^2{l}_{x,i}{l}_{y,i}\mathrm{\Gamma }({\alpha }_{x,i}+1/2)\mathrm{\Gamma }({\alpha }_{y,i}+1/2)}{\mathrm{\Gamma }\left({\alpha }_{x,i}\right)\mathrm{\Gamma }\left({\alpha }_{y,i}\right)}}\times {\displaystyle \frac{1}{{\left[1+{\left(K_x{l}_{x,i}\right)}^2\right]}^{{\alpha }_{x,i}+1/2}{\left[1+{\left(K_y{l}_{y,i}\right)}^2\right]}^{{\alpha }_{y,i}+1/2}}},$$

(4)

where $s_i$ denotes the surface RMS height for the i-th scale, and $l_{x,i}$ and $l_{y,i}$ are the correlation lengths along the x and y directions, respectively, for the i-th scale. The ratio between the RMS height and correlation length is correlated to the RMS slope. The scaling factors ${\alpha }_{x,i}$ and ${\alpha }_{y,i}$ control the power-law decaying rates in the x and y directions, respectively. The $\mathrm{\Gamma }(·)$ represents the Gamma function, and $(K_x,K_y)$ indicates the wavenumber domain. This model incorporates anisotropy and multiscale properties, providing a more accurate representation of the lunar surface roughness. The parameters $s_i$, $l_{x,i}$, $l_{y,i}$, ${\alpha }_{x,i}$, and ${\alpha }_{y,i}$ can be adjusted to capture the roughness characteristics across different scales, from fine to coarse features. The corresponding correlation function and roughness spectrum indexes cover a range of surfaces from Gaussian to exponentially correlated by adjusting these parameters. Figure 2 shows the correlation function and the corresponding roughness spectrum for different power indexes, illustrating the transition between Gaussian and exponentially correlated surfaces by adjusting the power index ${\alpha }_x$ and ${\alpha }_y$.

2.2. Dielectric Properties

The THz electromagnetic response of lunar regolith is determined by its dielectric properties. In general, the complex dielectric constant (also known as permittivity) can be expressed as:

$$\begin{array}{cc}\hfill {ϵ}_r\left(f,z\right)& ={ϵ}_r^{\prime }\left(f,z\right)-j{ϵ}_r^{″}\left(f,z\right)\hfill \\ \hfill & ={ϵ}_r^{\prime }\left(z\right)-j\left[tan\delta \left(f,z\right)\right]{ϵ}_r^{\prime }\left(z\right),\hfill \end{array}$$

(5)

where ${ϵ}_r^{\prime }$ and ${ϵ}_r^{″}$ denote the real and imaginary parts of the complex dielectric constant, respectively, and $tan\delta$ indicates the loss tangent as a function of depth z and frequency f.

According to the experimental analyses of the Apollo lunar returned samples reported in the Lunar Sourcebook [1], the relative dielectric constant of lunar regolith is density-dependent [1,50], and is expressed by:

$${ϵ}^{\prime }\left(z\right)=1.{919}^{\rho \left(z\right)},$$

(6)

where $\rho \left(z\right)$ represents the density at depth z in ${\mathrm{gcm}}^{-3}$.

The Apollo density profile model is given by:

$$\rho \left(z\right)=1.92{\displaystyle \frac{z+12.2}{z+18}},$$

(7)

The porosity of lunar regolith is interrelated as:

$$n\left(z\right)=1-{\displaystyle \frac{\rho \left(z\right)}{G{\rho }_w}},$$

(8)

where ${\rho }_w$ is the density of water at 4 °C and 1 ${\mathrm{gcm}}^{-3}$, and G is the specific gravity.

The loss tangent is related to both bulk density and the abundance S of the mineralogical composition (i.e., $\mathrm{F}\mathrm{e}\mathrm{O}+\mathrm{T}\mathrm{i}{\mathrm{O}}_2$) of the lunar regolith, as given by:

$$\tan \delta \left(z\right)={10}^{0.038S+0.312\rho \left(z\right)-3.26},$$

(9)

where S is in the range of 0∼30% [1].

The dielectric profiles of lunar regolith based on the experimental analyses of the Apollo returned sample are presented in Figure 3. In a recent study [43], a comparison was made between the density profile of Chang’E-4 and the Apollo density profile model. The results reveal that the density profile of Chang’E-4 is consistent with the Apollo curve at deeper regolith levels, while smaller than the Apollo model prediction at shallower depths, suggesting an inhomogeneous density distribution of the lunar regolith. Hence, it is essential to consider this inhomogeneity as it results in volume scattering within lunar regolith and can affect the accuracy of electromagnetic characteristics predictions. In addition, a reanalysis of the Apollo lunar returned samples by [51] found that the loss tangent may primarily relate to $\mathrm{T}\mathrm{i}{\mathrm{O}}_2$, rather than $\mathrm{F}\mathrm{e}\mathrm{O}$. Moreover, studies on the dielectric properties of lunar regolith, such as those by [52,53,54], have found that the loss tangent is also frequency-dependent.

In addition to dielectric inhomogeneity, variations in physical temperature can also induce volume scattering within lunar regolith.

The variation in physical temperature of lunar regolith with depth z and time t is modeled using the one-dimensional thermal diffusion equation [29,55,56,57]:

$$\begin{array}{c}\hfill \rho c{\displaystyle \frac{\partial T\left(z,t\right)}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(K{\displaystyle \frac{\partial T\left(z,t\right)}{\partial z}}\right),\end{array}$$

(10)

where T is the physical temperature; $\rho$ is the bulk density; c is the specific heat capacity; and K is the thermal conductivity.

The density profile as a function of depth is expressed as:

$$\begin{array}{c}\hfill \rho \left(z\right)={\rho }_d-\left({\rho }_d-{\rho }_s\right)e^{-z/H},\end{array}$$

(11)

where H denotes a location-dependent scaling factor that controls how density increases with depth. ${\rho }_s$ and ${\rho }_d$ are bounded densities [57] at the surface and at a certain depth, with ${\rho }_s=1.1{\mathrm{gcm}}^{-3}$ and ${\rho }_d=1.8{\mathrm{gcm}}^{-3}$, respectively.

The thermal conductivity varies with several factors such as composition, density, and temperature, which can be described by:

$$\begin{array}{c}\hfill K\left(T,\rho \right)=K_c\left(\rho \right)\left[1+\chi {(T/350)}^3\right],\end{array}$$

(12)

where $\chi$ is a dimensionless parameter related to the radiative component of thermal conductivity and $K_c$ is contact conductivity, which is assumed to be linearly proportional to density as:

$$\begin{array}{c}\hfill K_c\left(\rho \right)=K_d-\left(K_d-K_s\right){\displaystyle \frac{{\rho }_d-\rho }{{\rho }_d-{\rho }_s}},\end{array}$$

(13)

where $K_s$ and $K_d$ are contact conductivity values at the surface and at a certain depth, respectively.

In Figure 4, the profiles of the minimum, average, and maximum temperatures of lunar regolith, denoted as $T_i$, are displayed.

3. Physics-Based Emission Modeling

In passive remote sensing and under the context of emission, we identify the following contributive sources of emission:

(i)

Emissions from bounded inhomogeneous layer $T_l^{+}$, $T_l^{-}$;

(ii)

Emission from the host medium $T_h$;

(iii)

Emission from layer–boundary interactions $T_{l\leftrightarrow ls}^{\pm }$.

The emission from these sources undergoes multiple scattering inside the lunar regolith layer. The free space–regolith layer and regolith–host medium interfaces establish boundary interactions, which are included in each effective emission source. The emission from the host medium can be determined by employing Rayleigh–Jean’s law:

$$T_h=e_{host}\left({\displaystyle \frac{K_B{T}_{host}}{{\lambda }^2}}\right)\left[{\displaystyle \frac{\mu ϵ}{{\mu }_0{ϵ}_0}}\right],$$

(14)

where $e_{host}$ is the emissivity of the host medium, $K_B$ is the Boltzmann constant, $T_{host}$ is the physical temperature of the host medium, $\mu$ and $ϵ$ are the permeability and permittivity of the host medium, respectively, ${\mu }_0$ and ${ϵ}_0$ are the permeability and permittivity of the free space, respectively, and $\lambda$ is the wavelength.

The emission from sources (i), (ii), and (iii) undergoes multiple scattering inside the regolith layer and the free space–regolith layer and regolith–host medium interfaces; an interaction factor, which represents these processes, should be added for each effective emission source. The total emission can then be expressed as:

$$\mathbf{e}={\mathbf{L}}_{\mathbf{u}}{\mathbf{e}}_{\mathbf{u}}+{\mathbf{L}}_{\mathbf{d}}{\mathbf{e}}_{\mathbf{d}}+{\mathbf{L}}_{\mathbf{h}}{\mathbf{e}}_{\mathbf{h}},$$

(15)

where ${\mathbf{L}}_{\mathbf{u}}$, ${\mathbf{L}}_{\mathbf{d}}$, and ${\mathbf{L}}_{\mathbf{h}}$ denote the multiple interaction factors of upward, downward, and host medium emission sources, respectively.

To solve the total emission in Equation (15), we apply the matrix doubling method [58]. In matrix doubling, the emission source in the regolith layer, including the upward and downward emission source, is obtained by integrating over the thickness of the regolith layer, except the host medium itself. For ease of reading, we detail the concept of matrix doubling in the Appendix A.

The multiple operator ${\mathbf{L}}_{\mathbf{u}}$ of the upward effective emission source $T_l^{+}$ is expressed as:

$${\mathbf{L}}_{\mathbf{u}}={\mathbf{Q}}_{\mathbf{la}}{(\mathbf{I}-{\mathbf{T}}^{*}{\mathbf{R}}_{\mathbf{lh}}\mathbf{T}{\mathbf{R}}_{\mathbf{la}})}^{-1},$$

(16)

where $\mathbf{T}$ and ${\mathbf{T}}^{*}$ denote the forward and backward scattering matrices of the layer; ${\mathbf{Q}}_{\mathbf{la}}$ and ${\mathbf{R}}_{\mathbf{la}}$ represent the transmission and reflection matrices of the layer–free space interface; ${\mathbf{Q}}_{\mathbf{hl}}$ is the transmission matrix of the host–layer interface and ${\mathbf{R}}_{\mathbf{lh}}$ is the reflection matrix of the layer–host interface; and $\mathbf{I}$ denotes the identity matrix.

For the downward emission source $T_l^{-}$ and host emission source $T_h$, the multiple operators ${\mathbf{L}}_{\mathbf{d}}$ and ${\mathbf{L}}_{\mathbf{s}}$ are given, respectively, by the following:

$${\mathbf{L}}_{\mathbf{d}}={\mathbf{Q}}_{\mathbf{la}}{(\mathbf{I}-{\mathbf{T}}^{*}{\mathbf{R}}_{\mathbf{lh}}\mathbf{T}{\mathbf{R}}_{\mathbf{la}})}^{-1}{\mathbf{T}}^{*}{\mathbf{R}}_{\mathbf{lh}},$$

(17)

$${\mathbf{L}}_{\mathbf{h}}={\mathbf{Q}}_{\mathbf{la}}{(\mathbf{I}-{\mathbf{T}}^{*}{\mathbf{R}}_{\mathbf{lh}}\mathbf{T}{\mathbf{R}}_{\mathbf{la}})}^{-1}{\mathbf{T}}^{*}{\mathbf{Q}}_{\mathbf{hl}}.$$

(18)

The free-space layer and the layer–host interfaces may be rough; thus, the reflection and transmission matrices cannot simply be derived from the Fresnel reflectivity but must account for the roughness effect. In doing so, we adopt the AIEM model, as will be detailed in the following subsection.

3.1. Phase Matrices of Rough Boundaries

The reflection phase matrix is given by:

$$\begin{array}{cc}\hfill \langle S_{qp}{S}_{rs}^{*}\rangle & ={\displaystyle \frac{k^{2}A}{8\pi }}\mathbf{exp}[-s^2({k_z}^2+k_{sz}^2)]\hfill \\ \hfill & {\displaystyle \sum _{n-1}^{\infty }}{s}^{2n}\left(I_{qp}^n{I}_{rs}^{n*}\right){\displaystyle \frac{{\mathbf{W}}^n(k_{sx}-k_x,k_{sy}-k_y)}{n!}},\hfill \end{array}$$

(19)

where A is the illuminated area and $I_{qp}^n$ is the field coefficient, as detailed in [59,60]. The subscript $pq$ represents polarization, with $p=h$ and $q=v$ for horizontal and vertical polarization, respectively. s is the surface RMS height. k is the wavenumber, with $k_{sx}=ksin{\theta }_{s}cos{\varphi }_s$, $k_{sy}=ksin{\theta }_{s}sin{\varphi }_s$, $k_z=kcos\theta$, and $k_{sz}=kcos{\theta }_s$.

The transmission phase matrix is expressed as follows:

$$\begin{array}{cc}\hfill \langle S_{qp}{S}_{rs}^{*}\rangle & ={\displaystyle \frac{k^{2}A}{8\pi {\eta }_r}}\mathbf{exp}[-s^2({k_z}^2+k_{tz}^2)]\hfill \\ \hfill & {\displaystyle \sum _{n-1}^{\infty }}{s}^{2n}\left(I_{qp}^n{I}_{rs}^{n*}\right){\displaystyle \frac{{\mathbf{W}}^n(k_{tx}-k_x,k_{ty}-k_y)}{n!}},\hfill \end{array}$$

(20)

where ${\eta }_r={\eta }_t/\eta$ represents the ratio of the intrinsic impedance of the lower medium to that of the upper medium.

3.2. Phase Matrix of Regolith Layer

Volume scattering, which arises from distributed scatterers within the regolith layer, depends on the size, shape, orientation, and dielectric properties of the scatterers. In modeling volume scattering from a collection of small particles at THz frequencies, the dense medium phase and amplitude correction theory (DM-PACT) based on Mie theory [61] is adopted as follows:

$$\begin{array}{cc}\hfill {\sigma }_{pq}^0\left(v\right)& =4\pi {\mu }_s\left[{\displaystyle \frac{{\mathbf{I}}^s\left({\theta }_s,{\varphi }_s\right)}{{\mathbf{I}}^i\left({\theta }_i,{\varphi }_i\right)}}\right]\hfill \\ \hfill & ={\mu }_{s}T\left({\theta }_s,{\theta }_0\right){\omega }_0\left\{1-exp\left[-{\displaystyle \frac{\tau \left({\mu }_0+{\mu }_t\right)}{{\mu }_0{\mu }_t}}\right]\right\}\hfill \\ \hfill & \times {\displaystyle \frac{{\mathbf{P}}_v}{{\mu }_0+{\mu }_t}}{\mu }_{t}T\left({\theta }_t,{\theta }_i\right),\hfill \end{array}$$

(21)

where ${\mu }_j$ is the directional cosine, defined as ${\mu }_j=cos{\theta }_j$, with $j=0,t,s$, calculated by Snell’s law. $T({\theta }_i,{\theta }_j)$ is the power Fresnel transmission coefficient.

The phase function based on Mie theory [62] is given by:

$$\begin{array}{c}\hfill \overline{S}={\displaystyle \frac{4\pi r^2\eta }{k_{vs}{\left|E_0\right|}^2}}\Re \left[\begin{array}{cc}{\left(E_v^s{H}_h^{s*}\right)}_{v-inc}& {\left(E_v^s{H}_h^{s*}\right)}_{h-inc}\\ -{\left(E_v^s{H}_h^{s*}\right)}_{v-inc}& -{\left(E_v^s{H}_h^{s*}\right)}_{h-inc}\end{array}\right],\end{array}$$

(22)

where $\eta =\sqrt{\mu /ϵ}$; ℜ stands for the real part. $E_p^s$ and $H_p^s$ denote the scattered electric and magnetic fields, respectively, in p-polarization due to the scatterer when the incident electric field is $E_0$. $v-inc$ and $h-inc$ represent the vertical and horizontal incident wave polarizations, respectively. $k_{vs}$ is the volume scattering coefficient. r is related to the volume fraction $v_f$ of the spherical scatterer as written by:

$$r={\left(v_0/v_f\right)}^{1/3},$$

(23)

where $v_0=4\pi a_0^3/3$ denotes the volume of a sphere; $a_0$ is the radius of a sphere.

We adopt the phase correction factor with Gaussian correlation to eliminate the conventional assumptions of uncorrelated scatterers and far-field interactions for a collection of closely packed scatterers [63]. The expression of the phase correction factor is given by:

$$\begin{array}{cc}\hfill {\langle {\left|\mathrm{\Phi }\right|}^2\rangle }_n& ={\displaystyle \frac{1-e^{-k_{si}^2{h}_s^2}}{{d_s}^3}}+{\displaystyle \frac{e^{-k_{si}^2{h}_s^2}}{{d_s}^3}}{\displaystyle \sum _{q=1}^{\infty }}{\displaystyle \frac{{\left(k_{si}^2{h}_s^2\right)}^q}{q!}}\hfill \\ \hfill & ·\left[{\left(\sqrt{{\displaystyle \frac{\pi }{q}}}\left({\displaystyle \frac{w_s}{d_s}}\right)\right)}^{3}exp\left({\displaystyle \frac{-k_{si}^2{h}_s^2}{4q}}\right)-a\left(k_x\right)a\left(k_y\right)a\left(k_z\right)\right],\hfill \end{array}$$

(24)

where $h_s$ is the standard deviation of the scatterer position deviation from its deterministic position in the scatterer array; $d_s$ is the mean distance between scatterers; and $w_s$ is the correlation length of one scatterer position deviation from its deterministic position to the other scatterer position deviation in the scatterer array. The term $a\left(k_r\right)$ is defined as:

$$\begin{array}{cc}\hfill a\left(k_r\right)& =\left[\left(\sqrt{{\displaystyle \frac{\pi }{q}}}\left({\displaystyle \frac{w_s}{d_s}}\right)\right)exp\left({\displaystyle \frac{-k_r^2{w}_s^2}{4q}}\right)\Re \left[\mathrm{erf}\left({\displaystyle \frac{q{d}_s/w_s+j{k}_r{w}_s}{2\sqrt{q}}}\right)\right]\right],\hfill \end{array}$$

(25)

where $h_s$, $w_s$, and $d_s$ are determined by the geometric size and volume fraction of the scatterer. $\mathrm{erf}$() is the error function. The vectors $k_{si}=|{\overrightarrow{k}}_s-{\overrightarrow{k}}_i|$, ${\overrightarrow{k}}_s=k(sin{\theta }_{s}cos{\varphi }_s\widehat{x}$$+sin{\theta }_{s}sin{\varphi }_s\widehat{y}+cos{\theta }_s\widehat{z})$, and ${\overrightarrow{k}}_i=k\left(sin{\theta }_{i}cos{\varphi }_i\widehat{x}+sin{\theta }_{i}sin{\varphi }_i\widehat{y}+cos{\theta }_i\widehat{z}\right)$ are the scattered and incident propagation vectors.

Hence, the modified volume scattering phase function of a collection of scatterers can be reconstructed as:

$${\mathbf{P}}_v={\langle {\left|\mathrm{\Phi }\right|}^2\rangle }_n·\overline{S}.$$

(26)

4. Results and Discussion

In this section, we present numerical simulations for sets of lunar regolith conducted based on the developed model. The goal is to evaluate the scattering coefficients and emissivity in response to variations in the following lunar regolith parameters: dielectric constant, optical depth, and albedo. This analysis helps determine the relative contribution of these parameters to the total volume scattering that affects surface emission.

4.1. Dielectric Effect

Figure 5 displays the scattering coefficients obtained using a set of predefined regolith layer parameters at 480 GHz. The Mie scatterer parameters, including albedo ${\omega }_0$ and optical depth $\tau$, are 0.42 and 0.75, respectively. The background dielectric constant is assumed to be ${ϵ}_b=2.65-j0.008$, while the layer dielectric constant is varied as ${ϵ}_l=2.1-j0.002$ and ${ϵ}_l=2.3-j0.002$, with top boundary roughness $kl=12.12$, $ks=0.52$ and bottom boundary roughness $k{l}_b=7.1$, $k{s}_b=0.81$. Note that Figure 5 is based on the same set of parameters, except for the layer dielectric constant. Upon examining the results, we observe that a slight change in the dielectric constant from ${ϵ}_l=2.1-j0.002$ to ${ϵ}_l=2.3-j0.002$ leads to a higher but not significant drop-off rate of the angular scattering at both small and large angles. In addition, the polarization difference between VV and HH does not change over the range of angles between 10 and 70 degrees. Compared to pure surface scattering, the layer effect narrows the polarization difference and flattens the angular trend. Moreover, the contribution of volume scattering from the layer is more pronounced at larger incident angles, especially for HH polarization.

We plot the scattering coefficients versus incident angle by setting the identical regolith layer parameters as in Figure 5, except the background dielectric constant ${ϵ}_b=3.65-j0.008$. Figure 6 presents the scattering coefficients for ${ϵ}_l=3.0-j0.07$ and ${ϵ}_l=3.0-j0.03$. It is expected that such a tiny change in the attenuation factor would not cause a visible change in the scattering coefficient. However, at the THz spectrum, detecting such small changes in the dielectric constant of the regolith layer from the scattering strength is difficult, if not impossible. One notable feature in the model is the loss of sensitivity in volume scattering to the specific shape of the scatterers when they are randomly oriented in the lunar regolith layer. Thus, for randomly oriented scatterers, volume scattering is not sensitive to the specific shape of the scatterer.

Figure 7 shows two sets of simulated emissivity curves. The first set has layer dielectric values of ${ϵ}_l=1.1,1.3,1.5,$ and $2.0$, and the second set has ${ϵ}_l=2.5,3.0,3.5,$ and $4.0$. Both cases have a background permittivity of ${ϵ}_b=5-0.5j$ at 480 GHz, with a top boundary roughness of $kl=6.23$, $ks=0.42$, and bottom boundary roughness of $k{l}_b=8.38$, $k{s}_b=0.63$. The albedo and optical depth values are both 0.3. From the simulation results in Figure 7, it is evident that surface emission is highly sensitive to changes in the dielectric constant, which can be influenced by variations in soil moisture. Additionally, H polarization is more sensitive to dielectric changes between 40 and 60 degrees, while V polarization is sensitive to dielectric changes at small observation angles.

4.2. Roughness Effect

THz emission depends on various physical parameters, including surface height, correlation length, dielectric constant, albedo, and optical depth. We examine how these factors influence the emission from lunar regolith. Figure 8 illustrates the emission from lunar regolith with identical surface parameters but different correlation functions—exponential and Gaussian—at 480 GHz. The emissivity curves with Gaussian correlation are slightly lower than those with exponential correlation for all three normalized surface RMS heights, $ks=0.21,0.83,$ and $1.47$, with a surface correlation length $kl=8.37$ and a layer permittivity ${ϵ}_l=3-j0.01$. However, the difference is insignificant, and the curves are similar. At nadir observation, there is a slight increase in emission with an RMS height of $ks=1.47$, accompanied by a faster drop-off at larger observation angles. Although discernible, the changes resulting from variations in RMS height are relatively small.

Next, we examine the effects of correlation length on emission for Gaussian and exponential correlation functions. We compute emission with three different correlation lengths $kl=4.19,8.37,$ and $16.76$ at a surface RMS height of $ks=0.42$ and layer permittivity of ${ϵ}_l=3-j0.01$ at 480 GHz. The simulation results, presented in Figure 9, demonstrate that changes in correlation length, whether longer or shorter, have a negligible effect on emission for both surface correlation functions. In addition, changing the correlation function and correlation length is expected to be negligible. Note that, combined with Figure 8, the surface RMS height exerts more impact on emission than correlation length.

4.3. Coupling Effect

We examine the effects of optical depth on scattering, as illustrated in Figure 10. The simulation parameters are as follows: frequency $f=480$ GHz, albedo ${\omega }_0=0.26$, layer permittivity ${ϵ}_l=2.0-j0.01$, and background permittivity ${ϵ}_b=3.0-j0.01$. The surface parameters include top boundary roughness $kl=3.14$, $ks=0.52$ and bottom boundary roughness $k{l}_b=5.24$, $k{s}_b=0.62$. When the regolith layer is optically thick, it acts as a half-space. Hence, an expected saturation effect becomes evident when the optical depth $\tau$ approaches unity. By calculating scattering up to an optical depth $\tau$ of 1.6, it becomes clear that there is an increase in the scattering level with optical depth, as depicted in Figure 10. As $\tau$ increases, a saturation effect is observed, characterized by a reduction in the spacing between the scattering curves for both HH and VV polarizations. The surface boundary effect diminishes for larger dielectric values. Similar to an increase in optical depth, an increase in albedo also enhances the scattering level.

Figure 11 illustrates that as the albedo ${\omega }_0$ varies from 0.1 to 0.8 for HH and VV polarizations, there is an almost identical change in the scattering level of approximately 2.5 dB. The rough boundary surfaces are exponentially correlated with top boundary roughness $kl=3.14$, $ks=0.52$ and bottom boundary roughness $k{l}_b=4.19$, $k{s}_b=0.52$. The model parameters are as follows: $f=480$ GHz, $\tau$ = 0.5, layer permittivity ${ϵ}_l=2.0-j0.01$, and background permittivity ${ϵ}_b=4.0-j0.01$. Doubling the albedo results in a nearly uniform increase in the scattering level across the angular range of 10 to 80 degrees by approximately 2.5 dB. When the angle of incidence is 10 degrees, the scattering strength is similar for both VV and HH polarizations. However, as the angle of incidence increases, the HH polarization decreases faster than the VV polarization. Specifically, at an angle of 70 degrees, the VV polarization is about 1 to 1.5 dB higher than HH polarization.

Figure 12 shows that as the optical depth increases to 0.1, 0.5, 0.9, and 2.0 for albedo values ${\omega }_0$ of 0.1 and 0.2, emissions also increase for both H and V polarizations. The two cases have the same background permittivity of ${ϵ}_b=5-0.5j$ at 480 GHz, with a top boundary roughness of $kl=8.38$, $ks=0.63$, and a bottom boundary roughness of $k{l}_b=6.28$, $k{s}_b=0.42$. However, the increase in emission is more significant when the albedo is 0.1 compared to 0.2, which indicates that if the albedo increases further, the emission curves may converge to a narrow range, and the emission trend may reverse. Moreover, an increase in emission levels leads to decreased polarization differences.

The albedo is a measure of the scattering strength within an inhomogeneous layer, calculated by dividing the volume scattering coefficient by the extinction coefficient. A high albedo indicates low absorption or emission, as the extinction coefficient is the sum of absorption and volume scattering coefficients. As presented in Figure 13, with the model parameters as the optical depth $\tau =0.35$ and the background permittivity ${ϵ}_b=5-0.5j$ at 480 GHz, with a top boundary roughness of $kl=8.38$, $ks=0.63$, and a bottom boundary roughness of $k{l}_b=6.28$, $k{s}_b=0.42$, the increase in albedo results in a decrease in emission. An inhomogeneous layer lessens the discontinuity between the half-space below and the free space above it. It can significantly narrow the polarization difference by increasing the H polarization and slowing its decrease with the observation angle. Regarding the surface emission, the Brewster angle region always exhibits high emissivity for V polarization, regardless of the model parameters. However, adding the inhomogeneous layer, a dielectric cover over a continuous surface, reduces the Brewster angle effect. This reduction in the rise of V polarization towards the Brewster angle region is enhanced with an increase in albedo, as seen in Figure 13. Therefore, a potentially sharp rise in V emission around the Brewster region is flattened.

5. Conclusions

We have developed a THz emission model that provides a coherent solution for the scattering and emission from lunar regolith. The lunar regolith is modeled as a random medium bounded by irregular top and bottom boundaries and a host medium underneath. The surface scattering phase matrix of the irregular boundary has been constructed using the AIEM with the generalized power law spectrum. In addition, the volume scattering phase matrix of subsurface scatterers has been incorporated by employing the Mie theory, which accounts for the phase correction factor of the discrete random medium. The matrix doubling method has been adopted to solve the radiative transfer equations including surface scattering, volume scattering, and boundary interactions.

Actual thermal and dielectric properties of lunar regolith, from laboratory experiments on lunar samples returned from the Apollo and Chang’E missions, were used at THz as input to the emission model. The model provides physical insights on the scattering and emission processes of the lunar regolith parameters. Although the vertical profile of particle size distribution within lunar regolith is relatively unknown, numerical simulations suggest that surface roughness is the major contributor to the total scattering coefficient. In contrast, the dielectric contrasts between the layer and the boundaries are the dominant factor in determining the total emission. The surface roughness serves as an impedance matcher between the regolith and surface-leaving radiation, and its characteristics cannot be ignored in modeling the surface radiation. The dielectric constant, particularly the imaginary part, is highly sensitive to emissivity at THz, making it useful for detecting water ice in the lunar polar regions. In order to investigate the dynamic range of the scattering strength and emissivity, more thorough simulations must be conducted by varying regolith parameters at THz frequency. The relative importance of these parameters on the angular scattering and emission response must be determined quantitatively. We plan to validate the modeled THz surface properties of lunar regolith in our laboratory that is currently under investigation. Preliminary experimental findings related to the surface roughness of the simulated lunar regolith were presented in [64].

The future research efforts may be devoted to the following subjects:

(i)

The geometric shape of scatterers in the regolith layer were modeled as spheres in this model. In future work, the models should incorporate other types of scatterers over the entire frequency region of interest.

(ii)

A thorough examination of multiscale roughness on the lunar surface emissivity is necessary.

(iii)

It is important to treat the lunar regolith layer as a dense and random medium and compare the model results with experimental measurements.

(iv)

Retrieving geophysical parameters from surface emission is pivotal in the passive sensing of planetary surfaces. Since the emissivity lies within the $\left(0,1\right)$ range, even a small change in it leads to a significant difference in brightness temperature. Several factors, including surface RMS height, correlation length, dielectric constant, albedo, and optical depth, contribute to such slight variations in emission. Making full use of this slight variation is promising in deciphering the geophysical parameters of planetary surfaces.

(v)

The new THz emission model should be made readily applicable to other bodies with no atmosphere in the solar system, such as Mercury, Ceres, Vesta, Eros, and Phobos, as well as to bodies with atmosphere, such as comets (e.g., 67P/Churyumov– Gerasimenko) and Jupiter’s icy moons, like Ganymede.