**4. Electrothermal Simulation Approach**

We chose to resort to a *circuit-based* approach [7,11,26,27], which is a good trade-off between computational burden and accuracy. Such an approach makes use of the *thermal equivalent of the Ohm's law (TEOL)* and can be summarized as follows.


voltage *VTH* and of the current factor *K*, the bias- and temperature-dependent drift resistance *Rdrift*, as well as the bias- and temperature-dependent II mechanism. The TEOL is adopted, namely, the temperature rise over ambient Δ *T* = *T* − *T*0 is actually a voltage, while the dissipated power *PD* is treated as a current; this allows (i) enabling the temperature sensitivity of the key physical parameters, and (ii) describing the power-temperature feedback with an electrical network. Besides the standard electrical terminals (gate, drain, source/body), the cell subcircuit is also equipped with an input node carrying the "voltage" Δ *T* (provided by the thermal feedback block introduced below) and with an output node o ffering the "current" *PD* (to be fed to the thermal feedback block).


$$
\Delta T = T\_0 \cdot \left[ m\_k + (1 - m\_k) \cdot \frac{\Delta T\_{lin} + T\_0}{T\_0} \right]^{\frac{1}{1 - m\_k}} - T\_0 \tag{18}
$$

The calibration of the (positive) parameter *mk* will be detailed in Section 4.2. Nonlinear voltage-controlled voltage sources emulating (18) are applied to the *N* temperature rise nodes of the equivalent thermal network; the resulting circuit is referred to as a thermal feedback block. It is worth noting that the FANTASTIC-based approach improves the strategy exploited in [7,26], where the thermal feedback block was based on Foster networks extracted in a rather long pre-processing stage from *N* onerous transient COMSOL simulations of the DUT (performed by activating one heat source at a time), and the thermal coupling between horizontally-far heat sources was roughly described or even neglected.

• The cell subcircuits are then connected to the thermal feedback block in a commercial circuit simulation tool (like PSPICE, LTSPICE, Eldo, ADS, SIMetrix); as a result, the whole domain under test, composed by the DUT soldered on a DBC substrate, is transformed into a *purely-electrical macrocircuit*, which suitably accounts for ET e ffects: the temperature, and thus the temperature-sensitive parameters, are allowed to vary during the simulation run. The solution of this macrocircuit under both static and dynamic conditions is demanded to the powerful and robust engine of the circuit simulation tool, with very low computational e ffort and minimized occurrence of convergence issues compared to other numerical methods.

### *4.1. SPICE Subcircuit and DUT Discretization into Cells*

A sketch of the SPICE-compatible subcircuit for the transistor cell is represented in Figure 9, where the standard MOSFET instance at reference temperature *T*0 is indicated with *Mn*. The additional input node feeding the "voltage" Δ*T* and the output node providing the "current" *PD* are also represented.

**Figure 9.** Sketch of the SPICE-compatible subcircuit for the transistor cell.

The negative temperature coefficient of the threshold voltage *VTH*(*T*) described by (3) is accounted for by using the following approach. The voltage source **A** in series with the gate adds *VTH*(*T*0) − *VTH*(*T*) to *VG*, *VTH*(*T*) being given by (3), thus biasing *Mn* with *VG'* = VG + [*VTH*(*T*0) − *VTH*(*T*)]; as a result, the *e*ff*ective* overdrive voltage becomes

$$V\_{G'S} - V\_{TH}(T\_0) = V\_{GS} + \left[V\_{TH}(T\_0) - V\_{TH}(T)\right] - V\_{TH}(T\_0) = V\_{GS} - V\_{TH}(T) \tag{19}$$

and *Mn* conducts the current *ID*(*Mn*) given by

$$\begin{aligned} I\_D(M\_\text{lt}) &= K(T\_0) \cdot \left\{ 2[V\_{GS} - V\_{TH}(T)]V\_{DS\text{ch}} - V\_{DS\text{ch}}^2 \right\} & \quad V\_{DS\text{ch}} &< V\_{GS} - V\_{TH}(T) \\ I\_D(M\_\text{lt}) &= K(T\_0) \cdot \left[ V\_{GS} - V\_{TH}(T) \right]^2 & \quad V\_{DS\text{ch}} \ge V\_{GS} - V\_{TH}(T) \end{aligned} \tag{20}$$

The temperature dependence of the current factor *K*(*T*) expressed by (4) with (5) is emulated by using the current source **B** to derive the current

$$I\_{\mu} = I\_{\text{D}}(M\_{\text{ll}}) \cdot \left[ \left( \frac{T}{T\_0} \right)^{-m(T)} - 1 \right] \tag{21}$$

Consequently, the II-unaffected current *IDnoII* is obtained as

$$I\_{DudI} = I\_D(M\_n) + I\_\mu = I\_D(M\_n) + I\_D(M\_n) \cdot \left[ \left(\frac{T}{T\_0}\right)^{-m(T)} - 1 \right] = I\_D(M\_n) \cdot \left(\frac{T}{T\_0}\right)^{-m(T)}\tag{22}$$

where *m*(*T*) is given by (5).

II effects are activated by adding an avalanche current *IDII* = ξ·(*Ileak* + *IDnoII*) to *IDnoII* through the current source **C**, ξ = *M* − 1 being bias- and temperature-dependent according to (6). Hence, the drain current *ID* flowing through the drift resistance *Rdrift* is given by (9).

The bias- and temperature-dependent *Rdrift* described by (10) with (11) is taken into account by making use of source **D**, which imposes the voltage drop

$$V\_{drift} = I\_D \cdot R\_{drift} \Big( V\_{GS\prime} \cdot V\_{drift\prime} \, T \Big) . \tag{23}$$

Only half of the DUT was represented and simulated by exploiting its inherent symmetry. The e ffective active region (≈5 mm2) was partitioned into *N* = 79 cells, each with a 250 × 250 μm<sup>2</sup> area (the procedure leading to the choice of the *N* value will be detailed in Section 4.4). As a consequence, compared to those reported in Table 1, the values of the area-dependent parameters were scaled by dividing the current factor and the capacitances by 2 × *N*, and multiplying the resistances by 2 × *N*. The popular commercially-available OrCAD PSPICE [30] was chosen to perform circuit simulations. The main schematic (resembling the actual layout) is represented in Figure 10, along with the corresponding cell numbering and a detail of the connections between adjacent cells. The drain current of the resulting macrocircuit (an outcome of the ET simulation) has to be multiplied by 2.

**Figure 10.** Discretization of the e ffective active area of half DUT. **Left**: PSPICE subcircuit with hierarchical blocks corresponding to the transistor cells; **center**: cell numbering, with evidenced cells of interest for the ET analysis presented in Section 5; **right**: connections between some adjacent cells.

### *4.2. FEM Representation of the Component under Test*

Exhaustive details on the geometry and materials of the DUT (the CPMF-1200-S080B VDMOS) were taken from the datasheet and from reports available on a reverse-engineering website. The DUT was assumed to be soldered on a DBC substrate by means of a 50 μm-thick tin-platinum alloy (SnPt) layer ensuring both mechanical joint and electrical connection with the drain [31]. Figures 11 and 12 show the top view and cross-section of the assembly, which enjoys the same symmetry of the DUT. The soldering process can be automated by means of thin SnPt films, which are preformed and match the size of the die; in addition, alignment masks are typically exploited to ease their positioning. The DBC is composed by two copper (Cu) sheets (namely, bottom and top plates) with a ceramic layer placed in between them. Such a layer provides dielectric insulation to the assembly and is realized with alumina (Al2O3) for the proposed case study; however, alternative materials like silicon nitride (Si3N4) [32] and aluminum nitride (AlN) [33] are also adopted in the DBC manufacturing process. The drain contact of the DUT is soldered on the top plate, while the bottom plate ensures a good thermal interface with the 3 mm-thick Cu baseplate.

**Figure 11.** Sketch of the top view (not to scale) of the domain to be electrothermally simulated. All dimensions are expressed in mm.

**Figure 12.** Sketch of the cross-section (not to scale) of the domain to be electrothermally simulated. All dimensions are expressed in μm.

The 3D domain was represented in the environment of the commercial FEM-based COMSOL Multiphysics software package [34] by means of the *in-house* routine detailed in [35], which allows *automatically* building an extremely accurate geometry (Figure 13) and performing a smart selective optimization of the tetrahedral mesh (Figure 14); this process conveniently avoids a painstakingly long and prone-to-errors manual procedure for geometry/mesh construction. The number of elements (tetrahedra) and degrees of freedom (DoFs) of the resulting grid are 3.8 × 10<sup>5</sup> and 5.2 × 105, respectively. Figure 14 plainly illustrates that the mesh is highly fine over the die, while becoming gradually coarser by moving far away from the active region.

**Figure 13.** 3D representation of the geometry of the domain under test in the COMSOL Multiphysics environment (*draw mode*): (**a**) whole structure and (**b**) magnification of the 4H-SiC VDMOS die.

**Figure 14.** 3D tetrahedral mesh of the domain under test in the COMSOL Multiphysics environment (*mesh mode*): (**a**) whole structure and (**b**) magnification of the 4H-SiC VDMOS die.

As previously mentioned, all transistor cells (Figure 10) were individually associated to heat sources located over the top surface of the 4H-SiC DUT. In this representation, the heat is assumed to be dissipated over the whole effective active area, and not over the tens of thousands of individual channels; conveniently, this *unavoidable* approximation was demonstrated to be reasonable under dc and transient conditions (at least for not-too-short times) through a simulation analysis performed with TCAD Sentaurus and COMSOL. An isothermal boundary condition at *TB* = *T*0 was applied at the bottom of the assembly baseplate, whereas all other surfaces were assumed adiabatic (i.e., with zero outgoing heat flux). The material parameters adopted for the thermal simulations are listed in Table 2. Nonlinear thermal effects were accounted for by including the temperature dependences of the thermal conductivities described by

$$k(T) = k(T\_0) \cdot \left(\frac{T}{T\_0}\right)^{-\alpha} \tag{24}$$

$$k(T) = k(T\_0) - \beta \cdot (T - T\_0) \tag{25}$$


**Table 2.** Properties of the materials composing the assembly.

It must be noted that (24) applies to semiconductors and insulators, while (25) is followed by some metals.

The pre-processing calibration of parameter *mk* to be used for the Kirchhoff's transformation (18) was carried out with the following procedure. The domain under test was thermally simulated with COMSOL by varying the power *PD* dissipated by the whole effective active area (described with only one heat source) over a wide range (0.1 to 573.5 W for half device, the latter value leading to a peak temperature of 1400 K); the thermal conductivity dependences upon temperature were either activated (*nonlinear* conditions) or deactivated (*linear* conditions). The average temperature rise over that area was determined for both cases. Then the Kirchhoff's transformation was applied to the linear temperature rise, and *mk* was adjusted so as to obtain a good agreemen<sup>t</sup> between the nonlinear temperature rise computed by the transformation and the realistic one evaluated by COMSOL; it was obtained that *mk* = 0.785. The whole process lasted less than 1 h, as each nonlinear simulation was run in about 2 mins and 24 *PD* values were applied.

### *4.3. FANTASTIC-Based Derivation of the Equivalent Thermal Network*

The FANTASTIC tool, originally introduced by some of the authors in [6], where it was applied to merely-thermal 3D simulations of a state-of-the-art gallium-arsenide (GaAs) HBT, is based on the truncated balance (TRB)-based moment matching (MM) approach to model-order reduction (MOR) developed by one of the authors in [42].

In this tool, the dynamic heat conduction problem in a semiconductor device, assumed to be *linear*, is imported from either commercial (e.g., COMSOL) or open-source codes (e.g., SALOME SMESH), comprising: the mesh discretizing the geometry, as well as the definitions of heat sources, materials, and boundary conditions. Both hexahedral and tetrahedral meshes can be used. It is possible to define arbitrary heat capacity and tensorial thermal conductivity distributions. Neumann's, Dirichlet's, or Robin's boundary conditions can be applied. Superficial (i.e., indefinitely thin) and volumetric heat sources can be taken into account.

The FEM model of the thermal problem is then assembled by FANTASTIC. In particular, the mass matrix **M** and the sti ffness matrix **K** are constructed. High-order basis functions can be used: most commonly, as a trade-o ff between e fficiency and accuracy, tetrahedral meshes and second-order basis functions are considered. The *M* DoFs of the temperature rise distribution, forming the *M*-row vector <sup>ϑ</sup>(*t*), are solutions of the discretized linear dynamic heat conduction problem

$$\mathbf{M}\frac{d\boldsymbol{\Theta}}{dt}(t) + \mathbf{K}\boldsymbol{\Theta}(t) = \mathbf{q}(t) \tag{26}$$

in which the power density distribution vector **q**(*t*) takes the form

$$\mathbf{q}(t) = \mathbf{Q}\mathbf{P}(t) \tag{27}$$

where **P**(*t*) is an *N*-row vector with the powers dissipated by *N* independent heat sources, and **Q** is an *M* × *N* matrix, the *n*-th column of which is the power density distribution vector of the *n*-th source, with *n* = 1, ... , *N*. The port temperature rises of the *N* sources form the *N*-row column vector Δ **T**(*t*) given by [42]

$$
\Delta \mathbf{T} = \mathbf{Q}^T \boldsymbol{\theta}(t) \tag{28}
$$

As typical for MOR approaches, an *M* × *M*<sup>ˆ</sup> matrix **V** with *M*<sup>ˆ</sup> *M* is defined, which allows approximating ϑ(*t*) by means of a reduced number *M*<sup>ˆ</sup> of DoFs forming the *M*<sup>ˆ</sup> -vector <sup>ϑ</sup><sup>ˆ</sup>(*t*), so that

$$\mathfrak{d}(t) = \mathbf{V}\mathfrak{d}(t) \tag{29}$$

The **V** matrix is used to project the heat conduction discretized problem (26)–(28) with the Galerkin's method, thus deriving a DCTM in the form

$$
\hat{\mathbf{M}}\frac{d\hat{\mathbf{g}}}{dt}(t) + \hat{\mathbf{K}}\hat{\mathbf{g}}(t) = \hat{\mathbf{q}}(t) \tag{30}
$$

being

$$
\hat{\mathbf{q}}(t) = \hat{\mathbf{G}} \mathbf{P}(t) \tag{31}
$$

$$
\Delta \mathbf{T}(t) = \hat{\mathbf{G}}^{\mathrm{T}} \hat{\boldsymbol{\Phi}}(t) \tag{32}
$$

$$\text{in which}$$

$$
\hat{\mathbf{M}} = \mathbf{V}^T \mathbf{M} \mathbf{V} \tag{33}
$$

$$\mathbf{\hat{K}} = \mathbf{V}^{T}\mathbf{K}\mathbf{V} \tag{34}$$

are *M* ˆ -order matrices, and

$$
\hat{\mathbf{G}} = \mathbf{V}^T \mathbf{Q} \tag{35}
$$

is an *M* ˆ × *N* matrix. The **V** matrix is determined by the algorithm reported below.


At line 1 of the algorithm, the values <sup>σ</sup>*p*, with *p* = 1, ... , *Pn*, are automatically determined as a function of the desired relative error parameter ε as follows. First, the real positive quantities λ*n* < Λ*n* are properly estimated for the heat conduction problem with respect to the power impulse thermal response of the current *p*-th heat source. Next, the value *Pn* is determined as the smallest integer such that

$$4\exp\left(-P\_n\pi^2/\log(4/k')\right)\le\varepsilon\tag{36}$$

being *k'* = λ*n*/Λ*<sup>n</sup>*. It is then set

$$
\sigma\_p = \Lambda\_n \text{dn} \left( \frac{2p - 1}{2P\_n} K, k \right),
$$

with *p* = 1, ... , *Pn*, in which *K* is the complete elliptic integral of the first kind of modulus *k*, and dn is the homonymous elliptic function of modulus *k* = 1 − *k*2.

At line 3, the temperature response to the *n*-th heat source is solved in the complex frequency domain at the frequency value <sup>σ</sup>*p*. Thus, equation

$$(\sigma\_p \mathbf{M} + \mathbf{K})\boldsymbol{\Theta}\_p = \mathbf{Q} \mathbf{e}\_n \tag{37}$$

is solved for **<sup>Θ</sup>***p*, **e**n being an *N*-row vector having all zero elements but one at the *n*-th row. Since the complex frequency values are real positive, the coefficient matrices of the resulting linear systems are symmetric positive definite, and the most efficient multigrid iterative solvers can be used for their solution.

At line 4, the **<sup>V</sup>***p* matrix is achieved by appending to the columns of **<sup>V</sup>***p*−<sup>1</sup> a vector derived orthogonalizing **<sup>Θ</sup>***p* with respect to the columns of **<sup>V</sup>***p*−1.

At line 5, the DCTM C*p* is determined proceeding as in (30)–(32), matrix **V** in (33)–(35) being substituted by matrix **<sup>V</sup>***p*.

At line 2, the temperature response **<sup>Θ</sup>***p* in the complex frequency domain due to the *n*-th independent heat source is estimated at <sup>σ</sup>*p* using the DCTM C*p*−1. To this aim, equation

$$(\sigma\_p \mathbf{\hat{M}}\_{p-1} + \mathbf{\hat{K}}\_{p-1}) \mathbf{\hat{O}}\_p = \mathbf{\hat{Q}}\_{p-1} \mathbf{e}\_n \tag{38}$$

is solved for **Θ** ˆ *p*. This vector is used to determine the approximation **Θ** *p* = *Vp*−1**<sup>Θ</sup>** ˆ *p*. At line 3 of the algorithm, such estimation of **<sup>Θ</sup>***p* is used to speed up the solution of (37) by setting the initial estimation of the solution.

It is noted that for the *n*-th heat source, with *n* = 1, ... , *N*, the complex frequency values <sup>σ</sup>*p*, with *p* = 1, ... , *Pn*, are sorted in decreasing order. In this way, the linear systems introduced at line 3 of the algorithm are solved from the least to the most onerous ones in terms of CPU time. The computational burden is drastically reduced by this proper ordering, since the iterative solver can start from an estimate of the solution found in the previous step, which becomes increasingly accurate.

The DCTM achieved by this algorithm has dimension *M* ˆ = *P*1 + *P*2 + ... + *PN*. Let **Z**(*t*) and **Z**(*t*) be the *N*-th order power impulse thermal response [K/W] matrices of the discretized heat conduction problem (26)–(28) and of the DCTM, respectively. As shown in [43], it results in

$$\|\|\mathbf{Z}(t) - \hat{\mathbf{Z}}(t)\|\|\_{\mathcal{H}\_2} \le 2\varepsilon \|\mathbf{Z}(t)\|\_{\mathcal{H}\_2}$$

being

$$\|\|\mathbf{Z}(t)\|\|\_{\mathcal{H}\_2} = \sqrt{\int\_0^{+\infty} \text{tr}\Big(\mathbf{Z}^T(t)\mathbf{Z}(t)\Big)dt}.$$

the **Z**(*t*) Hankel norm. Similar results can be extended from the time to the frequency domain and can be stated for the whole spatial-temporal temperature distributions due to power impulses [43]. All these results provide a strong theoretical guarantee for the convergence of the method. As a consequence of (36), this convergence is very fast, being exponential with respect to the numbers *P* of <sup>σ</sup>*p*, with *p* = 1, ... , *Pn*. As a result, in practice accurate DCTMs can always be obtained with small space-state dimensions, *M* ˆ .

It is now observed that by solving at limited cost the generalized eigenvalue problem

$$\mathbf{U}^T \mathbf{M} \mathbf{U} = \mathbf{I}\_{\hat{M}}$$

$$\mathbf{U}^T \mathbf{K} \mathbf{U} = \mathbf{A}$$

having as unknown the *M* ˆ -order matrix **U**ˆ , in which **Λ**ˆ is an *M*<sup>ˆ</sup> -order diagonal matrix, and introducing the change of variables

$$
\Phi(t) = \mathbb{O}\xi(t).
$$

the DCTM equations (30)–(32) are transformed into

$$\frac{d\hat{\xi}}{dt}(t) + \mathbf{\hat{A}}\hat{\xi}(t) = \mathbf{f}^{\prime}\mathbf{P}(t) \tag{39}$$

$$
\Delta \mathbf{T}(t) = \mathbf{f}^{\mathrm{T}} \hat{\mathbf{f}}(t) \tag{40}
$$

**^**

where **Γ** ˆ is the *M* ˆ × *N* matrix **V** ˆ *<sup>T</sup>***G**<sup>ˆ</sup> . The temperature rise distribution is then reconstructed as

$$
\mathfrak{G}(t) = \Xi \hat{\mathfrak{E}}(t) \tag{41}
$$

**Ξ** = **VU** ˆ being an *M* × *M*<sup>ˆ</sup> matrix like **V**.

Equations (39) and (40) can be interpreted as the equations ruling the equivalent network sketched in Figure 15.

Such a network, which can in principle describe the thermal behavior of *any* electronic component, is particularly suited to be solved by means of nodal analysis in SPICE-like simulators, since all elements are voltage-controlled current sources, and thus the number of variables is limited to *M* ˆ . The topology is general and can be implemented into *any* circuit simulator.

**Figure 15.** DCTM equivalent thermal network (after Codecasa et al. [6]).
