Automated Netlist Generation for 3D Electrothermal and Electromagnetic Field Problems

We present a method for the automatic generation of netlists describing general three-dimensional electrothermal and electromagnetic field problems. Using a pair of structured orthogonal grids as spatial discretisation, a one-to-one correspondence between grid objects and circuit elements is obtained by employing the finite integration technique. The resulting circuit can then be solved with any standard available circuit simulator, alleviating the need for the implementation of a custom time integrator. Additionally, the approach straightforwardly allows for field-circuit coupling simulations by appropriately stamping the circuit description of lumped devices. As the computational domain in wave propagation problems must be finite, stamps representing absorbing boundary conditions are developed as well. Representative numerical examples are used to validate the approach. The results obtained by circuit simulation on the generated netlists are compared with appropriate reference solutions.


Introduction
In order to analyse complex electromagnetic (EM) problems, two main routes can be identified. First, a numerical approach for solving the full set of Maxwell's field equations offers the advantage of capturing all relevant effects but can be very demanding in terms of computational resources. One of the first methods for EM analysis was the finite difference time domain (FDTD) scheme proposed by Yee [1] in the 1960s. This method became popular and is still a standard approach for high-frequency EM simulations [2]. About ten years later, in the 1970s, Weiland [3] introduced the finite integration technique (FIT) as an extension to the FDTD scheme using integral unknowns and allowing for non-Cartesian and unstructured grids [4]. Additionally, later developments of the method led to the usage of integral unknowns to obtain an exact implementation of Maxwell's equations [5]. The method proved to be more efficient in terms of memory requirements and computing time. While the finite element method (FEM) was mainly used in structural mechanics for many years, the introduction of edge elements made it applicable for EM problems as well [6]. Secondly, one may use compact models to obtain an efficient representation of a complex system. For example, electrical engineers employ circuits to model and describe the behaviour of complex devices. Nevertheless, the generation of such compact models can be a tedious task requiring empirical knowhow to apply the appropriate approximations. To generate such circuit models, different techniques are available. For instance, mathematical analysis and physical insight allows to construct circuits representing the problem at hand as is done by Choi et al [7]. A circuit's topology and the required component values can also be obtained from experimental results as is done by Moumouni and Baker [8] and other groups. Another approach, e.g. followed by Codecasa et al [9], Eller [10] and Wittig et al [11], is to apply model order reduction (MOR) techniques directly to the field formulation of the problem from which a circuit description can be found more easily [9]. However, the resulting elements may have non-physical values. Now, if one is able to represent an EM problem by means of an electric circuit, one can use any circuit simulator to obtain the solution. The most popular representatives are SPICE programs that were introduced in the 1970s [12] and are still used as a synonym for circuit solvers. Later, extensions to deal with electrothermal (ET) simulations using circuits were developed within the SPICE framework [13,14]. The mathematical tool employed by most SPICE-like programs is still the modified nodal analysis (MNA) presented by Ho et al [15] in the 1970s.
First approaches to combine numerical field simulation with circuit elements were proposed in the 1990s. These were based on the FDTD scheme [16,17,18,19] because of the topological similarities between circuits and the finite difference scheme. Later, the insertion of lumped elements into FEM schemes in time or frequency domain was developed by Guillouard et al [20,21]. These approaches became known as field-circuit coupling and have evolved into an important research topic [22,23,24,25]. One possible field-circuit coupling method is the direct insertion of lumped circuit elements into the field model by applying these to edges of the discretisation grid [26].
To extract compact circuits representing EM problems in a generic way, numerous approaches can be found in the literature. In the partial element equivalent circuit (PEEC) method presented by Ruehli [27,28], equivalent circuits are derived from integral equations, allowing for a combined EM-circuit solution both in frequency and time domain. However, PEEC requires empirical approximations in addition to the applied discretisation. For quasistatic approximations, automated circuit generation based on the boundary element method was presented by Milsom [29]. The method presented therein yields circuits whose size depends on the electrical dimensions of the problem. Many methods for application-specific circuit extraction based on device or system responses are also available [30,31,32]. A methodology for the generation of equivalent circuits based on the semi-discrete Maxwell's field equations was also proposed by Ramachadran et al [33]. Therein, Yee's discretisation scheme is employed to cast Maxwell's curl equations as the concatenation of interacting fundamental circuits in which voltages and currents model the sought electric and magnetic fields. In this manner, circuit stamps are required for both primal and dual edges. Their interaction is organised by voltage controlled voltage sources (VCVSs) and current controlled current sources (CCCSs), respectively.
In the design of electronic devices, enhancing their functionalities is always of high interest. The according volume shrinking may give rise to high power densities that can lead to thermal issues. As an example, the introduction of stacked 3D chips intensifies the heat issue since the heat can be trapped between the stacked layers. Therefore, ET modelling is of great importance for device engineers. To handle the ET coupling in a circuit simulation framework, two general approaches are mainly available: the relaxation method which consists of the iterative coupling of an electric circuit with an external thermal-only field simulator [34,35,36], and the monolithic approach which consists of the direct coupling of the electric circuit with a thermal circuit [37]. The latter allows to run the simulation directly on the full ET circuit without any software package and thereby avoids the weak coupling between solvers. To extract ET circuits from a given 3D problem, various methods were proposed. Some of them were based on existing EM simulation methods and have been extended in functionality to also cover the ET case, as done by Lombardi et al [38] for the PEEC. Generating compact models from the calculated or measured response function is another popular approach and has been followed by Evans et al and Bernardoni et al [39,40]. For thermal problems, methods that derive an equivalent circuit directly from the mesh can also be found [41,42], but none of them accounts for the ET coupling. Karagol and Bikdash [43] presented an approximate representation obtained from a graph-partitioning algorithm of an FEM mesh resulting in a medium-sized ET circuit. A lumped-element representation of every FEM element has also been proposed for ET simulations [44]. For an exact representation of a semi-discretised 3D ET field problem, Casper et al [45] developed an automatic netlist generation method based on the FIT.
In this paper, we present a method to automatically generate netlists representing general 3D ET and EM coupled field problems. In our approach, neighbouring cells in the primal grid interact via the parallel connection of circuit elements as illustrated in Figure 1a. Thus, the size of the resulting circuit depends on the geometrical size of the problem and the fineness of the discretisation grid. This allows to use any available circuit simulator. Hence, the need for a dedicated field solver with a custom time integrator is alleviated. To accomplish this, we employ the FIT for discretising the relevant continuous field equations. In contrast to Yee's finite difference scheme employed by Ramachandran et al [33], FIT is a structure-preserving discretisation strategy which does not require further approximations in dealing with the field and material quantities. Moreover, the concept of integral quantities used in the FIT translates naturally into the framework of circuit descriptions. In this manner, we obtain an exact grid representation of the field equations that we map transparently into circuit stamps. In these stamps, which are associated only with primal edges, lumped elements are directly taken from the entries of the material matrices. These entries are properly integrated constitutive parameters. In Figure 1b, we summarise our approach. In addition to previous work on ET problems [45], we also present a more elaborated description and implementation of boundary conditions and excitations. We also point out that the methodology pre- The left branch of the diagram shows the standard solution approach by using a field solver. The right branch illustrates the approach described herein to generate a netlist that is then fed to a circuit simulator.
sented herein allows for straightforward field-circuit coupling and is especially useful for an accurate representation of small devices in a larger circuit. In order to simulate wave propagation problems on a finite computational domain by means of circuit simulation, we present absorbing boundary conditions (ABCs) [46,47]. The outline of this paper is as follows. In Section 3, we provide the required basics of the FIT. The fundamentals of the MNA are summarised in Section 4. Then, the main part of the paper starts with the circuit representation of ET field problems in Section 5. How to extract circuit stamps for EM field problems is presented in Section 6. Finally, we show numerical examples in Section 7 and conclude the paper in Section 8.

Continuous Thermal and Electromagnetic Formulations
Let us consider a domain D with boundary ∂D and characterised by the constitutive parameters {ε, ν, σ}, where ε is the electric permittivity, ν is the magnetic reluctivity and σ is the electric conductivity. For every facet A, every volume V and with impressed electric sources J i , the EM field {E, H} in D is given by Maxwell's equations We call henceforth (1) the E-H formulation for conciseness. Above, D and B are the electric and magnetic flux density, respectively, J c is the electric conduction current and is the electric charge density. To guarantee the uniqueness of the solution, Maxwell's equations are supplemented with the constitutive relations, viz.
with suitable initial and boundary conditions (BCs) on ∂D.
We also deem it convenient to obtain Maxwell's equations involving the auxiliary magnetic vector potential A. To this end, we recall that E = −∇ϕ − ∂A/∂t and ∇ · (σ g A) = f , with an auxiliary scalar potential ϕ, a gauging material parameter σ g and an arbitrary scalar gauging function f . Substitution of these definitions in (1) and applying Stokes' theorem yields which we refer to as the E-A formulation. Whenever conducting materials are involved, electric currents result in Joule losses Q J = σ (∇ϕ) 2 that enter as source term into the heat equation which is given by where ρc is the volumetric heat capacity, T is the temperature, λ is the thermal conductivity and q i represent any impressed thermal flux density. In general, when thermal effects are considered, all constitutive parameters are also a function of the temperature.
To perform circuit extraction, we shall consider the above ET and EM formulations separately as described in Sections 5 and 6, respectively.

Discretising the Thermal and Electromagnetic Formulations
We discretise the domain D into a pair of orthogonal grids given by the primal grid G and its dual G. The grid G consists of primal points P i , i = 1, . . . , N P , primal edges (lines) L n , n = 1, . . . , N E , primal facets (areas) A n , n = 1, . . . , N F , and primal volumes V i , i = 1, . . . , N V . Similarly, the grid G consists of dual points P i , i = 1, . . . , N P , dual edges L n , n = 1, . . . , N E , dual facets A n , n = 1, . . . , N F , and dual volumes V i , i = 1, . . . , N V . The grids G and G are dual to each other in the sense that a primal edge L n intersects a dual facet A n and a primal point P i is located inside a dual volume V i and vice versa. For a regular hexahedral grid, the grid staggering is depicted in Figure 2. Due to this duality, the number of primal and dual grid objects fulfils As mentioned above, we write L n for the n-th edge of the primal grid and we emphasise the duality of edges and facets by using the same index. Thus, A n is the dual facet corresponding to the primal edge L n and L n is the dual edge corresponding to the primal facet A n . Let us further introduce a short (index) notation for geometric objects. If L n or A n are used as an index, we simply write n instead. Whether n refers to an edge or a facet should become clear from the context. For the dual objects, we useñ instead of L n or A n . The notation for points and volumes and their duals is done accordingly. In Table 1, we summarise this notation.

Description
Normal Index Table 1: Normal and index notation format for different entities of the grid.
The grid counterparts of the field quantities are allocated to points, edges, facets or volumes and collected in column vectors. Typical examples from electromagnetics (EM) are the discrete electric potentials ϕ, fields e, currents j and charges q that are allocated to primal points, edges, facets and volumes, respectively. The number of bows indicates the dimension of the corresponding geometric object. Nevertheless, we typically write q instead of q for conciseness. Defining all other grid quantities accordingly, their allocation used in this paper is shown in Figure 2. Furthermore, if we want to indicate a grid quantity, e.g. an electric field, to be allocated to an edge L n that is part of the boundary of a facet A k (L n ∩ ∂A k = L n ), we use the indexed notation e k;n . To relate grid quantities on either the primal or dual mesh, topological matrices, corresponding to the continuous topological operators, have to be defined. The incidence between grid points and edges is given by the discrete gradient matrix G. For an oriented edge L n , G ni = −1 if P i is the starting point of the edge, G ni = 1 if P i is the ending point of the edge and G ni = 0 if P i is neither starting nor ending point of the edge (P i ∩ ∂L n = ∅). To express the incidence between grid edges and facets, we use the discrete curl matrix C. Given a primal facet A k and its oriented boundary ∂A k , C kn = −1 if edge L n is oriented in opposite direction than ∂A k , C kn = 1 if edge L n is oriented in the same way as ∂A k and C kn = 0 if L n does not touch A k (L n ∩ ∂A k = ∅). Finally, the discrete divergence matrix S denotes the incidence between grid facets and volumes. For a primal volume V i and its oriented boundary ∂V k , the entries of S are defined in analogy to those of G and C. Additionally, we have the dual gradient, curl and divergence matrices G, C and S, respectively, defined accordingly. Useful relations between the topological matrices on the primal and dual grids are given by G = − S , G = −S and C = C [48].
To relate quantities on the primal grid to quantities on the dual grid and vice versa, material relations are employed. For the problem formulated in Section 2, the following three different kind of relations can be identified: • Quantities allocated to primal facets must be related to quantities allocated to dual edges.
• Quantities allocated to primal edges must be related to quantities allocated to dual facets.
• Quantities allocated to primal points must be related to quantities allocated to dual volumes.
As representatives for the above listed constitutive relations, we formulate three discrete material laws as where M ν is the magnetic reluctance matrix mapping the discrete magnetic flux b allocated to primal facets to the discrete magnetic field h allocated to dual edges, M ε is the electric capacitance matrix mapping the discrete electric field e allocated to primal edges to the discrete electric flux density d allocated to dual facets, and M ρc is the thermal capacitance matrix mapping the time derivative of the grid temperatureṪ allocated to primal points to the discrete heat power Q allocated to dual volumes. For other constitutive parameters, the material matrices are defined following these three cases.
Having established the grid constructs G and G and the corresponding topological and material matrices, the E-H formulation (1) of Maxwell's equations upon such a grid pair can be written as [5] − where M σ is the electric conductance matrix, b is the discrete magnetic flux density, j i is the discrete impressed electric current density and q is the discrete electric charge. Similarly, with the discrete vector potential a, the discrete E-A formulation is given by together with the discrete form of the gauging where F is the discrete counterpart of the gauging function f and M G is a gauging matrix. Since M G maps from quantities on primal edges to dual facets, it shares properties with the material matrices (e.g. M ε ) and thus can be interpreted as a material matrix with a material value equal to σ g . Casting also the heat equation (3) into a spatially discrete form, we obtain where M λ is the thermal conductance matrix and Q J is the discrete vector of the Joule losses. For details on the computation of Q J , we refer the reader to the work by Casper et al [45]. The impressed thermal fluxes are given by their discrete representative q i .

Finite Integration Technique and its Relation to Other Discretisation Schemes
For the circuit extraction from 3D field models as presented in this paper, we require diagonal, symmetric and positive definite material matrices. To fulfil this requirement, we choose to use the FIT and assume the material to coincide with the primal grid cells, see Figure 2. While the FIT has also been formulated for anisotropic materials [49], diagonal material matrices are obtained only in the isotropic case or in the case when the principal axes of anisotropy coincide with the coordinate axes. Then, the entries of the different material matrices are given by where εñ n , σñ n , λñ n , νñ n , and ρcĩ i are obtained by a suitable averaging scheme [48,50] and | · | represents the measure (i.e. area, length or volume) of the corresponding geometrical object. When using the FIT as the discretisation scheme with a canonical numbering of the grid nodes, the curl and material matrices exhibit the block structure where P ξ ∈ {0, ±1} N E ×N P , ξ ∈ {x, y, z} are the grid differential operators for the different coordinate directions. Here, M ε was used as an example whereas an equivalent block structure applies for all other material matrices. Within the theory of the FIT, the discrete quantities ϕ, e, j , q that were introduced in Section 3 are defined by means of integration with respect to their corresponding geometrical object, such that (ϕ) i := ϕ(P i ), e n := where an analogous definition applies for a, h, d and b. Due to the applied integration, one speaks of grid voltages instead of discrete fields, of grid currents (fluxes) instead of discrete current (flux) densities and of grid charges instead of discrete charge densities.
As an alternative to the FIT, equivalent approaches such as the cell method [51] or the FEM can be used as long as the orthogonality and one-to-one relation between primal and dual grid objects is guaranteed. When using FEM, linear basis functions together with an appropriate mass lumping for the material matrices must be used to obtain an equivalent scheme [52].

A Primer on Circuit Theory and the Modified Nodal Analysis
In this section, we briefly review some fundamentals about circuit theory and the MNA [15,53]. Kirchhoff's current and voltage laws are derived and form the basics for circuit analysis.
Any circuit can be understood as a directed graph consisting of interconnected nodes and branches. Let v i be one of the N n nodal potentials in a circuit and b n one of N b directed branches. With the incidence matrix A ∈ {−1, 0, 1} Nn×N b linking nodes and branches, the voltage-potential relation is given by The entries of A are defined such that a in = +1 if the branch b n is directed away from node n i and a in = −1 if b n is directed towards n i . If n i is neither starting nor ending point of b n , then a in = 0. With this definition, the exemplary voltage V n on the branch b n directed from n i to n j is given by For time invariant geometries, the current continuity equation reads for an arbitrary volume V . By considering a volume V i around an arbitrary circuit node n i and assuming that capacitive charges are located either fully inside or outside of V i , the total charge and also the charge's change rate in V i is zero. Therefore, (8) becomes If ∂ V i is composed by a finite number s of conductors with cross-sectional areas A n , Kirchhoff's current law (KCL) is obtained as where I n is the total current through the facet A n . This relation is also depicted in Figure 3a. To express (9) for all nodes in the circuit (cf. Figure 3b), the incidence matrix A can be used such that where I ∈ R N b is a vector of all currents allocated to the branches and 0 is a vector of zeros of suitable dimension. In a circuit, the basic branch elements are conductors, capacitors, inductors as well as voltage and current sources. As these elements are allocated to branches, sets of N G , N C , N L , N V and N I branches are defined. Hence, A can be arranged into a block matrix with sub-blocks for these elements [53], viz.
A circuit example consisting of a current source, two resistors and two capacitors together with the corresponding block structure of A is shown in Figure 3c. A similar subdivision is done for the current and voltage vectors, With these definitions, the voltages are given by and (10) becomes The relation between voltages and currents for the different branches is established by constitutive diagonal matrices that contain the element-wise material parameters. These are the conductance, capacitance, and inductance matrices G, C and L, respectively.
Expressing the corresponding source branches by means of the source voltages V s and source currents I s , this relation becomes Combining (11), (12) and (13), we obtain the MNA formulation Note that in contrast to the standard MNA theory, (14) still requires regularisation, typically done by the introduction of a reference (ground ) node.

Circuit Representation of Electrothermal Field Problems
In this section, we derive the circuit representation for transient ET field problems. We apply the electroquasistatic (EQS) approximation [54] to Maxwell's equations given by (1) and consider the coupling with the transient heat equation given by (3). Then, the bi-directionally coupled system in differential form reads with suitable initial and boundary conditions. The coupling is manifested by the Joule heating given by Q J = σ(∇ϕ) 2 in one way and by the temperature dependent conductivity σ(T ) in the opposite way. Due to the EQS approximation, this formulation does not account for inductive effects but does consider resistive and capacitive effects. For simplicity, we neglect the temperature dependency of the permittivity ε and of the volumetric heat capacity ρc. Applying FIT upon the ET system of (15), the semi-discrete formulation reads with initial and boundary conditions yet to be applied. Note that the system (16a) requires a regularisation which is typically done by choosing a reference (ground ) node ϕ gnd = 0, where ϕ gnd is one of the entries of ϕ.
To generate the netlist for formulation (16), the electric and thermal sub-problems are considered separately. The connection is subsequently established by the Joule losses and the temperature dependent electric conductivity. Next, in Section 5.1 and Section 5.2, the netlist generation for the EQS case and the thermal case is presented, respectively. Temperature dependent materials are discussed in Section 5.3 and finally, the implementation of initial and boundary conditions is described in Section 5.4. This allows us to formulate an algorithm for the ET netlist generation as presented in Section 5.5.

Electroquasistatic Circuit Representation
Let us now concentrate on the EQS sub-problem given by (15a). Since inductances are neglected in the EQS case, (14) simplifies to Thus, by inspection of (16a) and (17), we are led to the following equivalences: • The incidence matrices A C and A G coincide with the FIT divergence matrix S.
• The capacitance matrix C coincides with the FIT capacitance matrix M ε .
• The incidence matrices A V and A I coincide with the identity matrix. Figure 4: Equivalent electric circuit stamp for the exemplary edge L m between points P i and P j .
• The conductance matrix G coincides with the FIT conductance matrix M σ .
• The nodal voltages v correspond to the FIT degrees of freedom ϕ.
• The source currents I s are given by the divergence of the impressed currents S j i .
• The source voltages V s correspond to the Dirichlet potentials ϕ Dir , which are related to the reference node ϕ gnd .
• These equivalences also prevail themselves in the physical units.
Summarised, the field-circuit relations for EQS read where I is the identity matrix of corresponding size and ϕ Dir represents the potentials on the Dirichlet boundary nodes. To find the circuit stamp of each edge in the grid upon which (16a) holds, we employ the equivalences of (18) from which the circuit topology is derived. From (18a), conductors and capacitors are placed along the branches of the circuit. According to (18b), the values of the conductors and capacitors are directly taken from the corresponding FIT material matrices. Current and voltage sources are connected between a circuit node and ground as indicated by (18c). Furthermore, (18d) shows that the circuit's nodal potentials are equal to the potentials at the grid points. According to (18e), the current sources in the circuit represent the divergence of the FIT impressed currents. Finally, if Dirichlet BCs are imposed, voltage sources in the circuit represent FIT Dirichlet potentials as given by (18f). We further discuss BCs in Section 5.4. To summarise, if we consider an exemplary grid edge L m , we obtain a representative EQS circuit stamp as shown in Figure 4. Note that the temperature dependence of the materials is neglected for now and will be discussed in Section 5.3

Thermal Circuit Representation
In this section, we describe the circuit representation of the sub-problem described by (15b). By comparing (16b) to (16a), we observe a slightly different equation structure. Thermal capacities are not subject to spatial differences and thus do not link to neighbouring nodes. Instead, a thermal capacitance influences the change rate of the absolute temperature of a node. Thus, thermal capacitances are placed on branches connecting the nodes to a reference node at zero temperature. This reference node is an additional non-physical node that is introduced to obtain a consistent circuit representation. In the literature, this approach is also referred to as the Cauer model representing a discretised image of the heat flow [55,56]. An equivalent approach is the Foster model, in which the capacitances are placed between the circuit nodes and the parameters are adjusted accordingly. In the Foster model, the heat propagation is instantaneous and does not account for the fact that an object requires some delay before changing its temperature. The MNA formulation of (14) must be extended by this additional reference ground node such that with where 1 is a column vector of ones of appropriate size and v gnd = T gnd = 0. With these definitions, the equivalences between the MNA formulation of (19) and the thermal formulation of (16b) are readily obtained as where T Dir represents the temperatures at Dirichlet boundary nodes. To derive the circuit stamp of each edge in the grid upon which (16b) holds, we employ the equivalences in (20) from which the circuit topology is derived. From (20a), we infer that capacitances, voltage and current sources are directly connected between a circuit node and the thermal ground. As in the EQS case, conductors connect two neighbouring nodes in the grid as seen from (20b). The values of the conductors and capacitors are taken directly from the matrices M λ and M ρc according to (20c). In this manner, the nodal potentials represent the sought temperatures as seen from (20d). The current source is the sum of Joule losses, in fact represented by CCCSs, and the impressed heat flux according to (20e). Finally, if Dirichlet BCs are given, these are modelled by voltage sources in the circuit  as stated in (20f). We further comment on BCs in Section 5.4. To summarise, if we consider an exemplary grid edge L m , we obtain a representative thermal circuit stamp as shown in Figure 5a. Note that the temperature dependence of the materials is neglected until now and will be discussed in Section 5.3.
To further highlight the equivalences between electric and thermal circuits, we would like to briefly comment on this. Many quantities in electrical circuits can find their equivalent in thermal circuits. For example, electric potentials are equivalent to temperatures while electric currents are equivalent to heat fluxes. In Table 2, some of these equivalences are summarised.

Temperature Dependent Materials
Most materials exhibit temperature dependent behaviour. For the kind of considered materials, the temperature mainly influences the electric and thermal conductivity of the involved materials. In this section, we describe an approach to account for a temperaturedependent electric conductivity when generating the corresponding SPICE netlist. Nevertheless, the presented approach can be applied accordingly to other material's temperature dependencies as well, e.g. the thermal capacity.
Since the SPICE elements that represent the electric conductivities are conductors placed in circuit branches (cf. Section 5.1), we first need to define the temperature of a branch b m . To this end, we take the average temperature T m of the nodes interconnected by the branch b m . Assuming that the temperature dependence of the electric conductivity is known, the temperature-dependent electric conductance of branch b m is given by We implement (21) by means of behavioural sources in the SPICE language.

Initial Conditions and Boundary Conditions
When a coupled problem of more than one transient differential equation is considered, each sub-problem requires its own initial conditions and BCs. Therefore, we impose these conditions on the EQS and thermal sub-problems separately. However, certain equivalences allow to follow the same procedure for both sub-problems. For any kind of transient problem, initial conditions are required. As every SPICE dialect supports specifying initial conditions, these can be directly imposed by the corresponding syntax in the netlist. For electric problems, different types of BCs are of interest. In low-frequency problems as the EQS case, Dirichlet or Neumann conditions are typically used. Dirichlet BCs conditions correspond to a fixed potential enforced at the boundary, while Neumann conditions prescribe the electric current through the boundary. Similarly, in thermal problems, Dirichlet BCs correspond to a prescribed temperature at the boundary, while Neumann BCs prescribe thermal fluxes through the boundary. Additionally, thermal problems commonly also involve Robin BCs that describe convective and radiative boundaries.
Dirichlet BCs are represented in the circuit by voltage sources between the ground node and the Dirichlet nodes. Homogeneous Neumann conditions are automatically fulfilled since no edge or branch leaves the domain. For simplicity, we do not consider inhomogeneous Neumann conditions in this paper. A Robin BC can be understood as a conduction between a boundary node and an external node n ∞ representing the fixed ambient temperature T ∞ . Therefore, Robin BCs are represented in the circuit by conductors connected between the boundary nodes and n ∞ as shown in Figure 5b. We collect the relevant BCs in Table 3.

Boundary Condition Implementation
Dirichlet lumped voltage sources hom. Neumann no edges leaving the circuit Robin additional non-physical ground write Cthi niT gnd M ρc;ĩi ic = T 0 8: write BLossi gnd niT I=Q J,i (t) 9: if P i is electric Dirichlet node then 10: write VDirEli ni gnd V Dir,i (t) 11: if P i is thermal Dirichlet node then 13: write VDirThi niT gnd T Dir,i (t) 14: end if 15: if an impressed current flows out of P i then 16: write IimpEli ni gnd ( S j i (t)) i

17:
end if 18: if an impressed heat flux flows out of P i then 19: write IimpThi niT gnd ( S q i (t)) i 20: end if 21: end for

Electrothermal Netlist Generation
To finalise this section, we formulate Algorithm 1 to automatically generate ET SPICE netlists. For every grid edge L m that connects grid points P i and P j , where i < j, the thermal conductance, electric capacitance and the temperature dependent electric conductance (cf. Section 5.3) are written to the netlist connecting nodes n i and n j of the circuit (lines 2-4). Due to the possible non-linearity of the electric conductivity (cf. Section 5.3), a behavioural source is used, where V ij is the voltage between node n i and n j . Initial conditions (ic) for the electric part can be included using the corresponding syntax for the capacitors. Here, we use exemplary zero initial conditions. Additionally, for every grid point P i , the thermal capacitance and the current controlled current source (CCCS) representing the Joule losses are added to the netlist connecting node n i and the ground node (gnd) of the circuit (lines 5 and 6) 1 . Initial conditions (ic) for the thermal part are specified by pre-charging the thermal capacitors with the initial temperature T 0 . To specify a CCCS in the SPICE language, a behavioural source is used. If P i is specified as an electric (thermal) Dirichlet node, an additional voltage source connecting node n i and the ground node of the circuit is inserted (lines 9-14). Furthermore, if an impressed current (heat flux) flows out of P i , an additional current source is added to the netlist connecting node n i and the ground node of the circuit (lines 15-20).

Circuit Representation of Electromagnetic Field Problems
In this section, we neglect thermal effects and describe the circuit representation of general 3D EM field problems as given by (1)- (2). However, the thermo-EM coupling can be established analogously. First, in Section 6.1, we introduce an auxiliary set notion to collect specific edges and facets of the grid. In Section 6.2, the E-H formulation (4) and the E-A formulation (5)-(6) of the Maxwell grid equations (MGEs) are transparently mapped into an electric circuit that fully describes the problem at hand. Finally, in Section 6.3, we extend our analysis in order to realise ABCs as circuit stamps. These are typically needed to limit the computational domain while minimising unphysical reflections caused by the domain truncation.

Auxiliary Sets of Edges and Facets
In the following sections, we will take sums over specific edges or facets of the grid.

Circuit Representation of the Maxwell Grid Equations
In this section, circuit representations of the MGEs given by (4)-(6) are derived. First, we start with the E-H formulation of (4) to find a corresponding circuit description which is presented in Section 6.2.1. Subsequently, in Section 6.2.2 a circuit description based on the E-A formulation of (5)- (6) is presented. For this formulation, a tree-cotree decomposition is used.

Circuit Representation Based on the E-H Formulation
A first electric circuit representing the MGEs is realised from the E-H formulation of (4). The goal is to find an expression for the electric voltage on one primal edge such that a circuit stamp for each edge is obtained. Let L m be the edge of interest with its corresponding dual facet A m . We collect all relevant quantities on the grid objects in the neighbourhood of L m by using the notation introduced in Section 6.1. To find the voltage e m on L m , let us consider the pair of interlocked facets A k and A m as illustrated in Figure 8a. For this part of the grid, it suffices to consider only the k-th row of (4a) and them-th row of (4b) giving Thanks to the properties C = C , Cmk ∈ {−1, 1} and C km ∈ {−1, 1}, we have that Cmk = C km and thus their product equals to unity. However, the product CmkC kn can be either −1 or 1 as illustrated by Figure 8a. Thus, we write (24) Equation (26) is the Kirchhoff's current law (KCL) associated with the primal edge L m with e m representing the voltage drop along the edge. In fact, according to the definitions of the field quantities in (26), it is easy to realise that • e m has unit of voltage (V), • j i,m has unit of current (A), • M ε;mm is positive and has unit of capacitance (F), • M σ;mm is positive and has unit of conductance (S), • M Σ ν;mm is positive and has unit of reluctance (H -1 ), • j c;mkn has unit of current (A).
Thus, by using voltage controlled current sources (VCCSs) to model j c;mkn accounting for the contributions from neighbouring edges, we can directly represent (26) with the circuit stamp depicted in Figure 7, which preserves the voltage drop between the terminals of L m . There are as many of these stamps as primal edges in G, and all of them interact via VCCSs. The concatenation of these elementary stamps constitutes the electric circuit representing the electromagnetic problem at hand. In Algorithm 2, the steps to generate the netlist representing the EM circuit for the entire grid are listed. A circuit stamp, as shown in Figure 7, needs to be created for every edge L m in the grid which is realised by a loop in the code. In each iteration, a resistor, inductor and capacitor with the values taken from the material matrices is added (lines 2-4). If an impressed current source shall be placed on L m , an independent current source with a predefined value j i;m is used (lines 5-7). Finally, an inner double loop is required to insert the controlled current sources that model the influence of the edges in the neighbourhood. For this purpose, we choose to insert CCCSs 2 being controlled by end for 13: end for the current M Σ ν;ñn e k;n dt with a gain of g Imk n := CmkM ν;kk C kn (M Σ ν;ñn ) −1 (lines 8-12). By abuse of notation, we denote the controlling device by e k;n .
We remark that a similar analysis on (4) in which magnetic conductivities and sources are considered instead of electric ones can be done 3 . This approach would lead to a circuit stamp that is dual to the one in Figure 7. Namely, in this dual stamp, the Kirchhoff's voltage law (KVL) is guaranteed for each dual edge and hk represents the electric current in the circuit. Furthermore, the resulting lumped elements, stemming from the material matrices, are placed in series with the discrete impressed magnetic current m i;k that plays the role of an independent voltage source exciting the circuit. The interaction between the dual stamps is mediated via current controlled voltage sources (CCVSs). In the general case, when both electric and magnetic sources are present, the electric circuit will consist of both the primal and dual stamps, which interact via CCCSs and VCVSs, accordingly.

Circuit Representation Based on the E-A Formulation
In this section, we describe a circuit representation based on the E-A formulation of (5)- (6). To guarantee uniqueness of the solution, the magnetic potential a is gauged by means of a tree-cotree decomposition [59]. Therefore, although the inferred circuit stamps are gauge dependent, the solution obtained for e is unique. We start by considering one primal facet A k ∈ A m (cf. Figure 8a). For this facet, it suffices to consider the k-th row of (5a) and them-th row of (5b) such that the E-A formulation of the MGEs for a generic edge L m reads As in Section 6.2.1, we aim at finding a unique circuit representation of edge L m . We start by observing that the system matrix CM ν C of (5) is singular 4 . This singularity manifests itself in the non-uniqueness of a. In fact, only the curl of a is uniquely defined. As a remedy, we must explicitly impose the gauging (6) upon (5). To this end, let us assume that we have constructed a suitable tree G t and a cotree G c out of the primal grid G, as exemplified in Figure 9. Then, we symbolically introduce the orthogonal permutation matrix P G to partition a into its tree a t ∈ G t and cotree a c ∈ G c components with N t and N c entries, respectively. Symbolically, this reads Figure 9: Example of a spanning tree G t (solid) and cotree G c (dashed) in a generic graph G. The case studied here consists of an orthogonal pair of primal and dual grids.
Similarly, we also divide M G and S into their tree and cotree components by applying the P G matrix accordingly, viz.
The identities of (28) and (29) together with the orthogonality property of P G enable us to rewrite the gauging of (6) as Additionally, since S 1t M Gt is a square and invertible matrix 5 , we may finally express the tree component a t as The column vector F, which represents the grid counterpart of the scalar function f quantifying the divergence of A, is of free choice. Therefore, we conveniently choose F = 0 impressing the Coulomb gauge [60] to straightforwardly arrive at The matrix E tc ∈ R Nt×Nc is known as the essential incidence matrix [59] and establishes a direct relation between the tree and cotree components a t and a c , respectively 6 . By 5 These two properties come from two facts: the squareness is a consequence of removing the row associated with the ground node required in circuit analysis. The invertibility is due to removing the non-null kernel space of the matrix CMν C by means of the tree-cotree decomposition. 6 Owing to this property, we identify in the rows of Etc ∈ R N t ×Nc the collection of all fundamental cut-sets associated with the selected tree and cotree. We recall that a fundamental cut-set is a set formed by the union of a single tree edge and the unique set of adjoining cotree edges. In this manner, the fundamental cut-sets are used to express Kirchhoff's current law in the general form of (31).
inspection and expansion of (31), we express the m-th component of the column vector a t as a t;m = − n∈Gc E tc;mn a c;n , with E tc;mn := i∈G M −1 Gt;mm S −1 1t;mĩ S 1c;ĩñ M Gc;ñn and i spanning over the nodes in G. In this manner, (32) removes the redundancy associated with a and guarantees a unique solution of (27).
Let us now get back to the generic primal edge L m once again. As introduced in Section 6.1, we use specific sets to refer to edges in the tree and cotree denoted by L m k;t , L m,0 k;t and L m k;c , L m,0 k;c , respectively. Isolating the voltage e m along edge L m embedded in facet A k from (27a), we arrive at We may then take the sum over all facets A k ∈ A m (cf. Figure 8b). Then, by splitting { a k;n } into tree and cotree components, we arrive at We observe that (34) and (36) can be interpreted as the KVL and KCL of an arbitrary primal edge L m provided that e m and a m represent the sought voltage and current, respectively. As a matter of fact, we observe that • e m has unit of V and represents the voltage drop along L m that is either in the tree or cotree set.
• a m has unit of Weber (Wb) and represents the current along L m . If L m ∈ G c , this current is a degree of freedom. Otherwise, it is modelled by a CCCS as stated in (32). For compact notation, we label this CCCS as I Σ;m := n∈Gc I Σ;mn , where I Σ;mn := −E tc;mn a c;n .
• M Σ ν;mm has unit of H -1 and is expected to be positive since M ν;kk > 0 for all materials. This term scales the displacement, conduction and impressed current as seen in (35).
With these observations, we may depict (33) and (35) by means of the circuit stamps of Figure 10a and Figure 10b for a cotree and tree edge, respectively. As we can see therein, a m is regarded as the electric current along the primal edge L m , while the voltage drop between its terminals is established by the VCVS V e;m and the CCVSs V c;m and V t;m that mediate the interaction with neighbouring edges. When the edge L m belongs to the tree G t , then the current a m is modelled by the CCCS I Σ;m as demanded by (32). Otherwise, for cotree edges, a m is a degree of freedom. The current a m is then split into several branches where a resistance, capacitance, the impressed current source, and CCCSs are connected. The concatenation of these fundamental stamps forms the electric circuit representing the magnetic vector potential formulation of the electromagnetic problem.
In Algorithm 3, we show the pseudocode to generate the netlist of the circuit stamps depicted in Figures 10a and 10b. The iteration over all edges L m in the grid is done in the outermost loop. In each iteration, one resistor and capacitor connecting the node n m to ground (gnd) must be added to the netlist (lines 2 and 3). If an impressed current is present at L m , an independent current source with value I i;m is added between n m and gnd (lines 4-6). For every tree edge, the currents I Σ;mn from the cotree edges in the neighbourhood are modelled by a parallel connection of CCCSs between n e;m and n m (lines 7-11). Using a double loop, the inductive currents I t;mkn and I c;mkn from the tree and cotree branches, respectively, are added as a parallel connection of CCCSs between n m and gnd (lines [12][13][14][15][16][17][18][19]. The controlling current is given by the current through device FIcn11 and the gain is g I mkn = CmkM ν;kk C kn (M Σ ν;mm ) −1 . Finally, the voltages V e;mkn , V c;mkn and V t;mkn are added as a series connection between n e;m and gnd with intermediate nodes indexed by k, n andn = n + 1 (lines [20][21][22][23][24][25][26][27][28][29][30]. For the voltage V e;mkn , the voltage between node nn and gnd controls a VCVS with a gain of g V mkn = N −1 F;m C kn C −1 km On the other hand, the voltages V c;mkn and V t;mkn are added as behavioural sources using DDT as the SPICE syntax for time derivatives and a gain of g I mkn = N −1 F;m C kn C −1 km . Note that for a cotree edge, the node n e;m coincides with the node n m . We remark that a similar analysis on (1) by considering only magnetic conductivities and sources instead of electric ones is also possible. Thus, with the auxiliary electric potential D = ∇ × F and ∇ · F = a, where a is an arbitrary gauging function, circuit stamps which are dual to those shown in Figure 10a and Figure 10b can be found. In these dual stamps, h m represents an electric current while f m , namely the grid counterpart of F, would be regarded as a voltage drop. Finally, if both electric and magnetic sources were present, then the entire circuit would consist of the aggregate of interacting primal and dual stamps.

Absorbing Boundary Conditions
In many electromagnetic field simulation set-ups, the computational domain must be bounded. To simulate free wave propagation, one must impose according conditions for the fields at the boundaries of the domain. These conditions are known as ABCs [46] and aim at minimising (ideally cancelling) unphysical incoming reflections. At the boundary of the domain, a distinction is made between normal and tangential components and between longitudinal and transverse derivatives 7 .
For starters, let ∂ r , ∇ t and ∂ t denote the longitudinal, transversal and temporal derivative operators, respectively. Furthermore, within the context of this analysis, we define ∇ := ∂ 2 r + ∇ 2 t 1/2 . We start the realisation of circuit stamps for ABCs by considering 7 In this regard, the longitudinal (transverse) derivative coincides with the normal (tangent) derivative at the boundary.
Algorithm 3 Electromagnetic SPICE netlist generation based on the E-A formulation. if L m ∈ G t then 8: for edge L n ∈ G c do 9: write FIsummn nem nm FIcn11 −E tc;mn for edge L k ∈ L m m do 13: for edge L n ∈ L m,0 k;c do 14: write FIcmkn gnd nm FIcn11 g I write BVcmkn ncmkn ncmkn . . .

30:
V=g I mkn DDT( a k;n ) 31: end for 32: end for 33: end for the time-domain wave equation for the electric field E in a homogeneous and isotropic medium 8 , viz.
Equation (37) can be expanded in terms of the so-called travelling wave operators as Above, we observe that in an arbitrary point in space, the field E can in general be considered as the superposition of an inward and outward travelling wave E + and E − , respectively, viz.
which upon substitution in (38) straightforwardly leads to the following set of travelling wave equations, inasmuch as both wave components E ± are linearly independent. Owing to their definition, the wave components E ± satisfy independently and simultaneously the following 9 , which tells us explicitly that the outward (inward) travelling wave operator cancels out the inward (outward) travelling wave at any point of interest. Now, let us get back to (38) to further expand the operators therein to obtain  As we can see, the above wave operators entail the calculation of the inverse of 1 + ∇ 2 t /∂ 2 r , which generally translates into a global integral operator [61]. In principle, the expansion of this integral operator around the observation point yields the exact explicit representation of the travelling wave operators in (40). However, this approach is contrary to the idea of realising simple and efficient ABCs for circuit simulations.
As a remedy, we may expand 1/ 1 + ∇ 2 t /∂ 2 r as in a Taylor series to yield The above equation holds for any component of E. Furthermore, if the field E propagates in free space, we may assume that the operation ∇ 2 t is negligible in the neighbourhood of an observation point within the spherical wavefront. Hence, we may rewrite (41) in terms of simplified inward and outward travelling wave propagators as since the general expression of E in (42) admits the superposition of inward and outward travelling waves E + and E − . Then, to minimise unwanted incoming reflections at a certain boundary given by r = r B , we may impose on E the condition in agreement with (39). The condition given by (43) is a first-order ABC known as Engquist-Madja condition [46]. It is a local condition because it is evaluated pointwise taking only into account the wave component propagating perpendicularly to the boundary. Thereby, the condition in (43) states that a practical transparent boundary at r = r B can be realised if it is guaranteed that the phase-amplitude of the field on the boundary at a certain time t B is equal to that one the field had at some previous instant t B − ∆t at a point r = r B − ∆t/ √ µε. The practical relevance of (43) comes from its simplicity.
Let us now interpret (43) within the context of the Maxwell grid equations in order to realise circuit stamps associated with the ABCs of (43). To this end, let us consider the grid wave equation for the electric grid voltage e, which for time-independent constitutive parameters and non-conducting source-free regions can be obtained from the grid curl equations of (4), viz.
With the structure of the matrices C, M ε and M ν given by (7) and the subdivision e = ( e x , e y , e z ) , the expression P z M ν;y P z + P y M ν;z P y e x − P y M ν;z P x e y − P z M ν;y P x e z = −M ε;x d 2 dt 2 e x (44) for the component e x is obtained. Similar expressions can be also obtained for the other two components e y and e z .
The grid counterpart of the ABC in (43) is obtained from (44) by extracting the grid wave equation associated with a generic primal edge L x;m ∈ G oriented along the x-direction. This yields k n P z;mk M ν;y;kk P z;kn e x;n + k n P y;mk M ν;z;kk P y;kn e x;n − k n P y;mk M ν;z;kk P x;kn e y;n − k n P z;mk M ν;y;kk P x;kn e z;n = −M ε;x;mm d 2 dt 2 e x;m , with n andk spanning over the corresponding primal and dual edges oriented along the indicated directions. We make the following observations upon the above grid wave equation.
• Two main grid wave components contribute to the time variation of e x;m along L x;m . The first of these grid wave components propagates along the z-direction and stems from the spatial variation of e x along this direction as stated by the grid derivatives in k n P z;mk M ν;y;kk P z;kn e x;n . The second one propagates along the y-direction and stems from the spatial variation of e x along this direction as stated by k n P y;mk M ν;z;kk P y;kn e x;n .
• Two secondary grid wave components contribute to the time variation of e x;m along L x;m . The first one propagates along the y-direction and stems from the spatial variation of e y along the x-direction as stated by k n P y;mk M ν;z;kk P x;kn e y;n . The second one propagates along the z-direction and stems from the spatial variation of e z along the x-direction as stated by k n P z;mk M ν;y;kk P x;kn e z;n .
• Owing to both the structure of (45) and the grid-like embedding where the propagation takes place, we may apply superposition to treat each propagation direction separately. Therefore, each grid wave component satisfies its own one-dimensional grid wave equation.
We construct the grid counterpart of (43) by adhering to the same principle that led to it. Namely, we take only into account grid wave propagation perpendicular to the boundary of interest. Thereby, let us assume that our boundary is located for example at z = z B . Thus, by invoking the superposition principle, we may write the grid wave equation associated with the component e (z) x;m that propagates along the z-direction as k n P z;mk M ν;y;kk P z;kn e (z) x;n − k n P z;mk M ν;y;kk P x;kn e (x) z;n = −M ε;x;mm where we have used the labels (z) and (x) to explicitly indicate directions of propagation. By considering that for sufficiently well-refined grids G and G we may expect M ν;y;kk to not vary significantly for the two relevant facets A y;k (cf. Figure 11), we may define an Figure 11: The ABC is imposed at edge L x;m in the grid G. The impinging propagating grid wave coming from edge L x;m−Mz is then absorbed. The two facets A y;k−Mz and A y;k that contribute to the averaging of M ν;y;kk are shown.
average value M ν;y;kk that enables us to write the leftmost term on the left-hand side of (46) as k n P z;mk M ν;y;kk P z;kn e (z) x;n = M ν;y;kk k n P z;mk P z;kn e (z) x;n , with a similar result for the rightmost term on the left-hand side of (46). Thus, together with the property P ξ;nk = − P ξ;nk , we may proceed to write (46) as k n P z;mk P z;kn e (z) Above, the term k n P z;mk P x;kn e (x) z;n is the grid counterpart of the continuous operator in the wave equation of (37). Analogously, we may assume that in the vicinity of L x;m at the boundary z = z B , the wavefront of the impinging grid wave is plane. Therefore, we may neglect transverse variations and we arrive at We observe that k n P z;mk P z;kn e (z) x;m is formally the second-order grid derivative of the scalar field e (z) x;m along the z-direction. This directly leads us to the grid counterpart of (37) in the vicinity of L x;m , viz.
Again, this equation can be factorised in terms of grid wave propagators similar to (42) to yield the grid version of (43), viz. , and length |L z;m−Mz |, see Figure 12. The voltage that excites the line is given by the voltage on edge L x;m−Mz , namely e x;m−Mz . Therefore, if one wants to implement the ABC at edge L m , an impedance Z 0;m must be assigned to L m . A subsequent circuit extraction can be carried out by implementing either Algorithm 2 or 3.

Numerical Examples
In this section, the presented methodology to generate electric circuit stamps representing 3D field problems is applied to several representative numerical examples. In Section 7.1, we use our netlist extraction method as described in Section 5 on an ET problem. The considered ET problem is a 3D field problem corresponding to the series connection of a capacitor and a resistor. While applying an external voltage, the transient heating due to the resulting current is simulated using SPICE and then compared to a field solver reference solution. Additionally, a circuit representation for the ET field problem of a microelectronic chip package is obtained and used for circuit simulation. In Section 7.2, we apply our method of circuit extraction for EM field problems as described in Section 6 to compute the resonant frequencies of a rectangular cavity with perfect electric conducting (PEC) boundaries. This example is quite illustrative and easy to implement because of the required BCs on the cavity walls. It simply suffices not to print the circuit stamp associated with those edges on the wall, meaning that the associated stamps are short-circuited. Furthermore, the availability of an analytic formula for the resonant frequencies permits a direct error assessment. Finally, in Section 7.3, the implementation of ABCs as discussed in Section 6.3 is carried out to investigate reflections at the end of a rectangular coaxial waveguide. For all presented examples, the Matlab R code Figure 13: Geometry of the ET validation example. The series connection of a resistive part and a capacitive part is excited with a voltage source V app imposed as a Dirichlet condition.
to generate the corresponding netlists from the discretised 3D field problem is openly available [62].

Electrothermal Circuit Validation
To validate our netlist extraction method on an ET problem, we consider the Joule heating in a 3D field problem represented by a series connection of an electric resistor and a capacitor. The temperature dependence of the electric conductivity is manifested via the temperature coefficient α = 3.9 × 10 −3 1/K. The relevant configuration is realised by a brick of two different materials as shown in Figure 13.  Table 4, all material properties are summarised for a reference temperature of T 0 = 293 K. A spatial grid with 9 × 9 × 9 cells is employed and the field problem is solved by using an in-house implementation of the FIT method with a first-order implicit Euler scheme as time integrator. The simulation time amounts to t end = 13 µs. For the simulation of the extracted electric circuit, we use the freely available LTspice software 10 . LTspice uses adaptive refinement in time for which an initial time step of ∆t init = 0.13 µs is used. The resulting non-equidistant time axis is refined by a factor of three and then used for the FIT solver. A voltage V app = 1 kV(1 − exp(−t/τ )) with τ = 0.1t end is applied at the electrodes as shown in Figure 13. Using this setting, two simulations are run. The first neglects the temperature dependence of the conductivities and thus a linear setting ensues. The second neglects only the temperature dependence of the thermal conductivity but accounts for that of the electric conductivity via the temperature coefficient α entailing a non-linear setting. To observe the transient behaviour, we select the resistor-capacitor interface point x 0 = ( , 0, 0) as observation point and plot the results in Figure 14.  For a quantitative error assessment of the solution, we define the measures where ϕ cir , ϕ FIT , T cir and T FIT are the potential and temperature solution vectors obtained via circuit and FIT simulation, respectively. To calculate these errors appropriately, the circuit solution is interpolated to the time axis employed by the FIT solution using cubic spline interpolation. We want to remark that the quantities in (48)  For an industry-relevant example, the proposed method is applied to the 3D microelectronic chip package [45] as shown in Figure 15a. The field problem is discretised as described in Section 3 and Algorithm 1 is used to generate the corresponding netlist. This netlist uses 101 147 circuit elements to describe a field problem that has been discretised using a grid with 9660 nodes. Running a transient analysis on this netlist, an error of ∆ ϕ ≈ 0.23 % and ∆ T ≈ 0.17 % compared to the field simulation is achieved. Figure 15b  shows the temperature of the hottest point in the chip package obtained by FIT and circuit simulation. The temperature distribution in the chip package resulting from circuit simulation is shown in Figure 15c. Thus, a good agreement of circuit simulation results for a 3D ET problem is achieved when compared to the corresponding field solver results.

Electromagnetic Circuit Validation
In this section, we validate the method presented in Section 6 for the circuit representation of EM field problems. To this end, a lossless rectangular resonant cavity with PEC boundaries and outer dimensions of a × b × d = 0.1 × 0.2 × 0.2 m 3 is simulated and its resonant frequencies are computed. The homogeneous material within the cavity is specified by the relative permittivity ε r = 2 and the relative permeability µ r = 1 for which the resonant frequencies can also be calculated by means of the formula [63] f mnp where c 0 is the speed of light and {m, n, p} are the indices of the resonant modes and are given by natural numbers including zero. For these resonant frequencies, the longitudinal transverse electric (TE) and transverse magnetic (TM) field components are given by are the corresponding field amplitudes. For a TE (TM) mode mnp to exist, H z (E z ) must not become zero.
To simulate the excitation of modes within the cavity, we discretise the interior of the cavity using a regular grid of 10 cells in each direction and apply PECs upon all cavity walls. Then, we obtain the resonant frequencies by solving the generalised eigenvalue problem given by where a resonant frequency is denoted by f FIT r;E . Alternatively, we can use appropriate excitations to analyse the resulting field at a set of given frequencies. For example, we can use an electric current source oriented along the positive z-direction and attached to the central grid point to excite TM modes. In a similar manner, TE modes are excited by means of a looping electric current source located in the cavity centre. We then solve the discretised problem given by for a set of angular frequencies ω, where j i is the current source vector whose entries are all zero except at the corresponding source edges. The frequency axis from 0.5 to 3 GHz is discretised using 2000 points for TM excitation and 3000 points for TE excitation. The results are evaluated on one edge for each excitation type. For the TM case, an edge in positive z-direction connected to the central grid point is used while for the TE case, an edge in positive x-direction connected to the point (5,6,10)cm is used. In Figure 16, the voltages e x (ω) and e z (ω) along these edges are plotted. From the peaks in the plots, the corresponding resonant frequencies f FIT r;TE and f FIT r;TM are identified 11 . To validate our circuit extraction method for EM problems, we generate the netlist of the resonant cavity according to Algorithm 2 and simulate the resulting circuit in LTspice by performing an AC analysis in the same frequency range as before. We then identify the circuit voltages corresponding to e x (ω) and e z (ω) and plot them also directly in Figure 16 for a comparison.  Table 5: Analytic and computed resonant frequencies in GHz for several resonant modes and the corresponding relative errors. As mode degeneracy in the cavity is relevant, an exact identification of the mode indices from the plots in Figure 16 is not possible.
the first few modes in

Signal Transmission Using Absorbing Boundary Conditions
In this section, based on the method described in Section 6.3, we present a simple validation example for ABCs in the context of circuit simulation. To this end, let us consider a coaxial transmission line of rectangular cross section oriented along the z-direction as depicted in Figure 17. For the simulation of an infinitely long line using a finite computational domain, the implementation of ABCs is required to counteract unwanted incoming reflections. We use an excitation signal at port 1 and simulate its propagation in time until it has reached port 2. Thus, we generate two simulation results in time domain. The first one corresponds to the case when port 2 is terminated with a perfect magnetic wall, that is an open port (Z 2 → ∞). The second one corresponds to the case when port 2 is terminated with the characteristic line impedance (Z 2 = Z 0 ). For both cases, a perfect magnetic conducting (PMC) (Z 1 → ∞) at port 1 is applied 12 . The length of the coaxial line is l = 150 cm, the width and height of the outer conductor are w o = h o = 3 cm and of the inner conductor w i = h i = 1 cm. While the conductors are modelled as PEC, the material between them is vacuum. Due to the expected propagation in z-direction, the longitudinal direction requires a finer discretisation compared to the transversal direction. Thus, we choose a grid of 3 × 3 × 150 cells. According to (47) for such a discretisation grid, the characteristic impedance for the edges connecting the inner and outer conductor along the plane of port 2 should amount to 13 Z 0;m = (M ν;y;kk M ε;x;mm ) −1/2 ≈ 376.7 Ω. For the given grid, there are eight such edges giving eight parallel conductances such that the total resistance at port 2 equals Z 2 = Z 0 = 8Z 0;m ≈ 47.09 Ω, which is also the characteristic impedance of the line. To excite the signal at port 1, the edges connecting the inner and outer conductors along the plane of port 1 are impressed with a current such that the total current from inner to outer conductor is , withÎ = 1 A, which is a Gauss pulse with a maximal frequency 14 component of f max = 1 GHz. The used constants are given by σ G = ln(10)/(πf max ) and by t 0 = 6σ 2 G ln(10). The simulation time t end = 10 ns is chosen such that the excited pulse can reach port 2. Having defined the geometry, the excitation and the simulation time, we also generate the corresponding netlist using Algorithm 2.
The generated netlist is simulated by means of a transient analysis in LTspice. On the other hand, the Leapfrog scheme is used as a time integrator within the FIT framework to carry out the simulation directly on the 3D grid. In both cases, we use the time axis generated by the adaptive time stepping algorithm provided by LTspice, which satisfies the Courant-Friedrichs-Levy (CFL)-condition being a stability requirement for the explicit Leapfrog scheme [1]. In the following, we compare the voltage V oi (z, t) between outer and inner conductor and the voltage V 2 (t) = V oi (l, t) at port 2. For Z 2 → ∞ and Z 2 = Z 0 , Figure 18 shows V oi (z, t) at different times computed by means of FIT and circuit simulation. Figure 18a shows the case in which port 2 is terminated by a perfect magnetic boundary (standard homogeneous Neumann) condition while Figure 18b shows the case when a matching impedance Z 2 = Z 0 according to (47) is applied at port 2. As predicted by the theory in Section 6.3, we observe that the matching impedance at 13 Note that in the calculation of Z0;m, the value employed for Mε;x;mm is taken directly from the parallel edge just in front of the boundary edge in accordance with the impinging grid wave front speed. 14 Confining the excitation to this maximal frequency component, we assure that the TEM mode is the only propagating mode on the line  Figure 19: Comparison of the output voltage V 2 with respect to time computed by FIT and circuit simulation for the case of (a) Z 2 → ∞ to realise total reflection and (b) Z 2 = Z 0 to realise ABCs. port 2 counteracts incoming reflections effectively. In Figure 19, we show V 2 (t) computed by means of FIT and circuit simulation for the two already considered cases. We observe therein that incoming reflections at port 2 result in an undesired overshooting of the voltage. For a quantitative comparison, we define the relative error of V 2 (t) between FIT and circuit results as and obtain ∆ Z 0 V 2 ≈ 1.046 % and ∆ ∞ V 2 ≈ 1.136 %.

Conclusion and Future Work
A method for the automatic netlist generation of general 3D ET and EM problems has been presented. The topology of each circuit stamp associated with edges in the regular primal grid has been derived by using FIT for spatial discretisation. Using the MNA, the FIT-discretised ET formulation has been mapped into a circuit that can be solved by any SPICE-like program. It has been shown that initial conditions can be easily prescribed as initial potentials for the lumped capacitances in the SPICE language. Furthermore, the implementation of mixed boundary conditions of Dirichlet, homogeneous Neumann and Robin type has been discussed. We have also shown that temperature dependent material models result in non-linearities in the lumped resistances requiring the implementation of behavioural VCCSs in SPICE. From the standard E-H formulation and the E-A formulation, we have derived circuit stamps representing general EM problems. In both circuit representations, the integrated electric field models the voltage between the stamp terminals while the integrated magnetic vector potential models the electric current in the E-A formulation. To guarantee uniqueness of the solution in the latter, we have employed Coulomb's gauge on the magnetic vector potential, that has been implemented by means of a tree-cotree decomposition of the primal discretisation grid. Thereby, the electric current along edges in the cotree are degrees of freedom whereas those along edges in the tree are modelled by CCCSs being controlled by currents in the cotree. For both representations, a dual circuit formulation exists if magnetic sources instead of electrical sources are considered. In the dual case, an auxiliary electric potential would be used instead of the magnetic vector potential. To demonstrate the correctness of our formulations, several numerical examples have been shown for the primal circuits involving electric sources only.
The formulation of inhomogeneous Neumann BCs could be a further extension to the presented approach. Furthermore, the method can also be applied to extract circuits from FEM models. To account for thermal effects in EM problems, the methods for the extraction of ET and EM circuit stamps can be combined to generate a thermo-EM circuit stamp. Methods to account for non-linear material characteristics in the EM case are still to be developed. However, in principle one can follow similar ideas to those presented in the ET case. For large field models, the resulting circuit can become very large. Therefore, to efficiently simulate such circuits, dedicated MOR techniques for circuits can be applied. The first of these techniques is known as the asymptotic waveform evaluation (AWE) proposed by Pillage and Rohrer [64] and extensions developed afterwards. The most prominent ones are the matrix Padé via a Lanczos-type process (MPVL) by Feldmann and Freund [65] and the passive reduced-order interconnect macromodeling algorithm (PRIMA) [66]. More recent approaches are based on the proper orthogonal decomposition [67] and other well-known general MOR techniques.