Coupled Simulation of Transient Heat Flow and Electric Currents in Thin Wires: Application to Bond Wires in Microelectronic Chip Packaging

This work addresses the simulation of heat flow and electric currents in thin wires. An important application is the use of bond wires in microelectronic chip packaging. The heat distribution is modeled by an electrothermal coupled problem, which poses numerical challenges due to the presence of different geometric scales. The necessity of very fine grids is relaxed by solving and embedding a 1D sub-problem along the wire into the surrounding 3D geometry. The arising singularities are described using de Rham currents. It is shown that the problem is related to fluid flow in porous 3D media with 1D fractures [C. D'Angelo, SIAM Journal on Numerical Analysis 50.1, pp. 194-215, 2012]. A careful formulation of the 1D-3D coupling condition is essential to obtain a stable scheme that yields a physical solution. Elliptic model problems are used to investigate the numerical errors and the corresponding convergence rates. Additionally, the transient electrothermal simulation of a simplified microelectronic chip package as used in industrial applications is presented.


Introduction
The present study is motivated from microelectronic packaging in power electronic applications, see Figure 1. An accurate and detailed thermal design is becoming increasingly important due to smaller chip sizes and, hence, increasing power densities. Therefore, an accurate electrothermal simulation is a more and more indispensable tool during design. The simulation requires the coupling of the transient heat equation with an electrokinetic problem through Joule heating effects and through temperature dependent material parameters. A challenge arising in these types of applications is the presence of thin wires used within the bonding process, i.e., the different scales of the wire and the surrounding package. The approach adapted here, which is not restricted to microelectronics, consists of modeling wires as 1D structures and computing the potential, resp. temperature, distribution using an additional 1D equation. The electric current, resp. heat flow, to the surrounding is then achieved D Γ Dir,l D ins D con Figure 1: Computational domain D with boundary ∂D consisting of a microelectronic chip package with bond wires. The package is modeled with a homogeneous conducting domain D con (chip in the center and contact pads located around the chip) and a weak conducting background domain D ins . The boundary part of a contact pad with index l is denoted Γ Dir,l ⊂ ∂D.
through a coupling of the 1D substructure with the 3D geometry. Such singular substructures are frequently encountered and dealt with in electromagnetics [1,2,3]. Here, we present a more detailed mathematical modeling of the problem both on the continuous as well as on the discrete level. An appropriate tool to model singularities in electromagnetics are de Rham currents [4] and their discrete counterparts [5]. We further show that the problem is closely related to fluid flow in porous media with fractures as studied, e.g., in [6,7]. Hence, mathematical tools for elliptic problems with Dirac measures [8] can be applied. The study was motivated by the analysis of micro-and nanoelectronic problems involving bond wire simulations in the scope of the EU FP-7 project nanoCOPS [9,10,11,12]. We use the finite integration technique (FIT) on a pair of rectilinear grids as proposed in [13,14,15]. Its use is motivated by the underlying rectilinear structure of typical chip package geometries. In this paper, due to the 1D-3D coupling, the curved geometry of a wire can be incorporated without suffering from a staircase approximation. Staircase approximations are very common for thin wire approximations in electromagnetics and a discussion of the associated drawbacks can be found in [16]. From a finite element (FE) context, the FIT can be understood as an FE method on a hexahedral grid where an appropriate mass lumping is applied to the material matrices [17,18,19]. In addition to Whitney FE, strong relations exist to other numerical schemes, such as mimetic finite differences, see, e.g. [20,21], allowing for an enhanced grid element variety.
The paper is structured as follows: in Section 2, we present the electrothermal coupled problem with the additional 1D-3D coupling. The numerical discretization is addressed in Section 3 together with a detailed discussion of boundary conditions. In Section 4, we relate the scheme to a FE method for flow problems with one-dimensional fissures. Finally, a numerical convergence rate study for a simplified elliptic model problem and the transient electrothermal simulation of a microelectronic chip package is presented in Section 5.

Continuous Electrothermal Problem
The computational domain D of the considered application consists of a microelectronic chip package with N applied bond wires as depicted in Figure 1. The subdomain modeled by conducting parts, such as the chip in the middle of the domain and the contact pads, is denoted with D con , whereas D ins denotes the insulating subdomain. Electrical signals are imposed on contact electrodes Γ Dir,i ⊂ ∂D and transferred to the chip by bond wires. For a concise notation, we assume that only one wire is present and release the need for an additional index. Due to the typically very small radius of a wire compared to the surrounding geometry, we model a wire as the 3D curve Λ with zero radius. Furthermore, we establish a 1D domain Λ by using a wire parametrization that maps from Λ to the 3D curve Λ. To distinguish quantities associated to the 3D domain D from those associated to the 1D domain Λ, we denote the latter by an overline. In the following, a formulation for this coupling is presented and we consider the general case of a wire and refer to bond wires as a specific application.

Strong Formulation
Let I := (0, t 0 ] be the time interval of interest. We consider capacitive effects only for the thermal 3D part and examine the coupling of an electrokinetic problem with the transient heat equation. Due to the small wire radii, we assume that the resistive losses in the wires are dominating and we thus neglect the losses in the 3D domain. Here, we use the Dirac distribution δ Λ to formulate the coupled 1D-3D electrothermal problem as for the 3D and 1D electric potentials ϕ and ϕ as well as the temperatures T and T , respectively, where we omit the dependency of ϕ on T and vice versa. Furthermore, all materials are regarded as time-invariant, the 3D electric and thermal conductivities are given by σ and λ, respectively, ρ is the volumetric mass density and c is the specific heat capacity. On the other hand, the 1D wire material coefficients are given by σ = |A|σ w and λ = |A|λ w , where |A| is the wire's physical cross-sectional area and σ w and λ w are the electric and thermal conductivities of the wire's material. The electrothermal coupling is established by the temperature dependent materials and the 1D resistive losses Q w . While the conductivities may depend on temperature, the temperature dependence of ρ and c is neglected. The 1D-3D formulation requires a coupling for the currents (thermal flows) and for the electric potential (temperatures). The former is denoted by the contributions q σ (q λ ) from 1D to 3D domain and by q σ (q λ ) from 3D to 1D domain. The latter is given by the coupling operator Π that will be defined in Section 2.3. Note that no explicit boundary conditions for the 1D case are required due to the strong coupling by Π. On the other hand, the 3D boundary conditions are given by where the Dirichlet (Neumann) part of the boundary is denoted by Γ Dir,i (Γ Neu ), Φ i is the constant Dirichlet potential on Γ Dir,i and n D is the outer unit normal of D. Convective heat exchange with the environment is modeled by the Robin boundary condition (1i), with the ambient temperature T ∞ and r(T ) = hT , where h denotes the heat transfer coefficient. Although omitted throughout this paper, combined convection with radiation could be imposed by setting r(T ) = hT +εσ SB T 4 , where ε denotes the emissivity and σ SB the Stefan-Boltzmann constant. The initial conditions read where ϕ init is the initial potential and T init the initial temperature. Note that the initial values ϕ init and T init are used for both the 1D and the 3D case.

de Rham Currents
In electromagnetics, singular sources, such as line and surface currents, point and surface charges are appropriately represented using de Rham currents [4], i.e., distributions on differential forms. In particular, for the inclusion of wires, line currents are introduced here. We largely avoid the formalism of exterior calculus and instead, identify differential forms with their vector proxies. The reader is referred, e.g., to [4,22] for further details on forms and currents. Let D p 0 (D) refer to the space of smooth differential p-forms. In particular, we can identify D 0 0 (D) with C ∞ 0 (D) and D 1 0 (D) with (C ∞ 0 (D)) 3 , respectively. A de Rham p-current is a map from D p 0 (D) into the real numbers. It can be defined by both vector fields and oriented manifolds. For instance, electric currents such as the current density as well as surface and line currents give rise to 1-currents. On the other hand, quantities as the energy density are represented by 0-currents, see [23]. A power density Q in the domain D gives rise to a 0-current The current density J σ := −σ(·, T )∇ϕ and heat flux density J λ := −λ(·, T )∇T , more generally J α with α ∈ {σ, λ} from now on, give rise to a 1-current An oriented curve Λ with unit tangent vector t gives rise to a 1-current with a line current I α (x), x ∈ Λ. This current may vary along the curve Λ since we allow a current exchange between the curve and its surrounding. The divergence of any 1-current K is defined as [5] div K(v) = −K(∇v), ∀v ∈ C ∞ 0 (D).
To illustrate how (1a) and (1c) can be reformulated in the setting of de Rham currents, we multiply these equations with v ∈ C ∞ 0 (D), integrate over D and integrate by parts to obtain in the sense of distributions. Hence, the de Rham current reformulation of (1a) without the wire contribution reads div J σ (v) = 0, ∀v ∈ C ∞ 0 (D).
Furthermore, by the definition of (2), the de Rham 0-current associated to ρc ∂ t T readṡ Then, the de Rham current reformulation of (1c) without the wire contribution readṡ From now on, we account for the presence of N wires by an additional index i. The embedding of a 1D wire current into the 3D domain is achieved by adding a line current to (7) and (9). To this end, let the i-th wire curve be parametrized as Λ i = {x i (s), s ∈ Λ = (0, 1)}, such that all one-dimensional quantities can be defined on the interval [0, 1].
The embedding into the 3D domain is achieved by considering the image of the map x i of the 1D current I where (x i ) * denotes the pullback by the map x i . After applying (5) to (10), and the total electric current and thermal heat flow wire contribution is obtained from the single wire 1-currents as I α := N i=1 I i α . Using (2), the wire Joule losses can be expressed as where Q w,i := σ i ∂ s ϕ i 2 . Again, the overall wire contribution Q w is obtained by summing over all wires.
With these definitions at hand, by adding (11) and (12) to (7) and (9), the 1D-3D coupled problem in terms of de Rham currents is obtained as for all v ∈ C ∞ 0 (D). ϑ s Λ Figure 2: Cross-section of a wire with its radius confined to a 1D line Λ in s-direction. Further denoted are the coupling radius r cpl and the angle ϑ around the wire.

Coupling Operator
What is left is the definition of the coupling operator Π as used in (1e) and (1f). Due to the singularity of the 3D solution at Λ, a direct coupling of 3D and 1D solution at Λ is not possible. Instead, following [8], we use an averaging scheme given by where u ∈ {ϕ, T } and γ is a scaling coefficient. Note that (s, r cpl , ϑ) refer to cylindrical coordinates around Λ, where r cpl is the coupling radius, see Figure 2. Typically, r cpl r, where r is the radius of the cylindrical wire, is chosen to circumvent resolving the wire.
For a 3D diffusion problem with Dirac right hand side as given by (1a) and (1c), cylindrical coordinates, an infinite domain D and Λ coinciding with the z-axis, the solution has the form where r 0 is a reference radius. With this logarithmic solution, γ = log(r/r 0 )/ log(r cpl /r 0 ) scales from the averaged value to the physical value of the line source solution at the wire radius r. With this definition of γ, applying (14) to (15) yields the corresponding 1D solution given by Finally, for arbitrary curves Λ, the above solutions are valid in a sufficiently small neighborhood of Λ.

Discrete Electrothermal Problem
Discretization of the electrothermal problem is carried out using the FIT [13,14,15] motivated by its natural relation to discrete differential forms [24]. We emphasize that no conceptual difficulty would arise when a comparable discretization scheme as, e.g., FE, would be used instead. A more detailed discussion of this aspect and the numerical discretization errors are given in Section 4. In this section, we define discrete currents in analogy to the continuous de Rham currents. Then, we successively introduce the 3D and 1D discretizations. Further, the averaging of the materials requires a dual grid that is presented in Section 3.3. Finally, in Section 3.4, boundary conditions are formulated and included in the discrete system of equations.
For the discretization, we choose standard lowest order nodal functions on a rectilinear grid, also referred to as nodal Whitney functions. Whitney functions are the discrete counterpart of differential forms and hence, discrete currents are defined as maps from the space of Whitney functions into the real numbers. In this way we obtain the discrete counterpart of (13) (denoting discrete de Rham currents with the same symbols as their continuous counterparts) as where E k refers to the nodal Whitney function associated with node k. In the following, we will derive the expressions that are required for the implementation of (17).

3D Discretization
For the discretization, a rectilinear grid with N N nodes, N E edges, N F facets and N C cells is used. With k = 1, . . . , N N and l = 1, . . . , N E , we denote with P k and L l the nodes (points) and edges (lines) of the grid, respectively. In this subsection, we omit the wire contribution. The standard discrete gradient operator is introduced as the primal node to edge incidence The discrete counterpart of (5) at the grid level reads for any discrete 1-current K and possible boundary contributions BC, see [5]. This boundary contribution is omitted for the time being and is discussed in Section 3.4. Thus, by applying (19) and evaluating the Whitney functions on the grid, we obtain with the coefficient vector j α of the current J α , see again [5], and the unit basis vector e k . Applying (20) to (17), we obtain the discrete system whereQ ρc denotes the coefficient vector of the currentQ ρc . The potential and temperature evaluated at the grid points are the degrees of freedom given by ϕ ∈ R N N and T ∈ R N N , respectively. Figure 3: 1D discretization of wire i visualized using lumped elements for the example of N 1D,i = 6 points with thermal and electrical resistances. The wire's current is discretized using piecewise constant basis functions.

1D Discretization and Coupling to 3D Domain
The FIT-discretization of I α results in a piecewise constant discretization on every wire part (curved element) Λ i j for wire i = 1, . . . , N . We model each wire part as a 1D element giving a series connection of N 1D,i − 1 elements for wire i, cf. Figure 3. We obtain the aforementioned partitioning of Furthermore, let represent the electrical/thermal 1D wire mass matrix. Then, the FIT discretization of the 1D current and heat flow reads We proceed by discretizing (11) as where we have introduced the coupling matrix ). It should be noted that, in view of (25), no approximation of the (possibly curved) wire geometry by the underlying 3D grid is required, which is the main advantage of the 1D-3D coupling. To simplify notation, we further introduce Figure 4) and infer the vector representation for div I i σ , combining (24) and (25) as Remark 1 The coupling matrix R i N can be understood as the discrete counterpart of the pullback (x i ) * . For the simple case when the 1D nodes are obtained by pulling back the 3D grid nodes, R i N ∈ {0, 1} N 1D,i ×N N becomes an operator that restricts to those nodes of the grid which are connected to wire i and contains only one non-zero entry in each row. In the general case, each 1D node is allocated in a grid volume defined by 8 nodes, this operator becomes R i N ∈ R N 1D,i ×N N and contains 8 non-zero entries in each row, obtained by Whitney interpolation. In any case, the sum of all entries in one row is always equal to one.
To account for the 1D Joule losses, we introduce the vector Q w,i , allocated at 1D edges, with N 1D,i − 1 entries each given by representing the discretized heat power of wire i with I i σ,j being the current in element j. Starting from (12), the Joule losses of the wire part are discretized as where the entries of X i abs are given by the absolute values of the entries of X i . Let T i avg = 1 2 X i abs T denote the averaged temperatures for each element of wire i. We account for nonlinearities in the wire material parameters as α i j ((T i avg ) j ). The last ingredient for the discretization is the discretized version Π ∈ R N 1D ×N N of the coupling operator Π as introduced in Section 2.3. It is obtained by interpolating x * Πu at the nodes of the 1D grid to yield the relations for wire i. Finally, by introducing the lumped wire stiffness matrix K w α as we conclude

Dual Grid
In this section, we introduce the dual rectilinear grid containing N N = N C points, N E = N F edges, N F = N E facets and N C = N N cells. We denote by A l , V k the facets (areas) and cells (volumes) of the dual grid, respectively. The non-singular quantities (i.e. not related to wires) can be identified with elements on the dual grid, as commonly done in FIT (cf. Figure 5). For instance, the discrete de Rham current J α is defined by the coefficients with the unit normal n A l of facet A l . Also, the discrete de Rham current Q is represented by Applying the standard FIT material approximation to the unknown current J α , we obtain where | · | denotes the area or length of the specified geometrical object. The α avg l are averaged material values defined on primal edges (or dual facets) that are obtained by a suitable averaging scheme [24]. Furthermore, t l refers to the unit tangent of L l and ϕ l,0/1 to the potential at the start/end point of L l , respectively. This defines the electrical and thermal conductance matrices Equivalently, the term for the heat powers gives where the ρc avg k are average material values defined on primal nodes (or dual cells) and are obtained by a suitable averaging scheme [10]. The thermal capacitance matrix M ρc ∈ Note that M α and M ρc are diagonal matrices with strictly positive entries in the FIT approach and sparse positive definite matrices in general, e.g., in a Galerkin FE setting. Hence, we infer To lighten the notation, we introduce the stiffness matrix K α (T) := G M α (T)G and collect the 1D potentials (temperatures) of all wires in a matrix denoted by ϕ ∈ R N ×N 1D (T ∈ R N ×N 1D ). Then, combining (21) with (38), the discrete electrothermal problem with wire contribution given by (28) and (31) reads By setting K α (T) := K α (T) + K w α (T) and Q(ϕ, T) := Q w ϕ, T , we can write (39) in a more compact form given by System (40) is a differential algebraic equation (DAE) and can be integrated in time by using a suitable integration scheme. Here, we use the implicit Euler method together with a fractional step method [25], which reads at time step n with uniform step size ∆ t The fractional step method has the same order of accuracy as the implicit Euler method, however, system (41) is decoupled and hence, easier to solve.

Boundary Conditions
Boundary conditions are applied for the electrical and thermal subproblem separately. Starting with the electrical part, to impose Dirichlet boundary conditions, we follow [26] and introduce the notion of active nodes, see Figure 6a. More precisely, the N N,a -active nodes associated with Γ Dir are obtained by removing all nodes contained in the closure of Γ Dir . Then, with the restriction operators R a N ∈ {0, 1} N N,a ×N N and R Dir N ∈ {0, 1} N N,Dir ×N N (denoted as trace operator T in [26]), we decompose the electric potential into active nodal potentials ϕ a = R a N ϕ and prescribed nodal potentials contained in the closure of Γ Dir given by ϕ Dir = R Dir N ϕ. Based on these definitions, we can reduce (40a) as by incorporating Dirichlet conditions. To shorten notation, we introduce and we obtain the electrical part with boundary conditions For incorporating the Robin boundary condition in the thermal equation, it is convenient to consider the augmented dual grid [26], where N N,∂D = N F,∂D dual facets on the boundary are introduced that complete the boundaries of the dual cells, see Figure 6b. For the augmented dual grid, the boundary contribution BC in (19) does not vanish.
Following [26], the additional term is quantified as where R ∂D N ∈ {0, 1} N N,∂D ×N N and R ∂D F ∈ {−1, 0, 1} N F,∂D × N F denote the restriction operators to the boundary nodes and the boundary dual facets (of the augmented dual grid), respectively. Note that the signs in R ∂D F are chosen in such a way that the orientation of the restricted fluxes with respect to the outer normal is taken into account. Let us define where k i refers to the index of the i-th boundary node, i = 1, . . . , N N,∂D and let T ∞ be a constant vector containing entries T ∞ . Then, by using the discrete Robin boundary condition and the notation M ∂D := ( R ∂D N ) M ∂D h R ∂D N , we obtain the electrothermal system with boundary conditions

Relation to Flow Problems with Fissures
Let us focus on the 1D-3D coupling and consider a simplified elliptic model problem, which represents, for instance, the electrical subproblem. We show that this model problem is closely related to a 1D-3D coupled problem of blood flow through tissues with thin tubular structures for the vessels as analyzed in [8,27]. More precisely, we show that the simplified elliptic wire problem is obtained in the limit of infinite permeability of the vessel-tissue interface. We thereby rely on the well-known similarities between the FE and the present FIT approach [17]. Several papers consider a singular Dirac right hand side for elliptic problems in an FE context [28,29]. It should be emphasized that in the present approach, as well as in [8,27], the singular 1D contribution depends on the solution itself, which further complicates the problem.
In a first step, we introduce a weak formulation together with the associated solution spaces for the model problem, where, for simplicity, we assume a linear conductivity. Furthermore, we consider only one wire Λ = {x(s), s ∈ Λ = (0, 1)}, such that its index can be omitted. We denote with u and u the 3D and 1D component of A major difficulty in the analysis is the singularity of u at Λ. In particular, u exhibits a logarithmic singularity and does not belong to the standard space H 1 (D), see Section 2.3. Instead, the solution space can be identified with the weighted Sobolev space V δ = H 1 δ (D), see [27,Section 2], where the weight is given by the distance to Λ to the power of 2δ with δ ∈ (0, 1). On the other hand, the 1D-variables belong to the space V := H 1 (Λ). We consider homogeneous Dirichlet boundary conditions on D for simplicity. The weak form of (49a) reads, find u ∈ V δ , such that where ·, · D refers to the duality product of V δ and V −δ , whereas (·, ·) {Λ,D} refers to the L 2 -inner product on {Λ, D}. Also, tr Λ denotes the trace operator from V −δ to L 2 (Λ). In the setting of this paper, the left hand side of (50) represents the divergence of the 3D current (heat flow), whereas the first term on the right hand side represents the current (flow) transport density from the 1D subdomain into the 3D surrounding. The weak form of (49b) reads, find u ∈ V , subject to see [27] for details. In [8], the law for q α is chosen as q α = β((Πu)•x−u), where β = | d ds x|β •x and β refers to the permeability coefficient for the vessel-tissue blood transfer. Note that the densities q α and q α in (50) and (51) are related through q α = | d ds x|q α • x. The coupled formulation reads, find (u, u) ∈ V δ × V subject to representing standard piecewise linear, globally continuous polynomial basis functions. The 3D polynomials coincide with the Whitney basis functions of Section 2. The connection from FE to FIT is established by introducing grid dependent inner products as follows Here, Λ j is the 1D dual element given by half of Λ j−1 and Λ j . This FE-FIT connection can be viewed as applying the trapezoidal rule in each dimension, see [30] 1 for details. Then, we Problem (57) is equivalent to the matrix system of equations Multiplying (58b) with R N and adding it to (58a) we obtain which is the same as (39) without the transient part and the source term f . It should be noted at this point that our coupling does not require the implementation of the complicated boundary term on the right hand side of (53). Contrary to (39), the coupling condition (29) is not directly contained in (58), but recovered in the limit β → ∞. In our setting, β represents a contact conductivity and hence, we consider the limit of perfect conduction between 1D and 3D part. Indeed, (58b) is equivalent to and passing to the limit β → ∞ yields Πu − u = 0. A strategy to bound the total error of the numerical scheme could consist in using the relation to the FE method, established above, and the triangle inequality as D'Angelo [27] has shown that for the FE error of the 3D variable there holds with ε ∈ (0, δ), provided that the solution is sufficiently regular, β small enough and the mesh grading strong enough. The same estimate was established for the 1D variable in a suitable norm. The norm · appearing in (62) is defined as where in L 2 δ (D) the weight is again given by the distance to Λ to the power of 2δ. Moreover, V 2 1+ε refers to a Kondratiev-type weighted space, see [27] for details. The second term on the right hand side of (61), representing the difference between FE method and FIT, is of order O(h), which can be obtained by bounding the error associated to the trapezoidal rule. Yet, it remains to show that the system (59), (60) is well-posed. Additionally, the limit β → ∞ for the continuous problem (52) and its FE approximation (53) is not included in [27] and has to be analyzed separately, which is beyond the scope of this paper.

Grid Generation with Local and Global Grading
As introduced in Section 3, we apply the FIT on a pair of rectilinear 3D grids. Any 3D grid that fulfills these properties is denoted by G and is constructed using a Cartesian product of 1D grids. For all grids G, the average 3D edge length is denoted by h. On the other hand, a 1D grid is used for the discretization of the wires and is denoted by G. We apply the same 1D grid for all wires with respect to their parametrizations. The 1D grid G is chosen to be equidistant in terms of the wires' parametrization s with the step size denoted by h. Note that for the numerical implementation, the 1D grid points of G coincide with 3D grid points of G. Therefore, for straight wires and equidistant 3D grids, the 1D step size h is always a multiple of the 3D step size h.
Due to the singular solution of (1), a graded grid is required to recover the expected convergence rates. Let us consider a 3D domain D with a line on which the solution becomes singular. We assume that an initial, equidistant grid G 1 is given, see Figure 7a. The region on which the refinement is applied is denoted by D r ⊂ D. As the 3D grid is composed by the Cartesian product of individual 1D grids, the refinement strategy for a 1D grid is described here. Based on [32], the 1D refinement around a singular point x 0 is applied using the layers where N is the number of refinement layers, µ the grading and b the radius of the refinement. Note that µ = 1 results in an equidistant refinement whereas µ → 0 gives a stronger grading towards the singularity. Furthermore, b is chosen such that no grid points of G 1 fall in the refinement region D r . The only exception to this rule is a possible grid point at x 0 . Note that for a refinement with N refinement layers with the assumption that x 0 was part of G 1 , 2N 1D points are added. Due to the rectilinear 3D grid, this results in a propagation of the refinement along the coordinate directions. Thus, 2N additional 1D points result in (2N + 1) 2 − 1 additional points for a two-dimensional refinement and (2N + 1) 3 − 1 additional points for a three-dimensional refinement. For line sources that are the main subject of this paper, a two-dimensional refinement is required. After applying the tensor product on the individually refined 1D grids, a locally graded 3D grid is obtained and denoted by G b µ,N , see Figure 7b.  With the local grading as introduced above, a global grid grading can be defined as a specific choice of the local grading. Let an initial grid contain the singular point only. Then, by choosing the refinement radius b = d/2, where d is the width of D, the refinement given by (64) results in a refinement on the entire computational domain D such that D r = D. A grid obtained by this special choice is denoted by G µ , see Figure 7c. Note that an additional choice of µ = 1 results in an equidistant grid G 1 without any refinement nor grading and the same number of grid points as G µ .

Numerical Examples
The model problems consist of one wire embedded in a homogeneous cube such that an exact representation of the geometry using a rectilinear grid is possible. In particular, three different models are considered. First, the coupling is realized point-wise using an external circuitry resulting in a 0D-2D coupling approach. Secondly, a straight wire embedded in the cube is simulated in a 1D-3D coupling approach. Thirdly, the general case of a bent wire is considered. For the first two models, an analytical reference solution is available while a fine reference is used to study the convergence of the third model. For these model problems, the electrical problem given by (48a) is solved.
The considered cube D is of side length d = 1 m and conductivity σ = 1 S/m. In all examples, the wire's radius is r = 1 µm and its cross-sectional area is A = πr 2 . The conductivity of the wire is given by σ = 1 × 10 15 Aσ manifesting an electrical example opposing thermal examples for which σ ≈ 1 × 10 2 Aσ. As discussed in Section 2.3, the coupling coefficient γ is used with a reference radius r 0 that we choose as the radius of an equivalent 2 cylindrical (or circular in the 2D case) geometry such that r 0 = d 2 /π. As a reference, Table 1 summarizes the described parameters.
To quantify the errors of the method, we introduce the following error measures. Let ϕ h and ϕ h denote the 3D and 1D FIT solution vectors, respectively, and ϕ and ϕ the corresponding analytical solutions. Since we neglect the Joule losses in the 3D domain (cf. Section 2), we consider the error in the solution's derivative only for the 1D solution. First, we define where ϕ and ϕ refer to the analytical solutions evaluated on the same grid as ϕ h and ϕ h . 2 The term equivalent refers to a geometry of equal cross-sectional area as the cube Secondly, where the norm of the analytical solution is compared to the norm of the FIT solution. Lastly, if no analytical solution is available, we use where ϕ and ϕ are FIT solutions computed on a very fine grid. The discrete and continuous norms used in (65)-(67) are defined by where DṼ, D S and DS are diagonal matrices with the 3D dual volumes, the 1D primal lengths and the 1D dual lengths on the diagonal, respectively. These matrices coincide with M ρc , M β and M

0D-2D Coupling
As a first validation, we consider a brick-shaped resistor with parameters as given in Table 1. The perfect electric conducting (PEC) wire of radius r and the surrounding boundary serve as the resistor's inner and outer electrodes. The coupling is established by connecting the inner and outer electrode using a series connection of lumped resistor and voltage source, see Figure 8a. Applying the parameters as shown in Table 1, we use the thin wire assumption for the inner electrode and thus model it by a single point as shown in Figure 8b. Then, the coupling is identified as a 0D-2D coupling. The dashed circle depicts the coupling circle of radius r cpl (cf. coupling condition (14)) that is shown in a magnified view in Figure 8c. To avoid resolving the inner electrode, a coupling radius of r cpl = max l L xy l r is chosen, where L xy l iterates over the lengths of all edges perpendicular to the wire. Assuming the resistance per unit length of the rectangular resistor to be given by Figure 8: (a) shows the cross section of a brick-shaped resistor with wire-shaped inner PEC electrode connected to an external circuit. In (b), the inner electrode is replaced by a 0D representation and the coupling circle is shown as a dashed curve. (c) shows a magnification of the coupling circle while also depicting the actual wire radius r and the wire's current I.
the resistance of the external resistor to be R 0 = 1 Ωm and the applied voltage as V 0 = 1 V, the potential ϕ at the inner electrode is of interest. For the current per unit length I 0 = V 0 /(R 0 + R int ) and a homogeneous conductivity σ, the analytical 2D solution of Laplace's equation is given by where r 0 = d 2 /π is the distance from the origin to the reference potential. After applying the coupling condition (14) to (71), the 0D potential ϕ = Πϕ is obtained. Applying (39a) to the here considered 0D-2D coupling, it simplifies to where, by abuse of notation, the matrices are 2D modifications of the usual 3D matrices and G 0 = (R 0 ) −1 . We impose the reference solution ϕ on the boundary of the domain (method of manufactured solution), apply the boundary conditions to (72) as described in Section 3.4 and solve the resulting system to obtain the FIT solution ϕ h . In this section, h is the average edge length of all edges in x-and y-direction. Using the analytical solution of (71) as reference, the convergence of the relative errors ε 3D L 2 and ε 1D L 2 for different grids with respect to h is shown in Figure 9. Considering the results for ε 1D L 2 , we first observe that a uniform grid G 1 gives a convergence rate far lower than one. Secondly, if we use a globally graded grid, the convergence order can be improved to a value around one. Thirdly, a locally graded grid G 10 −4 1,10 gives only a slightly higher convergence rate than G 1 . The convergence rates for the 2D error ε 3D L 2 behave similarly.

3D-1D Straight Wire Coupling
In this section, we consider a 3D-1D coupling and apply the FIT to solve the electrical problem. The investigated model consists of a straight wire in z-direction positioned in the xy-center of a cube of side length d that we refer to as the computational domain D, see Figure 10a. We again use the parameters as summarized in Table 1. Additionally, for the 1,10 Figure 9: Convergence of (a) ε 1D L 2 and (b) ε 3D L 2 with respect to h and different grid choices for the 0D-2D coupling.
where we choose I (z) = I 0 z/d and I 0 = 1 A/m. The analytical solution (73) is impressed on the boundary ∂D. Then, (48a) is solved for the 3D solution ϕ h and the 1D solution ϕ h . In Figure 11, ϕ h is plotted on Λ while ϕ h is plotted on D \ Λ for G 0.5 , h ≈ 1.822 × 10 −2 m and h = 3.125 × 10 −2 comparing r cpl = 0 m and r cpl = max l (L xy l ). The convergence of the relative errors as defined by (65) and (66) is investigated. In Figures 12 and 13, the convergence of the errors ε 3D L 2 and δ 1D L 2 is shown. In Figure 12a, the convergence of ε 3D L 2 with respect to h is shown for h = 3.125 × 10 −2 , r cpl = max l (L xy l ) and different choices of the 3D grid. For the local gradings G b 1,10 and G b 0.5,10 , the asymptotic convergence order is around one, independent of µ and similar to the order    Figure 13: Convergence of ε 1D L 2 with respect to h for a fixed but different h for each 3D grid and for the straight wire 1D-3D coupling.
for a global grid grading G 1 . For a globally graded grid G 0.5 , we observe a significantly higher convergence order than for an equidistant grid G 1 . In Figure 12b, the influence of the coupling radius r cpl on the error and on the convergence order is investigated. For a zero coupling radius, the 1D solution is taken directly from the 3D solution in the coupling points where the 3D solution is singular. Therefore, as expected, the error is large and the convergence very slow. For r cpl > 0 m, a convergence order of around three is observed while the solution and the order is independent on the exact choice of r cpl . In Figure 13, the convergence of δ 1D L 2 is shown with respect to h for a fixed but different h for each 3D grid choice. For large h, a wire is modeled rather by point sources along Λ than by a line source. For the limit case of h → 0 however, the line source case is recovered and the error with respect to the line source reference solution becomes smaller. The convergence is of almost second order and independent of the 3D grid and/or grading choice. This indicates that the 1D discretization error is dominating for the considered error measure and the considered grids. On the other hand, for the here considered grids, the errors ε 1D L 2 , ε 1D H 1 and δ 1D H 1 are smaller than 1 × 10 −4 and do not show further improvement for a refinement of the 1D grid. This behavior is attributed to a dominance of the 3D discretization error for this setting.

3D-1D Bent Wire Coupling
We aim for simulations of problems in which thin wires follow arbitrary curves. Therefore, in this section, we investigate the case of a single bent wire. Again, we use the parameters of Table 1 on a cube D as computational domain. To set up a model of a bent wire with path Λ as introduced in Section 2.2, we apply a parametrization using a Bézier curve given by with the start and end coordinates length d PEC = 40 mm, see Figure 10b. This setup is in analogy to the typical case that a wire connects two PEC contact pads. For the construction of the grids, we start with an equidistant 1D grid for the parameter s giving a characteristic 1D step size h. From the 1D grid, the 3D wire points are determined by (74) and, together with the boundary nodes in each direction, form the basis of the 3D grid. Additional grid lines are inserted due to the error evaluation domain Ω, because of the PEC cubes and between the wire and the boundary y = d. The number of grid lines used for the latter is given by N 1D /4 . We apply 0 V at the PEC electrode at x 0 , 1 V at the PEC electrode at x 1 , solve (48a) and consider the 1D solution as quantity of interest. We use a coupling radius of r cpl = 1 × 10 −2 /κ, where κ ≈ 6.08 m −1 is the maximum Frenet-Serret curvature [33,34] along the curve. In the following, we refer to the 1D and 3D reference solutions ϕ and ϕ as the solutions computed using h ≈ 1.884 × 10 −2 m and h = 1.563 × 10 −2 . In Figure 14a, ϕ is shown with respect to the wire parametrization s and Figure 14b shows the solution ϕ h on Λ and ϕ h on D \ Λ using a 3D visualization. Investigating the convergence of the error, we plot ∆ 1D L 2 and ∆ 3D L 2 with respect to h in Figure 15 using ϕ and ϕ as the reference. For both ∆ 1D L 2 and ∆ 3D L 2 , we observe a convergence order of around two.

Industry-Relevant Microelectronic Chip Package Simulation
We consider a microelectronic chip package model based on [10], see Figure 1. The contact pads are centered in z-direction with a height of 100 µm. The radius of all wires is r = 1 µm with a circular cross-sectional area given by A = πr 2 . The curvature of the wires is given by (74) and a height of H = 0.1 mm, where x 0 and x 1 differ for the different wires.
For simplicity, we assume all material characteristics to be linear and model D con with an  Figure 15: Convergence of (a) ∆ 1D L 2 and (b) ∆ 3D L 2 with respect to h for the bent wire 1D-3D coupling and r cpl = 1 × 10 −2 /κ.

Parameter
Description Value   electric conductivity σ con → ∞ (PEC) and D ins with σ ins = 1 × 10 −4 S/m. On the other hand, the thermal conductivities are given by λ con = λ Cu and λ ins = 0.87 W/(Km), where λ Cu = 401 W/(Km) is the conductivity of copper. The volumetric mass densities and specific heat capacities are given by ρ con = ρ Cu , ρ ins = 1500 kg/m 3 , c con = c Cu and c ins = 882 J/(Kkg), respectively, where ρ Cu = 8930 kg/m 3 and c Cu = 390 J/(Kkg) are the values for copper. We choose the wires' conductivities to be σ = Aσ Cu and λ = Aλ Cu , where σ Cu = 5.96 × 10 7 S/m is the conductivity of copper. As initial conditions, the package is homogeneously set to ϕ init = 0 V and T init = 300 K. A voltage of V = 100 mV is applied over each wire, where the central chip region is used as the ground potential. In terms of electrical boundary conditions, this translates to Dirichlet conditions on Γ Dir,l and homogeneous Neumann conditions on ∂D \ Γ Dir,l . Thermal boundary conditions are chosen to be convective on ∂D with the heat transfer coefficient h = 25 W/(m 2 K) and the ambient temperature T ∞ = 300 K. For each of the wires, the 1D-3D coupling is carried out using r cpl = 1 × 10 −4 H 2 κ and N 1D = 4. All mentioned simulation parameters are summarized in Table 2. The applied grid results from a commercial meshing tool for the 3D part and a subsequent insertion of grid lines at the points given by (74) to yield h ≈ 8.51 × 10 −5 m. The time discretization is given by N t = 10 with t 0 = 1 s. For a choice of h = 1/3, we solve (48) and obtain the results for ϕ h (t 0 ), ϕ h (t 0 ), T h (t 0 ) and T h (t 0 ) as shown in Figure 16 and 17. The presented results demonstrate the capability to predict the temperature distribution in a microelectronic chip package including the temperature profile of the bond wires.

Conclusion
Motivated by the electrothermal simulation of microelectronic chip packages including thin bond wires, we presented an approach that alleviates the necessity of fine grids to resolve the geometry of the wires. In the literature, such a problem is known as a 1D-3D coupling and  Figure 1). (b) shows the 3D solution T h (t 0 ) at z = 150 µm in the xy-plane.
was, e.g., investigated by [8] in the context of fluid flow in porous media with fractures. The main challenge of this problem was identified to be the solution-dependent line source term in the 3D Laplace equation. Using the setting of de Rham currents to describe the arising singularities, we proposed a continuous formulation of the 1D-3D electrothermal coupling. The coupling condition follows the work in [8] to average the 3D solution around a wire for the definition of the 1D solution. Additionally, we introduced a scaling factor that accounts for the distance from the wire to yield a physical solution for small but arbitrary coupling radii.
The discrete system was set up using the discrete counterpart of de Rham currents in a FIT formulation. A detailed description of all the involved discretization steps was presented. The theory concludes by showing that the relation of the present thin wire problem to finite element method (FEM) problems for fluid flow problems with fractures is established by an infinite permeability of the vessel-tissue interface.
To investigate numerical errors and convergence rates, electric model problems for a 0D-2D coupling and a straight and bent wire 1D-3D coupling were considered. The convergence rates of the corresponding errors were analyzed for different choices of the grid, the grid grading and the coupling radius. It was shown that the grid grading can improve the convergence rate substantially while the solution is independent of the coupling radius. We could see in particular that, due to the singularity of the 3D solution at the 1D domain, a direct coupling of 1D-and 3D-points (r cpl = 0) results in a very slow convergence. Lastly, a transient electrothermal simulation of a chip package with 12 applied wires was presented.
The numerical analysis of the nonlinear, transient and coupled electrothermal problem is still an open research topic. In particular, existence and uniqueness need to be established for case of an infinite permeability, which is not included in the analysis presented in [27]. Additionally, the error analysis for the special case of the FIT discretization is of interest.