A fast inversion method for highly conductive submicron wires on a substrate

Mirza Karamehmedović mirza@iwt.uni-bremen.de Dept. of Process and Chemical Engineering, University of Bremen, Badgasteiner Str. 3, 28359 Bremen, Germany Poul-Erik Hansen peh@dfm.dtu.dk Danish Fundamental Metrology, Technical University of Denmark, Matematiktorvet 307, DK-2800 Kgs. Lyngby, Denmark Thomas Wriedt thw@iwt.uni-bremen.de Dept. of Process and Chemical Engineering, University of Bremen, Badgasteiner Str. 3, 28359 Bremen, Germany


INTRODUCTION
Optical characterisation of submicron structures, for example in the production and quality control of functional materials, typically involves the solution of inverse scattering problems.A commonly used scheme for the numerical solution of such problems is the Kirsch-Kress Method (KKM), see Kirsch and Kress [1] and Colton and Kress [2] [ Section 7.3].Here, essentially, the scattered near field is reconstructed from the given data using a forward scattering model, whereafter the position and shape of the scatterer are estimated by minimising a boundary condition error within a class of admissible boundaries.This can be performed in a single, Tikhonov-regularised optimisation, or by conducting two separate, consecutive optimisations, the first ill-posed and the second nonlinear.Because of the separation of the ill-posed and the nonlinear elements of the inversion, the KKM is a decomposition method.
The reconstructed scattered near field is typically expressed in terms of radiation integrals, and since this field is evaluated many times during the boundary-error minimisation in KKM, the efficiency of the computation of the radiation integrals greatly influences the overall efficiency of the inversion.In the inversion of structures on substrates, the efficiency of the forward model is especially important, since the scattered near field in this case is given by radiation integrals that involve a half-space (or half-plane) Green's function.Such Green's functions, in turn, involve Sommerfeld-type integrals (see, e.g., Felsen and Marcuvitz [3] or Rahmat-Samii et al. [4]), and are therefore generally computationally expensive.To speed up the inversion, the forward model can be cast in terms of the so-called Method of Auxiliary Sources (MAS), where the surface current distributions in the radiation integrals are approx-imated by Dirac delta sources (line currents, dipoles etc.)For a description of the MAS, see [5]- [10] and the references therein.
The MAS was already used in conjunction with KKM in the numerical solution of inverse problems.For example, Tal and Leviatan [11] described a procedure for the reconstruction of the position and shape of PEC cylinders in two-dimensional free space, based on the amplitude and phase of the scattered far field.In the same paper, Tal and Leviatan presented numerical results for circular and elliptic cross section PEC cylinders, with the scattered field available in all directions of observation.Later, Lin and Kiang [12] used a combination of MAS and a Tikhonov-regularised KKM to estimate the position and shape of more complicated PEC cylinders in twodimensional free space, based on the amplitude and phase of the scattered near field given in all directions of observation.Obelleiro et al. [13] reconstructed boundaries of twodimensional PEC objects using generalised multipoles and based on the scattered near field sampled at a linear array of points.More recently, Alves and Martins [14] found the shape and location of inclusions and cavities in elastic bodies from the Cauchy data given at an external boundary.
In Karamehmedović et al. [15], we made an extensive numerical study of an efficient MAS-based model for forward scattering by PEC and penetrable submicron wires on substrates (see Eremina et al. [16,17] for a related model for gold and silver spheroidal particles on glass.)Our purpose here is to use the model of [15] to construct an efficient inversion scheme of the Kirsch-Kress type for sizing of such PEC and highly conductive penetrable wires.As a test problem, we estimate the cross-section radius of a PEC or silver (Ag) penetrable wire positioned on a semi-infinite silicon (Si) substrate.The measured data are simulated using the commercially available FEM solver COMSOL Multiphysics® [18,19].Thus the synthesis of the data is done using a Finite Element Methodbased scheme, and the inversion is done using a layer potential method -the underlying formulations being entirely different, we avoid so-called inverse crimes.In practical characterisation of submicron structures on surfaces, the scattered far field is typically measured only over a small aperture, rather than being available in all directions of observation.For PEC wires, we therefore limit the available far-field data to a 45 • aperture (or, due to the symmetry of the problem, a 90 • aperture).At the moment, we require a larger aperture for Ag wires.
The inverse problem is described in Section 2. In Section 3, the numerical method used for the solution of the inverse problem is specified.The results are presented and discussed in Section 4. Finally, conclusions and suggestions for future work are given in Section 5.

THE INVERSE SCATTERING PROBLEM
Figure 1 illustrates the measurement setup behind the considered inverse scattering problem.A PEC or Ag submicron wire of circular cross section Σ with boundary σ of radius r σ is positioned on a planar, semi-infinite Si substrate and illuminated by a time-harmonic, uniform plane wave of transverse electric (TE) polarisation and unit amplitude.The plane wave propagates in the normal direction on the substrate, and the operating free-space wavelength λ 0 is 325 nm.The time-dependence factor, in the following suppressed for sim-plicity, is e jωt .According to Palik [20], the complex refractive indices of silver and silicon are n Ag = 0.571 − j0.636 and n Si = 5.05506 − j3.20418, respectively, at the operating wavelength.A bistatic measurement of the scattered far field is done in an angular range [0 • , Φ], as shown in Figure 1 (note the symmetry of the problem.)The maximal measurement angle Φ from the direction of normal incidence is fixed between 0 • (measurement in the forward direction only) and 90 • (measurement in the whole available range above the substrate.) The inverse problem is to estimate the radius r σ of the wire boundary based on the knowledge of the incident field and on the measured far field.It is assumed known where the wire is centered, and the cross-section radius r σ is assumed to be within some specified interval [a, b].In the following, we use the Cartesian coordinate system (x, y, z) shown in Figure 1 and with the origin at the contact point between the substrate and the wire.Note the unit vectors x, y and z pointing in the respective coordinate directions.

THE NUMERICAL METHOD
The inverse problem of Section 2 is here solved numerically in two steps.The first, ill-posed step, minimises a functional of the form where E s f ar is the electric far field produced by the forward scattering model, C ∈ C N is a vector of complex values used in the forward model, as explained shortly, and where is the usual L 2 -norm on the interval [0, Φ].Once the optimal vector C is found, the near field associated with the measured far field E measured f ar can be approximately reconstructed.
In the second, nonlinear step, the reconstructed near field is probed within the domain that is expected to be occupied by the wire cross section.Recall that the positive numbers a and b are the expected lower and upper bounds, respectively, on the wire cross-section radius.
To describe how the objective functional f is computed, we need to briefly recapitulate the forward model for PEC particles on substrates examined in Karamehmedović et al. [15].In Figure 1, the incident field propagates in the negative y direction, and it is given by where k 0 = 2π/λ 0 is its wavenumber and η 0 ≈ 377Ω is the free-space impedance.With no wire present, the incident field is reflected (and refracted) at the air-substrate interface γ, giving rise to a reflected uniform plane wave (E re f l , H re f l ) = ( zE re f l , xH re f l ) = Γ r ( ze −jk 0 y , xη − 0 e jk 0 y ) in the lower half-plane and propagating in the negative y direction.
Here, the Fresnel reflection coefficient Γ r and transmission coefficient Γ t for normal incidence upon the air-substrate interface at 325nm are given by respectively, with n 0 = 1 being the refractive index of air.
Write E tot = zE tot for the total electric field with the wire present.In the PEC wire case, the total tangential electric field vanishes at the wire cross section boundary σ.In the penetrable case (Ag wire), the total tangential electric and magnetic fields are continuous across the boundary σ.In both cases, the tangential components of the total electric and magnetic fields are continuous across the air-substrate interface γ.The incident field satisfies (∆ + k 2 0 )E i = 0 in R 2 and it is continuous with continuous normal derivative across the boundaries γ and σ, so the scattered field E s = E tot − E i solves the inhomogeneous Helmholtz problem lim Here k Si = n Si k 0 is the wavenumber in the half-plane filled with the substrate, the vector ν is the unit outward normal to the cross section Σ, and the subscripts σ+ and σ− signify limit values at σ taken from the exterior R 2 + \ Σ and from the interior Σ of the wire cross section, respectively.A major point now is that, for the sizing of penetrable submicron wires presented in this paper, we assume the wire is sufficiently conductive for the total electric field inside its cross section Σ to be negligible.Thus, in the penetrable case, we here assume Equation ( 6) is satisfied and replace the transmission conditions (8) and ( 9) by the PEC Dirichlet boundary condition (7).Also, we wish to work with a scattered field expressed in terms of a single layer potential over the wire boundary σ, and are therefore interested in reformulating the problem ( 4)- (11) in terms of a problem with no inhomogeneous terms in the Helmholtz equations for the lower and upper half-planes.Defining a new incident field by z E i = E i = E i + E re f l for y > 0 and E i = E trans for y < 0, the corresponding scattered field E s = E tot − E i satisfies the problem Si ) E i = 0 for y < 0, and since E i is continuous with continuous normal derivative across the interface γ (due to 1 + Γ r = Γ t ).Let Φ 1/2 (•, x ) be the outgoing fundamental solution of the Helmholtz operator in the air-substrate medium, with singularity at x and continuous with continuous normal derivative across the interface γ.Then a solution of the problem ( 12)-( 16) is given by the single layer potential provided the surface current distribution g on the crosssection boundary σ satisfies the integral equation Introduction of the incident field E i (and, in the penetrable wire case, the PEC approximation) thus simplifies the original problem ( 4)-( 11) to a single integral equation.In the upper half-plane, the Green's function Φ 1/2 is given, e.g., in Eqn. ( 54)-(54c), Section 5.5e, p. 527 of Felsen and Marcuvitz [3] (here reproduced with our time-dependence factor e jωt and in our coordinate system): for x = (x, y) and x = (x , y ), with x, x ∈ R and y, y > 0.
Here, for each real η, The half-plane Green's function Φ 1/2 is generally numerically expensive, due to the presence of the integral in (19) accounting for the interface effect.Equation ( 18) is now discretised in terms of the Method of Auxiliary Sources, and using a particular approximation of the half-plane Green's function Φ 1/2 .Assume the bounds a and b on the cross-section radius r σ are given, and refer to Figure 2. A so-called auxiliary surface κ , conformal to the smallest expected wire cross-section, is fixed.
In the present case, κ is simply a circle of radius r < a and centered at (x, y) = (0, a).In the numerical examples presented here, we always set r = 0.86a.This choice of the auxiliary surface is based on numerical experimentation, and we do not claim that it is optimal.The scattered electric field E s in the exterior of the yet unknown cross section Σ in the upper half-plane is approximated by a finite linear combination of fields radiated by z-directed electric line currents (so-called auxiliary sources) and their reflections in the silicon substrate, Here, C ν , ν = 1, . . ., N , are unknown complex amplitudes of the auxiliary sources, and H (2) 0 is the Hankel function of zero order and second kind.The points x 1 , . . ., x N are the locations of the auxiliary sources, and for each air air point x ν = (x ν , −y ν ) is the reflection of x ν with respect to the x axis.See Figure 3 for an illustration.The auxiliary sources are uniformly distributed along the auxiliary surface κ .Note that, since H (2) 0 (k 0 | • |) is proportional to the radial outgoing fundamental solution of the two-dimensional Helmholtz operator in free space, the auxiliary sources and their reflections in (21) radiate in free space, i.e., with wire and substrate absent.This approximation of the half-plane Green's function for the air-substrate system by reflections is the only approximation of the interaction between the metallic cylinder and the substrate in the model.The vector C = (C 1 , . . ., C N ) ∈ C N of auxiliary source amplitudes is found by minimising the discrepancy, with respect to the L 2 ([0 • , Φ])-norm, between the measured electric far field E measured f ar and the electric far field E s f ar (C) produced by the MAS finite linear combination (21).In the numerical implementation, the far fields are sampled at a finite number of testing points x 1 , . . ., x M located a distance R from the origin and uniformly distributed within an angular interval [0 • , Φ], see Figure 3.The functional f to be minimised is thus proportional to Using the large-argument asymptotic expressions for the Hankel function H (2) 0 , we obtain a numerically more efficient, far-field form of the functional f , namely Here, the angles φ µ and φ ν (both within the interval [0, π]) are the azimuthal coordinates of the testing points and of the auxiliary sources, respectively, in the cylindrical (ρ, φ, z) coordinate system associated with the Cartesian system of Figure 1.
Minimising (22) or (25) with respect to the auxiliary source amplitudes C ν is a linear least-squares problem, and we here solve it by solving the associated normal equations; for example, with T µν = H (2)   Because all fields considered here are TE-polarised, the electric intensity vector is parallel with the wire surface, and thus the PEC boundary condition error is simply the total electric field.The symmetry of the problem is used to halve the number of actual near-field observation points.In all cases the near field is computed on a 160 × 160 grid.The parameters of the forward model used to obtain figures 4-6 were chosen with near-field quality, and not computational efficiency, in mind.The near-field computation takes approximately 260 sec for each of the figures.A total of 40 auxiliary sources are used, and the far field is sampled at 20 testing points in the angular range [0 • , 75 • ] from normal incidence.The radius of the auxiliary surface κ is set to 0.86r σ .The locations of the auxiliary sources radiating the approximation of the scattered near field E s are clearly visible as 'bright spots' in the figures.
The quality of the reconstructed near field seems highest in the direction of forward scattering (that is, along the positive y axis,) which is unsurprising since the illumination (E i , H i ) is normally incident and the far-field sampling, due to the symmetry of the problem, is centered about the y axis.The loci of the boundary condition error minima, as well as the phase jumps in the boundary condition error, give a rough reproduction of the top half of the wire cross section.The picture seems less sharp in the penetrable case, where the wire material has only finite conductivity, and where the PEC approximation is less adequate, especially for electrically small scatterers.
Since it is well-known that the wires have circular cross section, a natural approach would now be to minimise, in the least squares sense, the boundary condition error sampled within a family of circular arcs with some fixed central angle β 2 − β 1 .In this paper, we consider the numerically inexpensive special case of sampling along the y-axis The position of the global minimum in the sampled boundary condition error is interpreted as the approximate height of the wire.The resulting boundary error minimisation turns out to be relatively well-posed, especially for PEC wires.An example is given in Figure 7.For sufficiently high quality of the reconstructed near field, we did find it possible to obtain good estimates of the wire cross section radius also when testing the boundary condition error along circular arcs with central angle greater than zero.However, we experienced that the approach with β 2 − β 1 > 0 does not improve the accuracy of the inversion.Indeed, we observed that when the far-field observation aperture is narrower, and when the parameters of the forward scattering model are optimised for speed, the problem of minimising the boundary condition error over circular arcs can be highly ill-posed, even for PEC wires.
The approximation of the half-plane Green's function by reflections, as given in (21), is crude and can be improved by implementing well-known, numerically inexpensive modifica- tions.First, following Section 5.5e of Felsen and Marcuvitz [3], we can arrive at a far-field approximation of the interaction term in Φ 1/2 (•, x ) (see Eqn. ( 19)) valid in the upper half-plane and given by Here, the angle-dependent Fresnel reflection coefficient is defined by and α is the angle between the normal to the air-substrate interface and the straight line segment connecting the image source point x and the observation point x (as already defined, the angles φ and φ are the azimuthal coordinates of the points x and x , respectively.)It is readily seen that Γ r = Γ(0), and that our approximation by reflections is only accurate for small values of the angle α.
Second, the exact image of a time-harmonic, TE electric line source of constant amplitude and located at x is an electric z-directed current located in complex space and radiating a field proportional to (see, e.g., Lindell and Alanen [21,22]) Since |n 2 Si /n 2 0 | ≈ 35.8 is relatively large, this image can be approximated (see Section III of Lindell and Alanen [21]) by an electric z-directed line source phase shifted π radians relative to the original source, and located at the complex depth −y + 2jk −1 0 (n 2 Si − 1) −1/2 ≈ −y +(-0.03+j0.04)λ0 .
A detailed analysis of the inversion using the above alternative expressions for the half-plane Green's functions is outside the scope of this paper.We note, however, that our numerical experiments have shown that -at least when the far-field observation points x µ are restricted to a narrow angle interval from normal -the use of these approximations results in significantly higher condition numbers of the matrix T in (26), that is, Consequently, the computed amplitudes of the auxiliary sources are typically large, and, as evident in figures 4 and 6, this distorts the reconstructed near field in a region around the sources.The distortion can result in poor accuracy of the overall inversion method.In contrast, the reconstruction in Figure 8 is obtained using the freely available Ipopt global optimisation software [23].The optimisation was constrained, with the set of permissible values of C ν and C ν being the interval [−1V/m, 1V/m].Even though the resulting inversion method turns out to have a much smaller dynamic range (while being numerically much less efficient than simply solving the normal equations), it is seen that keeping the amplitude of the sources low improves the overall quality of the reconstruction of the near field.Figure 8 is here included for illustration only, and we leave a thorough investigation of inversion based on a global nonlinear optimiser to future work.Returning to the alternative approximations of the half-plane Green's function, the reason for the increase in the condition number of T is probably that, for sources close to the air-substrate interface, the far-field phase of these approximations varies more slowly with the angle of observation than is the case for the simple approximation by reflection.The slow variation in the phase namely diminishes the variation between the components of the columns of the system matrix T that correspond to sources close to the interface.As an example, Figure 9 shows the far-field phase for the combined unit electric line source with support at (x , y ) = (±0.172,0.2)λ 0 .This combined source is situated at the auxiliary surface used by our forward scattering model in the PEC-case inversion presented in Section 4. The forward model also includes sources that are closer to the airsubstrate interface.The phases in Figure 9 are computed using the constant reflection coefficient approximation, the variable reflection coefficient approximation and the complex image approximation.While the three approximations produce almost the same amplitude of the far field in the angular range Φ ∈ [0 • , 60 • ], the far-field phase produced by the constant reflection coefficient approximation clearly varies the most.
Using the Singular Value Decomposition (SVD) can reduce the condition number of T significantly, but an increase in the smallest acceptable singular value also results in a poorer approximation of the given far field.Our numerical experiments have not yet shown an advantage in using the SVD in conjunction with the correct far-field expressions, as opposed to solving the normal equations for the constant reflection coefficient approximation.

NUMERICAL RESULTS
In the following, the total number of far-field testing points is formally set to M = 40, and the number of auxiliary sources is formally N = 10 in the PEC case and N = 20 in the Ag case.Due to the symmetry of the problem with respect to the y axis, this gives overdetermined systems of 20 equations with 5 and 10 unknown auxiliary source amplitudes C ν , respectively.The inversion is implemented in MATLAB® and run on a single 2.66GHz standard PC core, with 4Gb memory.All far-field quantities are here evaluated at a distance of 12.5λ 0 from the wire.The reconstructed near field is sampled at approximately 170 (PEC case) and 90 (Ag case) equidistantly distributed points along the straight line segment [0.86a, b] on the y axis.
Results of the inversion are shown in tables 1 and 2. The tables show actual wire radii, estimates after inversion and the associated relative errors and execution times.The actual crosssection radii r σ and estimated radii r are given in wavelengths λ 0 .The estimates r are rounded to four decimal places, while the relative error ε = 100(r − r σ )/r σ , the relative far-field error and the execution time are rounded to two decimal places.For the results of tables 1 and 2, the form ( 22) of the functional f is used.As expected, the error in the PEC estimates is lower than that in the Ag case.For the considered inversions, the absolute value |ε| of the relative inversion error is less than 1.7% in the PEC case and less than 7% in the Ag case.In all cases, the total execution time is less than 1 sec, and it is mainly spent on the second step of the inversion, i.e., on sampling of the reconstructed near field along a subset of the yaxis.In the PEC case, the inversion works well over the whole test range of wire radii; however, at the moment it proves generally inaccurate for penetrable (Ag) wires of radius less than approximately 0.5λ 0 , even if the number of auxiliary sources is increased.Also, it seems necessary to provide the far-field data in a larger angular range for penetrable wires than for PEC wires.
In the PEC wire case, using a total of 8 auxiliary sources and the far-field approximation (25) of the objective functional f in the first step of the inversion yields an average value of |ε| of 3.24%, and a total execution time of 0.64 sec per inversion.For the penetrable wires, using a total of 20 auxiliary sources and the above far-field approximation results in an average value of |ε| of 5.11%, with execution time 0.80 sec per inversion.
Finally, for inversion using the variable Fresnel reflection coefficient Γ(k Si cos α), we found a good choice of the SVD tolerance to be 10 −4 .Over the same dynamical ranges as above, this choice results in the average absolute value of the relative error ε of 6.94% in the PEC case and 5.96% in the Ag case.The execution time is approximately 0.96 sec for each inversion.

CONCLUSION AND FURTHER WORK
A numerical inversion method was proposed for sizing of perfectly electrically conducting and highly conductive penetrable submicron wires on substrates.The method is essentially a Kirsch-Kress inverse scheme with the forward scattering model formulated in terms of the Method of Auxiliary Sources.The proposed scheme was applied to the sizing of PEC and silver submicron wires of circular cross section on a silicon substrate, with numerically simulated scattered far field available in a 45 • aperture (PEC case) and 90 • aperture (Ag case).The operating free-space wavelength was λ 0 = 325nm, and the range of wire radii for which the inversion was accurate was found to be at least 0.2λ 0 −1.3λ 0 in the PEC case, and at least 0.8λ 0 −1.3λ 0 in the Ag case.In this range, relative errors of less than 1.7% and less than 7% were achieved in the PEC and the penetrable case, respectively.The average relative errors were 0.62% and 4.92%, respectively, and all inversion times were shorter than 1 sec.Thus, at least for perfectly electrically conducting wires, the inversion scheme can accomodate a wide range of scatterer sizes using a single set of parameters.
Futher work includes investigation of the performance of the proposed inversion method in the two-dimensional transverse magnetic (TM) case and in the three-dimensional case.The types of auxiliary sources that might prove well-suited in these cases are transverse magnetic line currents and crossed Hertzian dipoles, respectively.An analysis of the sensitivity of the method to noise in the measurement data is also needed.It will be highly relevant to adapt the method to applications with further reduced measurement data, e.g., where only the scattering intensity is given.The minimisation of the boundary condition error was here done for simplicity by 'brute force,' that is, by simply evaluating the quantity at a large number of points along the y axis.This most of the total execution time, so implementing a faster minimisation scheme could significantly improve the efficiency of the overall inversion.The plots of the reconstructed near fields show not only the position of the top of the wire, but also indicate the shape of the illuminated part of the scatterer.This motivates the use of the presented forward scattering model in shape reconstruction of particles (on substrates) with complex morphologies described by multiple parameters, where the use of look-up tables can be impractical due to the large size of the tables.We expect improved quality of reconstruction if scattered far-field data are given for more than one angle of illumination or operating frequency.Instead of using the PEC approximation for highly conductive penetrable scatterers on substrates, better inversion accuracy might be achieved if the original transmission problem is approximated by an exterior impedance boundary value problem (see Colton and Kress [2, Section 9.5]).Finally, our choice of the simple farfield approximation of the half-plane Green's function regularised the inverse problem somewhat and, with the forward model used here, improved the accuracy of the inversion.It is of great interest to further investigate regularisation schemes that would diminish the distortion in the reconstructed near field.Here, the concrete setup of the inverse problem as well as the forward model should be taken into account.

FIG. 1
FIG.1Cross section of a wire on a silicon substrate.

FIG. 3
FIG.3The MAS forward scattering model: reconstructing the scattered near field.
Figures 4-6 illustrate the quality of the resulting reconstruction of the scatterer.The figures show the absolute value and the phase of the PEC boundary condition error

FIG. 9
FIG.9The phase of the electric far field with the forward scattering model using three different approximations for the half-plane Green's function.The support of the combined source is (x , y ) = (±0.172,0.2)λ0 .