W A Iff

J-P Hugonin

C Sauvan

M Besbes

P Chavel

G Vienne

L Milord

D Alliata

E Herth

P Coste

A Bosseboeuf

Electromagnetic analysis for OCT-based Through Silicon Via (TSV) metrology

This article reports on progress in the analysis of Time-Domain Optical Coherence Tomography (OCT) applied to the dimensional metrology of Through Silicon Vias (TSV, vertical interconnect accesses in silicon, enabling 3D integration in micro-electronics) and estimates the deviations from earlier, simpler models. The considered TSV structures are 1D trenches and circular holes etched into silicon with a large aspect ratio. As a prerequisite for a realistic modelling, we work with spectra obtained from reference interferograms measured at a planar substrate, which fully includes the dispersion of the OCT apparatus. Applying a rigorous modal approach, we estimate the differences to a pure ray tracing technique. Accelerating our computations, we focus on the relevant fundamental modes and apply a Fabry-Perot model as an efficient approximation. Exploiting our results, we construct and present an iterative procedure based on the minimization of a merit function, which concludes TSV heights reliably, accurately and rapidly from measured interferograms.

Introduction

With sales of semiconductor equipment breaking new records, the development and production of more compact electronic chips requires progress in optical metrology. One promising way to increase the compactness is 3D integration based on Through Silicon Vias (TSVs, see fig. 1) since it offers superior integration density and reduces interconnect lengths [START_REF]International Technology Roadmap for Semiconductors (ITRS), Assembly and Packaging[END_REF][START_REF] Wu | High aspect ratio silicon etch: A review[END_REF]. Accurate measurements are needed to monitor the depth uniformity of etched TSVs [START_REF] Liu | A 300-mm wafer-level three-dimensional integration scheme using tungsten through-silicon via and hybrid Cu-adhesive bonding[END_REF]. When TSVs are filled with metal (copper), their geometrical variations can affect the coplanarity and can warp the wafer, resulting in a low stacking yield. The accurate non-destructive measurement of TSV depths is an increasingly challenging task since the TSV diameter has shrunk to only a few microns in many cases.

Several optical techniques can be used for measuring TSV depths, including spectral reflectometry [START_REF] Ku | Reflectometer-based metrology for high-aspect ratio via measurement[END_REF][START_REF] Ku | Characterization of high density through silicon vias with spectral reflectometry[END_REF][START_REF] Ku | In-line metrology of 3D interconnect processes[END_REF][START_REF] Ku | Spectral reflectometry for metrology of three-dimensional through-silicon vias[END_REF][START_REF] Fursenko | In-line through silicon vias etching depths inspection by spectroscopic reflectometry[END_REF][START_REF] Fursenko | Through silicon via profile metrology of Bosch etching process based on spectroscopic reflectometry[END_REF][START_REF] Fursenko | 3D Through Silicon Via profile metrology based on spectroscopic reflectometry for SOI applications[END_REF][START_REF] Bauer | Very high aspect ratio through silicon via reflectometry[END_REF], backside infrared spectral reflectometry [START_REF] Marx | Wafer Thickness Sensor (WTS) for etch depth measurement of TSV[END_REF][START_REF] Teh | Backside Infrared Interferometric patterned Wafer Thickness Sensing for Through-Silicon-Via (TSV) etch Metrology[END_REF][START_REF] Gambino | Through-silicon-via process control in manufacturing for SiGe power amplifiers[END_REF], ellipsometry (for rather shallow trenches) [START_REF] Diebold | Perspective: Optical measurement of feature dimensions and shapes by scatterometry[END_REF], confocal chromatic microscopy [START_REF] Vienne | A Comparison of TSV Etch Metrology Techniques[END_REF] and time domain OCT [START_REF] Vienne | A Comparison of TSV Etch Metrology Techniques[END_REF][START_REF] Drexler | Optical Coherence Tomography[END_REF][START_REF] Alliata | Progress in Time-Domain Optical Coherence Tomography for TSV / 3Di stacking Metrology[END_REF][START_REF] Devanciard | Through Silicon Via Process Characterization by Integrated Inspection/Metrology Solutions in Visible and Infrared Domain[END_REF], and hybrid systems [START_REF] Devanciard | Through Silicon Via Process Characterization by Integrated Inspection/Metrology Solutions in Visible and Infrared Domain[END_REF][START_REF] Ahn | A Hybrid Non-destructive Measuring Method of Three-dimensional Profile of Through Silicon Vias for Realization of Smart Devices[END_REF]. Time domain OCT, whose principle is sketched in fig. 2, has its strengths in the case of deep TSV of high aspect ratio. It allows scanning a very large range along the optical axis (z-axis), whereby the resolution in z-direction does not depend on the size of the scanned z-interval contrary to spectral reflectometry [START_REF] Ku | Reflectometer-based metrology for high-aspect ratio via measurement[END_REF][START_REF] Ku | Characterization of high density through silicon vias with spectral reflectometry[END_REF]. The interpretation of timedomain OCT interferograms is unambiguous: The order of the detected signals agrees with the order of the reflecting objects along the optical axis. Additionally, OCT allows an economic exploitation of the light source power since the reference mirror has a reflectivity of almost 1; by comparison, spectral reflectometry uses the silicon surface as a reference mirror, which has a reflectivity in the order of 0.3 (in the case of through glass vias this is reduced to a mere 0.04).

OCT has got a sufficient lateral resolution so that it can be utilized not only for the measurement of groups of TSVs but also for the measurement of a large part of the single TSVs of nowadays size [START_REF] Vienne | A Comparison of TSV Etch Metrology Techniques[END_REF][START_REF] Drexler | Optical Coherence Tomography[END_REF][START_REF] Alliata | Progress in Time-Domain Optical Coherence Tomography for TSV / 3Di stacking Metrology[END_REF][START_REF] Devanciard | Through Silicon Via Process Characterization by Integrated Inspection/Metrology Solutions in Visible and Infrared Domain[END_REF]. In this publication, we concentrate on single TSVs; the considered TSV structures are 1D trenches (rectangular grooves, which are translation invariant in one direction) and circular holes (cylindrical structures); the sidewalls will be vertical and the bottoms will be flat. The measurement technique in this paper is connected with a focus on the TSV top and an averaging of the measurement beam over the TSV since the most important fundamental mode also averages over the TSV when propagating between TSV top and bottom. For the future, a simultaneous measurement of a group of TSVs -while resolving the single TSVs -could be an option. This could be accomplished by full-field OCT (FFOCT), which allows also a high lateral resolution [START_REF] Safrani | Spatial coherence effect on layer thickness determination in narrowband fullfield optical coherence tomography[END_REF][START_REF] Safrani | Ultrahigh-resolution full-field optical coherence tomography using spatial coherence gating and quasi-monochromatic illumination[END_REF]. Fig. 1. Schematic of an etched TSV (white color) before filling with copper. The TSV is characterized by its shape (1D trench or circular hole in this paper, its height 𝐻 and its width ∆𝑥 or diameter 𝐷); the figure shows the TSV cross section. Aspect ratios of 𝐻 𝐷 = 10 ⁄ are typical nowadays. At later process steps, the TSV walls are passivated, the TSV is filled with copper and the remaining silicon below the TSV is etched away. Fig. 2. Top: Scheme of the studied time domain OCT device for TSV measurements [START_REF] Drexler | Optical Coherence Tomography[END_REF][START_REF] Alliata | Progress in Time-Domain Optical Coherence Tomography for TSV / 3Di stacking Metrology[END_REF]. The motion of the reference mirror allows longitudinal resolution and the objective determines the lateral resolution. The detector records the intensity versus the z-position of the reference mirror. The top and the bottom of the TSV have the coordinates 𝑧 𝑡 and 𝑧 b (subscript "t" for "top", "b" for "bottom"), its height is 𝐻 = 𝑧 b -𝑧 t ; the z-position of the top silicon interface in the neighborhood of the TSV, where a reference interferogram is recorded (see section 2), is 𝑧 s ≈ 𝑧 t (subscript "s" for "substrate"); the TSV will be embedded into a larger region 𝑝 along the 𝑧-direction later and periodized for a numeric treatment in Fourier space. One group of adjacent fringessymbolized by a red sinusoidal pattern in the picturecorresponds to one reflecting interface. The incident beam is Gaussian and various objectives are available so that the beam waist at the top interface covers the entire TSV plus its surrounding; the signals coming from top and bottom of the TSV should be of the same order of magnitude. The phase of the reference signal exp[2i𝑘𝑧 + i𝜑(𝑘)] takes into account interferometer dispersion [START_REF] Drexler | Optical Coherence Tomography[END_REF][START_REF] Iff | Simulation of a Time-Domain OCT device for Through Silicon Via (TSV) fabrication[END_REF]; 𝑘 is the modulus of the vacuum wave vector. Bottom: example of two real interferograms at TSVs. Left: constructive interference of two adjacent fringe patterns; right: destructive interference.

Works on TSV height measurement by any optical technique incorporate up to our knowledge either non-electromagnetic models for light propagation inside TSV [START_REF] Ku | Reflectometer-based metrology for high-aspect ratio via measurement[END_REF][START_REF] Ku | Characterization of high density through silicon vias with spectral reflectometry[END_REF][START_REF] Ku | In-line metrology of 3D interconnect processes[END_REF][START_REF] Ku | Spectral reflectometry for metrology of three-dimensional through-silicon vias[END_REF][START_REF] Hyun | Methods to Measure the Critical Dimension of the Bottoms of Through-Silicon Vias Using White-Light Scanning Interferometry[END_REF] or time-consuming rigorous computations with the Rigorous Coupled Wave Analysis (RCWA) [START_REF] Fursenko | In-line through silicon vias etching depths inspection by spectroscopic reflectometry[END_REF][START_REF] Fursenko | Through silicon via profile metrology of Bosch etching process based on spectroscopic reflectometry[END_REF][START_REF] Fursenko | 3D Through Silicon Via profile metrology based on spectroscopic reflectometry for SOI applications[END_REF][START_REF] Bauer | Very high aspect ratio through silicon via reflectometry[END_REF]. This paper aims at an estimate and improvement of the accuracy of TSV depth measurements by application of electromagnetic but rapid calculations. Although our work is focused on time domain OCT, any optical measurement technique such as spectral reflectometry [START_REF] Ku | Reflectometer-based metrology for high-aspect ratio via measurement[END_REF][START_REF] Ku | Characterization of high density through silicon vias with spectral reflectometry[END_REF][START_REF] Ku | In-line metrology of 3D interconnect processes[END_REF][START_REF] Ku | Spectral reflectometry for metrology of three-dimensional through-silicon vias[END_REF][START_REF] Fursenko | In-line through silicon vias etching depths inspection by spectroscopic reflectometry[END_REF][START_REF] Fursenko | Through silicon via profile metrology of Bosch etching process based on spectroscopic reflectometry[END_REF][START_REF] Fursenko | 3D Through Silicon Via profile metrology based on spectroscopic reflectometry for SOI applications[END_REF][START_REF] Bauer | Very high aspect ratio through silicon via reflectometry[END_REF][START_REF] Marx | Wafer Thickness Sensor (WTS) for etch depth measurement of TSV[END_REF][START_REF] Teh | Backside Infrared Interferometric patterned Wafer Thickness Sensing for Through-Silicon-Via (TSV) etch Metrology[END_REF][START_REF] Gambino | Through-silicon-via process control in manufacturing for SiGe power amplifiers[END_REF] could benefit from the results of our rigorous modal calculations at single TSVs based on Fourier-Bessel functions [START_REF] Popov | Factorization rules of discontinuous functions applied to Fourier-Bessel basis[END_REF][START_REF] Bonod | Differential theory of diffraction by finite cylindrical objects[END_REF] and our simplified yet accurate, more efficient models derived from it.

Firstly, we study and model a real, non-Gaussian spectrum and investigate the impact of dispersion in the interferometer arms [START_REF] Drexler | Optical Coherence Tomography[END_REF] on the OCT interferograms. Then, we study the differences between ray tracing and a rigorous calculation on 1D trenches and circular holes. With rigorous computations being time consuming, and scalar techniques being fast but inaccurate, we give the limitations of their applicability and calculate the effective indices of trenches and holes on the basis of the propagation constants of the relevant fundamental modes in order to improve the accuracy of ray tracing-based procedures. The remaining limitations of the accuracy of ray tracing for small trench widths and hole diameters motivates the consideration of the coupling coefficients between the fundamental modes and the application of the Fabry-Perot model [START_REF] Born | Principles of Optics[END_REF] to TSV: This model, which is nowadays widely used at microand nanostructures [START_REF] Lalanne | Photon confinement in photonic crystal nanocavities[END_REF][START_REF] Ditlbacher | Silver Nanowires as Surface Plasmon Resonators[END_REF], allows fast calculations for various trench and hole depths once the coupling coefficients at top and bottom are known from a rigorous computation. Having studied the accuracy and limitations of the Fabry-Perot model, we apply it for an estimate of the maximum measurable depth of trenches and holes. The results obtained in section 2 -6 are checked and exploited experimentally in section 6, where we explain and compare our novel, electromagnetic interferogram evaluation procedure based on least squares applied to complex Fourier amplitudes with the common envelope-based one in the case of OCT applied to TSV [START_REF] Courteville | Contact free on-axis metrology for the fabrication and testing of complex optical systems[END_REF][START_REF] Wilhelm | A novel, low coherence fibre optic interferometer for position and thickness measurements with unattained accuracy[END_REF][START_REF] Larkin | Efficient nonlinear algorithm for envelope detection in white light interferometry[END_REF].

This publication is structured as follows: Section 2 is devoted to the modelling of the light source spectrum and interferometer dispersion. Section 3 models the light propagation in trenches and holes accurately. Section 4 presents the effective indices in trenches for the fundamental modes, applies the Fabry-Perot model to TSV and estimates the errors connected with this approximation. Section 5 estimates the maximum measurable depth. Section 6 presents our electromagnetic interferogram evaluation procedure and exploits our results experimentally. Section 7 concludes the work.

Accurate modelling of the spectrum of the light source and the interferometer dispersion

A prerequisite for the accurate modelling of any OCT interferogram at a complex structure is the knowledge of an OCT interferogram at a planar substrate. Having recorded such an interferogram (fig. 3), we obtain the complex spectral amplitude 𝑠(𝑘) of the light source by a Fourier transform with 2𝑘 the frequency of the fringes of the intensity pattern versus z. An ideal interferogram would be axis-symmetric and its Fourier coefficients would only consist of cosine-terms. In practice, real interferograms include the effect of dispersion in the interferometer arms [START_REF] Drexler | Optical Coherence Tomography[END_REF][START_REF] Iff | Simulation of a Time-Domain OCT device for Through Silicon Via (TSV) fabrication[END_REF]. Thus, there are also sine-terms or rather complex Fourier coefficients. The presence of dispersion leads to a larger full width at half maximum (FWHM) of the envelopes of the interferometric signals (fig. 3)in agreement with [START_REF] Hillman | The effect of water dispersion and absorption on axial resolution in ultrahigh-resolution optical coherence tomography[END_REF]. The impact of dispersion is similar to a decrease of the source spectral width and an increase of the coherence length. Therefore, a dispersion compensatione.g., in the reference arm of the interferometer -would clearly be beneficial. The spectra recorded for one source but various objectives have strictly speaking not Gaussian but complicated shapes with FWHMs of some percent of 2𝑘 c and 2𝑘 c = 2 • 2𝜋/1.329 µm the center frequency.

Rigorous calculation of the interferogram from a single etched TSV

We present the rigorous calculation of the OCT signal from a single etched TSV (circular hole or 1D trench). A prerequisite is that the spectral amplitude s(𝑘) of the broadband source used for the OCT measurement is already known from the recording of a reference spectrum, as discussed in the previous section. The OCT interferogram results from the interference between the amplitude 𝑎 TSV (𝑘) of the field reflected from the sample (see fig. 2) and the amplitude exp[2i𝑘𝑧 + i𝜑(𝑘)] from the reference arm, with 𝑘 = 𝜔/𝑐 the vacuum wave vector, 𝜔 the angular frequency and 𝑐 the speed of light. Integrating over the source spectral width [𝑘 min , 𝑘 max ], weighted by the complex spectral amplitude s(𝑘), considering the square, and time averaging yields the intensity pattern of the interferogram [START_REF] Drexler | Optical Coherence Tomography[END_REF] 𝐼

(𝑧) = 〈 (∫ s(𝑘) ⋅ {𝑎 TSV (𝑘) exp(-i𝜔𝑡) + exp[2i𝑘𝑧 + i𝜑(𝑘) -i𝜔𝑡]} 𝑑𝑘 2𝜋 ⁄ + cc 𝑘 max 𝑘 min ) 2 〉 𝑡 , ( 1
)
where 𝑧 is the position of the reference mirror (see fig. 2), cc denotes the complex conjugated expression and 〈⋅ 〉 𝑡 denotes the time average. The phase term 𝜑(𝑘) represents the interferometric dispersion due to different dispersion of the objective or the fibre propagation constants between both arms of the interferometer. In practice, we evaluate the integral over 𝑘 numerically, using few hundred integration points and midpoint sampling.

The amplitude 𝑎 TSV (𝑘) corresponds to the reflection in the far field from a single etched TSV (circular hole or 1D trench). We calculate it numerically with rigorous electromagnetic computations. The calculation is done in two steps. First, we illuminate a single via with a linearly-polarized Gaussian beam whose waist is located on the top interface and centered on the via axis and we calculate the near field, inside and close to the via. Then, from the nearfield calculation, we extract 𝑎 TSV (𝑘) by using a near-to-far field transformation and an integration over the solid angle corresponding to the numerical aperture of the collection objective.

Due to the large aspect ratio of a single via (typically in the order of 10), we use modal methods that rely on an efficient analytical integration of Maxwell's equations along the long dimension of the via, i.e., the vertical 𝑧 direction.

For the calculation of rectangular 1D trenches, we use the aperiodic Fourier modal method (a-FMM) [START_REF] Silberstein | Use of grating theories in integrated optics[END_REF]. Unlike the usual FMM (also known as RCWA) that is dedicated to light diffraction by periodic objects, a-FMM allows the calculation of light scattering by single objects. It is based on an artificial periodization of the object in the transverse 𝑥-direction with the crucial addition of a complex coordinate transform between the periods. The latter acts as a Perfectly-Matched Layer (PML) and prevents aliasing effects [START_REF] Hugonin | Perfectly matched layers as nonlinear coordinate transforms: a generalized formalization[END_REF]. TE-and TM-polarization are considered separately. Modal expansions in the different regions (air, via and silicon substrate) are connected by the S-matrix algorithm [START_REF] Chateau | Algorithm for the rigorous coupled-wave analysis of grating diffraction[END_REF][START_REF] Li | Formulation and comparison of two recursive matrix algorithms for modeling layered diffraction gratings[END_REF].

For the calculation of circular holes, we apply a modal method analogue to RCWA but adapted to body-of-revolution objects thanks to the use of Fourier-Bessel functions [START_REF] Popov | Factorization rules of discontinuous functions applied to Fourier-Bessel basis[END_REF][START_REF] Bonod | Differential theory of diffraction by finite cylindrical objects[END_REF]. This means that electric and magnetic field components, as well as permittivities and permeabilities, are developed in terms of Bessel functions in the radial direction and in terms of Fourier harmonics in the azimuthal direction. The field is expressed as 𝐄(𝑟, 𝜃, 𝑧) = ∑ 𝐄 𝑙 (𝑟, 𝑧)exp(i𝑙𝜃) 𝑙 , 𝐇(𝑟, 𝜃, 𝑧) = ∑ 𝐇 𝑙 (𝑟, 𝑧)exp(i𝑙𝜃) 𝑙 , with 𝑙 the azimuthal number. In case of an incident Gaussian beam that is linearly or circularly polarized and centered on the axis, the field contains only the harmonics 𝑙 = ±1. In addition to the method in [START_REF] Popov | Factorization rules of discontinuous functions applied to Fourier-Bessel basis[END_REF][START_REF] Bonod | Differential theory of diffraction by finite cylindrical objects[END_REF], we add a circular PML. We connect the different modal expansions by the S-matrix algorithm. Figure 4 shows the calculated near field inside and around a 1D trench in TE-( 𝐸 𝑦 component) and TM-polarization (𝐻 𝑦 component). The amplitude 𝑎 TSV (𝑘) is calculated with a standard near-to-far field transformation. We calculate the Fourier transform along 𝑥 of a crosssection of the field in air at a constant 𝑧-position a few nanometers above the via. Two important conclusions can be drawn from fig. 4. First, in the case of a 1D trench, light propagation inside the via is polarization dependent. TE-polarization is well-guided inside the TSV whereas TMpolarization leaves the TSV by larger part, meaning higher losses. Therefore, in the OCT interferogram, the signal coming from the bottom of the via will be much weaker with a TMpolarized incident beam and TE-polarization is more suitable for the measurement of deep 1D trenches. Secondly, in both polarizations, the field inside the via exhibits a clear standing-wave pattern in the vertical 𝑧-direction. A similar standing-wave pattern can also be observed in the calculation of circular holes (not shown here). We conclude that light propagation inside a via etched into silicon mainly arises from the excitation of a single mode. In the next section, we calculate the fundamental mode responsible for the standing-wave pattern and show that the problem can be accurately described by an approximate Fabry-Perot model that amounts to keep a single mode in the rigorous modal expansion.

Fabry-Perot model

In case of sufficiently deep vias, only the fundamental mode contributes significantly to light propagation inside the structure since higher-order modes in the modal expansion are more strongly damped. This motivates the application of the Fabry-Perot (FP) model [START_REF] Born | Principles of Optics[END_REF][START_REF] Lalanne | Photon confinement in photonic crystal nanocavities[END_REF], which amounts to neglect the higher-order modes and keep a single mode in the modal expansion. This approximation is all the more accurate that the via is deep. Note that, in contrast to most uses of the FP model in micro-and nanophotonics [START_REF] Lalanne | Photon confinement in photonic crystal nanocavities[END_REF][START_REF] Ditlbacher | Silver Nanowires as Surface Plasmon Resonators[END_REF], the fundamental mode that we consider here is a leaky mode whose amplitude is exponentially damped during propagation due to radiative leakage in the silicon cladding.

We first calculate in section 4.1 the fundamental mode that contributes to light propagation in trenches and in circular holes. Then, section 4.2 presents the principle of the FP model. Section 4.3 explains the interpretation of OCT interferograms on the basis of the FP model. Finally, we study in section 4.4 the accuracy of the FP model as a function of the TSV height.

Fundamental mode in trenches and holes

A key parameter in the FP model is the effective index of the mode that mainly contributes to light propagation. The effective index is defined from the mode propagation constant as 𝑛 eff ≔ 𝑘 z,0 𝑘 ⁄ with 𝑘 z,0 the propagation constant of the least attenuated mode and 𝑘 the modulus of the vacuum wave vector. It is important to keep in mind that 𝑘 z,0 and 𝑛 eff are complex numbers because of the radiative leakage in silicon. In case of a 1D trench, two different effective indices 𝑛 eff,TE and 𝑛 eff,TM have to be considered. In case of a circular hole, the least attenuated mode that is excited by a linearly or circularly polarized incident beam is characterized by an azimuthal number 𝑙 = ±1. There is a strong dependence of Re(𝑛 eff,TE ) and Im(𝑛 eff,TE ) on ∆𝑥 for small ∆𝑥, where the deviation of 𝑛 eff,TE from 1 is the largest. For ∆𝑥 → ∞, this deviation tends asymptotically to zero. The spectral dependence of silicon is weak so that 𝑛 eff,TE depends only on ∆𝑥/𝜆 with ∆𝑥 the trench width. Therefore, one can deduce 𝑛 eff,TE (∆𝑥, 𝜆) at a wavelength different from the OCT center wavelength 𝜆 c = 1.329 µm from 𝑛 eff,TE (∆𝑥 • 𝜆 c 𝜆 ⁄ , 𝜆 c ), as long as |𝜆 -𝜆 c |/𝜆 c is sufficiently small and under the assumption that the refractive index of silicon is rather independent on 𝜆 in this narrow spectral range around 𝜆 c = 1.329 µm. Qualitatively, 𝑛 eff,TM (∆𝑥, 𝜆) behaves similarly to 𝑛 eff,TE (∆𝑥, 𝜆) . However, Im(𝑛 eff,TM ) is larger, meaning larger leakage through the sidewalls. Illustratively, there is a smaller reflection for TM-polarization at planar surfaces and a larger transmission into the substratelike in case of Fresnel reflection and transmission.

Concluding, to an accuracy of 1% or better, TSVs in form of 1D trenches of widths larger than 5 µm may be calculated by scalar ray tracing, setting 𝑛 eff (𝜆) = 1, but the modelling should be rigorous for smaller 1D trench widths. In particular, the attenuation (imaginary part of the effective index) should not be neglected in the OCT data interpretation. For deep 1D trenches of small diameters, there are large losses of light at the sidewalls so that the signal from the bottom is weak and difficult to identify in the presence of noise. In that case, TEpolarization is clearly preferable to TM-polarization due to the smaller attenuation.

Figure 6 shows the |𝐸 𝑥 | 2 distribution of the fundamental mode of a circular hole for xpolarization. It is a linear combination of the both degenerate modes with 𝑙 = ±1. Figure 7 shows Re(𝑛 eff ) and Im(𝑛 eff ) as a function of the diameter 𝐷 at the center wavelength 𝜆 c = 1.329 µm. Again, there is a strong dependence of Re(𝑛 eff ) and Im(𝑛 eff ) on 𝐷 for small 𝐷, where the deviation of 𝑛 eff from 1 is the largest. For 𝐷 → ∞, 𝑛 eff tends to unity. The effective index only weakly depends on 𝜆 so that 𝑛 eff depends again only on 𝐷/𝜆. Note that Im(𝑛 eff ) is larger than for 1D trenches. Concluding, TSV in form of circular holes of widths larger than 5 µm may be calculated by a scalar method at the utilized center wavelength 𝜆 c , using the approximation 𝑛 eff = 1, but the modelling should be electromagnetic for smaller diameters.

Principle of the Fabry-Perot model and expression of the interferogram

The FP model amounts to neglect all higher-order modes inside the via region and to keep a single mode in the modal expansion. Since an air hole in silicon does not support guided modes, all modes are leaky and the fundamental mode that we consider is the least attenuated. Once this mode has been calculated, we solve two different problems, as sketched in fig. 8. First, we consider a single interface between air and a semi-infinite air hole in silicon. We calculate the amplitude 𝑎 ↑,refl (𝐻 + ) reflected at the top of an infinitely deep TSV, the amplitude 𝑎 ↓,trans (𝐻 -) transmitted into an infinitely deep TSV, the reflection 𝑟 t of the mode into itself, and the transmission 𝑡 t from the mode to plane waves in air (the subscript "t" refers to the top of the TSV). In a second step, we solve a second problem, a single interface between a semi-infinite air hole in silicon and a silicon substrate. We calculate the reflection 𝑟 b of the mode into itself (the subscript "b" refers to the bottom of the TSV). Once these few coupling coefficients are computed rigorously, the light scattered at a TSV can be computed analytically for any TSV height 𝐻 from the expression

𝑎 ↑ (𝐻) = 𝑎 ↑,refl (𝐻 + ) + 𝑎 ↓,trans (𝐻 -) exp(2i𝑘𝑛 eff 𝐻)•𝑟 b 1-𝑟 b •𝑟 t •exp(2i𝑘𝑛 eff 𝐻) 𝑡 t . (2)
Here, 𝑎 ↑ (𝐻) is the total amplitude reflected at the TSV and 𝑛 eff is the effective index of the fundamental mode. The numerator 𝑟 b 𝑡 t exp(2i𝑘𝑛 eff 𝐻) models downward propagation, reflection at the bottom, upward propagation, and transmission into air at the top of the TSV.

The denominator 1 -𝑟 b 𝑟 t exp(2i𝑘𝑛 eff 𝐻) models multiple reflections between top and bottom. The amplitude 𝑎 TSV (𝑘) is then deduced from 𝑎 ↑ (𝐻) given by ( 2) by a near-to-far field transformation. Multiple reflections are weak for sufficiently deep TSV, meaning 1 -𝑟 b 𝑟 t exp(2i𝑘𝑛 eff 𝐻) ≈ 1 (see also section 4.4), and we neglect them to further simplify the expression of 𝑎 TSV (𝑘). With this approximation, the amplitude 𝑎 TSV (𝑘) takes the simple form

a TSV (𝑘) ≈ 𝑎 t (𝑘) exp(-2i𝑘𝑧 t ) + 𝑎 b (k) exp{2i𝑘[-𝑧 t -𝑛 eff (𝑘)𝐻]} = 𝑎 t (𝑘) exp(-2i𝑘𝑧 s -2i𝑘∆𝑧 t ) + 𝑎 b (𝑘) exp[-2i𝑘𝑧 s -2i𝑘∆𝑧 t -2i𝑘𝑛 eff (𝑘)𝐻] (3)
with ∆𝑧 t = 𝑧 t -𝑧 s and 𝑧 s the z-coordinate of the top plane interface during the recording of a reference interferogram close to a TSV. The complex amplitude 𝑎 t (𝑘) represents the reflection at the top of the TSV, and the complex amplitude 𝑎 b (𝑘) the reflection from the bottom. The factor exp(-2i𝑘∆𝑧 t ) propagates light in air from 𝑧 s to 𝑧 t . Both 𝑧-positions are very close to each other but not strictly identical since the wafer is sometimes not exactly perpendicular to the optical axis and since two separate OCT measurements can have slightly different 𝑧-offsets.

The factor exp[-2i𝑘𝑛 eff (𝑘)𝐻] propagates light inside the via from top to bottom and back. Equation ( 3) combines the FP model (single-mode propagation) and the single-reflection approximation. It evidences the importance of the effective index 𝑛 eff (𝑘) of the fundamental mode. Insertion of (3) into (1), calculation of the quadratic expression and averaging over the time results in

𝐼(𝑧) ≈ 𝐼 0 + ∫ |𝑠(𝑘)| 2 exp[i𝜑(𝑘)] exp (-2i𝑘𝑧 s ) ⏟ =:𝑓(𝑘) 𝑎 t (𝑘) exp[2i𝑘(𝑧 -∆𝑧 t )] 𝑑𝑘 2𝜋 ⁄ + cc 𝑘 max 𝑘 min + ∫ |𝑠(𝑘)| 2 exp[i𝜑(𝑘)] exp (-2i𝑘𝑧 s ) ⏟ =:𝑓(𝑘) 𝑎 b (𝑘)exp[-2𝑘 Im(𝑛 eff )𝐻] 𝑘 max 𝑘 min × exp{2i𝑘[𝑧 -∆𝑧 t -Re(𝑛 eff )𝐻]} 𝑑𝑘 2𝜋 ⁄ + cc , (4)
where 𝐼 0 is an offset, which has been removed in the presented examples. The parameters 𝑎 t , 𝑎 b , and 𝑛 eff are calculated with the FP model. The expression 𝑓(𝑘) is obtained from a reference measurement at a planar substrate with the same machine settings (same light source power, same objectives, same focal 𝑧 -position)preferably taken shortly before or after the measurement at the TSV in the region around the TSV (see section 2). The intensity pattern of the reference interferogram is given by

𝐼 s (𝑧) ≈ 𝐼 s,0 + ∫ |𝑠(𝑘)| 2 exp[i𝜑(𝑘)] 𝑟(𝑘) exp[2i𝑘(𝑧 -𝑧 s )] 𝑑𝑘 2𝜋 ⁄ + cc 𝑘 max 𝑘 min = = 𝐼 s,0 + ∫ |𝑠(𝑘)| 2 exp[i𝜑(𝑘)] exp(-2i𝑘𝑧 s ) 𝑟(𝑘) ⏟ =𝑓(𝑘) 𝑟(𝑘) exp(2i𝑘𝑧) 𝑑𝑘 2𝜋 ⁄ + cc 𝑘 max 𝑘 min , (5)
where 𝑟(𝑘) is the Fresnel reflection coefficient at a plane air/silicon interface and the subscript "S" means "planar substrate". A Fourier transform of the reference interferogram then yields

𝑓(𝑘) 𝑟(𝑘) = 1 𝑝 ∫ 𝐼 s (𝑧) exp(-2i𝑘𝑧) 𝑑𝑧 𝑧 0 +𝑝 𝑧 0 ≈ 1 𝑝 ∑ 𝐼 s (𝑧 𝑗 ) • ∫ exp(-2i𝑘𝑧) 𝑧 𝑗 + 1 2 ∆𝑧 𝑗 𝑧 𝑗 - 1 2 ∆𝑧 𝑗 𝑑𝑧 𝑀 𝑗=1 . (6)
As discussed in section 6 later and in fig. 2 on the right and fig. 14, 𝑝 is the size of the periodized region, which contains the TSV length plus several coherence lengths 𝑙 c on top and bottom. 𝑀 is the number of sampling points. The intensity data 𝐼 s (𝑧 𝑗 ) are given on a non-equidistant mesh and the integral in ( 6) is approximated. Equation ( 6) finally yields 𝑓(𝑘) needed for simulations and the solution of inverse problems at TSV. As an example, we have simulated the intensity pattern 𝐼(𝑧) -𝐼 0 in dependence on the height 𝐻 for a circular TSV with a diameter 𝐷 = 5 µm and a beam waist 𝑤 = 2.7 µm (supplementary visualization 1 and 2 as well as fig. 8); each slide of the supplementary movies shows the intensity 𝐼(𝑧) -𝐼 0 for a given TSV height 𝐻. For large heights 𝐻, there are separated signals from top and bottom. For small heights, the signals (fringe patterns) from top and bottom overlap and interfere; constructive and destructive interference of the fringes take turns as 𝐻 is varied with a periodicity ∆𝐻 = 𝜆 c /{2Re[𝑛 eff (𝜆 c )]}.

Using the effective index for the interpretation of an OCT interferogram

For large TSV depths, the signals from top and bottom are well-separated in the interferogram, see fig. 9. The FP model provides a physically intuitive yet quantitative interpretation. According to [START_REF] Ku | Reflectometer-based metrology for high-aspect ratio via measurement[END_REF], the distance between the two peaks of the envelope depends on the TSV height 𝐻 and the group velocity 𝑣 g = 𝑐/𝑛 g with 𝑐 the speed of light and 𝑛 g the group index. The latter can be calculated from the effective index by [START_REF] Born | Principles of Optics[END_REF]

Moreover, the intensity of the second peak, which results from the reflection of the mode at the bottom of the via, is driven by (i) the imaginary part of 𝑛 eff (𝑘) and (ii) the coupling coefficients of the FP model, 𝑡 t and 𝑟 b in (2). If the mode attenuation is too strong (as is the case for small diameters and widths), the intensity of the second peak will be extremely weak. We touch here an intrinsic limit of the OCT for the metrology of etched TSVs: it is not possible to measure too small and deep structures with illumination from the top of the TSV. We quantify in section 5 the maximum measurable depth as a function of the diameter. Up to now, when interpreting an OCT interferogram from a TSV, the most commonly used technique consists in applying scalar ray tracing, which amounts to consider 𝑛 eff (𝑘) ≈ 𝑛 air (𝑘) = 1 and 𝑛 g = 1 [START_REF] Ku | Reflectometer-based metrology for high-aspect ratio via measurement[END_REF][START_REF] Ku | Characterization of high density through silicon vias with spectral reflectometry[END_REF][START_REF] Ku | In-line metrology of 3D interconnect processes[END_REF][START_REF] Ku | Spectral reflectometry for metrology of three-dimensional through-silicon vias[END_REF][START_REF] Vienne | A Comparison of TSV Etch Metrology Techniques[END_REF][START_REF] Alliata | Progress in Time-Domain Optical Coherence Tomography for TSV / 3Di stacking Metrology[END_REF][START_REF] Devanciard | Through Silicon Via Process Characterization by Integrated Inspection/Metrology Solutions in Visible and Infrared Domain[END_REF]. For vias with lateral dimensions smaller than 5 µm, we recommend the use of 𝑛 eff (𝑘) instead, which we have calculated and presented in section 4.1 and tabulated in the supplementary data files 1 and 2. The values of 𝑛 g (𝑘) are shown in fig. 10. By using the correct values for 𝑛 eff (𝑘, 𝐷) and integrating over 𝑘 as in (4), we take correctly into account the deviations from 𝑛 eff = 1 and 𝑛 g = 1. Moreover, (4) also accounts for possible distortions of the signal.

Validity of the Fabry-Perot model

In this section, we discuss the accuracy of the FP model, which means the neglect of higherorder modes. Large TSV heights are connected with large losses of higher-order modes, meaning their neglect by the FP model introduces only small errors; the accuracy of the FP model increases therefore with the TSV height (fig. 11). Similarly, small diameters lead to large losses of higher-order modes as well and the FP model benefits from this. Concluding, the FP model is the better the larger the depth and the smaller the diametersmeaning the larger the aspect ratio. For most of the TSV aspect ratios in practice, it works sufficiently well.

Figure 12 shows the error in the reflected amplitude when considering merely single-fold instead of multiple-fold reflection, meaning the approximation

1 [1 -𝑟 b • 𝑟 t • exp(2i𝑘𝑛 eff 𝐻)] ⁄ ≈ 1 in (2)
; the resulting error is quite small and decreases with increasing TSV height 𝐻. 2), modelling the multiple reflections between TSV top and bottom. The term is small and vanishes when increasing the height 𝐻. Thus, considering a single reflection at the TSV bottom allows reaching a good accuracy in most cases.

Maximum measurable depth

The FP model is an accurate electromagnetic model. Therefore, using the results and experience gathered in the previous sections, we can give an estimate on the maximum measurable depth when measuring from the top of the wafer. Our estimate is based on the attenuation experienced by the signal during its propagation to the bottom of the sample and back to top:

-At the top, merely a small fraction of the light is not coupled into the TSV but is reflected here in the best case (very small beam waist). We consider that 99% may be transmitted into the TSV. -We assume that 99% of the light coupled into the hole may be lost at the sidewalls during the propagation down-and upwards. -We neglect the multiple reflections between top and bottom since the attenuation due to a single pass is already large.

Under these assumptions, the signal from the bottom has only 0.99 × 0.01 = 0.0099 of the amplitude of the reflection at a planar substrate. Indeed, weaker signals might be detected, but this would require the use of advanced techniques such as a band pass filter for noise reduction [START_REF] Drexler | Optical Coherence Tomography[END_REF], the conduction of several measurements and application of averaging techniques as well as the numeric analysis of signals. Alternatively, different physical approaches such as measurements or additional measurements from belowanalogously to [START_REF] Marx | Wafer Thickness Sensor (WTS) for etch depth measurement of TSV[END_REF][START_REF] Teh | Backside Infrared Interferometric patterned Wafer Thickness Sensing for Through-Silicon-Via (TSV) etch Metrology[END_REF][START_REF] Gambino | Through-silicon-via process control in manufacturing for SiGe power amplifiers[END_REF] could be chosen.

Subsequently, we consider the value of 99% of the light coupled into the hole as a meaningful estimate for measurements from the top of the wafer. So, the maximum measurable depth 𝐻 max is given by the approach

exp [ i • 2𝜋 𝜆 c • i • Im(𝑛 eff ) 2𝐻 max ] = 10 -2 . ( 8
)
𝐻 max = ln ( 10 2 )
4𝜋 Im(𝑛 eff) 𝜆 c [START_REF] Fursenko | Through silicon via profile metrology of Bosch etching process based on spectroscopic reflectometry[END_REF] with the center wavelength 𝜆 c = 1.329 µm here. Figure 13 displays the value of 𝐻 max for 1D trenches (TE-and TM-polarization) as well as circular holes. The losses at the sidewalls are largest in the case of circular holes and they are larger for TM-than for TE-polarization at 1D trenchessee also section 4.1 and supplementary data files 1 and 2. A larger attenuation directly results in a smaller 𝐻 max .

Fig. 13. Maximum measurable depth 𝐻 max in dependence on the trench width ∆𝑥/µm or TSV diameter 𝐷 in the case of 1D trenches, using TE-polarization (purple) or TM-polarization (green) as well as circular TSV, using linear polarization (blue).

Electromagnetic TSV analysis procedure and experimental verification

Experimental setup

We use our electromagnetic modelling for retrieving the TSV depth from OCT measurements. A silicon wafer containing high aspect ratio 1D trenches as well as circular holes with an intended etching depth of 20 μm has been fabricated by deep reactive ion etching in an Inductively Coupled Plasma (ICP) reactor with a Bosch-like process with SF 6 and C 4 F 8 gases [START_REF] Wu | High aspect ratio silicon etch: A review[END_REF]. The intended TSV height 𝐻 is around 20 μm in the case of large trench widths Δ𝑥 and diameters 𝐷 but, due to the aspect ratio dependent etching (ARDE) effect [START_REF] Wu | High aspect ratio silicon etch: A review[END_REF], the etching depth 𝐻 decreases when Δ𝑥 and D are reduced. The measurements have been performed with a TMAP, a metrology tool commercialized by the company UnitySC [START_REF]UnitySC, 611 Rue Aristide Berges, 38330 Montbonnot Saint Martin[END_REF], equipped with a LISE (a TD-OCT sensor commercialized by the company FOGALE Nanotech [START_REF] Nanotech | Rue de l'Hostellerie[END_REF]). The spectrum of the light source after passing through all optical components has a complex shape and ranges from 𝜆 = 1.20 µm to 𝜆 = 1.43 µm, see section 2. The center wavelength is 𝜆 c = 1.31 µm here.

We measure each structure with the focus of the collection objective at the center of the TSV top. There are 4 selectable objectives with respective magnifications and NA (5x; 0.14), (10x; 0.26), (20x; 0.40), (50x; 0.42). If possible, we always use the objective with which signals from top and bottom of the TSV are obtained with a similar amplitude. Unfortunately, such a procedure is not possible in the case of small ∆𝑥 or 𝐷, where the fraction of light coupled into the TSV is limited by the NA of the 50x objective or the beam diameter at the entrance pupil of the objective; the beam waist at the focus is then 𝑤 = 2.7 µm. Here, the signal from the bottom is much smaller than the one from the top and close to the noise floor; moreover, it overlaps with the signal from the top due to small TSV heights compared to the source coherence length.

After a TSV measurement, we shift the objective some µm in 𝑥or 𝑦-direction without changing its z-position or other machine settings in order to record the interferogram at the planar substrate yielding 𝑓(𝑘) needed in (4); 𝑓(𝑘) is then obtained from [START_REF] Ku | In-line metrology of 3D interconnect processes[END_REF]. Since the interferograms at the planar substrate depend on the objective due to dispersion, they have to be recorded for each objective separately (see figs. 16 -17). Also due to the impact of the dispersion at the objective, its NA and the focal z-position on the fringe pattern [START_REF] Bosseboeuf | Assessment of technology and (thermo)mechanical behavior of MEMS devices by interference microscopy[END_REF][START_REF] De Groot | Principles of interference microscopy for the measurement of surface topography[END_REF] as well as possible fluctuations of the light source, we record the reference spectrum at a planar substrate in the vicinity of the structure after each TSV measurement without changing the objective nor its z-position.

Processing of the OCT data

The main goal is to obtain a reliable, accurate and fast retrieval of the TSV height 𝐻 from the OCT data. As well, we would like to construct a procedure, which could be easily generalized to the retrieval of additional TSV parameters such as the diameter 𝐷 . That is an inverse problem; for a measured known spectral amplitude or rather 𝑓(𝑘) in ( 4), which we obtain from (6) from an interferogram 𝐼 s (𝑧) measured at points 𝑧 𝑗 at a planar substrate, we have to find 𝐻 so that the signal 𝐼 sim,TSV (∆𝑧 𝑡 , 𝐻, 𝑧) simulated with (4) best matches the measured signal 𝐼 TSV (𝑧). The parameters 𝑎 t , 𝑎 b , 𝐻 and ∆𝑧 t in (4) are now fit parameters used to match the measured signals.

A problem class related to the solution of the inverse problem at TSVs and trenches is the solution of the respective problem at thin films. Popular data evaluation techniques here operate in the frequency spacealso for efficiency reasons (limited bandwidth of the utilized light sources) -and include the Fourier amplitude method [START_REF] De Groot | Principles of interference microscopy for the measurement of surface topography[END_REF][START_REF] De Groot | Angle-resolved three-dimensional analysis of surface films by coherence scanning interferometry[END_REF][START_REF] De Groot | Scanning interferometry for thin film thickness and surface measurements[END_REF][START_REF] Wan | Measurement of thin films using Fourier amplitude[END_REF] and the spectral nonlinear phase method [START_REF] De Groot | Principles of interference microscopy for the measurement of surface topography[END_REF][START_REF] De Groot | Angle-resolved three-dimensional analysis of surface films by coherence scanning interferometry[END_REF][START_REF] De Groot | Scanning interferometry for thin film thickness and surface measurements[END_REF][START_REF] Kim | Thickness-profile measurement of transparent thin-film layers by white-light scanning interferometry[END_REF][START_REF] Ghim | Thin-film thickness profile measurement using a Mirau-type low-coherence interferometer[END_REF]. These techniques are based on least square approaches, which match the measured amplitude or phase information respectively with computed results obtained on the basis of a separately recorded reference spectrum. A combination of both methods in order to exploit both amplitude and phase information seems reasonable [START_REF] Dong | Sensitivity analysis of thin-film thickness measurement by vertical scanning whitelight interferometry[END_REF][START_REF] De Groot | Transparent film profiling and analysis by interference microscopy[END_REF]. This is a reason why we match complex Fourier coefficients here, incorporating all the available information by least squares as a recognized, quantitative measurement data evaluation procedure; the resulting sum of squares (merit function) here is then minimized by the Marquardt-Levenberg algorithm (MLA) [START_REF] Dahmen | Numerik für Ingenieure und Naturwissenschaftler[END_REF][START_REF] Riley | Mathematical methods for physics and engineering[END_REF] as one possible choice and a standard.

Before solving the inverse problem, we apply filters to 𝐼 s (𝑧) and 𝐼 TSV (𝑧) in order to reduce the noiseone in direct space and one bandpass filter in Fourier space.

Our procedure requires a Fourier transform (6) of the measured interferograms 𝐼 s (𝑧) as well as 𝐼 TSV (𝑧). Our Fourier transform here is discrete in view of a numeric treatment, yielding few hundred coefficients typically. For this, we define a periodized computational domain (fig. 14 In case of 𝐼 s (𝑧), we apply 𝑔 s (𝑧) with a much smaller region of transmission (fig. 14,bottom). This choice suppresses noise of 𝐼 s (𝑧) in regions far away from 𝑧 s , where 𝐼 s (𝑧) is not meaningful. So, the reference signal becomes noise-free in large regions of the period so that 𝐼 sim,TSV (∆𝑧 t , 𝐻, 𝑧) obtained later from (4) with 𝑓(𝑘) computed from 𝐼 s (𝑧) has less noise. , with 𝑐 sim,TSV,𝑗 (∆𝑧 t , 𝐻): = 𝑓(𝑘 𝑗 ) 𝑎 t (𝑘 𝑗 ) exp(-

2i𝑘 𝑗 ∆𝑧 t ) + +𝑓(𝑘 𝑗 ) 𝑎 b (𝑘 𝑗 ) exp(-2i𝑘 𝑗 ∆𝑧 t ) • exp{-2i𝑘 𝑗 Re[𝑛 eff (𝑘 𝑗 )]𝐻} • exp{-2𝑘 𝑗 Im[𝑛 eff (𝑘 𝑗 )]𝐻}, (12)
"cc" meaning "complex conjugated" so that the merit function is real.

With the spectrum being quite narrow, we may approximate 𝑎 t (𝑘 𝑗 ) by a t (𝑘 c ) and a b (𝑘 𝑗 ) by a b (𝑘 c ) in the given case with 𝑘 c the center frequency (the subscript "c" means "center"). An example of 𝑚(∆𝑧 t , 𝐻) in dependence on 𝐻 is shown in fig. 18 (upper part); the shape of 𝑚(∆𝑧 t , 𝐻) around its global minimum is given by an envelope and an oscillation with the frequency 2𝑘 c Re[𝑛 eff (𝑘 c )]. When the fringe patterns of the simulated and measured signal from the bottom coincide, 𝑚(∆𝑧 t , 𝐻) has a local minimum regarding H; when there is a phase shift of π between them, it has a local maximum regarding H; the same holds for the fringe patterns associated with the top with respect to ∆𝑧 t . With the MLA being merely locally convergent, it would be stuck in a local minimum with respect to ∆𝑧 t as well as 𝐻. A rather similar problem has been reported in the case of white light interferometry [START_REF] Kim | Thickness-profile measurement of transparent thin-film layers by white-light scanning interferometry[END_REF] and spectral reflectometry [START_REF] De Groot | Scanning interferometry for thin film thickness and surface measurements[END_REF] at thin films. In order to make the MLA globally convergent, we introduce the complex coefficients 𝑎 sim,t for 𝑎 t (𝑘 c ) and 𝑎 sim,b for 𝑎 b (𝑘 c ) as optimization parameters enabling phase optimization (fringe pattern optimization) and additionally, we aim at a mere propagation of the envelope of the signals from top and bottom -not of the associated optimized fringe positions having a frequency of 2𝑘 c Re[𝑛 eff (𝑘

) exp[-2i(𝑘 𝑗 -𝑘 c )∆𝑧 t ] + + 𝑎 sim,b 𝑓(𝑘 𝑗 ) exp[-2i(𝑘 𝑗 -𝑘 c )∆𝑧 t ] • exp{-2i𝑘 𝑗 Re[𝑛 eff (𝑘 𝑗 )]𝐻 + 2i𝑘 c Re[𝑛 eff (𝑘 c )]𝐻} • • exp{-2𝑘 𝑗 Im[𝑛 eff (𝑘 𝑗 )]𝐻} ; (13)
a wavelength dependence of reflection and transmission at TSV top and bottom in the case of a broader spectrum could be incorporated into (13) by wavelength-dependent coupling coefficients as additional factors, which could be read and interpolated from precalculated tables.

Concluding, the simulated fringes are hardly moved when changing ∆𝑧 t or 𝐻 from iteration to iteration and the phase of the simulated fringe patterns belonging to top and bottom is optimized in each iteration by the optimization of 𝑎 sim,t and 𝑎 sim,b . So, the oscillations are removed here in case of 𝑚 ̃(𝑎 sim,t , ∆𝑧 t , 𝑎 sim,b , 𝐻), which is slowly varying and connects the local minima of the oscillating merit function 𝑚(∆𝑧 t , 𝐻). Fig. 18 (lower part) shows an example of the resulting merit function 𝑚 ̃(𝑎 sim,t , ∆𝑧 t , 𝑎 sim,b , 𝐻) in dependence on 𝐻 around its global minimum. The merit function is smooth so that even numeric differentiation of it with respect to ∆𝑧 t and 𝐻 is possible.

Since 𝑚 ̃(𝑎 sim,t , ∆𝑧 t , 𝑎 sim,b , 𝐻) depends linearly on 𝑎 sim,t and 𝑎 sim,b , the convergence of the MLA is rapid with respect to them (about 2 -3 iterations). The convergence of the MLA with respect to ∆𝑧 t and 𝐻 depends on the estimate of 𝑘 c : an error of several per cent here deteriorates the convergence notably, but an error in the order of ±1% doesn´t. Knowing the spectral amplitude at the sample, our errors in 𝑘 c are in the order of ±1%. The convergence also depends on the initial values. According to our experience, the convergence is best when starting the MLA with 𝑎 sim,t = 𝑎 sim,b = 0. In this way, a mismatch of the phases of the simulated and the measured fringes cannot occur in the first iteration (and in later steps, it does not occur anyway). The most reasonable choice for the z-position of the TSV top is ∆𝑧 t = 0; for 𝐻 , the expected value or several values within the expected range can be taken (e.g. equidistant sampling of an interval [𝐻 min ; 𝐻 max ]), usually leading to the same result. In the case of various results, the one minimizing the merit function [START_REF] Teh | Backside Infrared Interferometric patterned Wafer Thickness Sensing for Through-Silicon-Via (TSV) etch Metrology[END_REF] has to be selected.

At this point, we do not yet exploit the information on the position of the fringes of the signal. Such information is still contained in the phase of a t (𝑘 c ) and a b (𝑘 c ) and the value of 𝑘 c . However, with the fringes having a distance in the order of λ c {2Re[𝑛 eff (𝑘 c )]} ⁄ ≈ 0.66 nm, and our noisy OCT signals leading rather to a standard deviation of σ(H) in the order of several hundred nanometers in the presented examples (shallow, mostly small TSV, yielding weak, noisy signals from the TSV bottoms and overlapping signals from TSV top and bottom as in figs. 19 -20), the decision whether 𝐻 is a fraction of λ c {2Re[𝑛 eff (𝑘 c )]} ⁄ larger or smaller so that it matches the measured fringe positions of 𝐼 TSV (𝑧 j ) is subject to future worktogether with further noise reductions.

Our progress here is the electromagnetic modelling: We apply the correct values for 𝑛 eff (𝑘, 𝐷) , model the spectral amplitude and dispersion accurately, and in particular, we consider the interference of the fringe patterns from TSV top (1 st term in [START_REF] Teh | Backside Infrared Interferometric patterned Wafer Thickness Sensing for Through-Silicon-Via (TSV) etch Metrology[END_REF]) and bottom (2 nd term in ( 13)) contrary to [START_REF] Vienne | A Comparison of TSV Etch Metrology Techniques[END_REF][START_REF] Alliata | Progress in Time-Domain Optical Coherence Tomography for TSV / 3Di stacking Metrology[END_REF][START_REF] Devanciard | Through Silicon Via Process Characterization by Integrated Inspection/Metrology Solutions in Visible and Infrared Domain[END_REF][START_REF] Courteville | Contact free on-axis metrology for the fabrication and testing of complex optical systems[END_REF][START_REF] Wilhelm | A novel, low coherence fibre optic interferometer for position and thickness measurements with unattained accuracy[END_REF]. These two intensity patterns can e.g. interfere constructively (coinciding fringes, leading to a large contrast 𝐶(𝑧) and envelope in the region of overlap as in fig. 19) or destructively (phase shift of 𝜋 between the fringes of both signals, leading to a low contrast and small envelope in the region of overlap as in fig. 20). In extreme cases, the contrast of the fringes of 𝐼 TSV (𝑧 𝑗 ) has more than 2 local maxima, so that it is at the first glance not clear which of these maxima to associate with the TSV top or bottom and which to associate with a region of constructive interference or a sidelobe (fig. 19). Our novel procedure gives a clear answer here at a low computational cost without the need of a very good initial guess.

Experimental verification

We study the impact of our electromagnetic modelling on the basis of 1D trenches and circular holes of intended etching depth 20 µm. We consider TSVs of various diameters -smaller diameters leading to smaller depths due to the etching physics.

We apply the MLA to [START_REF] Teh | Backside Infrared Interferometric patterned Wafer Thickness Sensing for Through-Silicon-Via (TSV) etch Metrology[END_REF]. From a short look at few interferograms of 1D trenches of large width ∆𝑥, it is clear that there cannot be any etching depth larger than 24 µm. Thus, we restrict the modelled TSV height 𝐻 to the interval [4.01 µm; 24.00 µm] and start the MLA with the initial values 𝐻 = 5 µm, 𝐻 = 7 µm, ..., 𝐻 = 23 µm. In addition, we restrict ∆𝑧 t to [-4.0 µm; +4.0 µm], which is clearly sufficient (∆𝑧 t is in the order of few hundred nm). With a convergence radius for 𝐻 in the order of 10 µm, this sampling is more than sufficient. In most cases, the values found for ∆𝑧 t and 𝐻 agree with an accuracy in the order of 1 nm; more than one minimum of 𝑚 ̃(𝑎 sim,t , ∆𝑧 t , 𝑎 sim,b , 𝐻) is found in few cases and the smallest one is taken. The iteration number in the considered examples is usually in the order of 5 -25; on average, it is in the order of 10. The number of considered complex Fourier coefficients is 175 here. Having implemented our procedure in C++ and running it on a laptop with a 2.6 GHz processor, one MLA iteration takes 0.9 ms. So, our iterative procedure is rapid.

We compare the presented electromagnetic analysis of 𝐼 TSV (𝑧) with the widely applied envelope-based one: here, equation ( 20) from [START_REF] Larkin | Efficient nonlinear algorithm for envelope detection in white light interferometry[END_REF],

𝐶 3 = √ (𝐼 2 -𝐼 4 ) 2 -(𝐼 1 -𝐼 3 )(𝐼 3 -𝐼 5 ) 4 sin 4 (𝜓)
, [START_REF] Gambino | Through-silicon-via process control in manufacturing for SiGe power amplifiers[END_REF] is applied in order to remove the fringes from 𝐼 TSV and to extract the envelopemeaning the local fringe contrast 𝐶; 𝜓 is the phase step between adjacent sampling points, and 𝐼 1 , … , 𝐼 5 are the adjacent intensity samples; all quantities here are z-dependent and discretized. The distance of the fringes is assumed to be 𝜆 c 2 ⁄ . A reference spectrum or an interference between TSV top and bottom is not incorporated here.

In the case of large 𝐻 or destructive interference, the signals from top and bottom are separated clearly so that 𝐶(𝑧) has 2 distinct maxima (fig. 20, top right)one associated with the signal from the top and one with the signal from the bottom. We fit 𝐶(𝑧) with a sum of two Gaussian functions, optimizing their positions 𝑧 t , 𝑧 b , amplitudes 𝑎 t , 𝑎 b and widths FWHM t , FWHM b and conclude 𝐻 = 𝑧 b -𝑧 t (figs. 19-20, top right). No information on the spectrum, dispersion or interferogram at a planar substrate is exploited here apart from the knowledge of 𝜆 c . In the case of shallow TSV, there may be only one maximum of 𝐶(𝑧) or more than two due to the coherent addition of the interference pattern from the top and the interference pattern from the bottom. Here, pre-knowledge on the intended etching depth and physics is required in order to start the Gaussian fits with initial values that are already very good. In contrast to this, our novel technique does not need such help. Table 1 lists the results obtained by the common as well as the novel technique for TEpolarization at 1D trenches, table 2 for TM-polarization at 1D trenches and table 3 for linear polarization at circular holes. For each OCT measurement, there are two interferograms 𝐼 TSV (𝑧), one from the motion of the reference mirror in +𝑧direction and one from -𝑧direction. We treat and list them as separate measurements here and calculate weighted averages

𝐻 ̅ = (∑ 1 𝑒𝑟𝑟 std,𝑛 2 𝜁 𝑛=1 ) -1 • ∑ 1 𝑒𝑟𝑟 std,𝑛 2 𝐻 𝑛 𝜁 𝑛=1 (15)
from them with 𝜁 = 2 the number of measurements in this paper. The "±" signs in the tables indicate asymptotic standard errors " 𝑒𝑟𝑟 std " [START_REF] Riley | Mathematical methods for physics and engineering[END_REF][START_REF] Williams | GNUPLOT[END_REF] calculated from the MLA; they are computed from the MLA residual and derivatives by

𝑒𝑟𝑟 std = sqrt ( 𝑚 ̃(𝑎 sim,t , ∆𝑧 𝑡 , 𝑎 sim,b , 𝐻) 𝑁-𝑄 {diag[( 𝐽 + • 𝐽 ) -1 ]} 4,4 ) (16)
with 𝑚(𝑎 sim,t , ∆𝑧 t , 𝑎 sim,b , 𝐻) the residual for the final set of parameters, 𝑁 the no. of coefficients (175 in this paper), 𝑄 = 4 the no. of degrees of freedom and 𝐽 the Jacobian matrix, 𝐽 𝑗,𝑞 = 𝜕𝑐 sim,TSV,𝑗 𝜕parameter 𝑞 ⁄ , where 𝑞 denotes the parameter number and parameter 1 = 𝑎 sim,t , parameter 2 = ∆𝑧 t , parameter 3 = 𝑎 sim,b , parameter 4 = 𝐻 . The diagonal matrix element (4, 4) is taken since the desired quantity 𝐻 is the 4 th optimization parameter here. Equation ( 16) is only an estimate for the standard deviation and moreover, it does not include all sources of errors. For example, it neglects the repeatability. The Gaussian fits have been done by "Gnuplot" [START_REF] Williams | GNUPLOT[END_REF]; the asymptotic standard errors here neglect additional errors introduced by the underlying physical models -e.g. the fit of a non-Gaussian peak of 𝐶(𝑧 𝑗 ) by a Gaussian function or missing interference effects between signals from top and bottom for the common envelope-based procedure. According to our experience, 𝑒𝑟𝑟 std is too optimistic especially for the common envelope-based procedure and a small value of 𝑒𝑟𝑟 std here does not at all mean that it delivers more accurate final results than our electromagnetic analysis. Instead, 𝑒𝑟𝑟 std is merely suited for comparisons within a certain analysis proceduremeaning within a single column in table 1 -3, allowing the calculation of a weighted average [START_REF] Diebold | Perspective: Optical measurement of feature dimensions and shapes by scatterometry[END_REF]. Moreover, 𝑒𝑟𝑟 std,n allows a classification of the quality of the nth measurement. For example, measurements with 𝑒𝑟𝑟 std,𝑛 ≥ 0.4 in the presented test cases, indicating a strong influence of noise or weak signals here, could be discarded automatically by a computer. The comparisons of the columns 3 and 5 of table 1 -3 indicate that our novel electromagnetic analysis yields more reliable and accurate results for 𝐻 ̅ : our novel analysis agrees better with the etching physics (smaller depths in case of smaller diameters or trench widths). In addition, the reproducibility is clearly better in case of our novel technique. In contrast to this, the common technique yields strongly varying heights 𝐻 ̅ particularly when measuring small circular TSV of the same size (𝐷 = 2.5 µm) or similar size (𝐷 = 2.5 µm and 𝐷 = 3.0 µm). The tables for the 1D trenches contain an additional column on the right, dividing the amplitude of the signal from the bottom by the one from the top. This confirms our findings that the signal from the bottom is weaker (and thus more difficult to detect) for TM-polarization. Concluding, TE-polarization is preferable for the analysis of narrow 1D trenches.

In the case of the circular hole of diameter 𝐷 = 2.5 µm, there are 𝜁 = 6 measurements in total, so that the calculation of a standard deviation 𝜎(𝐻 ̅ ) is reasonable: (

𝜎(𝐻 ̅ ) = √
) 17
Evaluating these 6 measurements on the basis of ( 15) and ( 17), we get 𝐻 ̅ = 16.2 µm, 𝜎(𝐻 ̅ ) = 0.5 µm for the common technique and 𝐻 ̅ = 13.1 µm, 𝜎(𝐻 ̅ ) = 0.4 µm for the novel technique. The result obtained by the novel technique is in better agreement with the etching physics and has a smaller standard deviation than the one obtained by the common technique. The minimum width ∆𝑥 or diameter 𝐷 for which we can measure 𝐻 reliably is 2.0 µm for the 1D trenches and 2.5 µm for the circular holes, which coincides rather well with our findings for the maximum measurable depth (section 5). However, further progress is possible here in our opinion by reducing the influence of noise and sampling 𝐼(𝑧) more densely. In the case of a small ∆𝑥 or 𝐷, it would be beneficial to use an objective of larger NA and magnification (e.g. NA = 0.7 and 100x magnification, presuming a complete filling of the entrance pupil of the objective).

Summarizing, the precision of H and the axial resolution depend on the physical modelling, the data evaluation, but also on the dispersion and the coherence length 𝑙 c of the light source (we have chosen the settings of the OCT device available to us so that 𝑙 c has been as small as possible). The present paper has been focussed on the modelling and data evaluation; a section of a future publication will be devoted to the utilized light source and its spectrum.

Conclusion

In this paper, we have rigorously simulated a Time-Domain OCT device for improving the accuracy of Through Silicon Via (TSV) height measurementsalso in the case of overlapping fringe patterns, and compared the results to those of previous publications on the subject. We have considered circular holes and 1D trenches as the most common shapes of TSVs. The simulations have been conducted using realistic, recorded spectra of the light source, and considering interferometer dispersion.

The most accurate but also computationally heaviest technique applied here is the rigorous modal calculation based on a Fourier-Bessel basis for circular TSVs and on the aperiodic Fourier modal method (a-FMM) for 1D trenches. This technique allows us to calculate the modal propagation constants and effective indices for the TSVs, which can be used for an accuracy enhancement of usual methods that are based on ray tracing. These effective indices have a real part impacting the optical path length as well as an imaginary part characterizing the signal attenuation due to radiative leakage at the sidewalls. For TSV diameters below 5 µm, the effective indices differ considerably from the ones in air so that we recommend to take them into account. Aiming at a combination of an accuracy similar to a modal calculation with the physical intuition of ray-tracing techniques, we have introduced the Fabry-Perot model. This model is a very good approximation in the case of TSVs with a sufficient aspect ratio (in the order of 10) when only the fundamental mode determines the interaction between top and bottom significantly, which is mostly the case in practice. The reflection and transmission coefficients can be taken from the rigorously calculated S-matrix entries (coupling coefficients) and the propagation between top and bottom takes place by application of effective indices. On the basis of this model, we have estimated the maximum measurable depth of TSVs for illumination and detection from the top, which depends on TSV width or diameter, TSV shape and the wavelength range. In the case of small diameters or very deep TSV, the use of smaller wavelengths or the inclusion of additional measurements from below is advisable.

Based on our electromagnetic analysis, we have developed our insight at a novel analysis procedure for time-domain OCT at TSV. This analysis comprises data preprocessing by smoothening filters, the accurate knowledge of the spectrum, dispersion, effective indices, amplitude as well as phase information and the interference of the fringe patterns associated with TSV top and bottom. The novel procedure is based on an iterative least squares approach, but globally convergent in a given parameter interval due to the presented enhancements and can be readily automated. It is more reliable and accurate than the popular envelope-based technique, but nevertheless rapid. We think that it is a good basis for future progress, such as further accuracy enhancements, analysis of smaller TSVs and the conclusion of additional TSV parameters such as diameter profiles.

Fig. 3 .

 3 Fig.3. Experimental interferogram 𝐼(𝑧) recorded at a planar substrate (black), intensity reconstructed taking into account non-matched dispersion in the interferometer arms (red) and intensity reconstructed ignoring such dispersion (blue). The red curve agrees well with the measurement (black) in contrast to the blue curve[START_REF] Iff | Simulation of a Time-Domain OCT device for Through Silicon Via (TSV) fabrication[END_REF].

Fig. 4 .

 4 Fig. 4. Tangential electric field (left, TE-polarization) and tangential magnetic field (right, TM-polarization) for a 1D trench of air in silicon at λ c = 1.329 µm for n Si = 3.5 (trench width: 4 µm, trench height: 35 µm). TE-polarization is well-guided inside the TSV whereas TM-polarization leaves the TSV by larger part, meaning higher losses. There are standing wave patterns in z-direction with the periodicity 𝜆 c {2Re[𝑛 eff,TE (𝑘 c )]} ⁄ and 𝜆 c {2Re[𝑛 eff,TM (𝑘 c )]} ⁄ , respectively.

Fig. 5 .

 5 Fig. 5. Real and imaginary parts of the effective index for TE-(purple) and TM-(green) polarization as a function of the trench width (see also supplementary data file 1). The wavelength is the center wavelength 𝜆 c =1,329 µm of the OCT device and the refractive index is 𝑛 Si = 3.5. Re(𝑛 eff ) refers to the phase and Im(𝑛 eff ) to the attenuation due to radiative leakage. The values at slightly different wavelengths can be deduced from the fact that 𝑛 eff depends approximately merely on ∆𝑥/𝜆.

Figure 5

 5 Figure5shows the real (Re) and the imaginary (Im) part of the effective index in a 1D trench for TE-and TM-polarization as a function of the trench width ∆𝑥 at the center wavelength 𝜆 c = 1.329 µm. There is a strong dependence of Re(𝑛 eff,TE ) and Im(𝑛 eff,TE ) on ∆𝑥 for small ∆𝑥, where the deviation of 𝑛 eff,TE from 1 is the largest. For ∆𝑥 → ∞, this deviation tends asymptotically to zero. The spectral dependence of silicon is weak so that 𝑛 eff,TE depends only on ∆𝑥/𝜆 with ∆𝑥 the trench width. Therefore, one can deduce 𝑛 eff,TE (∆𝑥, 𝜆) at a wavelength different from the OCT center wavelength 𝜆 c = 1.329 µm from 𝑛 eff,TE (∆𝑥 • 𝜆 c 𝜆 ⁄ , 𝜆 c ), as long as |𝜆 -𝜆 c |/𝜆 c is sufficiently small and under the assumption that the refractive index of silicon is rather independent on 𝜆 in this narrow spectral range around 𝜆 c = 1.329 µm. Qualitatively, 𝑛 eff,TM (∆𝑥, 𝜆) behaves similarly to 𝑛 eff,TE (∆𝑥, 𝜆) . However, Im(𝑛 eff,TM ) is larger, meaning larger leakage through the sidewalls. Illustratively, there is a smaller reflection for TM-polarization at planar surfaces and a larger transmission into the substratelike in case of Fresnel reflection and transmission.

Fig. 6 .

 6 Fig. 6. Intensity distribution |𝐸 𝑥 (𝑟, 𝜃)| 2 for a circular TSV at 𝜆 = 1.329 µm. The refractive index of silicon is 𝑛 Si = 3.5.For 𝐷 = 2 µm, the electric field significantly touches the sidewallrequiring electromagnetic modelling; for 𝐷 = 10 µm, the field is well-confined in the holeindicating that a scalar method is accurate enough. In 𝑧-direction (along the optical axis), the field is damped exponentially; the damping increases with decreasing diameter D.

Fig. 7 .

 7 Fig. 7. Real part (left) and imaginary part (right) of 𝑛 eff (𝐷) for a circular air hole in silicon at 𝜆 c = 1.329 µm for 𝑛 Si = 3.5 (see also supplementary data file 2) [32]. The values at slightly different wavelengths can be deduced; approximately, 𝑛 eff merely depends on 𝐷/𝜆.

Fig. 8 .

 8 Fig.8. On the left: Illustration of the propagation scheme within the Fabry-Perot model; black arrows: light path; horizontal dashed black lines: reflecting facets of the FP cavity. Only the fundamental (leaky) mode is considered inside the TSV. Light is incident from the top; it is partially reflected in air as well as partially transmitted into the mode. The latter propagates downwards, is partially reflected at the bottom with the reflection coefficient 𝑟 b , propagates back upwards, is partially reflected and partially transmitted at the top with the reflection and transmission coefficients 𝑟 t , 𝑡 t , etc. The calculation results at the top and bottom interfaces, obtained separately, allow together with the effective index of the mode the calculation of the reflected amplitude for any height 𝐻, see[START_REF] Wu | High aspect ratio silicon etch: A review[END_REF]. We can then compute the interferogram by (1) and (2) efficiently for any TSV height 𝐻. The colored picture shows a cross section of the field distribution in and around the TSV on the basis of 𝐸 𝑥 (𝑥, 0, 𝑧). Upper picture on the right: example of an interferogram 𝐼(𝑧) -𝐼 0 in case of dispersion for a height 𝐻 at which the signals from TSV top and bottom interfere destructively; lower picture on the right: same for a height at which they interfere constructively. The interferograms in dependence on 𝐻 with as well as without interferometer dispersion are provided in terms of supplementary movies (visualization 1 and 2).

  𝑛 g (𝑘, 𝐷) = Re[𝑛 eff (𝑘, 𝐷)] + 𝑘 d d𝑘 Re[𝑛 eff (𝑘, 𝐷)].

Fig. 9 .

 9 Fig. 9. Example of a simulated interferogram for a circular TSV of large depth (TSV height H = 70.8 µm, TSV diameter D = 5 µm, beam waist w = 2.7 µm) [32].

Fig. 10 .

 10 Fig. 10. Group refractive indices n g,1D,TE (D, λ c ) (1D trench, TE-polarization, black curve), n g,1D,TM (D, λ c ) (1D trench, TM-polarization, red) and n g,circ (D, λ c ) (circular holes, green) as a function of the hole diameter at λ c = 1.329 µm [32].

Fig. 11 .Fig. 12 .

 1112 Fig. 11. Example of the relative error √ |𝑎𝑚𝑝𝑙 TSV,FP (𝜆) -𝑎𝑚𝑝𝑙 TSV (𝜆)| 2 |𝑎𝑚𝑝𝑙 TSV (𝜆)| 2 ⁄ in the reflected amplitude when applying the Fabry-Perot model to a circular TSV of 3.0 µm diameter and integrating over a half aperture angle of 15° (purple) and a half aperture angle of 60° (green); the beam is centered at the top of the TSV and has a waist of 𝑤 = 2.7 µm. The relative error decreases quickly with the TSV height 𝐻.

Fig. 14 .

 14 Fig. 14. Top: Computational domain of size 𝑝 = 𝑧 p -𝑧 o and weight function 𝑔 TSV (𝑧) applied in direct space to 𝐼 TSV (𝑧). Bottom: Computational domain of size 𝑝 = 𝑧 p -𝑧 o and weight function 𝑔 S (𝑧) applied in direct space to 𝐼 s (𝑧). The choice of the parameters here is applied uniformly to all subsequent examples. Having applied 𝑔 TSV (𝑧) and 𝑔 s (𝑧) respectively, both signals 𝑔 TSV (𝑧) • 𝐼 TSV (z) and 𝑔 s (z) • 𝐼 s (𝑧) are then Fourier transformed, since the algorithm operates in Fourier space. The computational domain of size 𝑝 is periodized along the 𝑧-direction, allowing discrete Fourier transforms. The results of the Fourier transforms are then weighted with the function 𝑔 ̃(𝑘), acting as a bandpass filter (fig. 15).

11 )

 11 Subsequently, 𝑔 s (𝑧) • 𝐼 s (𝑧) and 𝑔 TSV (𝑧) • 𝐼 TSV (𝑧) are Fourier transformed by[START_REF] Ku | In-line metrology of 3D interconnect processes[END_REF] and a bandpass filter ℎ ̃(𝑘) is applied to the Fourier coefficients of both signals (fig.15): 𝑘 0 )] for 𝑘 0 ≤ 𝑘 < 𝑘 min , 1 for 𝑘 min ≤ 𝑘 ≤ 𝑘 max , 𝑘 stop )] for 𝑘 max < 𝑘 ≤ 𝑘 stop , 0 for 𝑘 > 𝑘 stop .(The effect of the bandpass filter is illustrated in figs.16 -17.

Fig. 15 .

 15 Fig. 15. Weight function ℎ ̃(𝑘) applied to the Fourier transform FT[𝑔 s (𝑧) • 𝐼 s (𝑧) ] and FT[𝑔 TSV (𝑧) • 𝐼 TSV (𝑧)]. The choice of the parameters here is applied uniformly to all subsequent examples. FT[𝑔 TSV (𝑧) • 𝐼 TSV (𝑧)] • 𝑔 ̃(𝑘) yields the Fourier coefficients 𝑐 TSV,𝑗 and FT[𝑔 s (𝑧) • 𝐼 s (𝑧)] • 𝑔 ̃(𝑘) provides 𝑓(𝑘 𝑗 ) by means of (6) -needed subsequently for the simulation of the TSV, which yields 𝑐 sim,TSV,𝑗 in (12).

Fig. 16 .

 16 Fig.[START_REF] Vienne | A Comparison of TSV Etch Metrology Techniques[END_REF]. Interferogram at the planar substrate recorded with the 5x objective (left) and smoothened interferogram after application of the filter in direct space as well as the filter in Fourier space (bandpass filter). The signal is non-axissymmetric due to dispersion in the interferometer arms.

Fig. 17 .

 17 Fig. 17. Same as fig. 16 but for the 50x objective. The interferogram is slightly different and thus objective-dependent. Having filtered the signals at TSV and substrate, we solve the inverse problem by the MLA in Fourier space -where we already have the Fourier coefficients 𝑐 TSV,𝑗 as well as 𝑓(k 𝑗 ) on the basis of measurements and can propagate the harmonics by exponential factors in each iteration when computing the Fourier coefficients 𝑐 sim,TSV,𝑗 (∆𝑧 t , 𝐻) of the simulated TSV as modelbased counterpart. Applying the precomputed effective indices here, which are stored in a table and interpolated, the computational cost is very low. The solution of the inverse problem requires the minimization of the 2𝑁 differences in the Fourier coefficients 𝑐 sim,TSV,𝑗 (∆𝑧 t , 𝐻) and 𝑐 TSV,𝑗 . This leads to the merit function 𝑚(∆𝑧 t , 𝐻): = ∑ [𝑐 sim,TSV,𝑗 (∆𝑧 t , 𝐻) -𝑐 TSV,𝑗 ] 2 + cc 𝑁 𝑗=1

Fig. 18 .

 18 Fig. 18. Example of merit functions (sum of squares) in dependence on the TSV height 𝐻. Top: Merit function (12) oscillating in dependence on 𝐻; at local minima, the fringe patterns of modelled and measured signal from the TSV bottom coincide; at local maxima, they are phase-shifted by 𝜋 ; distance between two local minima: λ c {2Re[𝑛 eff (𝑘 c )]} ⁄ ≈ 0.66 µm. Bottom: Merit function (13) in dependence on 𝐻 around its global minimum when working with optimized fringe positions and moving only the envelope of the signal from the TSV bottom. The computed values (purple crosses) coincide with the parabolic fit (green), meaning that the merit function is noise-free. With respect to ∆𝑧 t , the merit functions look similar as with respect to 𝐻 . The optimized parameters are here (𝑎 sim,t ; ∆𝑧 t ; 𝑎 sim,b ; 𝐻) = (-0.061 -0.29i; -0.53 µm; -0.036 -0.17i; 19.70 µm).

Fig. 19 .

 19 Fig. 19. Top left: measured interferogram before bandpass filtering. Bottom left, green: interferogram simulated with best matching parameter set for TE-polarization at a 1D trench of width ∆𝑥 = 2.5 µm; purple color in the background: measured signal after the smoothening bandpass filtering. Top right: Contrast 𝐶(𝑧) of the presented interferogram (purple); green: Gaussian fit with 𝐶 fit (𝑧): = 𝑎 t exp[-𝑏 t (𝑧 -𝑧 t ) 2 ] + 𝑎 b exp[-𝑏 b (𝑧 -𝑧 b ) 2 ] + 𝑐𝑜𝑛𝑠𝑡. The selected interferogram is an example of constructive interference between the signals reflected at TSV top and bottom. The common envelope-based technique leads to underestimating the distance 𝐻 = 𝑧 b -𝑧 t due to this constructive interference (see tab. 1).

Fig. 20 .

 20 Fig. 20. Top left: measured interferogram before bandpass filtering. Bottom left, green: interferogram simulated with best matching parameter set for linear polarization at a circular hole of diameter 𝐷 = 5 µm; purple: color in the background: measured signal after bandpass filtering. Top right: Contrast 𝐶(𝑧) of the presented interferogram (purple). Green: Gaussian fit as in fig. 19. The selected interferogram is an example of destructive interference between the signals reflected at TSV top and bottom. The common envelope-based technique concludes a too large distance 𝐻 = 𝑧 b -𝑧 t here due to this destructive interference (see tab. 3).

  -𝐻 ̅ ) 2 .

  c )]. The resulting merit function reads 𝑚 ̃(𝑎 sim,t , ∆𝑧 t , 𝑎 sim,b , 𝐻): = ∑ [𝑐 sim,TSV,𝑗 (𝑎 sim,t , ∆𝑧 t , 𝑎 sim,b , 𝐻) -𝑐 TSV,𝑗 ]with 𝑐 sim,TSV,𝑗 (𝑎 sim,t , ∆𝑧 t , 𝑎 sim,b , 𝐻) ∶= 𝑎 sim,t 𝑓(𝑘 𝑗

	𝑁 𝑗=1	2	+ cc,

Table 1 . Results for H obtained from TE-polarization at 1D trenches

 1

	Trench width	𝐻/µ𝑚 from envelope	𝐻 ̅ : weighted averaged	from novel	𝐻 ̅ : weighted	|	𝑎 b 𝑎 t	• exp(-2𝑘 c 𝑛 eff 𝐻)|
	∆𝑥/	detection plus	𝐻/µ𝑚 (15)	technique	averaged	(values from electromag.
	µ𝑚	Gaussian fits	from Gaussian	with correct	𝐻/µ𝑚 (15)			analysis)
			fits	𝑛 eff	from novel
					technique
	3.0	18.9 ± 0.1	18.9	17.1 ± 0.1	17.2			0.69
	3.0	18.9 ± 0.1		17.2 ± 0.1				0.71
	2.5	14.6 ± 0.1	15.2	16.8 ± 0.2	16.4			1.52
	2.5	15.7 ± 0.1		16.0 ± 0.2				1.47
	2.0	16.6 ± 0.1	16.9	15.8 ± 0.2	15.8			0.97
	2.0	17.1 ± 0.1		15.9 ± 0.2				0.91

Comparison of results for the TSV height 𝐻 obtained from TE-polarization at 1D trenches. The weighted averaged 𝐻 ̅ /µ𝑚 obtained by the novel technique is more accurate than 𝐻 ̅ /µ𝑚 obtained by Gaussian fits of the envelope of signals. The first sub-line refers to +𝑧direction of the reference mirror, the second one to -𝑧direction each.

Table 2 . Results for 𝑯 obtained from TM-polarization at 1D trenches

 2 Same as tab. 1 for TM-polarization for the sake of completeness. The signals from the bottom are much weaker here, leading to a larger influence of noise and errors in case of any technique. We recommend the use of TE-polarization at 1D trenches instead.

	Trench width ∆𝑥/µ𝑚	𝐻/µ𝑚 from envelope detection plus	𝐻 ̅ : weighted averaged 𝐻/µ𝑚 (15)	𝐻/µ𝑚 from novel technique	𝐻 ̅ : weighted averaged 𝐻/µ𝑚 (15)	| (values from electromag. 𝑎 b • exp(-2𝑘 c 𝑛 eff 𝐻)| 𝑎 t
		Gaussian fits	from	with correct	from novel	analysis)
			Gaussian fits	𝑛 eff	technique
	3.0	19.9 ± 0.1	20.0	16.3 ± 0.2	17.4	0.48
	3.0	20.1 ± 0.1		17.7 ± 0.1		0.45
	2.5	19.2 ± 0.1	19.1	18.2 ± 0.3	17.8	1.26
	2.5	18.9 ± 0.1		17.6 ± 0.2		1.34
	2.0	19.6 ± 0.1	19.7	17.6 ± 0.2	18.1	0.50
	2.0	19.8 ± 0.1		18.7 ± 0.2		0.50

Table 3 . Results for 𝑯 obtained from linear polarization at circular holes

 3 Same as tab. 1 for linear polarization at circular holes. The signals from the bottom are weak here for small diameters 𝐷 ≤ 3.0 µm. The weighted averaged 𝐻 ̅ /µ𝑚 obtained by the novel technique is more reliable and accurate than 𝐻 ̅ /µ𝑚 obtained by the conventional technique (envelope detection and subsequent fit of the envelope by a sum of two Gaussian functions). The conventional technique does not resolve the two signals from TSV top and bottom properly in two cases (𝐷 = 2.5 µm and 𝐷 = 3.0 µm) in contrast to the novel technique. One reason for this is the incorporation of less information into the analysis by the conventional technique (merely the fringe contrast) compared to the novel technique (fringe contrast, fringe positions, interference of fringe patterns, reference spectrum, etc.). So, the novel technique performs better in the case of weak signals close to the noise level as well as interfering signals from TSV top and bottom.

	Diameter	𝐻/µ𝑚 from	𝐻 ̅ : weighted	𝐻/µ𝑚 from	𝐻 ̅ : weighted
	𝐷/µ𝑚	envelope	averaged	novel	averaged
		detection plus	𝐻/µ𝑚 (15) from	technique with	𝐻/µ𝑚 (15) from
		Gaussian fits	Gaussian fits	correct 𝑛 eff	novel technique
	12.0	20.7 ± 0.1	20.9	19.7 ± 0.1	19.7
	12.0	21.0 ± 0.1		19.9 ± 0.2
	7.0	16.4 ± 0.1	16.4	17.0 ± 0.3	17.4
	7.0	16.4 ± 0.1		17.6 ± 0.2
	5.0	19.7 ± 0.2	19.7	15.9 ± 0.1	15.9

The spectrum of the light source of the time domain OCT device applied in this paper is in the range of 0.9 𝜆 c ≤ 𝜆 ≤ 1.1 𝜆 c with 𝜆 c = 1.329 µm the center wavelength, which has the advantage that silicon is transparent there, allowing also measurements from the backside of the wafer when needed[START_REF] Alliata | Progress in Time-Domain Optical Coherence Tomography for TSV / 3Di stacking Metrology[END_REF][START_REF] Devanciard | Through Silicon Via Process Characterization by Integrated Inspection/Metrology Solutions in Visible and Infrared Domain[END_REF][START_REF] Courteville | Contact free on-axis metrology for the fabrication and testing of complex optical systems[END_REF][START_REF] Wilhelm | A novel, low coherence fibre optic interferometer for position and thickness measurements with unattained accuracy[END_REF].The reliability and longitudinal accuracy of OCT devices depends strongly on the processing of the measured signal, which itself depends on the underlying physical models. Rather popular modelling approaches in the case of time-domain OCT applied to TSVs assume a Gaussian spectrum of the light source, the absence of interferometer dispersion and apply ray tracing at the diffractive structures[START_REF] Alliata | Progress in Time-Domain Optical Coherence Tomography for TSV / 3Di stacking Metrology[END_REF][START_REF] Courteville | Contact free on-axis metrology for the fabrication and testing of complex optical systems[END_REF][START_REF] Wilhelm | A novel, low coherence fibre optic interferometer for position and thickness measurements with unattained accuracy[END_REF]. Then, the signal processing is based on the calculation of the envelope positions of the signals (fringe patterns) associated respectively with the top and with the bottom of the TSV[START_REF] Courteville | Contact free on-axis metrology for the fabrication and testing of complex optical systems[END_REF][START_REF] Wilhelm | A novel, low coherence fibre optic interferometer for position and thickness measurements with unattained accuracy[END_REF] and neglects the phase information. Such envelope calculation applied to the fringe pattern (fig.2) utilizes the equations in[START_REF] Larkin | Efficient nonlinear algorithm for envelope detection in white light interferometry[END_REF]. While such approximations may be sufficient in the case of large structures, they reach their limits for objects of small widths or diameters, where diffraction occurs, as well as for shallow objects, where the signals reflected at top and bottom overlap and interfere (fig.2, bottom), requiring phase considerations. The z-positions of the reflecting substrate and TSV bottom associated with the fringe patterns can only be concluded accurately by an appropriate, comprehensive data processing procedure including the correct modal propagation constants at the TSV. A difference of 1 or 2 µm in the concluded TSV depth between a simple and a good evaluation of an OCT signal is often realistic according to our experience.

Acknowledgements

We acknowledge the technical support of Unity-SC as well as Fogale Nanotech, which have provided the OCT device for the measurements. In particular, we want to thank Alexandre Tarnowka from Unity-SC for support, as well as Eric Legros and J.-P. Piel from Fogale Nanotech for introducing us into the soft-and hardware of the OCT device.

Funding

This work has been funded in the framework of the OLOVIA Project ANR 15-CE24-0028 supported by the French National Agency of Research and Institut d'Optique Graduate School.