Robust Numerical Solver for Nonlinear Semiconductor Problems

© 2025. Personal use of this material is permitted. However, permission to reprint/republish this material for advertising or promotional purposes or for creating new collective works for resale or redistribution to servers or lists, or to use any copyrighted component of this work in other works must be obtained from the IEEE. Link to publisher version with DOI: 10.1109/TMTT.2025.3576061 Robust Numerical Solver for Nonlinear Semiconductor Problems Mario Kupresak, Bruno Eckmann, Johannes Hoffmann, Michael Baumann, Philippe Peter, Jasmin Smajic, Senior Member, IEEE, and Juerg Leuthold, Fellow, IEEE Abstract—In this work, we develop a numerical solver, efficiently and robustly treating highly nonlinear semiconductor device problems. Beyond the capabilities of commercial tools, the solver can compute the time-domain capacitance and the spectrum of the device current. The solver is based on the finite element method and employs the successive under-relaxation scheme. Its capability has been assessed and validated in a study of an axisymmetric metal-oxide-semiconductor structure, presenting an archetypal scanning microwave microscopy calibration sample, with both nand p-doped semiconductors, including different excitation sources. Excellent agreement was obtained, when testing the tool against features of a commercial tool. By computing the capacitance for the applied low-frequency bias, combined with a high-frequency probe signal, the spectrum of the current flowing in the structure was evaluated, revealing mixproduct components. This allowed us to verify the solver against measurements, resulting in a very good agreement. Index Terms—Drift-diffusion model, finite element method, scanning microwave microscopy, semiconductor. I. INTRODUCTION OLVING highly nonlinear semiconductor problems in a robust and efficient manner is a challenge to the very day. This is especially important, since semiconductor devices are used in a plethora of applications, in the fields of electronics [1], photonics [2], and information and communication technology [3]. To accurately describe the nonlinear response of these devices, several charge transport models have been proposed in literature [4]. One of the most common models, namely the driftdiffusion model (DDM) [5], has been successfully employed, to study the response of different semiconductor structures, such as p-nand n+-n-n+-junction diodes, field effect transistors, etc. [6], [7], [8], [9], [10]. Essentially, this model investigates the drift and diffusion currents, flowing in the This work was supported by the European Partnership on Metrology, cofinanced by the European Union’s Horizon Europe Research and Innovation Program and by the Participating States, under Grant 23IND03 RF 4 6G. (Corresponding author: Mario Kupresak.) Mario Kupresak, Michael Baumann, Philippe Peter, Jasmin Smajic, and Juerg Leuthold are with the Institute of Electromagnetic Fields (IEF), ETH Zurich, Zurich 8092, Switzerland (e-mail: [email protected]). Bruno Eckmann and Johannes Hoffmann are with the Federal Institute of Metrology METAS, Bern-Wabern 3003, Switzerland. This paper is an expanded version of a conference paper presented at IEEE MTT-S NEMO 2024, Montreal, QC, Canada. structures, based on the concepts of the mobility and diffusivity of the associated charge carriers. We note that, compared with DDM, several more advanced theoretical treatments, such as the hydrodynamic model [8], [9], [10], or the Monte Carlo method [9], have been introduced. These models, however, go beyond the scope of the present study and will be considered in our future works. Computationally, DDM has been previously realized, within the framework of the finite element method (FEM) [11], a well-established and widely used numerical technique for solving partial differential equations. Especially, several related in-house developed solvers have been reported over the past decade by IEF, ETH Zurich [6], [7], [8], [12]. These works mainly consider 1-D and 2-D diodes, mentioned above, excited by DC and AC voltage sources, or an external EM wave, while implementing the continuous Galerkin FEM [11], [13]. The same works also incorporate the successive underrelaxation routine, proven to be robust for solving highly nonlinear systems. Additionally, in [12], the generated results have been confirmed with a commercial tool, namely the Atlas device simulation software of SILVACO [14]. Ultimately, the capabilities of numerical solvers have to be tested against experiments. To this end, the scanning microwave microscopy (SMM) technique may be employed. This technique has been extensively studied, to measure and characterize the properties of semiconductor devices [15], [16], [17], [18], [19]. However, an accurate and efficient computational tool, capable of predicting SMM measurements, is largely missing. In this work, we develop a numerical solver, established based on FEM and the successive under-relaxation routine, that can treat highly nonlinear semiconductor problems. The solver is then applied to explain SMM experiments. In particular, the solver can compute the nonlinear time-domain capacitance and the spectrum of the device current, thus exceeding the capabilities of commercial tools. Using the above features, the solver has been compared with measurements. To the best of our knowledge, the implemented tool, bridging the aforementioned computational gap, has not yet been reported in literature. We have employed the solver to study an axisymmetric metal-oxide-semiconductor (MOS) structure, a commonly encountered SMM calibration sample [17], [18], [19], constituted of a doped semiconductor substrate, an insulator layer, and metallic contacts, located at the top and the bottom of the structure. The analysis has been S 2 performed for both nand p-doped substrates, and DC and AC voltage sources. The preliminary numerical results of this work have been presented in our previous conference paper [20]. Compared with the current work, we have investigated in [20]: 1) a MOS structure with no rotational symmetry and a relatively thick insulator layer, 2) an n-doped substrate, 3) a DC voltage source, 4) only the electron concentration distribution, and no other features. This work is organized as follows. Section II discusses the studied geometry. The implemented solver and the experimental setup are also explained. The numerical results are demonstrated in Section III and discussed in Section IV. II. ANALYSIS SETUP A. Studied Geometry The geometry of the 2-D axisymmetric MOS structure, investigated herein, is depicted in Fig. 1. The indicated size and material parameters have been specified, in accordance with the experimental configuration, elaborated in Section II-C. The axisymmetric cylindrical coordinates are denoted by 𝜌 and 𝑧, 𝜀𝑟 is the relative permittivity of the corresponding silicon and silica regions, and 𝑁𝐴 is the acceptor concentration. Moreover, the metallic contacts, represented by the gold lines, are placed at the top and the bottom of the structure. A bias voltage 𝑉𝐵 has been applied to the former contact, whereas the latter one is grounded. Fig. 1. Studied 2-D axisymmetric MOS geometry: a p-doped silicon substrate, covered by a silica layer, with metallic contacts located at the top and the bottom of the structure (denoted by gold lines). The height of the silica layer is only 50 nm, requiring special attention in the simulations. The studied structure represents a typical SMM calibration sample. B. Numerical Solver Here, we elaborate the fundamental equations and the implementation details of the in-house developed numerical solver for a p-doped substrate, used in the experiments. However, the solver is also applicable to an n-doped substrate, and the associated equations are given in the Appendix. For the studied geometry, the implemented solver is based on the strongly coupled equations, namely the Poisson equation and the continuity equation of holes, defined as [6], [7], [9]: 𝜀0 𝑞∇·[ε𝑟(𝐫)∇𝛷(𝐫,𝑡)]={𝑛(𝐫,𝑡)−𝑝(𝐫,𝑡)+𝑁𝐴(𝐫),𝐫∊𝛺𝑆𝑖 0,𝐫∊𝛺𝑆𝑖𝑂2 (1) 𝜕𝑝(𝐫,𝑡) 𝜕𝑡 =µ𝑝∇·[𝑝(𝐫,𝑡)∇𝛷(𝐫,𝑡)]+𝐷𝑝∇2𝑝(𝐫,𝑡),𝐫∊𝛺𝑆𝑖. (2) In (1) and (2), 𝜀0 is the vacuum permittivity, 𝑞 is the elementary charge, 𝛷 is the electrostatic potential, 𝑛 and 𝑝 are the electron and hole concentrations, respectively, µ𝑝 and 𝐷𝑝 are the hole mobility and diffusion constants, respectively; 𝐫 and 𝑡 denote the space and time variables, respectively; 𝛺𝑆𝑖 and 𝛺𝑆𝑖𝑂2 are the corresponding silicon and silica regions, respectively. In a more general case, the Poisson equation, within the semiconductor region, should also include the donor concentration 𝑁𝐷: 𝜀0 𝑞∇·[ε𝑟(𝐫)∇𝛷(𝐫,𝑡)]=𝑛(𝐫,𝑡)−𝑝(𝐫,𝑡)+𝑁𝐴(𝐫)−𝑁𝐷(𝐫). (3) Note that while (1) needs to be solved for both the silicon and silica regions, (2) is computed only within the former region, since the charge concentrations in the latter region are vanishing. For the studied extrinsic semiconductor substrate, the electron concentration may be obtained from the hole concentration, through the following relation: 𝑛=𝑛𝑖2/𝑝, where 𝑛𝑖 stands for the intrinsic concentration of silicon. As a next step, (1) and (2) have been expanded in terms of the aforementioned axisymmetric cylindrical coordinates 𝜌 and 𝑧. Then, the equations have been normalized, to eliminate possible convergence issues of the studied multiscale problem [6], [7], [21], yielding (the dependence on 𝐫=𝑓(𝜌,𝑧) and 𝑡 has been omitted): 1 𝜌∗𝜕 𝜕𝜌∗(𝜌∗ε𝑟𝜕𝛷∗ 𝜕𝜌∗)+ 𝜕 𝜕𝑧∗(ε𝑟𝜕𝛷∗ 𝜕𝑧∗) ={ε𝑟𝑆𝑖(𝑛∗−𝑝∗+𝑁𝐴∗),𝐫∗∊𝛺𝑆𝑖 0,𝐫∗∊𝛺𝑆𝑖𝑂2 (4) 𝜕𝑝∗ 𝜕𝑡∗=𝛾𝑝∗[1 𝜌∗𝜕 𝜕𝜌∗(𝜌∗𝑝∗𝜕𝛷∗ 𝜕𝜌∗+𝜌∗𝜕𝑝∗ 𝜕𝜌∗)+ 𝜕 𝜕𝑧∗(𝑝∗𝜕𝛷∗ 𝜕𝑧∗) +𝜕2𝑝∗ 𝜕(𝑧∗)2]. (5) In (4) and (5), the normalized variables read: 𝜌∗=𝜌/𝐿0, 𝑧∗= 𝑧/𝐿0, 𝛷∗=𝛷/𝑉𝑇, 𝑛∗=𝑛/𝑛𝑖, 𝑝∗=𝑝/𝑛𝑖, 𝑁𝐴∗=𝑁𝐴/𝑛𝑖, 𝑡∗= 𝑡𝐷0/𝐿0 2, 𝛾𝑝∗=µ𝑝𝑉𝑇/𝐷0=𝐷𝑝/𝐷0, with 𝐿0=√𝜀0𝜀𝑟𝑆𝑖𝑉𝑇/(𝑞𝑛𝑖), 𝐷0=µ𝑝𝑉𝑇, 𝑉𝑇 being the thermal voltage, and the variables without asterisk denoting the unnormalized quantities. For a complete description of the investigated problem, besides the partial differential equations (4) and (5), initial and boundary conditions are required [7], [21]. The initial normalized charge concentrations may be expressed as: 𝑝∗(𝐫∗,𝑡∗=0)=1 2(𝑁𝐴∗+√(𝑁𝐴∗)2+4),𝐫∗∊𝛺𝑆𝑖 (6) 𝑛∗(𝐫∗,𝑡∗=0)=1/𝑝∗(𝐫∗,𝑡∗=0),𝐫∗∊𝛺𝑆𝑖 (7) and the initial solution of the potential may be obtained from the Poisson equation. For the boundary conditions (BCs), at the locations of the metals, we employ the Dirichlet BCs for the normalized potential and charge concentrations distributions, corresponding to the Ohmic contacts, where we have assumed no difference in the work functions between the metal and the semiconductor: 3 𝛷∗(𝐫∗,𝑡∗)={𝑉𝐵(𝑡)/𝑉𝑇,𝐫∗∊𝑙𝑐,𝑡𝑜𝑝 0,𝐫∗∊𝑙𝑐,𝑏𝑡 (8) 𝑝∗(𝐫∗,𝑡∗)=1 2(𝑁𝐴∗+√(𝑁𝐴∗)2+4),𝐫∗∊𝑙𝑐,𝑏𝑡 (9) 𝑛∗(𝐫∗,𝑡∗)=1/𝑝∗(𝐫∗,𝑡∗),𝐫∗∊𝑙𝑐,𝑏𝑡. (10) In (8)–(10), the applied DC or AC bias voltage is denoted by 𝑉𝐵(𝑡), and the line segments of the top and bottom metallic contacts are denoted by 𝑙𝑐,𝑡𝑜𝑝 and 𝑙𝑐,𝑏𝑡, respectively. For the other boundaries, we apply the homogeneous Neumann BC for the normalized potential and the normal component of the normalized hole current density 𝐉𝐩 ∗=−𝛾𝑝∗(𝑝∗∇∗𝛷∗+∇∗𝑝∗): ∇∗𝛷∗·𝐫⊥ =0 (11) 𝐉𝐩 ∗·𝐫⊥ =0 (12) with 𝐫⊥  denoting the unit normal to the corresponding boundary. Note that at the interface between the silicon and silica regions, only the BC in (12) is still valid. Within the framework of the continuous Galerkin FEM [6], [7], to obtain the discretized weak forms of (4) and (5), we apply the conventional method of weighted residuals, expand 𝛷∗ and 𝑝∗ in terms of the linear shape functions, and divide the computational domain 𝛺 into triangular elements. The term 𝜕𝑝∗/𝜕𝑡∗ has been dealt with by the backward difference formula (BDF), for the reasons of numerical stability [11]: 𝜕𝑝∗ 𝜕𝑡∗=(𝑝∗)𝑘+1−(𝑝∗)𝑘 ∆𝑡∗ (13) where the indices of the current and next time moments are denoted by 𝑘 and 𝑘+1, respectively, 𝑘=0,1,2,…,𝑁𝑡−1, with 𝑁𝑡 standing for the number of samples, and ∆𝑡∗ is the normalized time step. Since BDF is an implicit time discretization scheme, all the other time dependent terms in (4) and (5) will be evaluated at the next time moment. Combining FEM and BDF, the following set of matrix equations is obtained: [𝐹(1)]{𝛷∗}𝑘+1=[𝐹(2)]{𝑏}𝑘+1 (14) ([𝐹(3)]+∆𝑡∗𝛾𝑝∗([𝐷]+[𝐹(4)])){𝑝∗}𝑘+1=[𝐹(3)]{𝑝∗}𝑘 (15) 𝐹𝑖𝑗(1)=∑∬ε𝑟(𝜕𝑁𝑖∗ 𝜕𝜌∗𝜕𝑁𝑗∗ 𝜕𝜌∗+𝜕𝑁𝑖∗ 𝜕𝑧∗𝜕𝑁𝑗∗ 𝜕𝑧∗) 𝛺𝑒𝜌∗𝑑𝜌∗𝑑𝑧∗ 𝑁𝑒 𝑒=1 (16) 𝐹𝑖𝑗(2)=∑∬𝜀𝑟𝑆𝑖𝑁𝑖∗𝑁𝑗∗ 𝛺𝑒𝜌∗𝑑𝜌∗𝑑𝑧∗ 𝑁𝑒 𝑒=1 (17) 𝐹𝑖𝑗(3)=∑∬(𝜕𝑁𝑖∗ 𝜕𝜌∗𝜕𝑁𝑗∗ 𝜕𝜌∗+𝜕𝑁𝑖∗ 𝜕𝑧∗𝜕𝑁𝑗∗ 𝜕𝑧∗) 𝛺𝑒𝜌∗𝑑𝜌∗𝑑𝑧∗ 𝑁𝑒𝑆𝑖 𝑒=1 (18) 𝐹𝑖𝑗(4)=∑∬𝑁𝑖∗𝑁𝑗∗ 𝛺𝑒𝜌∗𝑑𝜌∗𝑑𝑧∗ 𝑁𝑒𝑆𝑖 𝑒=1 (19) 𝐷𝑖𝑗=∑[(𝜕𝛷∗ 𝜕𝜌∗)𝑒∬𝜕𝑁𝑖∗ 𝜕𝜌∗𝑁𝑗∗ 𝛺𝑒𝜌∗𝑑𝜌∗𝑑𝑧∗ 𝑁𝑒𝑆𝑖 𝑒=1 +(𝜕𝛷∗ 𝜕𝑧∗)𝑒∬𝜕𝑁𝑖∗ 𝜕𝑧∗𝑁𝑗∗ 𝛺𝑒𝜌∗𝑑𝜌∗𝑑𝑧∗] (20) {𝑏}={𝑝∗}−{𝑛∗}−{𝑁𝐴∗},𝑣∈𝛺𝑆𝑖. (21) In (14)–(21), 𝑖 and 𝑗 are the indices of the corresponding nodes, 𝑒 and 𝑣 are the indices of the corresponding elements and nodes, respectively, 𝛺𝑒 stands for a single element of the computational domain, 𝑁𝑒 is the total number of elements, 𝑁𝑒𝑆𝑖 is the number of elements in the silicon region, and 𝑁𝑖 and 𝑁𝑗 are the shape functions. Furthermore, the size of [𝐹(1)] and [𝐹(2)] is 𝑁×𝑁, where 𝑁 denotes the total number of nodes, and the size of [𝐹(3)], [𝐹(4)], and [𝐷] is 𝑁𝑆𝑖×𝑁𝑆𝑖, with 𝑁𝑆𝑖 representing the number of nodes in the silicon region. Furthermore, the size of {𝛷∗} and {𝑏} is 𝑁×1, where the entries of {𝑏} for the nodes of 𝛺𝑆𝑖 are computed using (21), with the size of the right-hand side quantities being 𝑁𝑆𝑖×1, and the corresponding entries for the nodes of 𝛺𝑆𝑖𝑂2 are vanishing. To solve the nonlinear system of equations (14) and (15), we have employed a successive under-relaxation scheme [6], [7], defined in the following way: {𝑝∗}𝑘+1,𝑚+1=(1−𝛼){𝑝∗}𝑘+1,𝑚+𝛼{𝑝∗}𝑘+1,𝑒𝑠𝑡. (22) In (22), 𝛼 is the relaxation factor, ranging from 0 to 1, 𝑚 is the index of iterations, and {𝑝∗}𝑘+1,𝑒𝑠𝑡 is the estimated hole concentration at the next time moment, computed by (15). Note that {𝑝∗}𝑘+1,𝑚=0={𝑝∗}𝑘. For a given normalized time moment 𝑡∗=(𝑘+1)∆𝑡∗, starting from 𝑚=0, first of all, {𝛷∗}𝑘+1 is computed from (14). Based on the potential distribution, the aforementioned matrix [𝐷] may be obtained. Then, {𝑝∗}𝑘+1,𝑒𝑠𝑡 is calculated using (15). If the error δ= |{𝑝∗}𝑘+1,𝑒𝑠𝑡−{𝑝∗}𝑘+1,𝑚|/|{𝑝∗}𝑘+1,𝑒𝑠𝑡| is larger than a limiting value, taken to be 10−4, one calculates {𝑝∗}𝑘+1,𝑚+1 from (22) and uses the same value to replace {𝑝∗}𝑘+1 in (14), so that the updated {𝛷∗}𝑘+1, [𝐷], {𝑝∗}𝑘+1,𝑒𝑠𝑡, and δ can be computed. This iterative procedure continues, until the above error threshold condition has been satisfied. When this is the case, {𝑝∗}𝑘+1,𝑒𝑠𝑡 is the final solution for the hole concentration at the next time moment in (15), i.e., {𝑝∗}𝑘+1={𝑝∗}𝑘+1,𝑒𝑠𝑡. Lastly, the same procedure is done for all the time moments. C. Experimental Configuration We describe here an experimental technique, to verify the developed numerical tool. In SMM, a scanning microwave signal is coupled to the tip of an atomic force microscope (AFM), and the complex-valued scattering parameter 𝑆11, that is, the ratio of the reflected and incident waves, is recorded. The technique is particularly suitable to detect small variations in capacitance [22]. Furthermore, a specific SMM operating mode, referred to as the dS/dV (or dC/dV), is usually employed. Here, a low-frequency (LF) bias modulates the charge carrier density in a semiconductor. A high-frequency (HF) signal, applied to the tip, and the above LF one, are mixed inside the semiconductor, leading to sidebands at mixproduct frequencies. Typically, this mode is used to distinguish between nand p-doped regions in semiconductors, and sometimes it is also used to determine the corresponding doping concentrations. 4 The sample, investigated in this study, is a development version of a now commercially available calibration substrate Fig. 2. Experimental setup: the scanning microwave interface is connected to a conductive SMM probe of an AFM. A LF bias 𝑉𝐵 at 28.125 kHz and a RF signal at 4.835 GHz are applied. The sample, shown in the inset, is a MOS microcapacitor formed by a silica layer, sandwiched between an Au/Ti pad and a p-type silicon substrate. for SMM from MC2 Technologies [23]. It consists of 64 identical sets, each containing 48 MOS micro-capacitors formed by Au/Ti pads of different sizes, a silica layer, with thickness ranging from 50 nm to 200 nm, and a p-doped silicon substrate, with a resistivity of 1 Ω∙cm. The experimental setup is schematically depicted in Fig. 2. A scanning microwave interface (Nanosurf) is connected to a conductive probe (25Pt300D, Rocky Mountain Nanotechnology) of an AFM. As mentioned earlier, SMM is operated in the dS/dV mode, in which a LF bias at 𝑓𝐿𝐹 𝑒𝑥𝑝= 28.125 kHz, with the adjustable amplitude 𝑉𝐵, can be applied. The RF frequency is set to 𝑓𝑅𝐹 𝑒𝑥𝑝=4.835 GHz, at which the coaxial cable, the probe and the impedance match form a resonant circuit, where the sensitivity of the instrument is enhanced. The device records the magnitude and the phase of the mix-product at 𝑓𝑅𝐹 𝑒𝑥𝑝+𝑓𝐿𝐹 𝑒𝑥𝑝. In this study, we have measured a structure formed by an Au/Ti pad, with a diameter of 4 µm, on top of a 50 nm silica layer. The tip is in galvanic contact with the pad. The RF and LF signals are conducted towards the tip and the pad through the RF cable. Note that the return paths for the LF and RF signals are quite different. The RF signal is reflected by the sample and travels back to the scanning microwave interface through the cable. The mix-products are generated in the semiconductor and travel as well back through the cable. The LF one is, in turn, only weakly coupled back through this path to the cable, because the distance between the outer conductor and the sample is too large for the LF, or in other words, the capacitor formed between the outer and central conductors is too small. Furthermore, the same signal travels through the sample and is capacitively coupled to the sample holder through the native silica layer at the bottom of the sample, with a thickness of only a few nanometers. The sample holder is coupled to the outer conductor of the RF cable through a relatively long wire. The inductance of this wire prevents the RF signals to travel along the wire, but conducts the LF ones. As a consequence, the applied LF voltage drops across the MOS contact. This technique has the advantage of requiring minimum sample preparation. The tip is scanned over a very small region of the micro-capacitor, while manually increasing the bias amplitude 𝑉𝐵. Note that, for the simulated topology, see Fig. 1, the Au/Ti pad is simply represented by the top metallic contact, with a fixed voltage, and that the native silica layer, between the substrate and the bottom contact, has not been considered, due to the very large capacitance. III. RESULTS In the analysis, we have assumed a uniform doping profile of the silicon substrate. To estimate the electron and hole mobilities in silicon, we have employed a well-known analytical model [24], that for the studied doping concentration and the lattice temperature of 300 K yields the following values: µ𝑛=1207 cm2/(V∙s) and µ𝑝= 271 cm2/(V∙s). We note that several other mobility models have been proposed in literature, capturing more advanced effects, such as the lattice, impurity, or carrier-carrier scattering, the saturation velocity, the surface roughness, etc. For a detailed discussion of the associated models, we refer the readers to the Atlas user manual [14] (and references therein). Moreover, to ensure stable solutions of (14) and (15), we have specified the time step ∆𝑡 around 9 ps. A very fine discretization of the thin silica layer and regions in the vicinity of the silicon–silica interface, associated with the depletion of the charge carriers, are required, to obtain accurate results. For the successive under-relaxation scheme, 𝛼=0.1 has been used, substantiated by numerous simulations, to ensure convergence of the studied features, indicating a highly nonlinear system. First of all, to evaluate the performance of the in-house developed tool, the generated results have been compared with the ones of Atlas, for the applied DC bias voltage. The associated results are depicted in Fig. 3, demonstrating an excellent agreement between the studied solvers. The solutions of (14) and (15) have been computed with our tool for multiple time moments, reaching the steady-state, whereas the results of Atlas have been obtained without the transient analysis. For the comparison, three bias voltages, namely 𝑉𝐵= (0,1,2) V, have been applied to the top contact. For the other study case, associated with the dS/dV measurements of the mix-product harmonic, when applying a LF bias signal, the time-domain capacitance of the studied structure shows nonlinear oscillations, as illustrated in Fig. 4. This observation has been made by considering a LF bias voltage 𝑉𝐿𝐹(𝑡)=𝐴𝐿𝐹sin(2𝜋𝑓𝐿𝐹 𝑠𝑖𝑚𝑡) V, with the amplitude 𝐴𝐿𝐹=0.4 and 𝑓𝐿𝐹 𝑠𝑖𝑚=28 kHz. The total capacitance of the MOS structure 𝐶 has been computed, using the following relation: 𝐶(𝑡)=2𝑊𝑒(𝑡) 𝑉𝐿𝐹 2(𝑡),𝑊𝑒(𝑡)=ε0 2∬ε𝑟 𝛺|𝐄(𝑡)|2𝜌𝑑𝜌𝑑𝑧 (23) where 𝑊𝑒(𝑡) is the total electric energy stored in the structure, with 𝐄(𝑡)=−∇𝛷(𝑡) representing the electric field. We have computed the capacitance for several bias voltages at different time moments, in the steady-state, as for the previously conducted DC simulations. Furthermore, the corresponding capacitance values at all other studied moments have been obtained through the data interpolation, in this way 5 eliminating numerical artifacts that occur, as the bias voltage is close to 0. Fig. 3. Comparison of the numerical results generated by the in-house developed solver and Atlas: (a) the studied structure; (b) the potential and (c) hole concentration distributions for three different bias voltages: 𝑉𝐵=0 𝑉, 𝑉𝐵=1 𝑉, and 𝑉𝐵=2 𝑉. Notably, the in-house developed solver and Atlas provide identical results. The plots have been evaluated at a fixed horizontal position 𝜌=1 µm, while varying 𝑧, denoted by the purple dotted line in (a). To compare the results of the developed solver and Atlas more easily, a narrow range of 𝑧 is demonstrated. For 𝑧<9 µm (not displayed here), the potential and hole concentration distributions are constant, with the values of -0.36 V and 1.7·1016 cm-3, respectively, same as the ones of the constant parts of the plots in (b) and (c). The vertical black line at 𝑧=10 µm illustrates the silicon–silica interface, see (b) and (c). Fig. 4. The time-domain LF bias voltage 𝑉𝐿𝐹(𝑡)= 0.4𝑠𝑖𝑛(2𝜋𝑓𝐿𝐹 𝑠𝑖𝑚𝑡) V (blue) and the associated total capacitance (red). To emphasize the nonlinear behavior of the computed capacitance, the capacitance containing only one frequency component, 𝑓𝐿𝐹 𝑠𝑖𝑚, is plotted as a reference (black). Fig. 5. The magnitudes of the currents for different harmonics, obtained by mixing the capacitance in Fig. 4, and a RF probe signal 𝑉𝑅𝐹(𝑡)=0.2𝑠𝑖𝑛(2𝜋𝑓𝑅𝐹 𝑠𝑖𝑚𝑡)𝑉. The spectrum shows the mix products at 𝑓𝑅𝐹 𝑠𝑖𝑚±𝑓𝐿𝐹 𝑠𝑖𝑚 and 𝑓𝑅𝐹 𝑠𝑖𝑚±2∙𝑓𝐿𝐹 𝑠𝑖𝑚. Fig. 6. The magnitudes of the simulated currents (blue) and the dS/dV measurements (red) at the mix-product frequencies 𝑓𝑅𝐹 𝑠𝑖𝑚+𝑓𝐿𝐹 𝑠𝑖𝑚 and 𝑓𝑅𝐹 𝑒𝑥𝑝+𝑓𝐿𝐹 𝑒𝑥𝑝, respectively. The simulation and experimental results are in very good agreement, with NRMSE being 7.6%. Considering an HF signal, on top of the above LF signal, nonlinear mixing may be observed, revealing the mix-product components in the spectrum of the current, as depicted in Fig. 5. To this end, a RF probe signal 𝑉𝑅𝐹(𝑡)= 0.2sin(2𝜋𝑓𝑅𝐹 𝑠𝑖𝑚𝑡)V, with 𝑓𝑅𝐹 𝑠𝑖𝑚=4.9 GHz, has been applied. Further, the charge 𝑄(𝑡)=𝐶(𝑡)𝑉𝑅𝐹(𝑡), with 𝐶 being the previously calculated capacitance, has been evaluated. Then, the spectrum of the current 𝐼, flowing in the structure, may be obtained by taking the Fourier transform of the expression 𝐼(𝑡)=𝑑𝑄(𝑡)/𝑑𝑡, that is, 𝐼(𝑓)=−𝑗2𝜋𝑓𝑄(𝑓), where 𝑓 is the frequency. Here, several periods of the LF signal have been considered, to ensure accurate spectrum resolution, detecting different frequency components. Note that the aforementioned frequencies 𝑓𝐿𝐹 𝑠𝑖𝑚 and 𝑓𝑅𝐹 𝑠𝑖𝑚, used in the simulations, are slightly different than 𝑓𝐿𝐹 𝑒𝑥𝑝 and 𝑓𝑅𝐹 𝑒𝑥𝑝, respectively, explained 6 in Section II-C, to compute the Fourier transform more efficiently. To the best of our knowledge, the above analyses of the LF capacitance and the spectrum of the device current cannot be completely performed by Atlas, Synopsis TCAD [25], and Ansys Lumerical [26]. The magnitudes of the simulated currents and the dS/dV measurements at the mix-product frequencies are plotted in Fig. 6, for the extended range of LF bias voltages, where the amplitude 𝐴𝐿𝐹 varies from 0.1 to 2, with a step of 0.1. The associated results are in very good agreement. IV. DISCUSSION It can be readily seen that the results of the potential and hole concentration distributions, generated by our in-house developed tool and Atlas, plotted in Fig. 3(b) and (c), are nearly identical. The agreement between the two tools has been evaluated in terms of a normalized root mean square error (NRMSE) [27], being around 1% for both studied features. In the silicon substrate, the potential is mainly constant, with the value of -0.36 V, for all the bias voltages, as illustrated in Fig. 3(b). This value is obtained by considering the corresponding boundary condition for the potential at the bottom Ohmic contact (𝑧=0), implemented in Atlas, for the studied p-doped semiconductor [14]. While approaching the silicon–silica interface, the potential changes drastically, showing nonlinear behavior, due to the depletion of holes, discussed below. Lastly, in the very thin silica layer, with no charge carriers present, the potential increases linearly, reaching the value of the applied bias at the final point, i.e., 𝑧=10.05 µm, where the top contact is located. More interestingly, in the vicinity of the above interface, the steady-state hole concentrations deviate from the doping profile, caused by the depletion of the charge carriers, attracted by the lower potential at the bottom contact. As intuitively expected, while increasing the bias voltage, the depletion regions are extended, see Fig. 3(c). Nevertheless, the same regions are rather small, as a consequence of the relatively high doping concentration. In contrast, in our related conference paper [20], with a roughly two orders of magnitude lower concentration, a considerably larger depletion region has been demonstrated, despite using a much thicker silica layer. Besides the aforementioned DC analysis, we have also conducted the simulations with AC sources, in the context of the dS/dV measurements. Here, as a first step, for the applied LF bias voltage, the total capacitance of the MOS structure has been computed, revealing the following features. First of all, as depicted in Fig. 4, the obtained time-domain capacitance clearly shows nonlinear behavior, deviating from the reference capacitance, assumed to be a function of only one frequency, namely, the frequency of the bias voltage 𝑓𝐿𝐹 𝑠𝑖𝑚. Furthermore, the capacitance displays an opposite trend, with respect to the input voltage. More precisely, in the time intervals, where the voltage increases, the capacitance decreases, and vice versa. Such an observation may be explained by the fact that for a positive voltage at the top contact, while increasing the voltage, the holes are more attracted by the grounded bottom contact, as mentioned above. As a result, the depletion regions become larger, thus reducing the capacitance. For the decreasing voltage, fewer holes are depleted, so that the capacitance increases. Within the next time interval, where the bias voltage is negative and decreasing, the holes become more confined, close to the silicon–silica interface, causing the capacitance to further increase. Lastly, when the same voltage increases, the capacitance is obviously decreased. Next, taking the above LF capacitance and the RF probe signal into account, the spectrum of the current may be calculated, demonstrating several components in Fig. 5. Indeed, besides the dominant component at 𝑓𝑅𝐹 𝑠𝑖𝑚, as a consequence of the nonlinear mixing of the aforementioned signals, the harmonics at 𝑓𝑅𝐹 𝑠𝑖𝑚±𝑓𝐿𝐹 𝑠𝑖𝑚 are clearly visible. On top, one may also observe two other harmonics, appearing at 𝑓𝑅𝐹 𝑠𝑖𝑚±2·𝑓𝐿𝐹 𝑠𝑖𝑚, whose magnitudes are significantly lower. Additionally, the spectrum contains several other higher-order components (not displayed here), with even lower magnitudes. We emphasize that the above studies of the capacitance and the device current, as far as we know, cannot be fully performed using Atlas, or some other commercially available semiconductor solvers, such as Synopsys TCAD and Ansys Lumerical. Finally, the magnitudes of the mix-product at 𝑓𝑅𝐹 𝑠𝑖𝑚+𝑓𝐿𝐹 𝑠𝑖𝑚 have been computed for different LF bias voltages. On top, for the same voltages, the dS/dV measurements at 𝑓𝑅𝐹 𝑒𝑥𝑝+𝑓𝐿𝐹 𝑒𝑥𝑝 have been performed. As depicted in Fig. 6, both quantities display very similar trends, and they are, overall, in very good agreement, with NRMSE being 7.6%. Notably, the agreement between the results may be improved by tuning several parameters in the numerical implementation. For example, one may use a more detailed description of the contacts, and extend the current Ohmic contacts to allow for different work functions between the metal and the semiconductor. Additionally, including more complicated doping profiles or mobility models, discussed above, may be also considered. We will perform a thorough analysis of the above parameters and their impact on the simulation results in future studies. V. CONCLUSION This work demonstrated a versatile in-house developed computational tool, capable of solving highly nonlinear semiconductor problems. The tool was applied to study a typical SMM calibration sample, namely the MOS structure, incorporating the rotational symmetry, both nand p-doped semiconductors, and DC and AC voltage sources. The tool was implemented based on the continuous Galerkin FEM and the successive under-relaxation scheme, for solving a highly nonlinear systems. To substantiate the numerical study, SMM measurements were carried out. First, the analysis was conducted by investigating the steady-state response of the MOS structure, for different DC bias voltages. The generated potential and carrier concentration distributions were confirmed with the ones of a commercially available solver, namely Atlas. Then, surpassing the capabilities of commercial solvers, the tool was employed to compute the total time- 7 domain capacitance of the structure for a LF bias voltage. By combining the LF capacitance and a RF probe signal, the spectrum of the device current was obtained, displaying mixproduct components. A very good agreement between the simulation results of the magnitudes of the currents and the dS/dV measurements at the corresponding mix-product frequencies was reported. APPENDIX For the case of an n-doped semiconductor substrate, where the majority charge carriers are electrons, the Poisson equation and the continuity equation of electrons are defined as follows: 𝜀0 𝑞∇·[ε𝑟(𝐫)∇𝛷(𝐫,𝑡)]={𝑛(𝐫,𝑡)−𝑝(𝐫,𝑡)−𝑁𝐷(𝐫),𝐫∊𝛺𝑆𝑖 0,𝐫∊𝛺𝑆𝑖𝑂2 (24) 𝜕𝑛(𝐫,𝑡) 𝜕𝑡 =−µ𝑛∇·[𝑛(𝐫,𝑡)∇𝛷(𝐫,𝑡)]+𝐷𝑛∇2𝑛(𝐫,𝑡),𝐫∊𝛺𝑆𝑖. (25) In (24) and (25), µ𝑛 and 𝐷𝑛 are the electron mobility and diffusion constants, respectively. All other parameters are the same as the ones defined in Section II-B, in the main text. The hole concentration may be obtained from the electron concentration, through the relation 𝑝=𝑛𝑖2/𝑛. After expanding (24) and (25) in terms of 𝜌 and 𝑧, and normalizing the equations, the following expressions are obtained: 1 𝜌∗𝜕 𝜕𝜌∗(𝜌∗ε𝑟𝜕𝛷∗ 𝜕𝜌∗)+ 𝜕 𝜕𝑧∗(ε𝑟𝜕𝛷∗ 𝜕𝑧∗) ={ε𝑟𝑆𝑖(𝑛∗−𝑝∗−𝑁𝐷∗),𝐫∗∊𝛺𝑆𝑖 0,𝐫∗∊𝛺𝑆𝑖𝑂2 (26) 𝜕𝑛∗ 𝜕𝑡∗=𝛾𝑛∗[1 𝜌∗𝜕 𝜕𝜌∗(−𝜌∗𝑛∗𝜕𝛷∗ 𝜕𝜌∗+𝜌∗𝜕𝑛∗ 𝜕𝜌∗)− 𝜕 𝜕𝑧∗(𝑛∗𝜕𝛷∗ 𝜕𝑧∗) +𝜕2𝑛∗ 𝜕(𝑧∗)2] (27) where 𝑁𝐷∗=𝑁𝐷/𝑛𝑖 and 𝛾𝑛∗=µ𝑛𝑉𝑇/𝐷0=𝐷𝑛/𝐷0. The initial distributions of the normalized charge concentrations may be written as 𝑛∗(𝐫∗,𝑡∗=0)=1 2(𝑁𝐷∗+√(𝑁𝐷∗)2+4),𝐫∗∊𝛺𝑆𝑖 (28) 𝑝∗(𝐫∗,𝑡∗=0)=1/𝑛∗(𝐫∗,𝑡∗=0),𝐫∗∊𝛺𝑆𝑖. (29) At the locations of the metals, the Dirichlet BCs for the normalized charge concentrations have been applied: 𝑛∗(𝐫∗,𝑡∗)=1 2(𝑁𝐷∗+√(𝑁𝐷∗)2+4),𝐫∗∊𝑙𝑐,𝑏𝑡 (30) 𝑝∗(𝐫∗,𝑡∗)=1/𝑛∗(𝐫∗,𝑡∗),𝐫∗∊𝑙𝑐,𝑏𝑡. (31) For the other boundaries, besides the normalized potential, the homogeneous Neumann BC has been applied to the normal component of the normalized electron current density 𝐉𝐧 ∗= 𝛾𝑛∗(−𝑛∗∇∗𝛷∗+∇∗𝑛∗), namely 𝐉𝐧 ∗·𝐫⊥ =0. By applying the continuous Galerkin FEM and BDF for the space and time discretizations, respectively, the following set of matrix equations is obtained: [𝐹(1)]{𝛷∗}𝑘+1=[𝐹(2)]{𝑏}𝑘+1 (32) ([𝐹(3)]−∆𝑡∗𝛾𝑛∗([𝐷]−[𝐹(4)])){𝑛∗}𝑘+1=[𝐹(3)]{𝑛∗}𝑘 (33) {𝑏}={𝑝∗}−{𝑛∗}+{𝑁𝐷∗},𝑣∈𝛺𝑆𝑖. (34) Similar to the previously studied p-doped semiconductor, the potential and electron concentration distributions are computed by the successive under-relaxation method, explained in Section II-B, in the main text. REFERENCES [1] J. A. del Alamo, “Nanometre-scale electronics with III–V compound semiconductors,” Nature, vol. 479, no. 7373, pp. 317–323, Nov. 2011. [2] S. Y. Siew et al., “Review of silicon photonics technology and platform development,” J. Light. Technol., vol. 39, no. 13, pp. 4374–4389, July 2021. [3] J. Smajic and J. Leuthold, “Plasmonic electro-optic modulators – a review,” IEEE J. Sel. Top. Quantum Electron., vol. 30, no. 4, p. 3300113, Aug. 2024. [4] T. Grasser, T.-W. Tang, H. Kosina, and S. Selberherr, “A review of hydrodynamic and energy-transport models for semiconductor device simulation,” Proc. IEEE, vol. 91, no. 2, pp. 251–274, Feb. 2003. [5] H. K. Gummel, “A self-consistent iterative scheme for one-dimensional steady state transistor calculations,” IEEE Trans. Electron Devices, vol. 11, no. 10, pp. 455–465, Oct. 1964. [6] A. C. Gungor, M. Celuch, J. Smajic, M. Olszewska-Placha, and J. Leuthold, “Electromagnetic and semiconductor modeling of scanning microwave microscopy setups,” IEEE J. Multiscale Multiphysics Comput. Tech., vol. 5, pp. 209–216, Sep. 2020. [7] A. C. Gungor, J. Smajic, F. Moro, and J. Leuthold, “Time-domain coupled full Maxwelland drift-diffusion-solver for simulating scanning microwave microscopy of semiconductors,” in IEEE PIERS-Spring, Rome, Italy, June 2019, pp. 4071–4077. [8] A. C. Gungor, T. Ehrengruber, J. Smajic, and J. Leuthold, “Coupled electromagnetic and hydrodynamic modeling for semiconductors using DGTD,” IEEE Trans. Magn., vol. 57, no. 6, p. 2400205, June 2021. [9] D. Vasileska, S. M. Goodnick, and G. Klimeck, “Computational electronics: semiclassical and quantum device modeling and simulation,” Boca Raton, FL, USA, CRC Press, 2010. [10] R. Granzner et al., “Simulation of nanoscale MOSFETs using modified drift-diffusion and hydrodynamic models and comparison with Monte Carlo results,” Microelectron. Eng., vol. 83, no. 2, pp. 241–246, Feb. 2006. [11] J.-M. Jin, “The finite element method,” in Theory and Computation of Electromagnetic Fields, Hoboken, NJ, USA, Wiley, 2010, pp. 342–398. [12] M. Kupresak, J. Smajic, and J. Leuthold: “Modeling and simulations of semiconductor structures at highest frequencies,” in IEEE CEFC, Jeju, Republic of Korea, June 2024. [13] J. Smajic, “How to perform electromagnetic finite element analysis,” Hamilton, UK, The International Association for the Engineering Modelling, Analysis & Simulation Community (NAFEMS), NAFEMS Ltd, 2016. [14] SILVACO, https://silvaco.com/tcad/. [15] A. Imtiaz, T. M. Wallis, and P. Kabos, “Near-field scanning microwave microscopy: an emerging research tool for nanoscale metrology,” IEEE Microw. Mag., vol. 15, no. 1, pp. 52–64, Feb. 2014. [16] J. Hoffmann, G. Gramse, J. Niegemann, M. Zeier, and F. Kienberger, “Measuring low loss dielectric substrates with scanning probe microscopes,” Appl. Phys. Lett., vol. 105, no. 1, p. 013102, July 2014. [17] M. Kasper et al., “Metal-oxide-semiconductor capacitors and Schottky diodes studied with scanning microwave microscopy at 18 GHz,” J. Appl. Phys., vol. 116, no. 18, p. 184301, Nov. 2014. [18] A. Buchter et al., “Scanning microwave microscopy applied to semiconducting GaAs structures,” Rev. Sci. Instrum., vol. 89, no. 2, p. 023704, Feb. 2018. [19] J. Hoffmann et al., “Comparison of impedance matching networks for scanning microwave microscopy,” IEEE Trans. Instrum. Meas., vol. 73, p. 6006109, Mar. 2024. [20] M. Kupresak, J. Hoffman, J. Smajic, and J. Leuthold, “Numerical solver for characterization of semiconductor devices,” in IEEE MTT-S NEMO, Montreal, Canada, Aug. 2024. [21] Andre De Mari: “Accurate numerical steady-state and transient onedimensional solutions of semiconductor devices,” PhD thesis, California Institute of Technology, Pasadena, California, 1968. 8 [22] T. M. Wallis and P. Kabos, “Measurement techniques for radio frequency nanoelectronics,” Cambridge, UK, Cambridge University Press, 2017. [23] F. Piquemal, J. Morán-Meza, A. Delvallée, D. Richert, and K. Kaja, “Progress in traceable nanoscale capacitance measurements using scanning microwave microscopy,” Nanomaterials, vol. 11, no. 3, p. 820, Mar. 2021. [24] N. D. Arora, J. R. Hauser, and D. J. Roulston, “Electron and hole mobilities in silicon as a function of concentration and temperature,” IEEE Trans. Electron Devices, vol. 29, no. 2, pp. 292–295, Feb. 1982. [25] Synopsys TCAD, https://www.synopsys.com/manufacturing/tcad.html. [26] Ansys Lumerical, https://www.ansys.com/products/optics/fdtd. [27] M. Kupresak, X. Zheng, V. V. Moshchalkov, and G. A. E. Vandenbosch, “Benchmarking of software tools for the characterization of nanoparticles,” Opt. Express, vol. 25, no. 22, pp. 26760–26780, Oct. 2017. Mario Kupresak received the B.S. and M.S. degrees in electrical engineering and information technology from the University of Zagreb, Zagreb, Croatia, in 2013 and 2015, respectively, and the Ph.D. degree in electrical engineering from KU Leuven, Leuven, Belgium, in 2021. From 2021 to 2022, he was a Postdoctoral Researcher with the ESAT-WAVECORE Group, KU Leuven. Since 2023, he has been a Postdoctoral Researcher with the Institute of Electromagnetic Fields, ETH Zurich, Zurich, Switzerland. His research interests include the modeling of the interaction between the electromagnetic waves and metallic nanostructures, semiconductor devices, and bioelectromagnetics. Bruno Eckmann received the B.Sc. and Mech.Eng. degrees from the University of Applied Sciences Northwestern Switzerland, Windisch, Switzerland, in 2011, and the M.Sc. degree in physics from ETH Zurich, Zürich, Switzerland, in (2019, focusing on quantum optics and quantum information technology. He joined the High Frequency Laboratory, Swiss Federal Institute of Metrology (METAS), Bern, Switzerland, in 2020. His current research interests include scanning probe microscopy with a focus on microwave microscopy and quantum sensing, including hardware and software developments. Johannes Hoffmann received the Dipl.Ing. degree in electrical engineering from the University of Stuttgart, Stuttgart, Germany, in 2005, the Ing. degree from the École Nationale Supérieure des Télécommunications (ENST), Paris, France, in 2005, in the course of a double diploma program, and the Ph.D. degree from ETH Zürich, Zürich, Switzerland, in 2009. He is currently with the Laboratory for RF and MW, Swiss Federal Institute of Metrology (METAS), Bern, Switzerland. His research interests include scanning microwave microscope, and general RF and microwave measurements that involve electromagnetics, uncertainty calculation, and numerical modeling. Michael Baumann received the M.Sc. degree in electrical engineering and information technology from ETH Zürich, Zürich, Switzerland, in 2018. He joined the Institute of Electromagnetic Fields as a Ph.D. student in the same year to work on highspeed plasmonic detectors. His research interests include integrated optics and nonlinear optics in silicon photonics, microwave photonics, and THz technology. Philippe Peter received the M.Sc. degree in computational science and engineering from ETH Zürich, Zürich, Switzerland, in 2024. He is currently working toward the Ph.D. degree with the Institute of Electromagnetic Fields, ETH Zürich. His research interests include computational electromagnetics and semiconductor modeling. Jasmin Smajic (Senior Member, IEEE) received his M.Sc. and Ph.D. degree from the Faculty of Electrical Engineering and Computing in Zagreb (Croatia) in 1998 and 2001, respectively, on the topic of numerical computing and optimization of static and time-varying electromagnetic fields. After his postdoctoral research at ETH Zurich, from 2002 to 2004 on the topic of full-Maxwell electromagnetic simulations of photonic crystals, he took a position of Principal Scientist at the ABB Corporate Research Centre in Baden-Dättwil (Switzerland) where he stayed until 2011. His work in ABB covered a wide range of projects in the field of computational and applied electromagnetics. From 2011 until 2020, he was Professor of Electrical Engineering at the University of Applied Sciences in Rapperswil (Switzerland) where he was leading Computational and Applied Electromagnetics Group. Since 2020, Jasmin Smajic was Senior Scientist and Lecturer at Swiss Federal Institute of Technology (ETH) in Zurich (Switzerland), where he was appointed Professor in 2025. Juerg Leuthold (Fellow, IEEE) is the head of the department of Information Science and Electrical Engineering (D-ITET) and he is the head of the Institute of Electromagnetic Fields (IEF) of ETH Zurich, Switzerland. His interests are in the field of Photonics, Plasmonics and THz with an emphasis on applications in communications and sensing. He was formerly affiliated with the Karlsruhe Institute of Technology (KIT) in Germany, where he was the Head of the Institute of Photonics and Quantum Electronics (IPQ) and the Helmholtz Institute of Microtechnology (IMT) in the time from 2004 to 2013. From 1999 to 2004, he was with Bell Labs, Lucent Technologies, Holmdel, NJ, USA, where he performed device and system research with III/V semiconductors and silicon photonics for applications in high-speed telecommunications. Juerg Leuthold received the Ph.D. degree in physics from ETH Zurich in Switzerland for work in the field of integrated optics and all-optical communications in 1998. Juerg Leuthold is a fellow of the Optica, a fellow of the IEEE, a member of the Swiss Academy of Engineering Sciences (SATW), and a member of the Heidelberg Academy of Sciences. He served the community as member of the Helmholtz Association Think Tank, as a member of the board of directors of the OSA (now Optica), as a general chair, program chair and member of many committees.