Towards Printable Natural Dielectric Cloaks via Inverse Scattering Techniques

The synthesis of non-magnetic 2D dielectric cloaks as proper solutions of an inverse scattering problem is addressed in this paper. Adopting the relevant integral formulation governing the scattering phenomena, analytic and numerical approaches are exploited to provide new insights on how frequency and direction of arrival of the incoming wave may influence the cloaking mechanism in terms of permittivity distribution within the cover region. In quasi-static (subwavelength) regime a solution is analytically derived in terms of homogeneous artificial dielectric cover with $\varepsilon<\varepsilon_0$ which is found to be a necessary and sufficient condition for achieving omnidirectional cloaking. On the other hand, beyond quasi-static regime, the cloaking problem is addressed as an optimization task looking for only natural dielectric coatings with $\varepsilon>\varepsilon_0$ able to hide the object for a given number of directions of the incident field. Simulated results confirm the validity of both analytic and numerical methodologies and allow to estimate effective bandwidths both in terms of frequency range and direction of arrival of the impinging field.

topology optimization approaches [7][8][9] at optical frequencies in the case of metallic objects, requiring extremely variable covers to induce cloaking exclusively for a narrow angular range of DoAs. Other all-ordinary dielectric cloaks have been proposed via global optimization techniques for radially cylindrical and spherical cloaks of metallic targets [10][11][12][13][14] . These approaches, although very attractive in terms of simplicity of the geometrical architecture of the covers, require very large refractive indexes 11,12 , rarely available in nature, or even near-zero dielectric constants 13 . Moreover they achieve only nearly optimal performance in terms of residual scattering radiation 14 . Another optimization approach based on the phase field method has been recently proposed 15 for cloaking metallic cylinder, by considering six different angular directions of the incoming wave. Finally, graded refractive index structures have been proposed for surface cloaking 16,17 .
Differently from the above contributions, mostly concerned with the cloaking of metallic objects and surfaces, in this paper, we investigate, through recent analytic 18 and numerical methods 19 , the conditions to pursue non-magnetic volumetric dielectric cloaking with natural permittivity materials (ε > ε 0 ) for penetrable scatterers. In this respect, we adopt a design procedure based on the solution of an inverse scattering problem 1 exploiting a gradient based local optimization approach which allows to easily handle many of the design constraints concerned with the cloaking problem. A physical insight into artificial-natural and natural-natural dielectric cloaking systems is also given. The presented methodology is mainly focused on investigations about volumetric dielectric scatterers, with a possible extension towards composite structures such as metal-dielectric metasurfaces 20 .

Methods: Analytic and Numerical Cloak Synthesis Procedures
We consider a 2D domain Ω embedding one (or more) penetrable non-magnetic homogeneous object(s) of arbitrary shape with support Σ 1 , see Fig. 1. A cloak region Σ 2 such that 1 2 ∪ Σ = Σ Σ is considered, and, for the sake of simplicity, the intersection between Σ 1 and Σ 2 form a null set. However, the analysis developed in the following is valid even if this hypothesis is removed.
Let us assume the Ω domain in the xy plane, and one (or more) plane wave(s) with unitary amplitude and electric field polarized along the ẑ axis (TM polarization) impinging towards the center of Ω, wherein = r x y ( , ) denotes the vector position in Ω, ω = 2πf the angular frequency and θ ν the DoA. Throughout the paper, the time harmonic factor e jωt is assumed and dropped.
For the sake of clarity, the dielectric properties of the overall region are grouped in: with Γ ≡ Ω Σ / the observation region where no "visible" scattering effects by the cloaking system would be desired and where the constitutive parameters are homogeneous (i.e., ε b does not depend on r in Γ).
From now on, arrow signatures on fields are suppressed, tacitly assumed all vectors being directed along ẑ. The electromagnetic scattering from such a cylindrical (i.e., longitudinal invariant) structure is due to the equivalent volumetric sources with support Σ, defined as J r j J r ( , , ) ( , , ) i s θ ω χ θ ω θ ω = + ν ν ν Figure 1. Geometry of the problem for the synthesis of dielectric cloaking: Ω is the computational domain which includes arbitrarily shaped scattering system Σ (cloaking region) and / Γ = Ω Σ (observation region). The Σ region, divided into Σ 1 (bare object) and Σ 2 (cloak cover), is illuminated by an incident field E i impinging from an angular directions θ ν . 0 Γ is a subset of Γ where the scattered field is enforced to be zero in the design procedure based on inverse scattering technique.
Scientific RepoRts | 7: 3680 | DOI:10.1038/s41598-017-03749-y In Eq. (3), J(·) is the so called contrast source given by the product between the total internal field E r E r E r ( , ) ( , ) ( , ) t i s ⋅ ≡ ⋅ + ⋅ and the contrast function defined as: The total field E t , expressed in the whole region Ω as the sum of the incident (or primary) and scattered (or secondary) field, can be conveniently expressed via integral formulation 21 : is the 2D Green's function of the homogeneous background, whose kernel has exact analytic form in a homogeneous background, namely the Hankel function of zero order and second kind 21 . According to the partition in Eq. (2), by definition, the contrast function is zero in the observation domain / Γ ≡ Ω Σ and it is mostly non-zero elsewhere. Let us also notice that, by definition, if vacuum is assumed as background medium (i.e., ε b = ε 0 ), the contrast function is coincident with the electric susceptibility. Equation (5) states that the physical cause of the scattering phenomenon is the contrast source J induced in Σ. It may be convenient to express the scattered field through a more compact notation with respect to Σ and Γ domains, as: [ ] with (7) s  Eqs (6) and (7) can be identified as the object and data integral equations of the ISP, respectively. From a physical point of view, the operators Σ  : Σ → Σ L L ( ) ( ) 2 2 and Γ  : L L ( ) ( ) 2 2 Σ → Γ map the contrast source into the corresponding scattered field in Σ and in Γ, respectively. Adopting this formulation, the cloaking effect can be pursued by enforcing In the following, solutions of Eqs (6) and (7) are pursued to achieve cloaking in and beyond the quasi-static regime, using both analytic and numerical approaches.
Cloaking in quasi-static regime. When Eq. (8) is assumed as desired specification of any cloaked system, the data Eq. (7) can be explicited as: ν Σ J r G r r dr r ( , , ) ( , , ) 0 with (9) Considering the overall system enclosed in a circular cylinder of radius b, in quasi-static condition (i.e., k b b→0), Eq. (9) can be solved in a straightforward and simple manner since the overall system is so extremely compact in terms of λ that the Green's function of the homogeneous background becomes constant over the entire domain Σ, i.e., in the quasi static-limit. Expliciting the contrast source according to Eq. (3), also the total field can be considered to be constant in the Σ domain in the quasi-static approximation, i.e., the Rayleigh scattering regime, thus giving: ∬ Splitting Eq. (11), namely the Contrast Cloaking Equation (CCE) 18 over the object region, χ ∈ Σ 1 1 and the cloak region 2 2 χ ∈ Σ , a necessary and sufficient condition to achieve cloaking comes out as a proper mix of positive/ negative values of the local contrast function, i.e., It is interesting to notice that the CCE generalizes in a very compact fashion PC for scatterers of arbitrary shape 22 for a general background medium, even when ε ε ≠ b 0 . Therefore, in the quasi-static limit, the designer can turn-off the effect of volumetric sources by locally compensating the positive-induced source associated to a positive contrast (e.g., the object) with the negative-induced one associated to a negative contrast (e.g., the cloak) or viceversa, regardless shape of the cloak system and DoA of the incident field. As a matter of fact, since the total field does not play any role in this derivation (being factored out from the integral), such a kind of cloaking is expected to behave as an omni-directional cloaking, i.e., its performances do not depend from the DoA of the incident field.
Scientific RepoRts | 7: 3680 | DOI:10.1038/s41598-017-03749-y Cloaking beyond quasi-static regime. When the overall cloaking dimension is not in subwavelength condition, the two main hypothesis for deriving the CCE are no longer valid. Beyond quasi-static condition, distributed effects take place with two consequences: (i) the Green's function is no more constant and (ii) the total field changes from point to point, taking into account all the non-local contributions of the scatterers in the whole domain Σ 21 .
However, thanks to these considerations, as also depicted in Fig. 2, we figure out an important finding: the cloaking mechanism can be achieved with all-positive values of the contrast provided that a proper spatial organization of the contrast layout is pursued. As a result, the cancellation effects occurring between positive-negative values of the contrast, that are very close at subwavelength scale, can be synthesized also for all-positive contrast values which are not so close in terms of wavelength (e.g., crest and trough of the working wavelength). On the other hand, due to the need of specific spatial distribution of the contrast function, the cloaking effect is expected not to be broadband. As a matter of fact, when the ratio D/λ increases, D being the diameter of the minimum circle enclosing Σ, the architecture scheme of the coating plays a crucial role with the possibility to get rid of metamaterials in Σ 2 through a proper arrangement of all-positive dielectric values ε ε ≥ r ( ) 2 0 . On the other hand, as stressed above, according to PC 2 , just one homogeneous layer of metamaterial is sufficient to achieve cancellation effects in quasi-static regime.
Under the above reasonings, the synthesis of cloaking profiles can be conveniently addressed as the solution of an ISP without any approximation on the mathematical model in Eqs (6) and (7). This can be performed through the minimization of a proper cost function, under the constraint r ( ) 0 2 χ > . Obviously, solutions accounting for both positive and negative values of the contrast function can be anyway pursued, but their investigation is beyond the scope of the present paper. Once fixed 1 χ in Σ 1 , the adopted formulation of the cost function depends on both contrast χ r ( ) 2 ∀ r 2 ∈ Σ and contrast source J r ( ) ∀ ∈ Σ r . In particular, the cost function is obtained as joint minimization of the weighted object and data equations (6) and (7). To this aim, it is convenient to modify the object equation according to the contrast source formulation 23 , just multiplying both members of Eq. (6) by χ. Moreover, the part of functional (13) relevant to the data equation has to be properly modified to take into account zero scattered field enforced in the optimization task. A possible way to address such a problem is given by the minimization of the following weighted cost function: where the apex ν stands for the ν-th impinging DoA and ||·|| denotes the usual L 2 -norm. In particular, in Eq. (13), the first addendum enforces, for a given set of plane waves, the solution (in the least square sense) of the scattering equation in Σ, while, interestingly, the second term stands for the minimization of the radar cross section (RCS), or the echo width 21 in a subset 0 Γ ∈ Γ, see Fig. 1. In this respect, it is worth noting that the finite bandwidth of the scattered fields 24 allows to not enforce zero scattered field in the whole observation region Γ, since it is sufficient to enforce such a value only in a finite number of points M over a surface 0 Γ . For circular cylindrical geometry (a circumference of radius R for the particular case at hand), according to the degrees of freedom of the scattered field ref. 24, the minimum non redundant number of sampling points M ≈ 2k b R can be considered as angularly equispaced along the circumference 0 Γ enclosing the cloaking system. As a last comment, it is worth noting that the minimization problem in Eq. (13) entails a non-quadratic form 19,25 , so that the minimization procedure may get stuck into "local minima" 26 . The possibility to incur in local minima is strictly related to the functional shape (depending on the data and constraints of the problem) as well as on the initial guess adopted to start the minimization procedure. In this respect, while the first set of parameters are fixed in advance by the designer, the initial guess is a "degree of freedom" of the synthesis problem, which can be conveniently exploited to obtain several equivalent solutions (in terms of cloaking effect) provided that a satisfactory weighted residual error is reached in the minimization of (13).

Numerical Results and Analysis
To keep things simple, we have considered the cloaking of a circular scatterer in free-space background (i.e., ε b = ε 0 ), made up of lossless allumina (ε s = 10ε 0 ) with radius a = 0.42 cm, see Fig. 3(a). However, it is worth noticing that such an approach can be easily adopted for any arbitrary (more sophisticated) shape of the scattering system (both object and cover).
The cloaking effect will be demonstrated, both in and beyond the quasi static-limit, adopting circular coatings with different radii.
To achieve cloaking in the subwavelength condition, we set the cover's radius to b = 0.6 cm (slightly larger than that of the bare object). In particular, the coating diameter 2b corresponds to 0.16λ at 4 GHz (quasi-static regime).
In quasi-static condition the CCE is exploited and, following Eq. (12), the contrast in Σ 2 is computed as χ χ = − Σ Σ / ) and the overall cloaking system is reported in Fig. 3(b).
Beyond subwavelength condition, the architecture for the natural dielectric cloak is obtained solving the optimization problem (13) via conjugate gradient fast Fourier transform method (CG-FFT) 25 , adopting a standard pixel representation for the discretization of the state and data equations, for both the scatterer and the cover ε r ( ) 2 in Σ 2 . Two different dimensions for the cloak have been considered in order to take into account the performance of the cloak with respect to its electrical dimension (compared to the operational wavelength). The discretization of the analysis domain has been fixed according to rules of the integral equation method (MoM) 21 . Specifically,  Table 1.  the dimension of the pixel has been set to 0.5 mm, which is slightly larger than λ/10 of the minimum wavelength (as referred to the allumina). This choice has been established as trade-off between the accuracy required by the numerical procedure and the need to deal with reasonable resolution in realizing the cloak trough solid printing techniques.
In order to take into account the dependence of the cloak from the impinging directions of the incoming wave, we have considered in the optimization problem (13), different number of DoAs, i.e. N = 2, N = 4 and N = 8 plane waves angularly evenly spaced on a complete arc of 360°. On the other hand, since the cloaking effects is required all around the cloak we have enforced the scattered field to be zero in M = 24 or M = 36 equispaced observation points located on a circumference placed in the close proximity of the cloaking system. The circumference's radius has been set to  . R b 1 45 , b being the different coating radius, with b = 1.2 cm and b = 2.4 cm (1λ and 2λ at 25 GHz, respectively) adopted in the synthesis procedure. These different values have been considered in order to take into account the different number of degrees of freedom of the scattered field pertaining to scatterers of different extent 25 .
According to the above reasons, we have considered a discretization grid of 48 × 48 pixels for the cloak with radius b = 1.2 cm (1λ) and 96 × 96 pixels for the cloak with radius b = 2.4 cm (2λ). Moreover, in the minimization procedure, constraints about natural permittivity (i.e., ε ≥ ε 0 ) and lossless (i.e., σ = 0) nature of the cover have been enforced at each step of the minimization procedure. In this respect, the maximum value of the permittivity to be used in the cover region has been set to that of the allumina (i.e., 10ε 0 ). It is worth noticing this constrain simplifies the fabrication of the cloak in terms of number of materials to be employed, as well as, not less important, allows to take under control, i.e. to not violate, the spatial discretization of the integral equations involved in the synthesis strategy.
The synthesis procedures was stopped when the functional reached a value of 10 6 Φ ≤ − . This ensures that both data and state equation are solved with a satisfactory accuracy. On the other hand, when the threshold value was not reached, the procedure was stopped when no significant changes arose in the functional value between two subsequent steps of the minimization procedure.
The initial guess of the dielectric cover has been set as homogeneous in the whole Σ 2 with relative permittivity of ε 2 = 3.5ε 0 or ε 2 = 5.5ε 0 , for N = 2, 8 and N = 4, respectively, in order to escape from local minima wherein the residual value of the functional (13) is unsatisfactory for achieving cloaking effects. It is worth to note that these initial values have been chosen after some a posteriori checks.
The results of the synthesis strategy are shown in Fig. 3(c-h), after a proper "post-quantization" procedure, which slightly adjusts the dielectric profile of the cover to predefined permittivity levels fixed at possible values of ε r = 1, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 10. Such procedure entails a simplification of the cloak manufacturing in terms of piecewise constant dielectric profile taking into account the availability of solid materials in the considered frequency range, while not substantially modifying the cloaking performance, which results only slightly worsened. Figure 4(a,b) show the real part of the total field for the uncloaked case and the plasmonic coating at 4 GHz: in Fig. 4(c-h), the real part of the total field for all the natural dielectric cloaks at 25 GHz is reported. The analysis was performed in COMSOL 27 importing the synthesized profile and meshing it by means of triangular shaped elements.
At the lower frequency, the plasmonic cover is in cloaking operation with the flat phase fronts well recognizable behind the object, whereas, as expected, at the same frequency, no scattering cancellations occur for the natural dielectric cover (not shown). At frequency of 25 GHz, the dielectric cloaks are all able to achieve cloaking mechanism for the DoAs considered in the design procedure (for the sake of brevity the cloaking effect is shown only for one of the designed θ ν impinging directions). From these results it is possible to observe that the cloaking can be achieved both for the cloak with radius of 1λ and 2λ, and that the smaller the cloak and the larger the number of DoAs, the more complicated the architecture of the covers to be synthesized is. Moreover, note that not all the discrete permittivity's values reported above are required by the synthesized covers. For example, for the cloaks of Fig. 3(c,h) only four kind of dielectrics are required with a further simplification in the architecture of the cloak.
In order to quantify the cloaking effect as a function of the frequency, we have calculated the scattering cross section (SCS) of the synthesized cloak, which is the average value of the RCS calculated at 72 equispaced points on a circle of radius R = 1.45b for each given DoA (the same shown in Fig. 4). It is worth noticing the numerical values of the SCS when calculated for different DoAs does not substantially change. The overall function SCS (θ ν , ω) has been calculated for different frequency values at a fixed DoA and for different DoAs at a fixed (designed) frequency: this leads to quantify the bandwidth performance of the synthesized cloaks in terms of frequency and omnidirectionality issue.
As shown in Fig. 5(a), in subwavelength condition (below 4 GHz), the plasmonic cover (dotted line) drastically reduces the scattering levels with respect to the uncloaked case (continuous line).
For what concerns the ordinary dielectric cloaks with radius b = 1.2 cm (1λ at 25 GHz), in Fig. 5(a), all they share the minimum around the designed frequency 25 GHz, but they possess very different angular responses as reported in Fig. 5(c), due to the fact that they have been designed for different incoming waves. The cloaking effect is clearly affected by changing the DoA: for the 2-views ordinary dielectric cloak (upward-pointing triangle line), the performance becomes worse when non-optimal DoA is considered, whereas, as expected, two different minima (as the numbers of views for which it has been designed) appear at DoA = 0, π. The 4-views ordinary dielectric cloak (dash-dot line) improves its omnidirectional performance, ensuring an overall response which is always below the uncloaked case around −14 dB (continuous line), especially when the 4 minima occur. The 8-views ordinary dielectric cloak (dashed line) shows the best omnidirectional performance, ensuring an overall response between −27 dB and −32 dB.
In Fig. 5(b-d), the SCS (θ ν , ω) is shown for the case b = 2.4 cm (2λ at 25 GHz). As reported in Fig. 5(b), the scattering from the ordinary dielectric cloak is larger in the low frequency window with respect their compact counterpart in Fig. 5(a), due to the fact that now their size is increased: around 25 GHz, they show a minimum with slightly worse performance in terms of dB reduction compared to the cloaks with b = 1λ in Fig. 5(a). Even the angular swing between optimal and worst DoA value, as reported in Fig. 5(d), is increased. In this case, the cloaking effect is still affected by changing the DoA, but in a different way with respect to the previous case. The 2-views ordinary dielectric cloak (upward-pointing triangle line) has wide regions between the two minima points (DoA = 0, π) for which it scatters more than the bare object. In this case, also the behavior of the 4-views ordinary dielectric cloak (dash-dot line) becomes worse, with several values of the SCS at 25 GHz above the uncloaked case (continuous line) of about +5 dB. Even the 8-views ordinary dielectric cloak (dashed line) loses its omnidirectionality even if it reachs values of about -27dB in its 8 minima.

Conclusion
The synthesis of all dielectric cloaks has been tackled through the solution of an inverse scattering problem where zero scattered field is properly pursued with artificial and natural dielectric materials within and outside the quasi-static frequency regime, respectively. It has been found and discussed that it is not strictly necessary the use of ENZ or ENG materials for cloaking beyond quasi-static limit, completely changing the paradigm suggested by quasi-static formulas and exploring the potentiality of cloaking techniques via scattering cancellation. For this reason, natural dielectric cloaks can be synthesized in a relatively compact and easy fashion. Moreover, the behavior of the cloak is non-resonant and a non-negligible operational bandwidth (at 3 dB) can be achieved (about 15%). Interestingly, being made of only natural dielectrics substances, the cloak could be suitably manufactured by means of solid multi-filament printing techniques when possibly the local permittivity values can be achieved using alternative dielectric mixtures with differing volume fractions and particle sizes 16,17 . As a result, this would pave the way to cloaks easy to fabricate with respect to natural dielectrics available in the frequency range of interest and exploiting non-uniqueness of solution in the design procedure based on inverse scattering techniques. We are currently working on the realization of such non-homogeneous dielectric covers for experimental measurement of their performance. Further efforts will be addressed to take into account additional design parameters and constraints, as well as to tackle the full 3D electromagnetic cases.