Equivalent circuits for electromigration

Electromigration continues to be a major concern for integrated circuit design. The susceptibility to electromigration is assessed by tracking the stress in metal lines under the influence of applied currents, which can be computationally expensive for large chips. Over the last few years, an efficient approach for tracking stress in large interconnect networks has been developed, and well-received, in part because it provides a model for stress dynamics in the form of a standard linear time-invariant dynamic system. In the context of this model, we will show that the dynamic behavior of the stresses and fluxes in metal lines is exactly identical to the dynamic behavior of voltages and currents in certain RC circuits that can be easily constructed for the metal lines. Thus, electromigration assessment for any metal interconnect structure can be done by simply simulating the corresponding equivalent RC circuit. This opens the door for benefiting from well-known techniques for fast circuit simulation, as well as methods for macro-modelling and analysis of RC circuits, in order to improve the performance and capacity of practical methods for electromigration assessment in large circuits.


Introduction
In a metal line carrying significant current density, the free electrons push and move the metal atoms in the direction of the electron wind, i. e., towards the anode end of the line; hence the name electromigration (EM) for this effect. The resulting atomic flow increases compressive stress at the anode and tensile stress at the cathode, which creates a stress gradient that presents an opposing driving force that retards EM [1]. If the levels of stress become high enough, a void may be created due to high tensile stress near the cathode, or a hillock (extrusion of metal through cracks in the dielectric) may form due to high compressive stress near the anode, which can either way result in circuit failure. With the confinement of metal lines in today's metal technology, voids are much more likely than hillocks and so one is often more concerned with the buildup of tensile stress. We will follow the common convention that tensile stress is positive and compressive stress is negative. A void is created once the stress exceeds a certain level of stress, called the critical stress, denoted σ crit , which is an effective parameter that depends on a number of factors. Stress is measured in units of pascal (Pa), and a typical value of the critical stress today [2] is about 500 MPa.
We focus on the on-chip power distribution network (PDN) because it is generally more susceptible to EM due to the fact that it carries mostly uni-directional currents. The PDN consists of the power grid and the ground grid. Modern grids span multiple layers (often all the layers) of metallization and they consist of meshes of power and ground lines. Without loss of generality, we will focus on the power grid. Under the influence of EM, metal atoms can travel between different connected branches on the same layer. However, in the modern dual-damascene semiconductor process, they cannot travel through a via to other metal layers above and below, because of the metal liner under every via which acts as a barrier to atomic movement but allows electron movement. As a result, EM-induced metal transport within a layer remains within that layer, so that the overall analysis problem is decomposed into sub-problems on different layers. Within any given layer, one will typically find a large number of such physically disconnected portions of the power grid. The vast majority of these structures turn out to be trees, i.e., they have no cycles. So, it is typical in the field to simply use the term interconnect tree to refer to these metal islands.
For any given interconnect tree under the influence of EM, the time duration before any voids have nucleated (i.e., formed) is called the nucleation phase. In this phase, the resistance of a line remains the same as that of a fresh (undamaged) line. Once a void has nucleated, the void growth phase begins and the void starts to grow. The growth rate is initially fast and can in some cases, depending on layout geometry, quickly cause an open circuit leading to early failures. In other cases, the line may continue to conduct current after void nucleation, and the void continues to grow in the direction of the electron flow until it saturates at some steady-state size. Correspondingly, the line resistance also increases towards some steady-state value. In order to track the evolution of the line (and tree) towards voiding and beyond, one needs to simulate and track the values of stress in the lines over time.
We will see in this paper that the relationship between stress and atomic flux in metal lines is identical to that between voltage and current in electrical circuits. Indeed we will see that, for any interconnect tree, we can construct an RC circuit whose voltage evolution over time is identical to the stress evolution in the tree. The two systems of equations are exactly identical, so that the RC circuit is an equivalent circuit. This is not simply an interesting curiosity; it can be very useful to draw on the extensive research and knowledge in the art of simulation of electronic circuits in order to improve the capabilities for EM assessment. For one thing, this would allow easy simulation of stress in interconnect trees using standard simulators like SPICE, thus benefiting from the many techniques that circuit simulators use to speed up the analysis. Also, as we will see, the stress analysis of a metal line can be replaced by the voltage analysis of an RC-chain, which is a very well-studied circuit structure, with many theoretical and numerical techniques for analysis and modelling, including reduced-order modelling and macromodelling.

Background
We will give a detailed review of the one-dimensional (1D) EM model developed by M. A. Korhonen [3], which we will refer to as the Korhonen model. Consider a metal line carrying current density j, with uniform height h and width w, confined in a rigid dielectric material with line length along the x axis, as shown in Fig. 1. By necessity, a 1D model assumes that there is uniformity across the other two dimensions, i.e., across the metal line cross-section in this case. Specifically, we assume that stress, electric current density, material density, as well as atomic transport rates are all uniform across the wire cross-section. This assumption is central to the Korhonen model, but it is only valid for interior sections of uniform metal lines, which excludes certain layout features of interconnect trees. Specifically, three layout features of the metal interconnect merit attention [4]: (i) the case where one or more lines on the same metal layer are all connected to a via that provides electrical conductivity to another layer, above or below; (ii) the case where three or more lines on the same layer are connected at the same physical location, without necessarily a via to another layer; and (iii) the case where one metal line connects to another line of different width. Because voids are commonly found at vias, the first case is of special importance and merits special case treatment. The other two cases also require special case treatment because of the abrupt changes in electrical current and, therefore, in atomic flux that are exhibited in these structures. We refer to the layout features of these three special cases as a layout junctions, and we categorize them according to the number of lines involved, so-called the degree of the junction. Following the terminology of [4], as shown in Fig. 2, a degree-1 junction (i.e., a junction of degree 1) is called a diffusion barrier, a degree-2 junction is called a dotted-I junction, a degree-3 junction is a T-junction and a degree-4 junction a plus-junction.
Because of the metal liner and the capping layer around every metal line, mass is conserved inside an interconnect tree, so that any change of material density in a metal line results in a change of the stress (pressure) in that line. Specifically, with C(x, t) as the atomic concentration in the line and σ(x, t) as the stress, at distance x from one end of the line at time t, a relative change in concentration corresponds to an incremental change of stress, i.e., as shown by Korhonen [3] based on the early work of Eshelby [5], where B is the effective bulk modulus [6], which depends on metal line geometry and microstructure. It has a typical value [2] of about 30 GPa in modern technology. Let C 0 be the concentration under zero stress, so that C 0 = 1/ Ω, where Ω is the atomic volume for the metal in the line, which is ≈1.66E-29 m 3 for Copper. Integrating Eq. (1) from (0, C 0 ) to (σ, C) gives The Korhonen model is a combination of two equations, an atomic flux equation that keeps track of atomic mass transport, and a continuity equation that enforces mass conservation, as we will review and sketch in the following simplified description.

Atomic flux equation
The first equation is based on the consideration that migration of metal atoms is primarily the result of (i) atomic drift due to the bombardment by the electron wind arising from the electric field in the metal line, E = ρj, where ρ is the metal resistivity and j is the current density, and (ii) atomic diffusion due to the stress gradient ∂σ/∂x, resulting in an opposing flux. We will use ϕ(x, t) to denote the atomic flux in the line, whose units are number of atoms per second per unit crosssectional area. The distance variable x, the electric current in the line (of density j) and the atomic flux ϕ are all assumed to have the same reference direction. So, both current and flux are positive when they are in the positive x direction. Considering the two opposing fluxes, one arrives at 1 where D is the coefficient of atomic diffusion, also simply called the diffusivity. Its value depends on the material properties and on temperature as D = D 0 exp (− E a /k b T), where D 0 ≈ 5.20 × 10 − 5 m 2 /s in modern technology and E a is the activation energy for atomic diffusion, typically taken as E a ≈ 1 eV, so that D ≈ 1.31 × 10 − 17 m 2 /s is typical today. Other terms include T as the metal line temperature in Kelvin; k b = 1.380 649 × 10 − 23 J/K is Boltzmann's constant; q* = |e|Z is the effective charge, where |e| = 1.602 176 634 × 10 − 19 C is the absolute value of the electron charge, and Z is the effective valence, which is a scalar whose value  Typical interconnect tree with various junction types. 1 Eq. (2) in [3] for atomic flux is in terms of the divergence of the chemical potential function, but it's also provided in terms of stress, and because divergence in 1D is simply the derivative, one gets the expression shown here as Eq. (3).

Continuity equation
The second part of the Korhonen model is a continuity equation that describes the stress evolution over time in response to the material transport. It expresses a simple balance: the more material that flows out of a region, the less material there is inside it over time. This may be tracked using the mathematical concept of divergence of a vector field. Divergence measures the net rate of reduction of a certain quantity (in this case, the number of atoms) within an infinitesimal region around a point, per unit volume. Divergence has both a differential form and an integral form. In 3D space, with orthonormal unit vectors u x , u y and u z , the flux Φ is a vector quantity with three components, so that Φ = ϕ x u x + ϕ y u y + ϕ z u z , and the differential form of the divergence of Φ is the scalar The integral form of the divergence operator ∇ in 3D relates to an arbitrary closed surface that is the boundary of a spatial region R that includes a point of interest, so that the surface is denoted ∂R , in which case the divergence of Φ at that point is the limit of a surface integral, as where '⋅' denotes the vector dot-product, n is the outward unit normal vector at the surface element ds and |R | is the volume of the region. For our EM problem, for a certain net atomic flux leaving the region, we get a corresponding rate of reduction in the material in that region, which can be expressed in terms of the atomic concentration in that region as ∂C ∂t Based on Eq. (2), we have ∂C/∂t = − (C 0 /B)(∂σ/∂t), so that ∂σ ∂t In general, the divergence can be computed using either the differential or integral form. In the 1D case, the differential form is convenient, so that ∇⋅Φ = ∂ϕ ∂x (9) where ϕ is the scalar component of Φ in the 1D x-direction of interest, leading to the continuity equation in the form used by Korhonen, ∂σ ∂t = BΩ ∂ϕ ∂x . (10)

The full model
Putting together the two Eqs. (4) and (10) results in Korhonen's equation which, assuming the diffusivity is constant along the length of the line, as is typically assumed in the field, takes the more commonly used form ∂σ ∂t According to [8], "The Korhonen model has been successfully used to explain a wide range of experimental behavior," and they cite [9][10][11][12] for this. Throughout this paper, we will assume that the applied currents, and therefore the currents in all branches, are constant over time. For a uniform line under constant current, the model takes the simpler form ∂σ ∂t which is basically the heat equation.

Extended Korhonen model (EKM)
The Korhonen model is a 1D model that applies at interior points of metal lines, but not at junctions. Before one can extend that model to multi-line interconnect trees, a model is needed for junctions that can be combined with the 1D model of the metal lines. This has been done in the so-called extended Korhonen model (EKM) [4]. The method assumes that a junction is a zero-volume point, and enforces mass conservation at that point so that the incoming and outgoing mass transport rates balance out. Then, and even though it's not strictly applicable, Korhonen's equation is applied at the junction anyway, with the hope that the error would be small because the equation is applicable nearby, in the lines. The method performs well in practice in spite of this approximation. The scheme also requires tracking the stress at "ghost points" around the junctions, which disturbs the structure of the system matrix in ways that are not desirable for numerical work.
Nevertheless, the artificial step of applying Korhonen's equation right at the junction is a shortcoming of this approach. Furthermore, while junctions are often small in size, they are sometimes designed to have considerable volume so as to act as reservoirs of atoms that can fill nearby micro-voids, if they form. EKM cannot handle such layout features, so this is another shortcoming of the method. In this paper, we will propose a new junction model that eliminates both these shortcomings, which will also turn out to be useful for development of an equivalent circuit.
In EKM, the junction model is used to stitch up and combine the various instances of Korhonen's equation applied to every line into a set of dynamic (ODE) equations for the whole metal structure under consideration. A preliminary step is to create reference directions for all the branches, which are generated by a breadth-first search graph traversal algorithm. These directions can be arbitrary, but must remain fixed throughout the analysis.
In general, voids can nucleate in many different places in interconnect trees, very often at junctions, but also internal to metal lines. The data seem to suggest [8] that internal voids appear mostly because of pre-existing micro-voids internal to the line as a result of the fabrication process. The EKM framework [4] assumes that there are no pre-existing voids, and no internal voids ever, so that all voids are assumed to nucleate at junctions. Furthermore, the treatment of voids in EKM depends on whether the void is at a degree-1 junction, which we will also call a terminal junction or terminal node of the tree, or at a junction with degree of more than 1 (non-terminal junctions). If the void forms at a non-terminal junction, such as in Fig. 3, the tree is broken up into multiple sub-trees [13], as in Fig. 4, by (i) tearing out that node from the tree graph, which creates multiple disconnected sub-trees, (ii) creating new copies of the original node for each of these sub-trees and reinstating in each of them the edge that was originally connected to that node, and (iii) applying a boundary condition at each of these duplicate nodes that represents the presence of a void there. As a result, we end up with a situation where multiple voids have been created at these (new) terminal junctions in the separate sub-trees. Therefore, for void analysis in the context of EKM, voids (once they have formed and the tree has been partitioned) will exist only at terminal junctions of interconnect trees. This also means that a void belongs to only a single metal line in a given sub-tree. It also means that no single line can have more than one void at a time. We adopt this tree partitioning approach of the EKM framework for treatment of voids in this paper. We will also expand on that approach, as described below.

Void analysis
The treatment of voids in EKM has one shortcoming which we will briefly describe then replace with better void models from the literature, as we will review below.
Voids can be partial or full (blocking the full cross-section of the line), they can increase in size until they reach a saturated size, or they can reverse course and decrease in size (healing) if the current changes direction, and they can move along the length of a metal line [8,21]. A simplified model that is suitable for 1D analysis was given in [18] (pg. 53), by making use of a phase field model [22], and this was the basis for the EKM work [4]. However, the implementation of the void model in [4] was based on the assumption that the void growth phase is very short relative to the nucleation phase, and so only the saturated void volume was part of the implementation, with no tracking of the void growth over time. More recently, it has become clear that void growth rates nowadays are such that the time to reach void saturation is not negligible [2], and so the 1D void model as implemented in the EKM framework should be improved so as to incorporate void growth tracking. We will see how this is done, leading to an equivalent circuit model that can track void growth.
According to the 1D void model in [18], the stress profile along the length of the line near the void is as shown in Fig. 5 for a void at the near end of the line (x ≈ 0), and in Fig. 6 for a void at the far end (x ≈ l). Right after the void is formed, with the stress having reached σ crit at that point in order for the void to nucleate, the stress right at the void surface must clearly be zero, as there is no confinement there at all. But the stress inside the metal very close to the void surface cannot change instantaneously and so must still be at σ crit initially (right after the void has formed) and will then decrease over time as the void grows and relieves the tensile stress. The distance between these two points is called the effective thickness of the void interface [18] and is shown in the two figures as the distance between the two dashed vertical lines, denoted by δ s , and described [18] as being "infinitely small in comparison with all other involved lengths." One can think of this as the thickness of a very thin "skin layer" -it was set at δ s = 1 nm in [18], and we will use that value as well. The inner boundary of this skin layer is denoted by the x s (t) position in the figures, and the stress at that point is denoted σ The length of the void, denoted l vd (t) is indicated in the figures as the distance to the mid-point of the skin layer, but that is debatable as this is not a precise distance in the actual 3D void. For simplicity, in our work, it is estimated as l vd (t) ≈ x s (t), for a void at the near end, and l vd (t) ≈ l − x s (t) at the far end. The model assumes that the stress is linear throughout the skin layer, so that the stress gradient in the skin layer is given by where the '+' is for a void at the near end and the '− ' at the far end. This gradient increases the magnitude of the atomic flux flowing out of the void into the bulk of the metal (which is the reason for void growth), and that flux can now be written, based on Eq. (4), as

Void length
The rate of change of the void size is proportional to the rate at which material is transported off the void surface and into the bulk of the    metal. This leads to [18] (pg. 56), Thus, for a void at the near end of line, we get while for a void at the far end of the line, For terminal voids, and because the void length is very small relative to the line length (typically under 2%), the common approach in the field is to numerically solve for the stress in the line while (i) assuming the length of the metal in the line remains fixed at its original length, and , depending on the location of the void, and this constitutes the boundary condition at that junction, which is needed for the solution.

Resistance change
Once a void has formed, the metal line resistance is impacted and the new resistance value depends on the void length, so that it is timedependent and grows towards a saturated value. The easiest way to think of this is to consider the voided line to be the series connection of two wires. The first, corresponding to the un-voided length of wire l − l vd (t), has a resistance that depends only on the Copper (ignoring the liner, whose resistivity is much higher), so it is given by ρ(l − l vd )/wh, where ρ, w and h are as defined earlier. The second is due to the void, where electric current flows only through the liner, so that the resistance is ρ lin l vd /h lin (w + 2h), where ρ lin is the resistivity of the liner material, h lin is the thickness of the liner and (w + 2h) is the sum of line dimensions around the cross-section of the wire, as in Fig. 3 of [18]. As a result, the resistance of the voided line is where R 0 = ρl/wh is the resistance of the original (un-voided) line.
Typical values for the metal parameters are given in [18] as h = 120nm,

New junction model
We will introduce a new junction model that overcomes the shortcomings of the model used in EKM as described in Section 2.4. Recall the divergence theorem, also known as the Gauss-Ostrogradsky theorem, which states that the surface integral of the flux through a closed surface is equal to the volume integral of the divergence of that flux over the region inside that surface. We normalize this by the volume, to state the theorem in this form So, the average flux divergence over a given finite volume can be found by tallying up the net outgoing flux through the surface bounding it, normalized by the volume. Based on this, let R be the 3D spatial region of some junction, and V be the volume of that region, then using Eq. (8) with σ R (x, y, z; t) as the stress in the 3D region, we can write so that the average value of σ R throughout the junction volume is proportional to the surface integral of the fluxes flowing across the junction boundary. It is reasonable to choose this average value of σ R over the junction volume to be the value of σ(t) for this junction in the context of a 1D model, to be used along with the Korhonen model. Indeed, 1D models are often constructed based on average quantities. Therefore, our new junction model will be based on the requirement, Because the flux across the boundary is only non-zero where a branch connects to the junction, and assuming that the fluxes in the lines are perpendicular to the boundary where they cross the junction region boundary, which is a reasonable simplification, then Eq. (22) provides where a k and l k are the cross-sectional area and the length of branch k, respectively, while B out (B in ) is the set of branches incident on this junction whose reference direction is away from (towards) the junction. Note that if we set V = 0, this equation reduces to that used in EKM [4]. In many cases, the junction volume is negligible compared to the volume of metal in the lines, so that this setting becomes reasonable. However, we have hereby given a more rigorous justification for this more general 1D junction model that allows for junctions that have significant volume and does not artificially apply Korhonen's equation at junctions.

Distributed equivalent circuit
This section contains the bulk of our contribution. We will develop and describe a number of equivalent circuits that make use of distributed RC transmission lines, whose simulated solution can give us the desired stresses for the whole interconnect tree, including junctions as well as interior points of metal lines, covering both the void nucleation and void growth phases.

Equivalent circuit for the nucleation phase
Consider first an interconnect tree without any voids. We start by identifying the connection between the equations for stress in a metal line and those for voltage in a distributed RC transmission line. Recall that Korhonen's equation for a uniform metal line with fixed current density Eq. (13) is given by where τ l = Δ k b T/(BΩD) which, it's easy to verify using Table 2, has units of s/m 2 . This special case of the heat equation is also called the diffusion equation and it happens to also be the governing equation for a distributed RC transmission line, i.e., a transmission line with no series inductance and no shunt conductance. This can be seen by starting from first principles and considering an infinitesimal section of the line, as in [23]. Another, indirect approach starts with the general transmission line equations 2 [24] (pg. 438), ∂v ∂x = − ri − l ∂i ∂t (25) ∂i ∂x = − gv − c ∂v ∂t (26) where v(x, t) and i(x, t) are the voltage and current along the line (i is positive in the positive x direction), r and l are the series resistance and inductance per unit length, and c and g are the shunt capacitance and conductance per unit length. Then, we set l and g to zero, corresponding to a simple RC line with no inductance and no leakage, to get ∂v ∂x Differentiating the first equation with respect to x and substituting from the second equation provides where rc has units of (ohm/m)(F/m) = s/m 2 , just like τ l in Eq. (24). The two Eqs. (24) and (28) are virtually identical, so that it should be possible to use a transmission line to model and simulate the stress in a metal line, as we will see below.

Line equivalent circuit
For a uniform metal line, let this be the kth line in the interconnect tree, we now define the following per-unit-length parameters, where ψ= Δ 0.01 C 2 /m 3 is a scale factor that provides desirable unit conversion and helps avoid running into numerical errors due to very small or very large numbers. It is easy to verify, using Table 2, that r k has units of resistance per unit length (ohm/m), c k has units of capacitance per unit length (F/m) and which is the same τ l in Eq. (24).
Next, let ξ= Δ 1V/MPa be another scale factor, which we use to define as the voltage associated with the stress σ(x, t) at every point in the interconnect tree, then multiply both sides of Eq. (24) by ξ and use Eq. (30) to get Thus, indeed the governing equation for stress in a uniform metal line is identical to the governing equation for voltage in a distributed RC transmission line with the parameters r k and c k defined in Eq. (29).
Formally, with τ l = τ k from Eq. (30), Such a transmission line can be either solved analytically or simulated numerically in simulators like HSPICE® that comprehend transmission line models. The resulting voltages can be directly translated to stress by dividing by ξ.
It remains to consider the flux dynamics and whether they can be equally well represented in an equivalent circuit. To that end, notice that the first equation in Eq. (27) translates to ∂σ/∂x = − r k i/ξ from which, using Eq. (4), we have where I k = Δ a k j k is the current in the kth branch of the original interconnect tree corresponding to the current density j k . To simplify the notation, we will define the following current notation for the kth branch, where I Sk has units of current (A) with the same sign as j k , and we can write Thus, flux can be easily found from current in the transmission line, via an affine relationship Eq. (36) that includes the constant I Sk .
The above relations for flux and stress motivate the definition of the line equivalent circuit shown in Fig. 7, which includes both the RC transmission line, with parameters r k and c k , as well as an ideal current source carrying the current I S defined above. At any value of x, the stress in the metal line corresponds to the voltage in the transmission line at that point, while the flux is the (scaled) sum of the current in the transmission line at that x position and the current in the current source.

Junction equivalent circuit
Using Eq. (36) for the atomic flux, we can write the intensity of the atomic transport (number of atoms going through a cross-section of a metal line per second) as The subscript k will be dropped when there's no risk of confusion, such as for a single isolated line. Combining this with the junction model (23), we get and multiplying both sides by ξψ converts σ to v and gives ( This is the motivation for defining what we will call the junction which, it's easy to verify using Table 2, indeed has units of capacitance (F). Thus, the stress-flux model (23) for a junction is equivalent to the voltage-current model which we would like to reproduce at the corresponding node in the equivalent circuit. Because this is basically the result of the applying Kirchhoff's current law (KCL) at a node with capacitance C J , then the equivalent circuit for a junction is simply a node with capacitance C J to ground, fed by currents from its connected branches in accordance with Eq. (41). Note that, as defined, the junction capacitance is proportional to the junction volume, so that it represents the capacity of the junction to hold a certain number of atoms at a given stress value. As we develop the model further, we will see that C J will often be set to zero, because it will be negligible in relation to the total line capacitances tied to that node. But, in general, it may be specified in the circuit description file whenever it's available, based on the junction volume extracted from the chip layout database.

Tree equivalent circuit
We will now combine the equivalent circuits for lines with the equivalent circuits for junctions to get an equivalent circuit for the whole interconnect tree. To start, for the case of a single isolated metal line, terminated by two degree-1 junctions, we will see that the equivalent circuit is as shown in Fig. 8, where the line equivalent circuit has been simply connected to the two junction equivalent circuits at its terminals. From basic circuit theory, this circuit can be redrawn as in Fig. 9 for simplicity, without any impact on the circuit currents and voltages. We can easily show that this circuit realizes the combination of the line Eq. (32) and the junction Eq. (41). All internal points of the line obviously satisfy Eq. (32), so that equivalence is established due to Eq. (33). As for the two junctions, apply KCL at the node at x = 0, i.e., at junction 1 with capacitance C J1 , to get C J1v (0) = − (i(0) +I S ), while for the node at x = l, i.e., at junction 2 with capacitance C J2 , get C J2v (l) = (i(l) +I S ), both of which satisfy Eq. (41) for the case of a single isolated line, so that equivalence is established due to Eq. (23) being equivalent to Eq. (41).
For interconnect trees consisting of multiple lines, multiple instances of the above line equivalent circuit can be connected in the same way as the actual metal lines are connected in the interconnect tree. Then, any junctions whose volumes are deemed significant would be represented by additional junction capacitances at the relevant line terminations, based on Eq. (40). An example is shown in Fig. 10 for the case of two metal lines connected at a dotted-I junction.
Furthermore, the equivalent circuits for multiple trees that form a whole power grid can be jointly simulated in order to find the stresses everywhere, in all trees. Very little extra work is needed for this extension, for simulation of the nucleation phase. Things get a bit more complicated when the void growth phase simulation is to be extended to whole power grids, as we will briefly describe at the end of the next section.

Equivalent circuit for the void growth phase
Consider next the case of an interconnect tree in which a void has formed, say at the far end of the line, at x = l. Recall that, as in the EKM framework, this means that this is a terminal junction of the tree, i.e., a junction of degree-1. We will define a few scaling factors and parameters that are needed for developing an equivalent circuit for the void end of the line. First let η= Δ 1V/nm, which we will use as a conversion factor in order to track void length by means of a node voltage in the equivalent circuit. Then, define the following conductance G s and capacitance C v , where δ s is the void-metal interface thickness defined earlier. It is easy to verify, using Tables 1 and 2, that G s and C v indeed have units of conductance and capacitance, respectively. For a void at the far end of the line, consider now the circuit in Fig. 11 which includes a voltage-controlled current source (VCCS) carrying the current G s v(l). As we will see, the new (relative to Fig. 7) circuit portion on the far right will provide (through G s ) enforcement of the void boundary condition (14) and will give us the ability (through C v ) to track the void length evolution over time, based on Eq. (16), i.e., l vd ′ (t) = − Ωϕ(l). Recall, using Eq. (37) that Because KCL at the node marked v(l) provides i(l) = G s v(l), then Eqs. (43) and (4) lead to Using the expressions for G s Eq. (42) and I S Eq. (35), we get so that the boundary condition (14) is enforced. As for the void length, KCL at the capacitor node provides where we have again used the fact that i(l) = G s v(l). Substituting for C v Fig. 7. The equivalent circuit for a metal line. Fig. 8. The distributed equivalent circuit for the case of a single isolated line, redrawn.
so that we can indeed track the void length by monitoring the voltage v v in the equivalent circuit, as This equivalent circuit can be simplified somewhat using simple circuit transformations to get the circuit shown in Fig. 12, and then again into the final form in Fig. 13. Finally, in case the void has nucleated at the near end of the line, i.e., at x = 0, then the same analysis provides the equivalent circuit diagram in Fig. 14.
Once these equivalent circuits for individual lines are used in the context of a simulation of a multi-line tree or for multiple trees in a larger power grid, the impact of void growth on line resistance, via Eq. (19), and therefore on branch currents must be taken into account. This can be done by customization of the overall circuit simulation, in order to recompute the currents (j k and therefore I Sk ) for every line whenever a significant change of resistance has occurred. These details are beyond the scope of this paper.

Lumped equivalent circuit
As mentioned earlier, transmission lines can be solved analytically, but that is just as complex as solving the original Korhonen's equation analytically. Transmission lines can also be simulated using existing tools like HSPICE®, but this can be expensive, and the detailed highfrequency analysis performed by the simulator is not needed for the slow EM system. Instead, and as is often done in practice, RC transmission lines can be approximated, with very good accuracy, by a lumped RC line, as shown in Fig. 15. The lumped line is composed of N segments, each of which is a π-RC approximation of a short segment of the original line. A segment of the line of length δ k , corresponding to a total resistance of R k = δ k r k and total capacitance of C k = δ k c k , is approximated by the combination of a single (lumped) series resistor R k and two (lumped) shunt capacitors to ground of value C k /2 each. Internal nodes of the line end up with a capacitor C k each, while the two line terminal nodes get only C k /2 each. Thus, if line k has length l k , then δ k = l k /N and the line is characterized by the lumped-element values Note again that capacitance is representative of the volume of a metal region (in this case the segment volume a k δ k ) and represents the capacity of that region to contain metal atoms. The HSPICE® user guide [25] (pg. 203) recommends as default a number of N = 20 segments for lumped element RC approximations of transmission lines. Beyond this, it says, one gets negligible improvement for the increased simulation time. It is interesting that this guidance is independent of the line length. We have found N = 20 to work very well in practice, as also observed in [13] where an even smaller N of 16 was used. With this, the equivalent circuit for a single isolated metal branch in the nucleation phase becomes as shown in Fig. 16. Similarly, the corresponding circuits for lines with voids are shown in Figs. 17 and 18. These lumped equivalent circuits can then be interconnected in the same way as the original metal branches, possibly combined with junction capacitors (at non-voided junctions), in order to represent the whole interconnect tree.

Validation
We will give a number of simulation results that demonstrate the validity of this approach, covering both the void nucleation phase and the void growth phase.

Nucleation phase
We will describe the implementation of the above approach and give comparisons to exact solutions from [3,17]. We have developed a generic SPICE sub-circuit, an example of which is shown in Fig. 21. The figure shows a sub-circuit type consisting of 20 π-RC segments, called PIRC20, which makes use of another sub-circuit for a single π-RC link, called pirc_seg and included at the bottom of Fig. 21. The PIRC20 subcircuit captures the contributions of a single metal line to the equivalent circuit of an interconnect tree, based on the equivalent circuit shown in Fig. 16. This sub-circuit can be used to construct the equivalent circuit for any given interconnect tree, by connecting multiple instances of the PIRC20 sub-circuit in the same way as the metal branches are connected. Additional capacitors to ground may then be added for any junctions whose volumes are deemed to be significant. An example is shown in the SPICE circuit description given in Fig. 19, for the 4-line interconnect tree with a plus-junction shown in Fig. 20. Circuit simulation can then be     The exact solution can be found analytically in certain special cases. For the simple case of a single isolated line, the exact solution is available from both [3,17], as             where s n = (2n + 1)π, and κ = DBΩ/(k b T), which has units of m 2 /s, while l is the line length as usual. This specific form of the solution follows the expression 3 used in [17]. For the case of a Copper line of width 1 μm, height 1 μm and length 250 μm, discretized into 20 segments as in PIRC20, at T = 400 K and carrying a current density of 1 × 10 9 A/m 2 , the circuit parameters turn out to be (in round numbers), R k = 528 kΩ, C k = 251 mF and I S = 69 μA, based on the physical parameters in Table 1. Setting the junction volumes to zero, the comparison of the exact solution at the cathode junction to that from the SPICE solution using PIRC20 is shown in Fig. 22, and the agreement is excellent. This test was repeated for line lengths of 150 μm, 80 μm, 40 μm and 10 μm, and the agreement is excellent in every case, as in the 10 μm case shown in Fig. 23.
Throughout the following, all test runs will be based on zero junction volumes.
We also tested the case of an interconnect tree composed of two identical lines in a dotted-I arrangement, against our implementation of the exact solution based on the analysis 4 in [17]. Each line is L = 250 μm long and 1 μm wide, with both reference directions to the right, carrying j 1 = 1 × 10 9 A/m 2 in line 1 and j 2 = 3 × 10 9 A/m 2 in line 2, at 400 K. The comparison at the three junctions is shown in Fig. 24, showing excellent  π-RC segments. 3 The sign difference (after the σ 0 term) relative to Eq. (3) in [17] is because they've assumed that the reference direction for the current is opposite to that of the x distance variable, as can be seen by comparing their Eq. (1) to our Eq. (12).
accuracy. We also tested the case of unequal widths, keeping w 1 at 1 μm while setting w 2 = 2w 1 and keeping everything else the same. In this case, in the absence of an exact solution, we compared to the EKM implementation [4] and the results are again excellent as shown in              Another interesting case is the 3-line T-junction arrangement, given in the SPICE description in Fig. 28, for which the comparison to EKM is shown in Fig. 29. Finally, for the plus-junction described in the netlist in Fig. 19 and shown in the layout diagram in Fig. 20, the comparison to EKM is shown in Fig. 30.

Void growth phase
We will describe the implementation of the above approach and give comparisons to exact solutions from [26]. We have developed a generic SPICE sub-circuit, an example of which is shown in Fig. 31 for a void at x = 0 and Fig. 32 for a void at x = l. The figures show two sub-circuit types consisting of 20 π-RC segments each, called PIRC20V0 and PIRC20VN, which make use of another sub-circuit for a single π-RC link, called pirc_seg and included in the figures. These sub-circuits capture the contributions of a single (voided) metal line to the equivalent circuit of an interconnect tree, based on the equivalent circuits shown in Figs. 17 and 18. One can construct the equivalent circuit for any given interconnect tree by connecting multiple instances of these PIRC20V0 and PIRC20VN sub-circuits, along with the PIRC20 sub-circuits shown earlier for any un-voided line, in the same way as the metal branches are connected. Additional capacitors to ground may then be added for any junctions whose volumes are deemed to be significant. Circuit simulation can then be applied to the circuit to find all voltages, therefore all stresses and void lengths.
We will explore the transient stress response in the line as well as the void growth evolution over time. When possible, we will compare to exact solutions, which can be found analytically in certain special cases [26]. For the simple case of a single isolated line, and assuming that a void has just nucleated at x = 0 at time zero, the condition σ(0, t) = 0, ∀t is enforced as a boundary condition in [26]. They also assume that the initial stress throughout the line is negligible, i.e., σ(x, 0) = 0, ∀ x, an assumption that is presumably made in order to be able to get an analytical solution. With this setup, the exact solution is given by 5 where c n = (2n − 1)π/2, with κ = DBΩ/(k b T) as before, while l is the line length. For our test case, we specified a Copper line of width 1 μm and length 250 μm, discretized into 20 segments as in PIRC20V0, at T = 400 K and carrying a current density of − 2 × 10 9 A/m 2 , as shown in the SPICE circuit description in Fig. 33. Notice that there are additional physical constants (DELTA for δ s and ETA for η) that have to be specified up-front. Also, a detailed. IC statement is needed to set up the correct initial profile of stress in the line at the time of voiding, which in this case was set to zero so as to allow a fair comparison with [26]. Having set the junction volume (at x = l) to zero, we compared the time evolution of σ(x, t) for the SPICE solution v.s. the exact solution, at every node of the RC chain, i.e., at x = 0, δ k , 2δ k , …, l k , where δ k = l k /20. The comparison is shown in Fig. 34, where the rectangular area in the top-left corner is expanded and shown in Fig. 35.
The SPICE solution is shown as the solid red curves while the exact solution is the dashed black curves, and the agreement is clearly excellent. Throughout the following, all test runs will be based on zero junction volumes. Using our approach, we can handle more realistic cases, such as when the initial stress in the line is non-zero, and where a non-zero residual thermal stress can be included. For example, for the same 250μm line, we first simulated the line before voiding, saved the state at the time of voiding and then applied that state as the initial stress profile during the subsequent simulation of the voided line. The result is shown in Fig. 36, where we have shown the profile of the stress along the length of the line at different points in time. Notice that the stress at x = 0 drops very quickly from σ(0, 0) = 500 MPa to a constant σ(0, t) = 0 over a period of seconds, while the rest of the line (at 5% of the line length away from the void, and beyond) takes months or years to reach its steady state. This fast transient behavior of the stress at x = 0 is shown in Fig. 37 and testifies to the very high initial flux at the void surface due to the high initial stress gradient σ crit /δ s there. Fig. 38 shows the void growth for a single line with δ s set to 10 nm, 1 nm and 0.1 nm. It is clear that the solution is not sensitive to the choice of δ s in this range, as one would like to see in fact. We have found that this property of this model [18] remains valid over the wide range of 1pm ≤ δ s ≤ 1μm.
Looking next at the void length over time, we will compare to the exact solution provided in [26], as   5 The sign difference (before the first term) relative to Eq. (8) in [26] is because they've assumed that the reference direction for the current is opposite to that of the x distance variable, as can be seen by comparing their Eq. (1) to our Eq. (12).
where l vd,sat = q*ρ|j|l 2 /(2BΩ) from [26]. The comparison to SPICE using PIRC20V0 for the same 250 μm line with j = − 2E9A/m 2 is shown in Fig. 39. Another comparison is also shown for a 100μm line carrying j = 5E9A/m 2 , with the void at the far end of the line, so we use PIRC20VN, and the results are in Fig. 40. The agreement is excellent in both cases. Finally, the resistance increase over time for these two lines is shown in Figs. 41 and 42, based on Eq. (19).

Reduced line
Given the PIRC-20 lumped model, we can apply the PACT [27] model-order reduction approach, as described in [28] to get the generic scalable compact model shown in Fig. 43, which we use to replace PIRC20, resulting in what we call the PACTN4 sub-circuit, given in Fig. 44. For the same 250 μm line considered in Fig. 22, the comparison to the PACTN4 reduced line is shown in Fig. 45. Comparisons are also given for the dotted-I arrangement, with equal widths, in Fig. 46, showing some degradation at the fast transients, but otherwise quite acceptable. The major advantage of using the PACT transformation is to reduce the number of nodes in the system, thus allowing one to handle very large interconnect trees. If the interconnect tree has n junctions and m branches (recall that in a tree, m = n − 1) then using PIRC20 would generate an equivalent circuit with (20n − 19) nodes while, for PACTN4 the resulting circuit has only (2n − 1) nodes, which is a ≈10× reduction. In other words, this reduction allows the analysis of 10× larger interconnect trees. This clearly shows the benefit of using the equivalent circuits approach! It would have been very difficult to imagine, looking only at the stress equations, that this kind of reduction of the problem size is at all possible.

Conclusion
We have demonstrated a deep connection between the EM-induced stress-flux dynamics in metal lines and the voltage-current dynamics in RC circuits. Effectively, the atomic transport in an interconnect network behaves in the same way as the electronic current in the metal lines; they obey the same circuit laws. We have done this for both the void nucleation phase and the void growth phase, and provided sample SPICE sub-circuits that can be used to generate a large RC network for any given metal interconnect network. Simulation results demonstrate excellent agreement between the two "worlds" of stress and voltage.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.