1. Introduction

A generic problem in optics and electromagnetics is to design objects in a manner that allows their scattering response to be engineered. That is, for a set of incident wave fronts we want to obtain a specified scattered field amplitude. A prominent example is diffractive optics, which is mainly concerned with designing and fabricating thin optical elements that modulate the incident wave front in a single plane [1]. The restriction to thin elements is highly suitable for lithographic fabrication methods, but prohibits angular selectivity of the response function. While a number of ad hoc solutions for achieving angular selectivity have been proposed [2–5], the extension to volume elements remains challenging. In principle, multiplex volume holograms can be used to implement essentially any desired system response. In practice, however, problems such as low diffraction efficiency and the experimental burden of the recording process limit the range of applications. More recent efforts to extend the variety of custom designed system functions are based on photonic crystals and meta-materials [6,7], where the three-dimensional (3D) permittivity distribution is sculptured by manipulating the diffracting structure on a nanoscopic scale. While comprehensive models of the fabrication techniques used to create nano-photonic structures reveal a variety of challenges [8,9], experimental work as reached the point, where practical applications can be considered [10].

While the ultimate goal of designing structures that transform any incident wave into any desired scattered field may not be completely possible, control of scattering for a subset of incident waves remains a realistic focus. This paper is concerned primarily with electromagnetic cloaks which disguise or render objects invisible to an observer by encapsulating it inside a special material that sculpts the scattered field in an appropriate manner to perform the cloaking function.

Research efforts to address these design problems include a proposal for smart obstacles [11]. Based on control theory, objects are designed which change their boundary properties as a function of the incident field. In effect, this determines a prescribed nonlinear material response which realizes the desired scattering properties.

Other work based on ray-optical models suggests the design of passive artificial materials acting as cloaks for embedded objects [12–15]. Rigorous simulations have confirmed the ray optical analysis at least qualitatively [16,17]. In fact, perfect cloaking may not be achievable for rather fundamental reasons [18]. We also note the strong relationship between the design of electromagnetic cloaks and the classic problem of non-radiating sources [19].

In our work, we address the generalized problem of designing scattering structures as a problem complementary to inverse scattering. This allows us to classify design problems and select suitable methods based on concepts developed for imaging applications. Most importantly, classifying the permittivity contrast function into weakly scattering and strongly scattering requires rather different strategies for solving the inverse problem.

For structures with low permittivity contrast, design can be addressed by specifying the scattered field amplitudes in the Fourier transform domain. In essence, it is possible to apply the concept of linear diffraction tomography to specify the desired object distribution [20].

For high permittivity contrast functions, the inverse problem is known to be very difficult. While advances have been made using iterative techniques, long computational times are necessary and convergence is not guaranteed. To address this class of inverse problems, we developed a method which is simple to implement and which can provide a good estimate of the scattering structure [21]. More recently, we adapted this method to address the synthesis problem as well [22]. This was motivated by the fact that the structure synthesis problem is less constrained than the imaging problem. In principle, any structure which provides the desired scattering characteristics is a solution to the synthesis problem. The remaining concern is how easily the predicted scattering structure can be fabricated, i.e. the range of permittivity and the scale of structural features.

The inversion and synthesis method we use is grounded in diffraction tomography. It is a Fourier inversion procedure followed by a nonlinear signal processing step that involves filtering the cepstrum of the secondary source distribution. This obtains a relationship between the permittivity distribution and the associated k-space of scattered field amplitudes.

In this paper, we apply cepstral filtering to the problem of designing electromagnetic cloaks, i.e. covers for objects which are encapsulated by a suitable dielectric function. While this addresses a problem of widespread interest, it also allows us to expand the scope of our method as compared to the work presented in Ref. 22. In order to design cloaks with an isotropic response, we need to control essentially the entire plane wave spectrum of the scattered field rather than only a small number of discrete scattering angles. In addition, control over the scattered field has to be accomplished for any incident plane wave. We demonstrate how the latter can be accomplished by exploiting the symmetry of the scattered field data used as the starting point for our design algorithm. A further key point is the use of rigorously simulated field far field data, rather than experimental data. This provides greater fidelity and flexibility to our method and results.

We do not claim this method provides a rigorous solution to the design problem. Instead, it provides us with a fast and efficient method to arrive at an approximate solution that can be further improved with only a few iterations of a suitable iterative method. However, validation of the design with a rigorous diffraction model shows good correspondence between the simulated scattered far field and the desired scattered field distribution.

2. Nonlinear inverse scattering algorithm

Consider a penetrable scattering object, V(r) in free space. The target V(r) is related to the permittivity by V(r) = k ²[ε(r) -1]. The scattered field Ψ_s(r,k r^₀), due to the interaction of an incident wave Ψ₀(r,k r^₀) with the target is given in two dimensions by the integral equation

where r^₀ is a unit vector in the direction of the incident field. The solution to this equation requires knowledge of the total field Ψ(r,k r^₀) within the object volume D, but this is not possible when V(r) is unknown. Using the first Born approximation [20], one assumes the total field Ψ(r,k r^₀) in the integral can be replaced by the known incident field Ψ₀(r,k r^₀). This makes the inversion problem linear and permits one to find a solution. The integral in Eq. (1) reduces to a Fourier transform relation between the scattered field Ψ_s(r,k r^₀) and the target V(r). Consequently, each measurement sample taken of the scattered far field can be related to one sample in k-space of the Fourier transform of V(r) by applying the data mapping illustrated in Fig. 1. Physically this requires that the scattering from the object be extremely weak in order for the total field everywhere within the object to be well approximated by the incident field. When this requirement is not valid, the same data inversion step yields information about the secondary sources, or the so-called contrast source function V _B(r). This is

This first Born reconstruction is modulated by the field pattern within D, which will be different for every illumination direction r ₀. In diffraction tomography, the scattered field data for all incident field directions is combined in k-space and a Fourier inversion of that data provides an estimate for V(r). This is possible because Ψ(r,k r^₀)≈Ψ₀(r,k r^₀) and V_B(r)≈V(r) to at least within low-pass spatial filtering limits resulting from the available k-space coverage. When the first Born approximation is not valid, then V_B(r)≈V(r)〈Ψ(r)〉 where 〈Ψ(r)〉 is a complex noise-like term with a characteristic range of spatial frequencies determined by the bandwidth of the source.

 

Fig. 1. k-space interpretation of scattering in the first order Born approximation: (a) geometry for plane wave scattering of a permittivity distribution; (b) k-space representation of the incident wave and one scattered plane wave component. The scattering amplitude is proportional to the object spectrum at the Ewald circle.

Download Full Size | PDF

The image of the product V∙〈Ψ〉 exhibits spatial fluctuations characteristic of the wavelengths being employed as well as the spatial fluctuations of the permittivity. For incremental illumination wavelength changes, the 〈Ψ〉 term will change quite considerably but V(r) need not. The first step in the nonlinear or homomorphic filtering procedure is to take the logarithm of the product of V∙〈Ψ〉 and perform a spatial filtering operation in the Fourier domain to remove the field component. Since

there are numerical problems when the magnitude of the V∙〈Ψ〉 is close to zero where its logarithm becomes singular. In addition, when the phase of V∙〈Ψ〉 has a range that exceeds 2π, the resulting wrapped phase introduces spurious spatial frequencies in the log-Fourier space, or cepstrum. Phase unwrapping is exceedingly difficult, especially in two and higher dimensional problems, because zeros in the field are associated with phase or wavefront dislocations. Phase discontinuities correspond to broad spatial frequency features which make successful homomorphic filtering impossible. We mitigate these problems by preprocessing the data to satisfy a minimum phase condition.

3. Minimum-phase based homomorphic filtering

The concept of minimum phase is well understood for one-dimensional problems. A one-dimensional signal is minimum phase if, and only if, F(x + iy) has a zero-free upper half-plane. That is, it has no zeros for y > 0. The meaning of a zero-free half-plane in two dimensions is less easily conveyed but an important feature of a minimum phase function is that the phase is continuous and bounded between -π and π. The function has “minimum phase” in the sense that phase is not wrapped.

It is possible to enforce the minimum phase condition on a function by applying Rouche’s theorem [14]. This states that if a band limited function H(r) has N zeros within some contour, another band limited function F(r) has M zeros in the same contour, and |H| > |F| on that contour, then H + F will have N zeros within the contour. In other words, the sum of two functions will have the number of zeros equal to the number of zeros of the larger magnitude function. Consequently, adding a sufficiently large minimum phase background or reference wave H to a band limited function F, guarantees that the sum will be a minimum phase function. It is therefore possible to preprocess V∙〈Ψ〉 by introducing a reference point in k-space around where the scattered field data resides to make a minimum phase function prior to calculating its logarithm.

As a practical first step, the data in k-space is made causal by moving it to one quadrant of k-space and adding a reference point at the origin. This is equivalent to adding a reference wave to V∙〈Ψ〉 which needs an amplitude just large enough to ensure that the phase of V∙〈Ψ〉 is continuous and lies within the bounds of -π and +π. This is readily determined by inspecting the phase of the modified V ∙〈Ψ〉 calculated by inverse Fourier transforming the k-space data. Second, implementation of the homomorphic filtering algorithm requires a low pass filter to be applied in the log-Fourier, or cepstral, domain to suppress the wavelike features in the resulting image associated with 〈Ψ〉. This spatial filtering is successful to the extent to which the field internal to the scattering structure has spatial frequencies that are distinct from those of log (V).

4. Structure synthesis by k-space manipulation

The method can now be applied to structure synthesis. For weakly scattering permittivity distributions we could engineer the scattering behavior by filling k-space with the associated k-space amplitude distribution and obtain the permittivity, at least in an approximate sense from Fourier inversion. Numerical optimization may help to meet additional constraints imposed on the spatial distribution and its spectrum.

For objects with a strong permittivity contrast we use the inverse scattering method as outlined. To ensure a physically realizable solution, we illustrate our method by modifying the scattered field of a known scattering object. For example, in Fig. 2(a), we show the exact scattered field for a cylinder of refractive index n=√ε=2.0 and diameter 4λ when illuminated from below. This data lies on a single Ewald circle in k-space. When many illumination directions circumscribing the cylinder are used to fill k-space, as shown in Fig. 2(b), the data can be inverted to recover an image of the cylinder. This is a straightforward inverse Fourier transform for a weakly scattering cylinder, but in our case, the inverse Fourier transform from this strongly scattering cylinder produces the field image shown in Fig. 3. After cepstral filtering, a much improved Fig. 4(a) is obtained. This can be compared with the ideal result from spectral filtering of a similarly sized cylinder, modified by the same low-pass filter in Fig. 4(b).

 

Fig. 2. (a) Rigorous solution of plane wave scattering off a homogeneous cylinder, (b) k-space constructed from the far field for a set of different incident field directions.

Download Full Size | PDF

Starting with the k-space data depicted in Fig. 2(b), scattering in specific directions can be suppressed by reducing the amplitude of corresponding points in k-space. For example, a disc centered at the k-space origin with greatly reduced magnitude corresponds to constructing a scattering distribution of a cylindrical object modified to no longer scatter in any direction as long as the incident wave vector satisfies |k|/2<K. Reintroducing a non-zero amplitude at k=0 permits incident waves to scatter in forward direction only, thereby approximating an invisible object instead of a completely non-scattering object.

 

Fig. 3. Magnitude of the inverse Fourier transform of the k-space distribution shown in Fig. 2(b).

Download Full Size | PDF

We show in Fig. 5, an example of modifying k-space as described above. To design an invisible object where scattering is observed only in the forward direction, the k-space signal was removed within a circle centered at the k-space origin out to a radius of k ₀/2, where k ₀ is the wave number of the input data set. All other spatial frequencies of the k-space data set were left unchanged as well as the point at the k-space origin. The later represents the direction of the incident plane wave. The reason for zeroing k-space out to k ₀/2 was to allow illumination by a plane wave with at least twice the wavelength of the input data set for which no scattering should occur since the corresponding Ewald sphere lies entirely within the zeroed region. We note that a clean and precise nulling of the forward scattered field was not imposed because of the highly nonphysical nature of our k-space constraint. It is reasonable, however, to impose very small amplitudes in the chosen k-space regions to avoid errors and non-physical attributes. Quantitative information about the permittivity is lost in our cepstral inverse scattering method. However, the refractive index modulation can be calibrated from the knowledge of the refractive index of the original cylinder which was used to calculate the scattering data that serve as input to our algorithm. In this case, the refractive index was n=2.0.

 

Fig. 4. (a) Reconstruction of object permittivity contrast based on (a) cepstral filtering, and (b) spectral filtering.

Download Full Size | PDF

Verification that the synthesized structure has the anticipated scattering properties is performed by computing the scattered fields in the forward direction using a rigorous method such as finite-difference time-domain (FDTD) [23,24] or finite-difference frequency-domain (FDFD) [25] algorithm. We would expect the scattering objects predicted by our cepstral filtering approach to at least provide a good initial structure that can be further optimized using a more rigorous method.

 

Fig. 5. Reconstruction of object permittivity contrast after k-space engineering based on (a) cepstral filtering, and (b) spectral filtering.

Download Full Size | PDF

5. Rigorous simulation results

To verify the nonlinear signal processing technique described above, rigorous simulations were performed using a FDFD method [25]. This method was chosen over the more popular FDTD method to more efficiently accommodate scattering objects that may be highly resonant or incorporate materials with refractive index near zero.

Figure 6 summarizes scattering from a perfect dielectric cylinder (n=1.5) in free space. The object itself is depicted in Fig. 6(a) where the size scale has been normalized to the free space wavelength. Scattering was simulated using FDFD where a plane wave was incident from the left. The total-field is depicted in Fig. 6(b), while Fig. 6(c) shows the scattered field. From the scattered-field data, the pattern in Fig. 6(d) was computed to illustrate the preferred directions of scattered energy. Each lobe corresponds to a preferred direction.

 

Fig. 6. (a) Refractive index distribution of scattering object. Here the object is a perfect cylinder with n=2.0. (b) (Movie 2403kb) Total-field computed by FDFD simulation. [Media 1] (c) Scattered-field computed by FDFD simulation. (d) Pattern of energy scattered from object. Lobes correspond to preferred directions of scattering.

Download Full Size | PDF

Using cepstral filtering, the homogenous cylinder was then modified to suppress scattering in all but the forward direction. The new object and its scattering behavior is depicted in Fig. 7. Direct comparison between Fig. 6(d) and Fig. 7(d) demonstrates the performance of our design method. Forward scattering clearly dominates the angular spectrum and scattering in all other directions is suppressed by over 95%. The circular symmetry, which we have chosen for the design problem further ensures that the response of the structure does not depend on the direction of the incident field.

 

Fig. 7. (a) Refractive index distribution of filtered scattering object. (b) (Movie 2219kb) Total-field computed by FDFD simulation. [Media 2] (c) Scattered-field computed by FDFD simulation. (d) Pattern of energy scattered from object. Lobes correspond to preferred directions of scattering.

Download Full Size | PDF

The performance of the design we obtained with the cepstral method is further highlighted by comparing the angular spread in Fig. 7(d) with that of Fig. 8(b) which shows the spread computed for the case of a perfect plane wave. The finite width of the lobe in Fig. 8 is caused by the finite window size used to compute the scattering pattern. Window size was limited by the memory requirements of the rigorous diffraction model. The similar width of the lobes in Fig. 7 and 8 thus allows us to conclude that within the bounds of accuracy provided by the rigorous diffraction model, the design in Fig. 7(d) indeed shows only the desired behavior, i.e. a single scattered field mode in forward direction.

We emphasis again, however, that the cepstral algorithm will provide only an approximate solution to the design problem. The comparison between Fig. 6(c) and Fig. 7(c) clearly illustrates the improvement with respect to the desired response. In Fig. 6(c) the scattered near field clearly separates into the three dominant components that can be identified in the far field pattern. In contrast, Fig. 7(c) can be interpreted as the scattered field distribution merging slowly into a single plane wave mode in the forward direction. Nevertheless, the near field in 7(c) clearly shows aberrations from the perfect shape of a plane wave. In other words, if a sufficiently large simulation window is used to compute the far field patterns, the aberrations would manifest as small side lobes very closely spaced to the main forward side lobe.

We also note that while our investigation is aimed at designing electromagnetic cloaks, the corresponding design problem is different from the approach discussed in Refs. [12–17]. The latter proposes an electromagnetic structure which is guiding the electromagnetic wave around the objects to be phase matched with the incident wave in the near field. In contrast, our method is only concerned with the magnitude of the scattered far field. This includes as acceptable solutions scattered plane waves which propagate in forward direction, but which are phase shifted relative to the incident field.

 

Fig. 8. Angular spread of the far field for plane wave. The finite simulation window for the near field in (a) results in a finite angular spread of the far field pattern predicted by the numerical diffraction model.

Download Full Size | PDF

6. Conclusions

In this paper we presented the reconstruction of strongly scattering objects, designed to exhibit scattering behavior that generates the incident field (a plane wave) in the forward scatter direction. To solve the synthesis problem, an inverse scattering method based on homomorphic filtering was adopted. This nonlinear filtering technique allows one to recover an estimate of the scattering object V(r) from the image of V ∙ 〈Ψ〉 obtained using conventional diffraction tomography. The quality of the reconstruction depends on a number of factors, including the quantity and quality of the data made available. Note that simple low-pass filtering of the spectrum of V ∙ 〈Ψ〉 as opposed to the cepstrum, will do nothing.

In principle, when solving the scattering-object synthesis problem, one can specify the entire volume of k-space provided that choice corresponds to a physically meaningful set of scattered field data. In the examples presented here, annuli of zeros were imposed in the k-space of a known scatterer corresponding to zero scattering being required over a range of scattering angles and wavelengths. Inversion of the resulting k-space data set generates object structures which are modifications of the original object from which the data was originally acquired. The synthesized structures were validated using the FDFD method to rigorously compute their scattered-fields.

Our interest for this investigation was particularly focused on the design of invisible scatterers, which could act as electromagnetic cloaks. The design method does not result in a structure which suppresses the scattered field, but which scatters only in the forward direction. Rigorous simulation of the designed structure confirms the desired response within the accuracy provided by the simulation model. An isotropic response is achieved by exploiting the cylindrical symmetry of the scattered field data used as a starting point of the algorithm and the spatial filter in the cepstral domain. While it is clear that our results for designing invisible structures of circular symmetry cannot easily be realized, we note that this constitutes the first attempt where our method was used to design structures with a response that is specified over the entire plane wave spectrum, rather than for a small number of discrete scattered field angles. This suggests a large variety of potential applications for cepstral filtering as a design method. As one further example we mention objects, which emulate scattering from a quite different shape, i.e. shape-shifters. While the fabrication and experimental demonstration of both electromagnetic cloaking as well as of shape-shifting objects will require additional progress, cepstral filtering provides us with a unique framework for addressing the nonlinear relationship between the desired response and the structure of the scatterer.