Analysis of Electrothermal Effects in Devices and Arrays in InGaP/GaAs HBT Technology

: In this paper, the dc electrothermal behavior of InGaP/GaAs HBT test devices and arrays for power ampliﬁer output stages is extensively analyzed through an efﬁcient simulation approach. The approach relies on a full circuit representation of the domains, which accounts for electrothermal effects through the thermal equivalent of the Ohm’s law and can be solved in any commercial circuit simulator. In particular, the power-temperature feedback is described through an equivalent thermal network automatically obtained by (i) generating a realistic 3-D geometry/mesh of the domain in the environment of a numerical tool with the aid of an in - house routine; (ii) feeding the geometry/mesh to FANTASTIC, which extracts the network without performing simulations. Nonlinear thermal effects adversely affecting the behavior of devices/arrays at high temperatures are included through a calibrated Kirchhoff’s transformation. For the test devices, the thermally-induced distortion in I – V curves is explained, and the limits of the safe operating regions are identiﬁed for a wide range of bias conditions. A deep insight into the electrothermal behavior of the arrays is then provided, with particular emphasis on the detrimental nonuniform operation. Useful guidelines are offered to designers in terms of layout and choice of the ballasting strategy.


Introduction
Gallium arsenide (GaAs) heterojunction bipolar transistors (HBTs) are the dominant technology for handset power amplifier (PA) design by virtue of their power density, cut-off frequency, and efficiency [1]. However, designing rugged circuits with GaAs HBT devices requires extra care because of strong electrothermal (ET) effects arising from (i) low thermal conductivity of the substrate, (ii) lateral heat confinement by mesa isolation, and (iii) high operating currents. Thermal-aware design methodologies relying on suitable ET simulation tools are highly desired to alleviate or avoid performance and reliability degradation. Unfortunately, the choice of the simulation approach is challenging. Full 3-D ET device simulations based on the finite-element method (FEM), e.g., with Atlas from Silvaco or Sentaurus from Synopsys, are computationally onerous or even unviable when dealing with complex structures like practical transistor arrays and packages.
A more efficient, yet accurate enough, approach is based on the generation and solution of a SPICE-compatible purely-electrical macrocircuit that accounts for ET effects by means of the thermal equivalent of the Ohm's law (TEOL). In this macrocircuit, the power-temperature feedback is included through an equivalent thermal network (ETN, an electrical circuit relying on the TEOL), the components of which can be optimized in a pre-processing stage

•
More details of the individual transistor model and its subcircuit implementation are provided. • Differently from [15], where the arrays were assumed to lie on an unthinned GaAs substrate (as typical for known-good-die identification), here they are considered in a realistic phone-board environment, i.e., the substrate is thinned and attached on a laminate, the bottom of which is at T B = 358 K. • Similar to [21], in this work the linear power-temperature feedback is described by invoking FANTASTIC [22,23], which is fed with the COMSOL geometry/mesh accompanied with additional information (on position/shape of heat sources, boundary conditions, and thermal conductivities), and rapidly extracts an ETN based on the R TH matrix without performing simulations. Contrary to conventional ETNs (like the one used in [15]), the FANTASTIC network allows reconstructing the overall temperature field for selected bias conditions in a post-processing step. • In [15], the Kirchhoff's transformation was applied by assuming that all materials share the same nonlinear thermal behavior as GaAs. Unfortunately, this was found to lead to a perceptible overestimation of ET effects. Here, more realistic results are achieved by carrying out a suitable preliminary calibration procedure, similar to that made in [21].
By virtue of the above points, this work can be reviewed as an improved, and selfconsistent, version of [15].
The remainder of the paper is outlined as follows. In Section 2, the technology/layout details of test devices and arrays are given. Section 3 probes into the simulation approach; in particular, the transistor model, its subcircuit representation, the COMSOL geometry/mesh generation, the FANTASTIC extraction of the R TH -based ETN, and the assembling of the final macrocircuit are described. Section 4 reports and discusses the dc ET simulation results carried out by PSPICE. Conclusions are then drawn in Section 5.

Devices and Arrays
The structures investigated in this paper are mesa-isolated InGaP/GaAs NPN HBTs manufactured by Qorvo using an HBT-only process (referred to as HBT8) with two metal layers (e.g., [23,24]), the key features of which are reported in Table 1. The individual transistor, hereinafter also denoted as unit cell (Figure 1), is composed by four emitter fingers, each with a 2 × 20.5 µm 2 area (the total emitter area is then equal to 164 µm 2 ). The emitter is composed by a stack of four layers, namely (from the top): • an In 0.5 Ga 0.5 As cap to reduce the contact resistance with the gold (Au)-based emitter metallization; • a grading In x Ga 1-x As layer (with x spanning from 0.5 to 0) used to ensure a good lattice continuity with the underneath layer; • a GaAs layer acting as a set-back for an easier manufacturing process; • an n-doped In 0.49 Ga 0.51 P emitter layer, the bottom surface of which corresponds to the metallurgical base-emitter junction. Table 1. Key features of the investigated HBT technology.

Parameter Value
Common-emitter current gain at 300 K and medium current levels β F0 135 Open-emitter breakdown voltage BV CBO 27 V Open-base breakdown voltage BV CEO 17 V Peak cut-off frequency f T for V CE = 3 V 40 GHz Collector current density J C at peak f T for V CE = 3 V 0.2 mA/µm 2 Maximum oscillation frequency f MAX 82 GHz Figure 1. Schematic cross-section of the HBT unit cell highlighting the materials, the base-emitter junction, and the base-collector space-charge region (SCR).
The base is a thin p-doped GaAs layer lying on a thick n-type GaAs mesa. Both the base and collector contacts are designed with materials compatible with the manufacturing process in order to reduce the parasitic resistances. More specifically, the base contact is composed by the series of titanium (Ti) and platinum (Pt) layers that alloy through the leftover InGaP emitter on the base mesa, while the collector contact is located at the bottom of this mesa on the highly-doped GaAs subcollector layer and consists of a stack of Au-germanium (Ge)-nickel (Ni) layers. A metallization with two Au levels is exploited to electrically connect the unit cells and plays a role in the thermal exchange by favoring heat shunt and spread [24].
The electrical isolation from potential neighboring components is ensured by ionimplant damage (referred to as GaAs-ISO) resulting in an amorphous GaAs layer. Far from the active regions, the metal layers are separated from the GaAs substrate through a thin silicon nitride (Si 3 N 4 ) layer providing some shunt effect. A 10-µm-thick polybenzoxazole (PBO) layer covers the whole structure.
The analysis first focuses on devices fabricated on a 620-µm-thick GaAs substrate (thickness normally used for testing) and provided with 65 × 65 µm 2 pads (in a groundsignal-ground configuration) for bare-die experimental characterization through RF probes. The devices have been designed with a single unit cell, as well as with two and three paralleled cells, respectively (Figure 2, top). In the analysis reported in Section 4, the temperature T 0 = 300 K (hereinafter also denoted as reference temperature) is assumed to be applied to the substrate backside, which can be practically done with a thermochuck. Subsequently, the investigation is conducted on transistor arrays for PA output stages, where the unit cells are arranged in columns to (i) ensure an almost uniform distance from them to the through-wafer via providing the ground; (ii) meet the very stringent die size requirements. The arrays are assumed to be operating in a typical phone-board environment. Architectures with even and odd number of unit cells per column are examined to offer helpful guidelines to designers. In particular, arrays comprising 24 and 28 unit cells are considered, which are both arranged in four columns composed by six and seven cells each, respectively (Figure 2, bottom). In contrast to the simplified analysis performed in [15], (i) the arrays lie on a 100 µm-thick die, obtained by thinning the original 620-µm-thick wafer; (ii) the die is attached with epoxy to an 830 × 830 µm 2 -wide and 270 µm-thick laminate, which includes eight 600 × 600 µm 2 -wide 12 µm-thick Cu plates connected by nine circular Cu vias, all embedded in a dielectric; (iii) as previously mentioned, the laminate bottom is held at T B = 358 K, which is typical in the phone-board environment.

General Description
Similar to our previous papers [5][6][7][8][9][10][11][12][13][14][15]21], we chose to employ a circuit-based approach, which provides a good trade-off between computational overhead and accuracy. Such an approach makes use of the TEOL and can be summarized as follows: • Each four-finger unit cell is represented with one SPICE-compatible subcircuit (Section 3.3) implementing a simple analytical transistor model (Section 3.2). This assumption relies on the following considerations: (i) the two metals uniformly distribute the temperature over the base-emitter junction; (ii) the electron currents emerging from the four closely-spaced individual emitters are expected to spread and give rise to only one heat source. The subcircuit uses (i) a basic/standard bipolar transistor at reference (and unchangeable) temperature T 0 as a core component, and (ii) linear and nonlinear controlled sources to account for the variation of the temperature-sensitive parameters during the simulation run, as well as for other specific mechanisms. According to the TEOL, the temperature rise ∆T j = T j − T 0 averaged over the base-emitter junction (which mainly influences the ET device behavior) is actually a voltage, while the dissipated power P D is treated as a current. In addition to the standard transistor terminals (emitter, base, and collector), the unit-cell subcircuit is also equipped with an input node carrying the "voltage" ∆T j and with an output node offering the "current" P D .

•
The power-temperature feedback is described with a SPICE-compatible thermal feedback block (TFB), the construction of which is carried out in a pre-processing stage. The TFB contains an ETN including the matrix of self-heating (SH) R TH s of the unit cells and mutual R TH s among them. The inputs of the ETN are the powers P D dissipated by the cells (represented with currents), and the outputs are their temperature rises ∆T jlin = T jlin − T 0 for the test devices or ∆T jlinB = T jlin − T B for the arrays (all emulated with voltages) under linear thermal conditions. • The ETN is automatically determined through the following procedure. First, an accurate 3-D geometry/mesh of the domain is built in the COMSOL environment using an in-house routine; then, the geometry/mesh, along with additional information concerning position/shape of heat sources, boundary conditions, and thermal conductivities, is fed to FANTASTIC, which extracts the ETN in a really short time without the need of user's intervention/expertise or onerous FEM simulations (Section 3.5). Generally, the whole process is very fast and error-free. The adoption of FANTASTIC is an improvement over our prior contribution [15], where (i) the R TH matrix was calculated by performing N purely-thermal static COMSOL simulations by activating only one heat source at a time, and (ii) the simple ETN adopted did not allow a post-processing reconstruction of the whole temperature map in the domain.

•
As mentioned above, the ETN only accounts for linear thermal conditions. However, nonlinear thermal effects can be significant when particularly high temperatures are reached. Such effects are taken into account by making use of the Kirchhoff's transformation, which converts the linear temperature rises ∆T jlin (test devices) or ∆T jlinB (arrays) offered by the ETN into the nonlinear counterparts ∆T j = T j − T 0 and ∆T jB = T j − T B , respectively. Contrary to [15], here the transformation was properly calibrated (Section 3.4) to improve the ET simulation accuracy.

•
Besides the ETN, the TFB also includes N voltage-controlled voltage sources that apply the calibrated Kirchhoff's transformation to the ETN linear outcomes; as a result, the nonlinear temperature rises ∆T j (test devices) and ∆T jB (arrays) are computed; only for the arrays, the increment T B − T 0 is added to ∆T jB to get the N nonlinear ∆T j = T j − T 0 to be provided to the unit-cell subcircuits.

•
The subcircuits are then connected to the TFB in the environment of a commercial circuit simulation tool. As a result, the whole domain is transformed into a purelyelectrical macrocircuit, which inherently accounts for ET effects (Section 3.6): the temperature, and thus the temperature-sensitive parameters, are allowed to vary during the simulation run. The task of solving this macrocircuit is given to the powerful and robust engine of the circuit simulation tool, with very low computational effort and minimized occurrence of convergence issues compared to other numerical methods.

Bipolar Transistor Model
The analytical model used to describe the dc behavior of the unit cell is simple, accurate enough, and associated to a low-effort parameter extraction process; this provides high flexibility to the overall approach. The collector current in forward active mode is given by [13] where: is the reverse saturation current density at the reference temperature T 0 = 300 K; is an I E -dependent dimensionless term introduced to describe the attenuation dictated by high-injection (HI) effects, i.e., the Kirk-induced gain roll-off.
In this approach, the temperature dependence of I C is taken into account with a V BEj shift, while J S0 , η, and V T0 are kept at their T 0 values (e.g., [25]). Coefficient φ [V/K] is assumed to vary with emitter current I E according to the following logarithmic law [11,13,[26][27][28][29]: where k [eV/K] is the Boltzmann's constant, and q [C] is the (absolute value of the) electron charge; parameter φ 0 can be extracted by comparing (1) with I C -V BE characteristics measured at various backside temperatures at low current levels for V CB = 0 V [29]. As far as factor M is concerned, any model can in principle be adopted. For the GaAs HBTs under analysis, we chose the classic Miller model given by [30] BV CBO [V] and n AV being the open-emitter breakdown voltage and a fitting power factor, respectively. The HI attenuation term B HI is modeled as [5,13] where α F is the common-base (CB) forward current gain, while J HI [A/µm 2 ] and n HI are fitting parameters. The common-emitter (CE) forward current gain β F is modeled as where β F0 is the gain at T 0 , at medium current levels (i.e., before the Kirk-induced fall-off), and in the absence of Early effect, while ∆E G [eV] is the positive difference between the bandgaps of emitter and base yielding a negative temperature coefficient (NTC). The CB gain α F in (4) is related to the CE counterpart by The base current is given by (e.g., [13,31]) where use has been made of (1).
We note that replacing the simple M formulation (3) with a more complex one (e.g., [31]), the model can be adopted for bipolar transistors fabricated in any technology if a proper parameter calibration procedure is performed (in silicon devices, ∆E G is negative, since the emitter bandgap is narrower than that of the base due to high doping levels).
The parameter extraction procedure is straightforward (details are reported in [29]). The values adopted for the simulations described in Section 4 are reported in Table 2.

Parameter
Value

SPICE Unit-Cell Subcircuit
A sketch of the subcircuit is reported in Figure 3, which evidences the standard bipolar transistor model as a core component at temperature T 0 , as well as the additional ∆T j (input) and P D (output) terminals. Linear/nonlinear controlled voltage/current sources are added to enable the variation of the temperature-sensitive parameters during the simulation run, as well as to account for physical mechanisms not included in the standard transistor (for this reason, the resulting subcircuit is also popularly referred to as wrapper model). The II-free collector current I CnoAV accounting for the temperature dependence of V BEj , HI, and Early effects is computed with the nonlinear source denoted with A. The II-unaffected base current I BnoAV is calculated by source B as I CnoAV /β F , β F being evaluated by the nonlinear source C from V CB , I E , and ∆T j according to (5). The II (avalanche) current I AV is determined as (M − 1) × I CnoAV by source D, where factor M described by (3) is provided by the nonlinear source E; the collector current I C is obtained by adding I AV to I CnoAV , while the base current I B is given by I BnoAV -I AV . The lumped resistances R B and R E emulate the parasitic base and emitter resistances, respectively. The subcircuit also includes further nonlinear sources (omitted in Figure 3) to describe the cell behavior in saturation mode. The dissipated power P D under dc conditions (to be given to the TFB) is determined as

Construction of the Geometry/Mesh in COMSOL and Calibration of the Kirchhoff's Transformation
The 3-D geometry/mesh of each domain was constructed in the environment of the COMSOL software package [16] by making use of an in-house routine that relies on the MATLAB-COMSOL Livelink and the Toolbox for files in GDSII format provided by U. Griesmann [17]. The procedure is described as follows: first, the GDS layout file (i.e., the masks used for the technological process) is input to the routine along with the thicknesses of the layers, mask biases, and material parameters; then, the routine automatically draws the 3-D geometry, and finally generates and optimizes the mesh in the COMSOL environment. Such a process is error-free and requires a much shorter time compared to a painstakinglylong manual approach.
A single heat source (geometrically coinciding with the base-collector SCR [17]) and a single subcircuit were assigned to each unit cell ( Figure 1); as mentioned earlier (Section 3.1), this intrinsically assumes that the fingers within the cell are tightly thermally coupled, and thus not prone to thermal hogging. Figure 4 depicts the geometry of the 1-cell test device, while Figure 5 shows the mesh of the 2-cell one.  For the arrays, the horizontal symmetry was exploited to construct the geometry/mesh of only half of the domains ( Figure 2, bottom), thus alleviating the computational burden; the missing portions were virtually restored by applying an adiabatic boundary condition to the plane of symmetry. Figure 6 illustrates the geometry of (half of) the 24-cell array, while Figure 7 shows the mesh of (half of) the 28-cell array, both in COMSOL.  Besides the geometry/mesh of the domains and the position of the heat sources, FANTASTIC [22,23] requires information on the boundary conditions and thermal conductivities to extract the ETN including the R TH matrix (Section 3.5). As mentioned in Section 2, the backside of the GaAs substrate was assumed to be at T 0 for the test devices, while the laminate bottom was set to T B for the arrays. All the other surfaces were considered adiabatic (i.e., with zero outgoing heat flux), such an assumption being justified by the particular scratch-protection coating employed in the technology (thick PBO).
The thermal conductivities associated to the materials of the test devices (at T 0 ) and arrays (at T B ) are listed in Table 3. Table 3. Thermal conductivities of the materials composing the devices and arrays. Material 19.6 × 10 −6 (8), α = −0.33 [33] In 0.5 Ga 0.5 As 0.048 × 10 −4 [33] 3.9 × 10 −6 (8), α = 1.175 [33] In x Ga 1-x As (0 < x< 0.5) 0.092 × 10 −4 [33] average in the layer 7.4 × 10 −6 (8 The ETN is generated by FANTASTIC under linear thermal conditions (temperatureinsensitive thermal conductivities). However, as very high temperatures are reached, nonlinear thermal effects can no longer be neglected. More specifically, the thermal conductivities of the materials vary with temperature T according to the following laws: where (8) applies to semiconductors and insulators, and (9) to some metals; the accepted values of the α and β coefficients are also reported in Table 3. In order to account for the temperature dependences described by (8) and (9), we resorted to the Kirchhoff's transformation [18,19], which converts the linear junction temperature rises (∆T jlin and ∆T jlinB ) into the nonlinear counterparts (∆T j and ∆T jB ) through [36] ∆T for the test devices, and for the arrays. Contrary to the simplified analysis in [15], where m k was chosen equal to 1.25 (α value for GaAs) for all the domains under investigation, here a preliminary calibration procedure was performed for this parameter. Let us refer to the 1-cell test device for the sake of simplicity. This HBT was simulated with COMSOL over a wide range of dissipated powers P D by activating the thermal conductivity dependences upon temperature (8), (9) (nonlinear thermal conditions), and the average temperature rise over T 0 at the base-emitter junction ∆T j was computed for each P D . The linear temperature rise ∆T jlin over the same P D span was evaluated by multiplying P D by the R TH obtained by COMSOL under linear conditions. Then the Kirchhoff's transformation was applied to ∆T jlin , and m k was tuned so as to ensure the best agreement between the nonlinear temperature rise ∆T j calculated by the transformation and the realistic one evaluated by COMSOL; the optimum m k value was 0.8333 ( Figure 8). The same operation was repeated by activating one cell belonging to an array and considering the average junction temperature rise over T B ; in that case, the optimum m k value was found to be 0.8932. As a result, m k = 0.8333 was used in (10) for the ET simulations of the test devices (Section 4.1), whereas m k = 0.8932 was adopted in (11) for the ET simulations of the arrays (Section 4.2).

FANTASTIC
The FAst Novel Thermal Analysis Simulation Tool for Integrated Circuits (FANTAS-TIC), originally presented in [22,23], was conceived and developed to approximate a FEM model of heat conduction in an electronic device, having typically millions of DoFs, with a dynamic compact thermal model (DCTM), having only tens of DoFs, and the corresponding ETN [21]. The extraction of the DCTM and ETN does not require to solve the FEM model, as it is performed in short times through a refinement of the truncated balance-based moment matching approach to model-order reduction (MOR) introduced in [37]. Moreover, it also allows the extraction of a static compact thermal model (SCTM) and the related ETN relying on the R TH matrix (hereinafter indicated with R TH ); for this simplified case, FANTASTIC is even quicker, as briefly discussed below.
The heat conduction problem in the structure, assumed to be static and linear, is imported from either commercial (e.g., COMSOL) or open-source numerical tools in the form of: mesh discretizing the geometry, position/shape of the heat sources, boundary conditions, and thermal conductivities (also mass densities and specific heats are required for the dynamic case). Both hexahedral and tetrahedral meshes can be used. Arbitrary tensorial thermal conductivity distributions can be defined. Neumann's, Dirichlet's, or Robin's boundary conditions can be applied.
A FEM model of the problem is then assembled by FANTASTIC. In particular, the stiffness matrix K is constructed. High-order basis functions can be adopted; the typical choice is to select tetrahedral meshes and 2nd-order basis functions, as a good trade-off between accuracy and efficiency. The M DoFs of the temperature rise distribution, forming the M-row vector ϑ, are solution of the discretized heat conduction problem Kϑ = q (12) in which the power density distribution vector q takes the form q = QP D (13) In (13), P D is an N-row vector with the powers dissipated by the N 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 heat source, with n = 1, . . . , N. The port temperature rises of the N heat sources form the N-row column vector ∆T jlin (∆T jlinB for the arrays) defined as [37] ∆T jlin = Q T ϑ (14) In order to extract a compact thermal model, an M ×M matrix V withM M is defined, which allows expressing ϑ by means of a reduced numberM of DoFs, thus forming theM-vectorθ so that ϑ = Vθ The V matrix is used for projecting the discretized heat conduction problem (12)-(14) by the Galerkin's method, deriving an SCTM in the form Kθ =q (16) whereq =ĜP D ∆T jlin =Ĝ T ϑ (18) in whichK is anMth-order matrix andĜ is anM × N matrix. The V matrix is determined by the Algorithm 1 reported below.

Algorithm 1: SCTM extraction
Set V:=0 for each heat source n=1, . . . , N do 1 Solve (21) for Θ n 2 Update matrix V by appending Θ n 3 Generate a SCTM projecting (12)-(14) onto V At line 1, the temperature response to the n-th heat source is solved for the static heat conduction problem. Thus, equation is solved for Θ n , e n being the vector selecting the n-th column of Q. Since the coefficient matrices of these linear systems are symmetric positive definite, the most efficient multigrid iterative solvers can be used for their solution. At line 2, the V matrix is updated by appending vector Θ n to its columns if it is linearly independent with respect to them.
At line 3, the SCTM is determined proceeding as in (16)- (18), and has dimensionM ≤ N M. A much smaller dimension is thus obtained with respect to the dynamic case; moreover, since a much smaller number of discretized heat conduction problems are solved, a large speedup of the algorithm is also achieved.
It is now observed that by solving at limited cost the eigenvalue problem having as unknown theMth-order orthogonal matrixÛ, in whichΛ is anMth-order diagonal matrix, and introducing the change of variableŝ ϑ =Ûξ (23) the SCTM Equations (16)- (18) are transformed into the equivalent form ∆T jlin =Γ Tξ (25) whereΓ is theM × N matrixV TĜ . The spatial distribution of the temperature rise is then Being Ξ = VÛ an M ×M matrix like V. Theξ vector encompasses the DoFs of the thermal field. These DoFs are the node temperature rises in the extracted ETN sketched in Figure 9, which is governed by (24) and (25) and represents a simplified version of the dynamic counterpart [21][22][23], as it benefits from a much lower (by one order of magnitude) number of nodes, resistances, and controlled sources. The ETN transforms the port powers P D into the port temperature rises ∆T jlin by inherently accounting for the boundary conditions initially applied to the FEM model. The port response of this network defines the N th -order resistance matrixR TH [K/W] of the SCTM given byQ . Since by construction theM columns of the V matrix span the N columns of matrix K −1 Q, it is straightforward to prove thatR TH = R TH , where R TH is the thermal resistance matrix of the discretized heat conduction problem (12)- (14) given by Q T K −1 Q. This further strengthens the general result relating the thermal impedance matrices of the compact thermal model and of the discretized heat conduction problem in the dynamic case [38]. The resulting ETN is particularly well-suited to be solved by means of modified nodal analysis in SPICE-like circuit simulators, since all circuit elements are voltage-controlled, and thus the number of variables added by the SCTM is limited toM. The topology is general and can be implemented into any circuit simulation program.
It is worth noting that, after the circuit simulation, as a post-processing stage, the whole spatial distribution of temperature rise in the examined domain can be reconstructed at negligible computational cost and memory storage using (26), for both thermal-only and ET simulations. The temperature map can be plotted by any proper tool, like e.g., Paraview. This option is not allowed using conventional ETNs like in our former paper [15].

Construction of the Macrocircuit
The ETN derived by FANTASTIC was enriched with N nonlinear voltage-controlled voltage sources to account for the Kirchhoff's transformation, and the TFB was then obtained. The ∆T j and P D nodes of the subcircuits (the unit cells) were connected to the TFB, thereby giving rise to the whole TEOL-based SPICE-compatible macrocircuit, a simplified scheme of which is reported in Figure 10. As previously mentioned, the solution of the macrocircuit was delegated to PSPICE [20], although any other commercial circuit simulation software (e.g., LTSPICE, Eldo, Keysight ADS [39], and SIMetrix) could in principle be used. Figure 10. Sketch of the merely-electrical macrocircuit including ET effects through the TEOL. The paralleled (sharing the same base, emitter, collector contacts) subcircuits describing the unit cells are connected to the TFB containing: the ETN preliminarily evaluated by FANTASTIC, the calibrated Kirchhoff's transformation, and a block adding 58 K (for the arrays only).

Extension to the Dynamic Case
Unlike advanced bipolar transistor models equipped with temperature-dependent parameters and a thermal node (like HiCUM [40], VBIC, AHBT, and Mextram504, all available in ADS), the proposed unit-cell model/subcircuit is only suited to fairly well describe forward active and saturation modes under dc conditions. RF simulations including ET effects can be enabled by resorting to a variant of our approach based on one of the above models in the ADS environment; the strategy can be described as follows.

•
The selected transistor model must be provided with a power (output) node, and the internal one-or two-pair thermal network has to be deactivated.

•
All parameters of the model must be extracted from experimental data, which is a nontrivial task. • As mentioned in Section 3.5, if FANTASTIC is also fed with the mass density and specific heat for all materials, it can be enabled to extract a DCTM of the domain and the associated ETN accounting for the dynamic heat propagation [21][22][23]. Such an ETN, together with the Kirchhoff's transformation sources, will constitute the SPICEand ADS-compatible TFB. • Lastly, the macrocircuit has to be built in ADS by connecting the model instances among them and with the TFB.
It is worth noting that the adoption of the one-or two-pair thermal networks embedded in the transistor models instead of our TFB (i) would lead to a significant inaccuracy in terms of SH of the single cell (the typically-used single pair is not enough for the transient SH response [23]), and (ii) would exclude the mutual thermal interactions among unit cells, which however play a relevant role, as demonstrated in Section 4.
As an alternative, one could resort to an ET solver available in recent ADS releases, which couples a quasi−3-D layout-based numerical thermal tool to the circuit simulator (the thermal networks embedded in the model instances being disabled) [41]. However, (i) this strategy is only applicable to device models equipped with a thermal node; (ii) using such a solver for layout optimization is very labor-intensive; (iii) the iterative process leading to convergence is resource-hungry.

Test Devices
The analysis of the test devices is two-fold: it is intended (i) to offer an overview of the ET-and II-induced positive-feedback mechanisms limiting the safe operating region, and (ii) to explore the thermally-stabilizing effect ensured by the base ballasting. In the latter case, an integrated resistor R Bext = 400 Ω was used per each cell, as this is the typical ballasting strategy chosen for the arrays investigated in Section 4.2.
Let us first consider the single-cell device, which does not require an ETN, but only a linear SH R TH to which the Kirchhoff's transformation is applied. The R TH was computed to be about 440 K/W, which is in good agreement with the experimental value extracted by means of the classic method proposed in [42]. Figure 11 shows the V BE -constant I C -V CE and ∆T j -V CE characteristics. It can be inferred that the curves of the unballasted device are affected by a flyback (also denoted as snapback or turnover) mechanism followed by a negative-differential-resistance (NDR) branch [6,13,31,[43][44][45][46], which can be simulated/measured by (i) incrementing I C or (ii) connecting a resistor R Cext to the collector terminal, sweeping V Cext on the available resistor node (Cext), and evaluating V CE as V CextE − R Cext ·I C . It is worth noting that increasing V CE beyond the value corresponding to the flyback would have instead led to a thermal runaway (shown for V BE = 1.25 V) and sudden device failure. Adding a base-ballasting resistor (as suggested in [31,47,48]) with R Bext = 400 Ω prevents the flyback (and thus the runaway); however, it reduces the internal (junction) V BEj at medium/high current levels, decreasing and limiting the collector current.  Figure 12 shows the behavior of the 2-cell test device under I BTOT -constant conditions. As can be seen, a bifurcation phenomenon is triggered: beyond a critical V CE , cell #2 tends to conduct the whole current, while #1 gradually turns off [26,[49][50][51]. For I BTOT = 0.5 mA, the critical V CE is equal to 8 V for the unballasted case, and reduces with increasing I BTOT . It is not possible to identify a priori which cell will take all the current, since it is determined by random (and unavoidable) technology and layout fluctuations. In PSPICE the cells are assumed electrically identical, and the current hogging in cell #2 is favored by a marginally higher SH R TH due to a slight layout asymmetry. As far as the total collector current I CTOT is concerned, a smooth NDR region due to the β F NTC is observed in the V CE range where I CTOT is equally shared by the two cells. The NDR mechanism is then replaced by a marked I CTOT reduction within the bifurcation region due to the 'faster' temperature growth with V CE in the hotter cell, which implies a significant β F decrease; such a behavior is also denoted as collapse of collector current [49] or collapse of current gain [26,50,51]. As V CE exceeds 13 V, the II current I AV2 (≈90 µA at V CE = 15 V) can no longer be neglected with respect to I BTOT (=500 µA); as a result, the avalanche-less current I BnoAV2 , almost equal to I BTOT + I AV2 , perceptibly grows due to the I AV2 increase, and in turn raises I CnoAV2 ≈ I C2 ≈ I CTOT . The inclusion of base ballasting with R Bext = 400 Ω per cell pushes the critical V CE to 16.6 V, restoring a uniform behavior over a large voltage range [47].  Figure 13 illustrates the behavior of the 2-cell test device under V BE -constant conditions. Let us first examine the unballasted case. If I CTOT is swept, the current is equally divided between the unit cells until a flyback takes place, followed by an NDR branch still showing uniform operation; at low V BE (e.g., 1.2 V), a bifurcation will also occur for relatively low cell temperatures [31]. A similar behavior was observed for more paralleled BJTs in SOG technology [7]; in that case, an uneven current distribution was found to arise beyond the uniform NDR branch under I CTOT -controlled conditions. By increasing V CE , the whole device (or at least one of the two cells) would have blown up beyond the value corresponding to the flyback; this means that our simulations do not confirm the observation of a 'safe' collapse encountered in [50,51]. Let us next consider the device with base-ballasted cells subject to the same biasing conditions ( Figure 14). As far as the V BE = 1.2 V case is concerned, the flyback point moves to the left, thus shrinking the V CE -controlled safe operating region; this can be ascribed to the increased avalanche current flowing in the unit cells, in turn induced by R Bext = 400 Ω, which favors the positive feedback action related to II [31,48]. On the other hand, the thermally-induced bifurcation disappears; it was found than the bifurcation-triggered discrepancy between I C1 and I C2 reduces with increasing R Bext , and R Bext = 400 Ω is the threshold value for which the uniform behavior is fully restored. For higher V BE s, the ballasting makes the flyback mechanism vanish, thus improving the thermal stability of the device; however, the current capability significantly plummets at high current levels. It is also important to analyze the ET behavior of the 2-cell device under I ETOTcontrolled conditions, typical for differential pairs and comparators [52,53]. Results obtained by increasing the total emitter current I ETOT for V CB = 8 V are shown in Figure 15 for both the unballasted and ballasted cases. It is shown that without ballasting a bifurcation phenomenon is triggered for a critical I ETOT (33 mA), and the whole I CTOT eventually flows in cell #2, rapidly increasing its temperature; a similar behavior would be obtained by fixing I ETOT and sweeping V CB [6,31]. Applying the resistors R Bext = 400 Ω, thermal effects are weakened and the bifurcation disappears, limiting the maximum temperature reached by the device. We now consider the ET behavior of the test device composed by three unit cells, assumed ideally identical in PSPICE. The first simulation was performed under CE conditions by increasing V CE with a constant I BTOT = 0.7 mA (Figure 16). Let us focus on the unballasted case. It is found that cell #2 starts conducting more current due to the thermal coupling with both the adjacent (outer) cells #1 and #3. For higher V CE values, a counterintuitive behavior takes place: a bifurcation mechanism involving cells #1 and #3 occurs at V CE = 7 V; in particular, cell #3 eventually bears more current, whereas #1 turns off, as induced by the slightly higher SH R TH of #3 with respect to #1. By further increasing V CE , cell #3 prevails over #2 since the SH R TH of #3 (396.2 K/W under linear thermal conditions) is perceptibly higher than that of cell #2 (361.7 K/W), which experiences a more effective heat flow through metal 2 (top metal). The strongly uneven current distribution for V CE > 7 V turns into a collapse in the I CTOT -V CE curve. For V CE > 13 V, cell #3 conducts the whole current, and therefore I B3 = I BnoAV3 − I AV3 is almost equal to I BTOT = 0.7 mA; the increase in I AV3 with V CE due to the enhanced II leads to a growth in I BnoAV3 , which dominates over the NTC of β F3 , thus driving an I C3 ≈ I CTOT increase. Adopting R Bext = 400 Ω for all cells pushes the uneven current distribution to V CE > 15.7 V, thus leading to a much wider uniform operating region. It is worth noting that increasing the emitter area of cell #1 by only 1 µm 2 , the onset of the bifurcation takes place at the same critical V CE , but the behavior of cells #1 and #3 reverses: cell #1 prevails over #3, and eventually sinks also the current of #2. Nevertheless, the I CTOT -V CE curve would coincide with that obtained for cells sharing ideally identical areas. This means that it is impossible to foresee which of the outer cells will dominate, since it depends on unavoidable technological/layout discrepancies. The practical implication is that the failure analysis should look at planes of symmetry, and not at the specific failed cells.
An overview of the I CTOT -V CE characteristics for various I BTOT values is illustrated in Figure 17 for both the unballasted and base-ballasted cases; it is apparent how the collapse locus (which can be reviewed as the boundary of the safe operating area) significantly moves rightward by virtue of the external base resistances [47]. In the above analysis, it was stated that the 3-cell test device benefits from a lower SH R TH of the inner cell #2 compared to the outer cells. We decided to quantify the related thermally-stabilizing effect by performing a test simulation where the R TH of cell #2 was considered identical to those of the lateral cells (396 K/W under linear thermal conditions) and the mutual R TH s between adjacent cells were assumed to coincide (115 K/W). Figure 18 shows that the classic collapse with #2 conducting the whole current is obtained. Unfortunately, the collapse onset in the I CTOT -V CE characteristic takes place at the same V CE as in the real test device with lower R TH for cell #2; this means that the bifurcation mechanism occurring for the outer cells in the real case is as deleterious as the strong thermal coupling affecting #2 in the ideal device with uniform SH R TH s. Figure 18. Unballasted 3-cell test device with unit cells sharing ideally identical SH R TH s: simulated collector currents against collector-emitter voltage V CE . Also shown is the I CTOT -V CE characteristic corresponding to the real structure (solid black line).

Transistor Arrays
As already mentioned in Sections 2 and 3.4, only half of each array for PA output stages was drawn, meshed, and simulated by exploiting the horizontal symmetry. This means that only 12 (14) cells were taken into account for the 24-cell (28-cell) arrays. Our simulation approach allows monitoring the dc ET behavior of the arrays at cell-level, that is, besides the junction temperature, all the key parameters, voltages, and currents are available for each cell; this would have been unviable by performing experiments since the collector layers of the individual cells are all shorted to one another, by being in a single isolation tub.
Hereinafter, I CTOT and I BTOT are the total collector and base currents of the semi-arrays. The analysis focuses on I BTOT -constant CE conditions, since GSM PAs are more or less biased with a constant I BTOT . A base ballasting with the nominal value R Bext = 400 Ω per unit cell was applied following a strategy denoted as segmented or split, in which such resistances only appear in the dc path and there is no RF performance penalty. The bottom of the laminate was held at T B .
Let us first consider the halved 24-cell array. Figure 19 reports the PSPICE results corresponding to I BTOT = 2 mA; the simulation lasted less than 100 s on a normal PC, in spite of the very small V CE step used. A significant current/temperature nonuniformity is observed as V CE exceeds 9 V, leading to an I CTOT collapse; more specifically, the right (internal) column (#7 to #12) conducts the entire current, while the cells belonging to the left (external) one (#1 to #6) tend to turn off. Below V CE = 11 V, the symmetric cell pairs of the right column (namely, #9 and #10, #8 and #11, and #7 and #12) share the same current. For V CE higher than 11 V, a bifurcation mechanism occurs for all these pairs, as dictated by slight layout asymmetries: in particular, the top cells (#7, #8, and #9) prevail over the bottom counterparts (#10, #11, and #12, respectively). The uneven behavior is thus enhanced, and I CTOT decreases more steeply with V CE . Over the V CE span from 11 to 12 V, cell #9 is the hottest cell since it suffers from the stronger thermal coupling with the surrounding cells. As V CE exceeds 12 V, only the adjacent cells #7 to #9 are conducting, whereas cells #10 to #12 run dry; as a consequence, mutual thermal interactions among cells play a minor role, while the behavior is dominated by the SH R TH s of the "active" cells, namely, 541, 533.3, and 531.5 K/W for #7, #8, and #9, respectively, under linear conditions (the closer the cell to the die border, the higher the R TH ). For a higher V CE , the total current first focuses over cells #7 and #8 (at V CE = 12.6 V, ∆T j7 = 500 K, and ∆T j8 = 700 K), and then, if #8 survives, over cell #7 only. It must be remarked that introducing an intentional (very small) technological discrepancy between #9 and #10 to favor a slightly higher conduction of #10, for V CE > 11 V the bottom cells of the right column dominate over the top ones; by further increasing V CE , cells #11 and #12 sink the whole current, which eventually focuses in the outer #12 suffering from the highest R TH . The I CTOT -V CE curve, including the collapse region, remains instead unaltered.
If the applied I BTOT is higher, the current nonuniformity (i.e., the collapse onset) and the subsequent bifurcation phenomenon involving symmetric cells occur for slightly lower critical V CE s, as ET effects are exacerbated by the higher currents (and dissipated powers). Clearly, the temperatures reached by the right-column cells at the bifurcation are higher than in the lower-I BTOT case. Figure 20 illustrates the scenario corresponding to I BTOT = 3 mA; here it is shown that, as the bifurcation takes place, the ∆T j shared by the inner cells #9 and #10 has already exceeded 600 K; this suggests that in a real array the metallization and surrounding interlevel dielectrics of these cells is likely to melt before the bifurcation arises. In the total absence of ballasting, the V CE range enjoying uniform current and temperature distribution dramatically reduces; beyond a critical (and low) V CE , the pair #9, #10 starts taking more current, which translates in the collapse onset in the I CTOT -V CE characteristic; then, the bifurcation between symmetric cells in the right column arises, and, for a slightly higher V CE , the current flows only in cell #9. An example is reported in Figure 21, which corresponds to I BTOT = 3 mA. It is observed that, as #9 dominates, the cells symmetric with respect to #9 (i.e., #8 and #10, and #7 and #11) tend to exhibit a similar behavior (i.e., they conduct almost the same current and are at the same temperature). Again, by intentionally applying a technological discrepancy leading to a slightly higher current capability for cell #10, the current will focus only on this cell in the deep collapse region, without distorting the I CTOT -V CE characteristic. Finally, the case of emitter ballasting with R Eext = 4 Ω per cell is examined, as it is considered approximately equivalent to R Bext = 400 Ω by circuit designers; however, in this case a significant degradation in terms of f T and f MAX (both by some GHz), as well as of RF gain, is induced. Figure 22 depicts the individual collector currents of the unit cells and the associated temperature rises above T 0 with I BTOT = 3 mA applied. By virtue of the reduced ET feedback with respect to the unballasted array, (i) the collapse locus is shifted ahead of about 3.5 V and (ii) the bifurcation involving the right-column cells occurs for much higher temperature values; (iii) in the thermally-unstable nonuniform region, the behavior of the individual currents vs. V CE resembles that observed in the unballasted case: beyond the bifurcation onset, a current hogging over cell #9 takes place. Another interesting result is that the stabilizing effect is much less effective than the base ballasting with R Bext = 400 Ω per cell: for this specific I BTOT , the critical V CE triggering the collapse is about 6.5 V, whereas R Bext = 400 Ω allows extending this value to 9 V.  , and no ballasting, as well as the corresponding temperature rises over T 0 affecting cell #9. A wide I BTOT range was selected so as to cover the typical operating current densities [54]. The considerable reduction in the safe operating region for the unballasted case and the inferior aid ensured by R Eext = 4 Ω with respect to R Bext = 400 Ω are apparent. From the overall analysis, the following relevant findings emerge:

•
The collapse onset in an I BTOT -constant I CTOT -V CE curve is associated to a uneven current/temperature distribution, wherein the right-column cells (in particular, the inner ones) bear almost all the current.

•
The steeper I CTOT drop in the collapse region is induced by a bifurcation mechanism involving the symmetric cells belonging to the right column. In the base-ballasted case (with R Bext = 400 Ω) this leads to three adjacent cells conducting all the current (either the top or the bottom ones, depending on small technological/layout discrepancies); a further V CE increase makes the current flow in two cells, and eventually in only one cell (the outer one). In the unballasted and emitter-ballasted case (with R Eext = 4 Ω), for V CE slightly higher than that entailing the bifurcation, the current flows in only one of the inner cells (#9 or #10). This leads to a very sharp and linear temperature increment vs. V CE of this cell (plainly illustrated for cell #9 in Figure 23b). • Such a linear nature of the ∆T j9 -V CE behavior can be straightforwardly explained as follows. Neglecting II effects, which is reasonable in both the unballasted and emitter-ballasted cases, ∆T j9 is approximately equal to R TH99 (T j9 )·β F (T j9 )·I BTOT ·V CE + 58 K (I B9 ≈ I BTOT ), where R TH99 is the SH thermal resistance of cell #9, and 58 K is the difference between T B and T 0 ; as V CE increases, there is a compensation between the NTC of β F and the increase in R TH99 with temperature due to nonlinear thermal effects.
We shall now focus on half of the 28-cell array; all cells are base-ballasted with R Bext = 400 Ω. Figure 24 depicts the PSPICE results obtained under CE conditions for I BTOT = 3.5 mA. It is found that beyond V CE = 8.5 V the current distribution becomes uneven, leading to the onset of collapse in the I CTOT -V CE curve. The right-column cells (#8 to #14) carry all the current, while the left-column counterparts (#1 to #7) turn off. In particular, the current is mostly conducted by the inner cells: the highest by the innermost #11, a slightly lower one by each of the symmetric cells #10, #12, an even lower by #9, #13, and so on. At V CE = 10.8 V, a bifurcation occurs for these symmetric pairs, and #10, #9, and #8 prevail over #12, #13, and #14, so that only the top group composed by the adjacent #8 to #11 cells conducts. Further increasing V CE , the current would first focus over the 3-cell group #8 to #10, then over the pair #8, #9, and would eventually flow only in the outermost cell #8, which is affected by the highest R TH . On the other hand, the reduction in the number of conducting top cells is found to occur after the temperature rise of cells #11, #10, and #9 has well exceeded 600 K; consequently, in the practical case, the metallization and dielectrics over such cells are expected to lose their integrity when still the whole group #8 to #11 is conducting. As mentioned in Section 3.5, unlike conventional ETNs, the one derived by FAN-TASTIC allows the reconstruction of the whole linear temperature rise field ∆T lin (x, y, z) = T lin (x, y, z) − T 0 , which can be easily processed through the calibrated Kirchhoff's transformation to get the nonlinear ∆T(x, y, z), for selected biasing conditions. Figure 25 shows the ∆T(x, y, z) map for half of the 28-cell array biased with I BTOT = 3.5 mA and V CE equal to 6 V (as representative of the uniform operating region), 10 V (collapse region: all the right-column cells are conducting), and 13.3 V (deep/marked collapse: after the bifurcation mechanism, the current flows only through the group #8 to #11). Lastly, an interesting and fair comparison is performed between the halved portions of the base-ballasted 24-and 28-cell arrays by applying various I BTOT values ensuring the same total base current density J BTOT = I BTOT /A ETOT in both cases. Results are reported in Figure 26; more specifically, Figure 26a shows the J CTOT -V CE characteristic (J CTOT = I CTOT /A ETOT ), and Figure 26b depicts the junction temperature rise of cell #9 for the 24-cell array, and of cell #11 for the 28-cell one. The analysis demonstrates that the 28cell array featuring an odd number of cells (seven) per column is less thermally robust than the 24-cell counterpart with an even number of cells (six) per column (the collapse locus is shifted leftward). This is attributed to the fact that the nonuniform current/temperature distribution triggering the collapse in the 28-cell structure implies that the innermost cell #11 takes more current than the others, while two cells (#9 and #10) concurrently bear this task in the 24-cell array. It is worth noting that, although this seems to be an expected results, there are still many designers using arrays with an odd number of cells per column.

Conclusions
In this paper, an efficient simulation approach has been used to analyze the dc electrothermal behavior of test devices and arrays for output stages of power amplifiers in InGaP/GaAs HBT technology. The approach relies on a full circuit representation of the structures under investigation (also referred to as macrocircuit), the solution of which can be evaluated by any circuit simulation software and requires tens of seconds at most, despite the complexity of the analysis. The macrocircuit is based on the thermal equivalent of the Ohm's law, and includes subcircuits to describe the unit cells, as well as a thermal feedback block to account for the power-temperature feedback. The thermal feedback block is obtained by combining (i) a FEM thermal tool aided by an in-house routine to generate an exceptionally accurate 3-D geometry/mesh of the structure, (ii) the FANTASTIC code to automatically get an equivalent thermal network without the need of simulations or user's experience, (iii) a preliminarily-calibrated Kirchhoff's transformation to include nonlinear thermal effects. The test devices have been analyzed under all the bias conditions of interest. An overview has been given on the thermally-and avalanche-induced distortion in the I-V characteristics, the limits of the safe behavior have been identified, and the beneficial effect of base ballasting has been explored. As far as the HBT arrays are concerned, the approach better shows its potential, as the individual currents, temperatures, and key parameters of all unit cells can be monitored. Some of the most important findings are: the base ballasting has been found to be more effective than emitter ballasting, which is typically suggested for breakdown-limited bipolar transistors; the arrays are more thermally-robust if arranged in columns with an even number of unit cells, where two central cells concurrently bear the heat coming from the outer ones as a nonuniform operating condition occurs. The approach is well suited to support engineers and designers in making choices oriented to develop more thermally-rugged and reliable circuits through, e.g., layout and/or nonuniform ballasting optimization.