Full waveﬁeld simulation versus measurement of microwave scattering by a complex 3D-printed asteroid analogue

Context. The small bodies of the Solar System, and especially their internal structures, are still not well-known. Studies of the interior of comets and asteroids could provide important information about their formation and also about the early Solar System. Aims. In this paper, we investigate the possibility of obtaining information about their inner structure from their response to an incident electromagnetic ﬁeld in preparation for future space radar missions. Our focus is on experimental measurements concerning two analog models with the shape of 25143 Itokawa, a small rubble pile asteroid monitored by the Japanese space agency’s (JAXA) Hayabusa mission in 2005. Methods. The analog models prepared for this study are based on the a priori knowledge of asteroid interiors of the time. The exper- imental data were obtained by performing microwave-range laboratory measurements. Two advanced in-house, full-wave modelling packages – one performing the calculations in the frequency domain and the other one in the time domain – were applied to calculate the wave interaction within the analog models. Results. The electric ﬁelds calculated via both the frequency and time domain approach are found to match the measurements appropriately. Conclusions. The present comparisons between the calculated results and laboratory measurements suggest that a high-enough cor- respondence between the measurement and numerical simulation can be achieved for the most signiﬁcant part of the scattered signal, such that the inner structure of the analog can be observed based on these ﬁelds. Full-wave modeling that predicts direct and higher order scattering effects has been proven essential for this application.


Introduction
This study concerns radar measurement as a potential technique for detecting scattering from inside of an asteroid of a complex shape. The internal structure of the small bodies of the Solar System remains an unanswered scientific question that is central to theories on the formation and evolution of the Solar System Tsiganis et al. 2005;Morbidelli et al. 2005;Desch 2007;Herique et al. 2018) as well as to the development of planetary defense strategies and asteroid deflection (Cheng et al. 2018;Michel et al. 2018). Small bodies and their interior structures have been investigated in several recent planetary space missions, most prominently, ESA's Rosetta mission which, among other things, was aimed at reconstructing the nucleus of the comet 67P/Churyumov-Gerasimenko by its tomographic radar instrument known as CONSERT (Kofman et al. 2007. Among the most significant future missions planned in the area of low-frequency radar investigation is the HERA mission , the European component of Asteroid Impact and Deflection Assesment (AIDA), which is a shared effort between ESA and NASA (Hérique et al. 2019).
A radio frequency radar measurement is one of several potential approaches to obtaining direct measurement data The measurement data obtained for the asteroid analog models are freely available on this web page https://www.fresnel.fr/ 3Ddirect/index.php following a basic registration. ) from the subsurface structures of asteroid. Due to the absence of liquid water, a low-frequency and low-power radio wave can penetrate a significant distance from several hundred meters to a few kilometers inside an asteroid regolith (Kofman 2012). Distinguishing surface and interior scattering is, however, a challenging task which necessitates advanced mathematical modeling of the electromagnetic field.
Here, our goal is to validate this distinguishability by measuring and simulating the full-field wave propagation and to show via analog models that our approach has the potential to be used in in situ radar investigations. We compare the numerical results to experimental microwave radar data measured for a 3D-printed analog object with an exterior surface shape that follows that of the asteroid 25143 Itokawa (Fujiwara et al. 2006) and an interior structure built based on current a priori knowledge of the asteroid interiors (Jutzi & Benz 2017;Carry 2012). Itokawa is a small rubble pile asteroid, an ordinary chondrite S-type assemblage, which is supposed to contain mainly pyroxene and olivine (Abell et al. 2006;Nakamura et al. 2011). Its exact shape was obtained when it rendezvoused with the Japanese space agency's (JAXA) Hayabusa mission in 2005 (Kawaguchi et al. 2008).
We performed analyses in both the frequency and time domain. In the former case, we used the dyadic Green's function recursively in order to model the full measured spectrum and calculate the field in the time domain using the inverse Fourier transform. In the latter case, a pulse with a limited bandwidth is propagated in the time domain using the finite element time-domain (FETD) method (Takala et al. 2018a;Sorsa et al. 2019). Our ultimate aim is to advance the development of both frequency-and time-domain full-wave tomography techniques for complex-structured targets, in particular, small bodies of the Solar System (Eyraud et al. 2018(Eyraud et al. , 2019Sorsa et al. 2020).
The results obtained show that the applied full-wave modeling approach is valid for a target with a realistic shape and a priori 3D structure and that the computational resources are sufficient. Together with our earlier numerical simulation results (Sorsa et al. 2019;Takala et al. 2018a), the present observations of the peak signal-to-noise ratio (peak S/N) between the measured and modeled signal suggest that the internal details of a target with a realistic shape and a priori 3D structure can be detected via radar observations. Also, higher order scattering wavefronts were observed to be predictable. Based on these results, the full-wave modeling approach seems necessary to distinguish the subsurface structures in a situation where the measurement point distribution is sparse, which is a likely scenario for a space mission.
This study is structured as follows: In Sect. 2, we describe our asteroid analog model and the experimental radar laboratory measurements. In Sect. 3, the frequency and time-domain full-waveform modeling approaches are described. Section 4 presents the results of the laboratory measurements and the comparisons. Finally, in Sect. 5, we discuss the results and in Sect. 6, we present our conclusions.

Asteroid analog model
In this study, we investigate the interior structure of two 3D-printed analog models ( Fig. 1) with the largest dimension of 20.5 cm and the exterior surface shape of the 535 m diameter asteroid 25143 Itokawa (Hayabusa Project Science Data Archive 2007). These analogs are referred to as the homogenous model (HM) and the detailed model (DM) based on their relative electric permittivity distribution, ε r , which is constant in HM and has a detailed structure in DM. Following Sorsa et al. (2019, and in prep.), we know that the DM is composed of a mantle, an interior, and an ellipsoidal detail compartment. The compartment-wise constant permittivity was modeled on the basis of the a priori knowledge of typical asteroid minerals (Herique et al. 2002 as well as on studies of asteroids' internal porosity distribution, in particular, the binary system and impact studies (Jutzi & Benz 2017;Carry 2012). The impact studies (Jutzi & Benz 2017) imply that the porosity is likely to increase towards the surface, due, for instance, to a granular mantle. The density versus volume estimates obtained for binary systems (Carry 2012) suggest that the interior structure can include inhomogeneities with a considerably higher porosity compared to its average value. The relative permittivity of HM and of the interior structure of the DM was thus assumed to be ε r ≈ 4 at a few dozen of MHz, matching roughly with that of solid or fragmented rocky minerals, for example Kaolinite or Dunite (Hérique et al. 2016). The mantle permittivity was approximated as ε r ≈ 3, modeling a granular regolith with the same mineral decomposition, and the ellipsoidal detail was assumed to be a void with (ε r = 1).
A tetrahedral mesh consisting of 21 125 nodes and 109 433 tetrahedra was applied as a metamaterial structure. The inflated edges of the mesh were 3D-printed out of a prototyping  ABS plastic filament (PREPERM ABS450, Premix Oy) with a relative permittivity of ε r ≈ 4.5 at 2.4 GHz 1 . To obtain the sought a priori permittivity values for the interior and mantle compartments, the edge width was adjusted (Saleh et al. 2020) so that the mixture between the plastic and air would roughly match the volumetric filling levels predicted by the classical Maxwell Garnett (MG), exponential (EXP), and complex refractive index model (CRIM). Denoting by η the volume fraction of the plastic filament, the MG mixing formula is given by The exponential law (Birchak et al. 1974) takes the form: where a is a case-specific exponential constant. As the exponential model reference, we use the value of a = 0.4, which has been recognized to be well-suited for mixtures of snow, air, and liquid water (Sihvola et al. 1985), where both the real and imaginary part of the permittivity vary. When a = 1/2, as, for example, in Birchak et al. (1974), the exponential model coincides with CRIM, that is, Based on a priori estimates obtained from the mixing laws, the volume fraction of the filament was set to be 66 and 90% for the mantle and interior compartments, respectively. The median structural parameters of the 3D-printed mesh constituting the analog object are found in Table 1. Notes. The volume fraction has been evaluated with respect to a 35 mm diameter sphere.

Permittivity measurements
To obtain refined estimates for the permittivity in the different parts of each analog object, we investigated the bistatic far-field scattering patterns of three 3D-printed test spheres (Eyraud et al. 2015) with volume fraction levels of 66, 90, and 100%. These spheres were based on a 35 mm diameter tetrahedral mesh, whose effective diameter (due to the edge inflation) was found to be 35.5 mm with respect to the asteroid analog (Appendix A). Table 2 shows a priori and a posteriori estimates given by the mixing laws (Sect. 2.1) together with the measured values for the mantle and interior permittivity. The estimates given by the mixing laws are somewhat higher compared to the measured values, which might be partly due to the effect of the fragmented microstructures of the 3D-printed filament resulting from the printing process. The closest match is obtained with MG with a difference of 6% and <1% to the measured mantle and interior permittivity. In relating the measurements to the effective test sphere diameter 35.5 mm, the permittivity was observed to be 1.5% lower compared to the value obtained for the original 35 mm diameter. Since the effective diameter depending on the modeled detail grows along with the size, with the largest diameter being 35.5 mm, details larger than the test spheres up to the size of the asteroid analog can have 0-1.5% lower permittivity compared to the values in Table 2; the larger the size, the greater the difference (Sorsa et al., in prep.) (Appendix A).

Experimental procedure
The electromagnetic fields scattered by the asteroid analog model were measured using the anechoic chamber of the Centre Commun de Recherches en Microondes (CCRM) in Marseille, France (Fig. 2), which makes it possible to perform measurements at a realistic distance with respect to the diameter of the analog model. Many electromagnetic scattering problems (nm to m wavelengths) can be experimentally simulated with microwaves (cm wavelengths) on a scale of a few centimeters (Geffrin et al. 2012;Vaillon & Geffrin 2014;Barreda et al. 2017;Hettak et al. 2019). Microwave experiments allow us to make accurate measurements of the quantitative complex scattered fields in controlled conditions in order to extract relevant information about the target. Thanks to the centimeter-sized analog targets, such experimental simulations are possible, on the one hand for large objects (tens or hundreds of wavelengths) that involve complex scattering effects and, on the other hand, for nanoscale analog models, which are very difficult to handle, characterize, and control in their original size.

Experiment setup
The measurement distance in the laboratory was set to be 1.85 m, which scales to a 4.83 km distance (orbit), considering the actual size of 25143 Itokawa. To detect the mantle and void inside the analog object, backscattering data were recorded from a direction in which the void is closest to the surface in the horizontal plane intersecting the centre of mass (Fig. 2). The measurements were performed in the quasi-monostatic configuration described in Fig. 2. The source and the receiver were moved together over a given range with respect to both the polar angle φ and atsimuth θ of the spherical coordinate system with the center of mass fixed to the origin (Zwillinger 2002), maintaining a constant 12 • separation angle between the transmitter and receiver. In both directions, the angular range is 30 • and the angular step is 3 • , satisfying the Nyquist sampling criterion with respect to the information content of the scattered field (Bucci & Franceschetti 1987). Consequently, both transmitter and receiver positions form a regular 11 × 11 point grid distributed over the surface of the spherical measurement coordinate system. The polarization of the field is linear along the e φ unit vector, with φφ (vertical) polarization (Fig. 2c).
The fields were measured with continuous waves between 2 and 18 GHz, with a 0.05 GHz frequency step and they were calibrated at each frequency so that the incident electric field at the origin has an amplitude of 1 S m −1 and a null phase. Consequently, the measurements at a given location correspond to a flat-spectrum point source, meaning that any time domain pulse shape can be synthesized based on the measurement data. In this study, we synthesize a Ricker window pulse with the center frequency of 6.00 GHz and two quadrature amplitude modulated Blackman-Harris (BH) window pulses with center frequencies of 10.1 and 12.9 GHz and bandwidths of 5.45 and 5.70 GHz, respectively. Altogether, these cover the investigated frequency range regarding the intensity range from −10 to 0 dB with respect to the pulse amplitude. The limit −10 dB has been chosen to ensure that the signal bandwidths are appropriately contained within the measured frequency range. Both the Ricker and BH windows are frequently used in the processing of the groundpenetrating radar measurements (Priska et al. 2019;Daniels & Institution of Electrical Engineers 2004).
The most significant factors limiting the accuracy of our experimental setup are the ≤1 mm positioning and ≤1 degree orientation error allowed by a specifically designed 3D-printed support plate (Fig. 2) and the >20 dB S /N of the radar measurement, which was obtained here by evaluating the difference between the measurement and analytically calculated field for a metallic reference sphere. In addition, the accuracy of the final experiment is affected by the modeling approaches utilized in analog preparation and numerical wavefield simulations. Figure 3 shows the 2D cuts of the radargrams obtained from the measurements considering the entire measurement domain described in Fig. 2 for the two analogs, DM and HM. The radargrams for the DM (top in the Fig. 3) and for the HM (bottom) are rather similar, but for the DM, the signal is denser than in the case of the HM, suggesting the presence of an internal structure. Based on this radargram, we deem that there are no significant hot spots in the data and thereby we choose to investigate the centermost measurement position (Fig. 2) as the primary point of interest. A68, page 3 of 12 A&A 643, A68 (2020) Notes. The a priori estimates refer to the nominal permittivity value of the ABS450 filament which has been measured by the manufacturer at 2.4 GHz frequency. The a posteriori estimates have been obtained by performing bistatic far-field scattering measurements (Eyraud et al. 2015) with the 3D-printed spheres as targets and, in the case of the mixing models, referring to the permittivity of the 100% filled sphere. The best match between a priori and a posteriori values is obtained with the MG mixing model. Notes. Reflections (1) to (7) concern the DM; the HM is only relevant to reflections (1), (6), and (7). The times indicated are two-way traveltimes.

Data
In order to enable a comparison between the measured and the simulated signal, we divided the temporal domain into mantle I (reflections due to air-mantle and mantle-interior interfaces), void (reflections due to interior-void and voidinterior interfaces), mantle II (reflections due to interior-mantle and mantle-air interfaces), and a higher order scattering zone (Fig. 4). These zones, summarized in Table 3, were determined on the basis of a plane wave propagation, as predicted by geometrical optics and considering the different paths involving a single or dual reflection. They contain the effects originating from the subdomains depicted in Fig. 4. The higher order scattering zone involves mainly multipath and multiple scattering effects. As the specific time points of interest, we consider the two-way traveltimes for the scattering boundaries (1)-(7) described in Table 3. The widths of the zones determined are based on these arrival times, taking into account the duration of the pulse.

Modeling of the wave interaction within the asteroid
To obtain the electromagnetic field scattered by an asteroid analog, we model the electromagnetic interaction with two in-house modeling packages based on (A) frequency and (B) time domain calculations. These numerical simulations are compared with the laboratory measurements in an ideally controlled environment.

Frequency domain
The interaction between the wave and the asteroid can be modeled in a harmonic domain using the integral formulation of the Maxwell equations. The electric field scattered by an inhomogeneous structure, contained in the spatial domain Ω, is written in the integral form. The scattered field, E s , on the receiver positions, r, included in the domain, Γ, is thus obtained with the observation equation: Here, E represents the electric field and G the free space dyadic Green's function between a point r in the Ω domain and a point r in the Γ domain, χ(r ; f ) = k 2 (r ; f ) − k 2 o is the contrast term where k(r ; f ) is the wavenumber at the point r from Ω and k o the wave number in the vacuum. The field E satisfies the following coupling equation: where E i (r ; f ) is the incident field at point r in Ω domain and G(r , r ; f ) the Green's dyade of the free space between the points r and r of the Ω area. Equation (5) is solved numerically using the method of moments (Harrington 1987) and a biconjugated gradient stabilized method (van der Vorst 2003); for more details, see Merchiers et al. (2010). The Toeplitz structure of the dyadic Green function is exploited in this resolution (Barrowes et al. 2001). This improves the computation speed and reduces memory requirements, which is particularly necessary for large objects in relation to the wavelength. Once the scattered field is calculated in the frequency domain, the field in the time domain is obtained in response to a Ricker pulse p(t) using an inverse Fourier transform F −1 as follows: Here, f c is the center frequency of the pulse and C = π 3/2 f c exp (1)

Time domain
Given a complex-valued permittivity distribution ε r = ε r + jε r , the wave can be modeled in the time domain by solving the wave equation for the total field, E, numerically via the following  Fig. 2 for the two analogs, DM and HM. The distance corresponds to the two-way travel after a time-distance conversion considering a propagation in vacuum. The entire frequency band (2-18 GHz) was used for this radargram. The magnitude of the scattered field is given in decibels.  Table 3. The gray and black arrows indicate the directions of the transmitter and the receiver, respectively. The wavefront (7) corresponds to a round trip inside the analog, starting from reflection number (6). High-order scattering zone includes multi-path and multiple scattering contributions. first-order system (Takala et al. 2018a;Sorsa et al. 2020): for i = 1, 2, 3 in the spatio-temporal set [0, T ] × Ω ∪ Γ with Tr denoting the trace operator which evaluates the sum of the diagonal entries of its argument matrix and Denoting the signal frequency as f , the absorption parameter (conductivity) takes the form σ = 2π f ε r . The coordinates are scaled by a spatial scaling factor of s (meters), setting the A68, page 5 of 12 A&A 643, A68 (2020) velocity of the wave in vacuum to one. With this scaling, the time, t, position, r = (x 1 , x 2 , x 3 ), absorption σ, and frequency, f, in SI-units can be obtained as (µ 0 ε 0 ) 1/2 st, sx, (ε 0 /µ 0 ) 1/2 s −1 σ, and (µ 0 ε 0 ) −1/2 s −1 f , respectively, with ε 0 = 8.85 × 10 −12 F m −1 denoting the electric permittivity of vacuum and µ 0 = 4π × 10 −7 H m −1 denoting the magnetic permeability, which is assumed to be constant in Ω. The right-hand side of Eq. (8) denotes the antenna current density given by J(t, r ) = δ p (r)J(t)e A , transmitted at the point, p, in the far-field domain Γ. The timedependence of the current is J(t). Its position, p, is determined by the Dirac's delta function, δ p (r), and orientation by the vector, e A . By finding the numerical solution, Green's function, G, can be approximated by [0, T ] × Ω ∪ Γ.
The numerical solution of the system (Eq. (8)) within the near-field subdomain [0, T ] × Ω can be obtained through the finite element time-domain (FETD) method (Sorsa et al. 2020). The incident and scattered field can be simulated in the far-field subdomain [0, T ] × Γ via surface integrals defined on the boundary of Ω (Takala et al. 2018a). An iterative time-domain solver following from the FETD discretization involves sparse matrices which necessitate effective parallel processing for their large size and the need to perform a high number of iteration steps. A state-of-the-art high-performance computing cluster equipped with high-end graphics processing units (GPUs) provides an efficient platform for running such an iteration (Takala et al. 2018a). In this study, the FETD computations were performed using the 32 GB RAM GPUs of the Puhti supercomputer 2 (CSC -IT Center for Science Ltd., Finland).

Cube as a reference case
Since this is the first time that a complex-shaped and comparably large (w.r.t. wavelength) target is measured in this quasimonostatic configuration, we first validate the present full-wave approach with a cube made of polyethylene. The side length of this cube is equal to 119 mm, corresponding to 1.2 λ and 11 λ with respect to the longest and shortest wavelength (inside the cube) within the investigated frequency band. The permittivity of the cube was found to be equal to ε r = 2.35, based on the bistatic far-field scattering patterns of a sphere made of the same material (Eyraud et al. 2015). The comparison between the laboratory measurements and numerical simulations in the frequency domain, presented in Fig. 5, and the corresponding time-domain results (after a Ricker window centered at 6.00 GHz applied to the fields and an inverse Fourier transform) in Fig. 5c validate the configuration for large objects.

Asteroid analog
The numerical simulations concerning the DM and HM analog models were performed using the measured a posteriori permittivity values (Table 2). Table 4 describes the signal properties applied in these simulations for multiple different scales.

Results with the full-wave solver in the frequency domain
The numerical simulations in the frequency domain were performed for the HM and DM analogs. Figure 6 shows a comparison 2 https://research.csc.fi/csc-s-servers In the case of (c), the signal in the time domain was obtained using the inverse Fourier transform (Eq. (7)) and a Ricker window centered at 6.00 GHz.
in the frequency domain between the measurement and the numerical simulations in the case of DM and the centermost point. The simulated field is very close to the experimental one below the frequency of 6.00 GHz. The differences that appear at around 6 GHz can be due to local inhomogeneities in the permittivity distribution or possible resonance effects related to the homogenization limit of the analog object. This point is detailed in the discussion (Sect. 5). Once the scattered field has been calculated in the frequency domain, it can be obtained in the time domain via the inverse Fourier transform (Eq. (7)). The scattered fields in the time domain, in response to a Ricker pulse with a center frequency of 6.00 GHz, are shown in Figs. 7a and b for DM and HM, respectively. Figure 7c presents the difference between the magnitude of the temporal field in these two cases. Figure 3 shows that disturbances due to echoes other than the analog echo are low as the signal has a significant value only in the time zone corresponding to the target's response. In this response zone, the calculated and measured fields are close to each other for both the DM and HM. In order to analyze these in more detail, we used the temporal arrival times (Table 3) expected for a single reflection on each interface (1)-(6) and for a dual reflection on the back interface (7) (in relation to Fig. 4). Notes. In the case of (A), the time domain data was obtained using the Ricker window and, in (B), the first time derivative of a quadrature amplitude modulated Blackman-Harris (BH) window. The dimensions are given in term of wavelengths, with λ, the wavelength in the medium having a permittivity of ε r = 3.4 + j0.04 (interior part).  Fig. 6. Scattered field for the DM analogs in the configuration described in Fig. 2 for φφ polarization and the central transmitter-receiver position (Fig. 2). Measurement is in (-) and calculation in (--). The phase shown is the phase of the field after compensating the propagation distance.
In Fig. 7, the first peak in both the measurement and the simulation correspond to the reflection between the air-mantle (1) and mantle-interior (2) interface. The ambiguity between these two responses is due to the pulse duration. The peaks corresponding to the reflections by the back of the asteroid, that is, by the interior-mantle (5) and mantle-air (6) interface, are separate and fully visible in the two curves. The distinguishability of the reflections from the void boundaries (3) and (4) is comparably weak. This somewhat counterintuitive finding suggests that the main part of the energy coming from the void does not directly propagate to the receiver. The final peak (7) corresponds to the second reflection due to the mantle-air interface,  . Scattered fields for the DM and the HM analogs in the configuration described in Fig. 2 for φφ polarization and the central transmitter-receiver position (Fig. 2). Measurement is in (-) and calculation in (--). A Ricker window centered at 6.00 GHz is used. The vertical lines (1)-(6) shown in (a) and (b) correspond to the time points of interest, whose spatial correspondence is shown in Fig. 4 and Table 3. that is, a doubled two-way path inside the analog; the time interval between (1) and (7) is roughly double the one between (1) and (6) for each analog and for both the measurement and simulation.
In the case of the DM (Fig. 7a), the measured and the simulated A68, page 7 of 12 A&A 643, A68 (2020) fields are relatively intense in the higher order scattering zone compared to the case of HM, where the fields are almost zero in this zone proving that the effects of multipath scattering are likely to be significant with the presence of internal heterogeneity(ies) (Fig. 7b). This is obvious when looking at the magnitude of the difference between the fields obtained in these two cases (Fig. 7c). This is expected as the target has a high contrast and has a relatively large diameter corresponding to the wavelength, suggesting that the wave interaction within this kind of target necessitates a full-wave approach akin the one used here, which takes into account all the multiple interactions.

Results with the full-wave solver in time domain
Figures 8 and 9 illustrate the measurement and simulation results for the two center frequencies 10.1 and 12.9 GHz of the Blackman-Harris pulse ( Table 4). The results are shown for the following three cases: (i) HM data; (ii) DM data; and (iii) the difference between DM and HM data. The approximate of the maximum observed measurement error given by the reference measurement (Sect. 2.3.2) is indicated as a shadowed region around the measurement data, reflecting the credibility of the measurement. As previously, the vertical lines (1)-(7) in the figures indicate the a priori estimated travel-times of the wavefronts scattering from the external and internal boundaries. The measured signal (solid red line) exhibits distinct peaks which localize at the areas where scattering is expected. This effect is clearly visible with the 12.9 GHz center frequency data (the right panel in the Fig. 8), while being somewhat less prominent in the case of the lower center frequency 10.1 GHz (the left panel). Of the two center frequencies, the superior distinguishability of the interior details is obtained with 12.9 GHz which yields clear electrical field amplitude peaks observed in the scattering interfaces in both the measured and the simulated data. For the DM data (the right middle panel), these peaks are localized appropriately in the air-mantle and interior-void-interior interfaces, as expected based on the estimated travel-times (the right panel in the Fig. 8), suggesting that the detail structure inside the asteroid has been detected. The difference curves (the right bottom panel) indicate that the simulated and measured amplitude peaks coincide approximately providing evidence on the detail detection via both the experimental and computational methods. The comparison between the simulated and measured HM data shows that there are no interior details in HM as the curve lacks most of the peaks observed in the case of the DM. The deviations between the measurement and simulated data can be observed to increase along with the time, which is highlighted in the higher order scattering (strongly non-linear) zone, where the best match corresponds to the DM and 12.9 GHz center frequency. Figure 9 presents a moving peak S/N between the measurement and FETD simulation for the DM, HM, and DM-HM difference signals. This peak S/N includes the effect of both the measurement and modeling accuracy of which the former is determined by the radar instrument performance (S /N > 20 dB), positioning and orientation errors (≤1 mm and ≤1 degree, respectively), and the latter one by the accuracy of the numerical FETD simulation and of the modeled permittivity, which can involve deviations of ≤1.5% due to the variation of the edge inflation effect for different surface curvatures (see Appendix A).
Of the two center frequencies, the superior peak S/N is obtained in the case of 12.9 GHz, where the peak S/N for the DM is maintained at or above 10 dB on an interval 11.95-13.4 ns, starting from the mantle I zone, covering the void zone, and the major part of the mantle II zone. Additionally, this level is reached and surpassed in the beginning of the higher order scattering zone. The corresponding peak S/N for the HM is above this level in the mantle II zone, which is obviously due to the backscattering from the exterior surface opposed to the antenna position, but lower in the mantle I and void zones, which can be understood as the mantle and void structure being absent in the HM and, thus, the scattering signal has a lower amplitude over those zones. The peak S/N of the DM-HM difference corresponding to a 12.9 GHz center frequency is above 0 dB otherwise except in the higher order scattering zone. That is, the simulated difference peaks correspond roughly to the measured ones. Because of a lower amplitude, the peak S/N of the difference is lower overall than for either of the DM and HM signals alone.

Discussion
This study concerns two asteroid analog models, a homogeneous model (HM) and a detailed model (DM), with the surface shape of the asteroid 25143 Itokawa. DM and HM were investigated through experimental radar measurements and numerical modeling. The permittivity of the interior and mantle were evaluated via an electromagnetic measurement of spherical sample objects. The results we obtained are in a good agreement with the estimates given by the classical mixing models, in particular, the Maxwell Garnett formula. The frequency modeling results show that there is an appropriate match between the measured and simulated fields electric fields in the frequency domain, in both amplitude and phase, up to a frequency of 6.00 GHz (Fig. 6). Above that, the differences observed between the measured and simulated fields may be due to the local inhomogeneities of the 3D-printed mesh and to the resonances related to the metamaterial structure. Based on a visual inspection of the analog object during the 3D printing process, the tetrahedral mesh structure might deviate slightly in vicinity of the compartment boundaries, which might be reflected in the laboratory measurements. At a 6 GHz frequency, the wavelength in the metamaterial and air (2.7 and 5.0 cm, respectively) correspond to about ten times the average edge width and edge length in the interior compartment of the 3D-printed mesh (2.4 and 4.4 mm, respectively), which might be a potential cause of a resonance. These effects are, nevertheless, not very significant in the analysis of time signals, which confirms that the permittivity values measured for the test spheres match with the direct-path travel-times of the wave within the analog. Namely, by setting the permittivity to the measured value, the measured two-way travel-time of the wave, that is, the arrival time of the signal peaks, can be observed to match with the outcome of the numerical simulation. Similarly, when comparing the simulated and measured signals, the observed total attenuation rate, that is, the total effect of absorption and multiple scattering, matches roughly with the expectation given by the permittivity measurement. Moreover, the loss tangent (ε r /ε r ) of the 3D-printed structures is about 0.01 which, based on the present knowledge of asteroid materials and minerals, is close to the expected range for a real asteroid corresponding to around 10-20 dB km −1 attenuation at 10 MHz center frequency (Kofman 2012).
The simulated fields are appropriately matched by the laboratory measurements, verifying that it is possible to distinguish the structural details as separate reflections dominated by the ones from the surface facing and opposing the radar antennae. With respect to frequency-domain modeling (6 GHz center The vertical lines (1)-(7) depict the a priori estimated travel-times (Table 3) of the wavefronts scattering from the surface, mantle, and void matching roughly with the peaks of the measured field and the simulated DM data. These peaks are emphasized in the difference data. Of the center frequencies 10.1 and 12.9 GHz, the latter one yields a superior match between the measured and simulated data.
A68, page 9 of 12 A&A 643, A68 (2020) Center frequency 10.1 GHz Center frequency 12.9 GHz Fig. 9. Moving peak S/N between the measurement and FETD simulation for the DM (dashed blue), HM (solid purple), and DM-HM difference (solid brown) signal. Here, the mean power of the simulated signal is normalized to that of the measured one and the length of the moving window is 0.5 ns. Of the center frequencies, 10.1 and 12.9 GHz, the superior overall level is obtained in the latter case, where the peak S/N for DM is maintained over 10 dB (horizontal dashed line) on the interval 11.95-13.4 ns starting in the Mantle I zone, continuing through the void zone almost to the end of the mantle II zone. frequency), the time domain results match very well with the laboratory measurements even in the high-order scattering zone. We find a clear difference between the DM and HM analogs in both numerical simulations and laboratory measurements in this high-order scattering zone and also in the mantle I and mantle II zones. With regard to the modeling in the time domain and with a higher center frequency (10.1 and 12.9 GHz), this difference also appears in the void zone, where direct scattering by the void is observed in the case of the DM. It is obvious that the temporal resolution increases with the bandwidth but also that the scattered field depends on the parameters of the applied signal. We regard this as an expected result because of the complexity and non-linearity of the scattering phenomena. For the complexity of the current wave propagation problem, the present frequency and time domain approach constitute two alternative strategies to simulating the signal. We leave a more profound comparison of these methods to a future work.
Full-wave approaches complement other methods developed to calculate the interaction of waves with a target of this size and which are based on optical physics, such as the ray-tracing technique (Ciarletti et al. 2015) and solvers based on physical optics (Berquin et al. 2015). The latter methods do not require the solution of the rigorous diffraction problem and are generally less expensive in terms of computational requirements but, on the other hand, they do not allow for a full reconstruction of the signal. The present strategy for modeling the full wavefield is all the more important since the higher order scattering effects, that is, multi-path and multiple scattering effects are likely to be significant (see for example Fig. 7). In addition, depending on the target geometry and pulse duration, the scattered wavefronts are likely to overlap and involve indirect propagation effects, implying that the full-wave approaches, which take these effects into account appropriately, are fully justified and may even be necessary in order to optimize the outcome of the structural analysis, in particular, as the interior details were observed to affect the higher order scattering wavefronts in this study.
In our earlier studies which focused on simulated measurements (Sorsa et al. 2019;Takala et al. 2018a), a peak S/N of around 10 dB has been found sufficient for finding a reconstruction of the interior structures. The peak S/N analysis of the present experimental study shows that in the case of the DM analog with internal details, this accuracy can be reached and even surpassed in mantle I, void and mantle II zone which correspond to the direct scattering originating from the mantle and void compartment. Together these results suggest that the structural details inside the DM analog object, that is, the surface layer and void, can be distinguished based on the experimental measurement data. Achieving a comparable accuracy with in situ measurements might be possible when assuming that the positioning accuracy is similar, that is, with a ≤1 degree angle error (for the potential effect of the orientation inaccuracy, see Takala et al. 2018b) and ≤1 mm position error in the analog scale (≤2.6 m scaled to the actual size of 25143 Itokawa). Namely, comparing the S/N of the laboratory measurement (Sect. 2.3.2) to the peak S/Ns observed, it is obvious that the measurement noise is low compared to the errors related to the experimental setup and modeling. Furthermore, our earlier estimates for 1 AU distance from the Sun suggest that the possible measurement errors due to the Galactic noise and the radiation by the Sun will remain at a sufficiently low level with respect to observing the interior structures (Takala et al. 2018a;Sorsa et al. 2019). In particular, their effect can be shrunk via repetitive measurements. Regarding the literature about the signal quality and noise, our present assumptions match also with the general knowledge of the ground-penetrating radar (Daniels & Institution of Electrical Engineers 2004), and by the successful CONSERT measurement .
Since the required density of the spatial discretization and, thereby, the computational cost of the simulation increases along with the signal frequency, the applicability of actual target sizes and signal frequencies might still be a computational challenge regarding real in-situ planetary radar measurements. With the current model parameters and asteroid model, we are able to simulate the wave propagation at least up to a target diameter of 264 m at 10 MHz center frequency a bandwidth of 4.42 MHz, which might be sufficient for a real mission configuration (Binzel & Kofman 2005;Kofman 2012;Herique et al. 2018). Furthermore, based on our preliminary numerical experiments, the computing platforms utilized in this study, that is, the 32 GB RAM GPUs of the Puhti supercomputer and the 128 2 TB RAM symmetrical multiprocessing CPU cores of the Mesocentre AMU, will allow us to expand the target diameter to at least to 480 m or, alternatively, the center frequency to 20 MHz. Considering the potential mission relevance of these numbers, for example, the 163-m moon of the asteroid 65803 Didymos which is the target of the HERA mission falls into these limits . Also, larger asteroids can be considered, the maximum size being dependent on the computing platform's performance. Targets that do not fit into the memory available can be modeled by restricting the wave simulation process into a subdomain. If needed, such a restriction could be done based on the present principles to divide the scattering processes into lower and higher order scattering zones.

Conclusion
The results of this study, obtained in microwave-range laboratory experiments and computer simulations with an analog model based on the shape of the asteroid 25143 Itokawa as a target, suggest that full-wave modeling is a valid approach for investigating the interior structure of asteroids based on radar measurements in both the frequency and time domains. Comparing the present results to earlier numerical simulation studies suggests that a high enough peak S/N can be obtained for the direct scattering from the mantle and void structures, that is, in mantle I, void and mantle II zones. Our results underline that due to the complexity of the scattered wavefield, full-wave modeling might be essential for optimizing the outcome of in situ measurements which are potentially sparse in space missions. Namely, it distinguishes the overlapping direct and indirect wavefronts scattering from different parts of the target and, thereby, allows for the detection of otherwise weakly detectable internal details. For this capability, we are also able to predict multiple scattered and multi-path wavefronts (i.e., data corresponding to higher order scattering zone) which can potentially be utilized in the observation of the interior structures. As the present experiment setup and permittivity modeling approach are applicable with various measurement configurations, asteroid shapes, and interior structures, they can potentially be applied in the examination of subsequent asteroid models in the future.