Adrien Poupardin
email: apoupardin@estp-paris.eu

Eric Calais

Philippe Heinrich

Hélène Hébert

Mathieu Rodriguez

Sylvie Leroy

Hideo Aochi

Roby Douilly

Deep submarine landslide contribution to the 2010 Haiti earthquake tsunami

The devastating M w 7.1 Haiti earthquake in 2010 was accompanied by local tsunamis that caused fatalities and damage to coastal infrastructure. Some were triggered by slope failures of river deltas in the close vicinity of the epicenter, while others, 30 to 50 km to the north across the Bay of Gonâve, are well explained by the reverse component of coseismic ground motion that accompanied this mostly strike-slip event. However, observations of run-up heights up to 2 m along the southern coast of the island at distances up to 100 km from the epicenter, as well as tide gauge and DART buoy records at distances up to 600 km from the epicenter, have not yet received an explanation. Here we demonstrate that these observations require a secondary source, most likely a submarine landslide. We identify a landslide scar 30 km from the epicenter off the southern coast of Haiti at a depth of 3500 m, where ground acceleration would have been sufficient to trigger slope failure in soft sediments. This candidate source, 2 km 3 in volume, matches observations remarkably well assuming that the sediment collapse obeys a viscous flow with an initial apparent viscosity of 2×10 5 Pa s. Although that particular source cannot be proven to have been activated in 2010, our results add to a line of evidence that earthquake-triggered submarine landslides can cause significant tsunamis in areas of strike-slip tectonic regime.

Introduction

The devastating M w 7.1 Haiti earthquake in 2010 occurred within the Caribbean-North America plate boundary, where oblique relative motion at 19 mm yr -1 is partitioned between shortening on the North Hispaniola fault to the north and strike-slip motion on the E-W-striking Septentrional and Enriquillo faults throughout the island (Fig. 1; Symithe et al., 2015). Although the earthquake was first assumed to have occurred on the Enriquillo fault, several independent studies later showed that it actually ruptured a previously unmapped fault -the Léogâne fault -with a source mechanism combining strike-slip and reverse faulting [START_REF] Calais | Transpressional rupture of an unmapped fault during the 2010 Haiti earthquake[END_REF][START_REF] Hayes | Complex rupture during the 12 January 2010 Haiti earthquake[END_REF][START_REF] Hashimoto | Fan-delta uplift and mountain subsidence during the Haiti 2010 earthquake[END_REF][START_REF] Symithe | Coseismic Slip Distribution of the 2010 M7.0 Haiti Earthquake and Resulting Stress Changes on Regional Faults[END_REF], in a setting resembling the 1989 Loma Prieta earthquake in California (e.g., [START_REF] Dietz | The October 17, 1989, Loma Prieta, California, Earthquake and its aftershocks: Geometry of the sequence from high-resolution locations[END_REF].

In spite of its moderate size and primarily strike-slip mechanism, the earthquake was accompanied by local tsunamis that caused at least three fatalities in Grand Goâve [START_REF] Fritz | Twin tsunamis triggered by the 12 January 2010 Haiti earthquake[END_REF] and some damage to coastal infrastructure. Thanks to the rapid deployment of an International Tsunami Survey Team (ITST), run-up height, flow depth and inundation were documented at more than 20 sites along the Haitian coastline (Fig. 1; [START_REF] Fritz | Twin tsunamis triggered by the 12 January 2010 Haiti earthquake[END_REF]. Maximum tsunami heights reached 3 m in the close vicinity of the earthquake epicenter in the bays of Petit Goâve and Grand Goâve along the northern coast of the southern peninsula, where the tsunami was triggered by slope failures of river deltas [START_REF] Hornbach | High tsunami frequency as a result of combined strike-slip faulting and coastal landslides[END_REF][START_REF] Fritz | Twin tsunamis triggered by the 12 January 2010 Haiti earthquake[END_REF]. A total of 30 to 50 km to the north across the Bay of Gonâve, eyewitnesses reported drawdown near the locality of Luly, well explained by the reverse component of the coseismic ground motion [START_REF] Hornbach | High tsunami frequency as a result of combined strike-slip faulting and coastal landslides[END_REF].

However, observations of run-up heights up to 2 m along the southern coast of Haiti, at distances up to 100 km from the epicenter, as well as tide gauge and DART buoy records at distances of 300 and 600 km from the epicenter, respectively, have not yet received a satisfactory explanation. [START_REF] Fritz | Twin tsunamis triggered by the 12 January 2010 Haiti earthquake[END_REF], using the early source model from the National Earthquake Information Center (NEIC), were able to roughly match these observations by scaling coseismic slip by a factor of 4 while keeping a constant seismic moment, hence downscaling the regional rigidity coefficient. A similar ad hoc procedure was also proposed by Newman et al. (2011) for the 2010 Mentawai, Indonesia, earthquake in order to account for the discrepancy between coseismic slip inverted from teleseismic waves and actual slip derived from geodetic measurements. This artificially lowered rigidity can indeed account for slow ruptures at subductions -"tsunami earthquakes" [START_REF] Kanamori | Mechanism of tsunami earthquakes[END_REF] -such as in Mentawai in 2010 or in Java in 2006 [START_REF] Hébert | The 2006 July 17 Java (Indonesia) tsunami from satellite imagery and numerical modelling: a single or complex source?[END_REF], but the Haiti earthquake fault ruptured at 2.6 km s -1 , a regular rupture velocity (de Lépinay et al., 2011).

Here we revisit this issue by (1) testing most of the finitesource models proposed for the 2010 Haiti earthquake, (2) investigating the role that an earthquake-triggered landslide off the southern coast of Haiti may have played in tsunami generation. We show that this latter hypothesis is very likely, in line with the growing set of evidence linking tsunamis along strike-slip faults to submarine landslides (e.g., [START_REF] Ma | The origin of the tsunami excited by the 1989 Loma Prieta earthquake faulting or slumping?[END_REF][START_REF] Yalçıner | Tsunamis in the Sea of Marmara: Historical documents for the past, models for the future[END_REF][START_REF] Rodriguez | Tsunamigenic submarine landslides along the Xauen-Tofiño banks in the Alboran Sea (Western Mediterranean Sea)[END_REF].

Tsunami calculation method

The numerical method involves modeling of the initiation, propagation and run-up of the tsunami waves. For the earthquake source, we calculate the initial seafloor perturbationassumed to be instantaneously and fully transmitted to the water column -using coseismic static ground displacement for a dislocation in an elastic half-space [START_REF] Okada | Internal deformation due to shear and tensile faults in a half-space[END_REF]. For the landslide source, we assumed that the sediment collapse obeys a viscous flow, consistent with the mud or shale nature of the sediments involved [START_REF] Assier-Rzadkiewicz | Numerical modelling of a landslide-generated tsunami: the 1979 Nice event[END_REF]. We keep the dynamic viscosity constant during the whole simulation, which may not be realistic in terms of landslide propagation but has little effect on the free surface, which is mostly sensitive to the initial landslide acceleration (e.g., Løvholt et al., 2015;[START_REF] Poupardin | The 1979 Submarine Landslide-Generated Tsunami in Mururoa, French Polynesia[END_REF].

The tsunami propagation calculation accounts for dispersive terms in the deep ocean (Boussinesq equations) and neglects them near the coast (St Venant equations), retaining only the lowest-order nonlinear terms [START_REF] Poupardin | Traveltime delay relative to the maximum energy of the wave train for dispersive tsunamis propagating across the Pacific Ocean: the case of 2010 and 2015 Chilean Tsunamis[END_REF]. In the latter case, we use a shock-capturing method to propagate strong nonlinear waves. In the former one, we use an iterative Crank-Nicholson method to solve the Boussinesq equations (details in [START_REF] Poupardin | The 1979 Submarine Landslide-Generated Tsunami in Mururoa, French Polynesia[END_REF]. Propagation equations are solved using a finite-difference scheme with a spatial discretization that uses centered differences for linear terms and forward differences for advection terms [START_REF] Hébert | Far-field simulation of tsunami propagation in the Pacific Ocean: impact on the Marquesas Islands (French Polynesia)[END_REF].

We solve the numerical problem using four successive levels of nested bathymetric grids from the deep ocean to the local coastal areas of interest. The two lower-resolution grids use the GEBCO World Bathymetry with a resolution from 1600 to 400 m. We digitized local bathymetric charts and included existing digital bathymetric data in the bays of Portau-Prince, Jacmel and Santo-Domingo to build local grids with a resolution of 25 m in order to accurately model nearshore resonance and amplification. We used the 1600 m resolution grid to propagate the tsunami to the DART Buoy and the 25 m one for the Santo Domingo tide gauge in order to account for resonance due to coastal bathymetry.

Earthquake-generated tsunami

We tested a uniform slip of 1 m using the NEIC fault geometry and rupture mechanism, as in [START_REF] Fritz | Twin tsunamis triggered by the 12 January 2010 Haiti earthquake[END_REF] hereafter denoted the F model), as well as four finite-source models with variable coseismic slip derived from combinations of GPS, InSAR, coastal uplift and teleseismic data [START_REF] Hayes | Complex rupture during the 12 January 2010 Haiti earthquake[END_REF]Meng et al., 2012;[START_REF] Symithe | Coseismic Slip Distribution of the 2010 M7.0 Haiti Earthquake and Resulting Stress Changes on Regional Faults[END_REF][START_REF] Saint Fleur | Seismotectonics of southern Haiti: A new faulting model for the 12 January 2010 M7.0 earthquake[END_REF] hereafter denoted the H, M, S, and SF models, respectively). All finite-source models show that rupture occurred on a north-dipping blind fault (Léogâne fault) with two-thirds of the moment released by strike-slip motion and one-third by reverse motion, consistent with up to 60 cm of coastal uplift observed in the epicentral area [START_REF] Hayes | Complex rupture during the 12 January 2010 Haiti earthquake[END_REF]. The NEIC source model, based on teleseismic data only, shows the same coseismic motion partitioning, but on a south-dipping fault with a north-verging reverse component.

Static coseismic displacements from the finite-source models show up to 0.5-1 m of seafloor uplift north of the rupture (Fig. 2) and less than 0.01 m of subsidence along the southern coastline. The F model shows an opposite pattern because of its fault dip opposite to the finite-source models and larger coseismic displacement because of the ad hoc coseismic slip scaling factor applied (see above). This model is inconsistent with the coseismic geodetic observations and the coseismic coastal uplift observed in the epicentral area. It matches observations at the Santo Domingo tide gauge and DART buoy 42407 reasonably well in amplitude and period at low frequencies, but its coseismic slip was arbitrarily tuned to match these observations (Fig. 3).

All finite-source models predict a small tsunami in the epicentral area, with the M and SF models predicting wave heights of up to 1 m along more than 10 km of coastline length, inconsistent with the observations reported by the ITST team [START_REF] Fritz | Twin tsunamis triggered by the 12 January 2010 Haiti earthquake[END_REF]. Indeed, coseismic coastal uplift in the M and SF models is larger and more spread out than the observations reported in [START_REF] Hayes | Complex rupture during the 12 January 2010 Haiti earthquake[END_REF]. All finitesource models predict a local tsunami at Luly, with 0.5-1 m run-up height (Fig. 2), consistent with observations at that location [START_REF] Hornbach | High tsunami frequency as a result of combined strike-slip faulting and coastal landslides[END_REF]. This local tsunami likely occurred as a result of the intersection of the coastline with the sharp edge of the bathymetric low that occupies most of the Bay of Port-au-Prince. Only the F model predicts significant wave heights along the southern coast (Fig. 4), but it also predicts run-up heights of 0.5-1 m in the highly populated and low-lying coastal region surrounding the capital city of Port-au-Prince that were not observed after the earthquake (Fig. 2).

That none of the finite-source models come close to the observations in the near field in Jacmel and Pedernales or in the far field at two locations 450 km apart likely indicates the presence of a secondary tsunami source. An offshore landslide comes to mind, as they are well known to trigger tsunamis even at large distances (e.g., [START_REF] Heinrich | Numerical modeling of tsunami generation and propagation from submarine slumps: the 1998 Papua New Guinea event[END_REF][START_REF] Okal | Source discriminants for near-field tsunamis[END_REF][START_REF] Fritz | Extreme runup from the 17 July 2006 Java tsunami[END_REF][START_REF] Hébert | The 2006 July 17 Java (Indonesia) tsunami from satellite imagery and numerical modelling: a single or complex source?[END_REF]. Fritz el al. (2013) indeed already suspected that this type of source may be required to explain field observations near Jacmel and Pedernales.

Landslide-generated tsunami

We use a high-resolution bathymetry recently acquired offshore of the southern coast of Hispaniola [START_REF] Leroy | Segmentation and kinematics of the North America-Caribbean plate boundary offshore Hispaniola[END_REF] to search for scars in the seafloor morphology that could be indicative of large submarine landslides within the area significantly shaken by the earthquake. The fact that the earthquake source approximately matches observed arrival times at the tide gauge and DART buoy is an indication that a possible landslide source should be located fairly close to the epicentral area, which is a requirement for seafloor acceleration to be sufficient to trigger slope failure.

We find the clearest and largest landslide signature 30 km due south of the tip of the cape that bounds the bay of Jacmel to the west, at a water depth of 3500 m (Fig. 5). The landslide is expressed on the seafloor by two arcuate scars at the slope break between a flat plateau (Haiti Plateau) marked by a sediment wave field and a narrow canyon likely feeding the northeastern corner of the 5000 m deep Haiti sub-basin [START_REF] Mauffret | Neogene intraplate deformation of the Caribbean plate at the Beata Ridge[END_REF]. The scar shows two horseshoeshaped lobes with sharp edges, opening to the east onto the narrow canyon that likely acted as an evacuation pathway for the sediments. Unfortunately, no chirp or high-resolution seismic data are available crossing that feature, which has not been cored either, so that one cannot determine the age of the sediment collapse. We try to identify the bathymetric scar using GEBCO grids, which contain data older than the 2010 Haiti earthquake, but we could not find any evidence for it. This could mean that the scarp is post-earthquake, but it could also mean that the resolution of the GEBCO data is insufficient to detect such a small-scale feature.

We test that landslide as the plausible secondary source necessary to explain the observations in the near field (southern coastline) and far field (tide gauge and DART buoy) of the earthquake. We use the detailed bathymetry (Fig. 5) to re-construct the sedimentary volume involved in the landslide, as in ten [START_REF] Brink | Size distribution of submarine landslides and its implication to tsunami hazard in Puerto Rico[END_REF]. It consists of filling in the failed area according to the adjacent scar height, to reconstruct the seafloor morphology prior to the slide. By using this method, we find a volume of 2 km 3 , a maximum thickness of 150 m, and a slope of 10 • . Once triggered, the failed volume flows with maximum velocities in the 10-65 m s -1 range depending on the assumed viscosity (Fig. 6). The initial acceleration remains approximately the same at ∼ 0.5 m s -2 for all tested dynamic viscosities but it quickly decreases for the dynamic viscosity of 2 × 10 5 Pa s, resulting in smaller amplitudes within the model domain. The influence of the viscosity is evaluated by considering water heights calculated at a synthetic gauge located offshore, a few kilometers north of the landslide (gauge G 1 in Fig. 8). As expected, the larger the viscosity, the smaller the amplitudes of water waves, whereas periods are unchanged (Fig. 7). The largest viscosity results in maximum wave amplitudes below 50 cm at the location of gauge G 1 .

In our preferred model, the southern coast of Haiti is reached by water waves within 5 to 10 min of the landslide triggering. Maximum water heights after 1 h of simulation (Fig. 8) show that water waves are amplified along an arc segment extending from the landslide area towards the Bay of Jacmel. The maximum wave heights reach 2 to 3 m along a coastal segment of 10 km for a viscosity of 2 × 10 5 Pa s (Fig. 8) consistent with ITST observations [START_REF] Fritz | Twin tsunamis triggered by the 12 January 2010 Haiti earthquake[END_REF].

Model results at synthetic gauge G 2 located in the bay of Jacmel show that maximum wave amplitudes range from 3 to 6 m depending on the viscosity of the sliding material (Fig. 9). ITST measurements at Jacmel (see Fig. 5c in [START_REF] Fritz | Twin tsunamis triggered by the 12 January 2010 Haiti earthquake[END_REF], around 3 m, favor a sliding material with a viscosity in the high range of the ones tested here. Due to the local amplification of the bay, the maximum water height at Jacmel occurs 11 min after the underwater slope failure consistent with observations. Further along the coast, the comparison of calculated and observed time series at the Santo Domingo tide gauge is shown in Fig. 10. Data are available with a 1 min sampling rate and the tide gauge is located at the end of a small harbor that is poorly described in the simulation. Nevertheless, the arrival time of the first wave (about 45 min) and the periods and amplitudes of water waves are approximately reproduced by the model. A time lag of 8 min exists between observations and simulations (Fig. 10), which could be explained by imprecisions of the tide gauge clock, a delay of the collapse triggering or imprecisions of the simulation due to errors in bathymetric data.

In the offshore direction, wave propagation is subject to frequency dispersion, which the model takes into account by solving Boussinesq equations. In Fig. 10, we compare the calculated time series with recorded data at DART buoy 42407 with a sampling rate of 1 min. Due to frequency dispersion, the highest water height in the calculated wave train is not the first one but the third, with an amplitude of 0.5 cm. The arrival times as well as amplitudes of the first two waves match the observed signal well. Nevertheless, the oscillations of the following waves (around 3 min) have periods that are shorter than the observed ones (around 4 min), which could indicate that the simulated landslide volume is insufficient to explain the observations.

Discussion

As a first verification, we checked that the landslide tested here, located 70 km from the epicenter, satisfied the magnitude-distance relationships proposed by [START_REF] Salamon | Empirical constraint son magnitudedistance relationships for seismically induced submarine tsunamigenic landslides[END_REF]. Then, we asked whether ground acceleration during the 2010 Haiti earthquake would be sufficient to trigger slope failure at the location of the earthquake described above. We estimate the resonance frequency of the landslide using f = V s /4H , where V s is the shear-wave velocity and H the landslide thickness [START_REF] Parolai | New Relationships between Vs, Thickness of Sediments, and Resonance Frequency Calculated by the H /V Ratio of Seismic Noise for the Cologne Area (Germany)[END_REF]. For V s = 1000 m s -1 and H = 100 m f is 2.5 Hz, a high-enough frequency to justify using peak ground acceleration (PGA) as a proxy for the acceleration that determines slope failure. It is generally recognized that 0.1-0.2 g PGA is the threshold of stability for the triggering of landslides in benthic sediments (e.g., [START_REF] Keefer | Landslides caused by earthquakes[END_REF][START_REF] Meunier | Regional patterns of earthquake-triggered landslides and their relation to ground motion[END_REF]. [START_REF] Tanyaş | Presentation and Analysis of a Worldwide Database of Earthquake-Induced Landslide Inventories[END_REF] show that 80 % of earthquake-induced terrestrial landslides are observed in the PGA interval 0.1-0.8 g. Offshore, sedimentation rates, the state of sediment consolidation and the occurrence of weak layers determine the potential of seismic shaking to increase pore pressure up to the failure threshold. As a result, large PGA does not always trigger mass flow, whereas low PGA can sometimes be sufficient [START_REF] Hampton | Submarine landslides[END_REF][START_REF] Viesca | Nucleation of slip-weakening rupture instability in landslides by localized increase of pore pressure[END_REF][START_REF] Brink | Seismicity and sedimentation rate effects on submarine slope stability[END_REF][START_REF] Pope | Which earthquakes trigger damaging submarine mass movements: Insights from a global record of submarine cable breaks?[END_REF]. For example, the study of Janin et al. (2019) on Holocene sediment collapses along the North Anatolian Fault in the Marmara Sea shows that water circulation induced by seismic shaking in the material can have a significant effect on landslide triggering.

There was no direct measurement of regional ground acceleration during the 2010 Haiti earthquake, but indirect estimates from rigid-body displacements and structural damage infer a PGA value of 0.2-0.4 g in Port-au-Prince [START_REF] Olson | Documenting liquefaction and lateral spreading triggered by the 12 January 2010 Haiti earthquake[END_REF][START_REF] Goodno | Effects of the January 2010 Haitian earthquake on selected electrical equipment[END_REF][START_REF] Hough | Estimation of peak ground acceleration from horizontal rigid body displacement: A case study in Port-au-Prince, Haiti[END_REF], ∼ 30 km east of the earthquake epicenter. A numerical ground motion study based on the coseismic slip distribution from [START_REF] Hayes | Complex rupture during the 12 January 2010 Haiti earthquake[END_REF] finds the mean PGA in Port-au-Prince to be 0.20-0.33 g [START_REF] Mavroeidis | Finite-Fault Simulation of Broadband Strong Ground Motion from the 2010 M w 7.0 Haiti Earthquake[END_REF]. We computed a first-order estimate of the horizontal ground acceleration at the location of the offshore landslide described above using the "Next Generation Attenuation Ground Motions code" [START_REF] Boore | NGA08_GM_TMR, Next generation attenuation ground motions for specified period (T ), magnitude, and distance (R) code, revised version[END_REF] for a M7 earthquake (Fig. 11). We find PGA values on the order of 0.1 g at a distance of 30 km to the earthquake epicenter on soft rock (shear wave velocity near surface Vs30 of 530 m s -1 ). We also used the more precise, but more complex, approach based on a dynamic rupture simulation of the 2010 earthquake on a fault of realistic ge- 2017). This combines the deterministic low-frequency part from the dynamic rupture simulation of [START_REF] Douilly | Three-dimensional dynamic rupture simulations across interacting faults: The M w 7.0, 2010, Haiti earthquake[END_REF] with the stochastic high-frequency one calculated through the specific barrier model [START_REF] Papageorgiou | A specific barrier model for the quantitative description of inhomogeneous faulting and the prediction of strong ground motion. Part I. Description of the model[END_REF] for a crossover frequency of 1 Hz. ometry, as described in [START_REF] Douilly | Three-dimensional dynamic rupture simulations across interacting faults: The M w 7.0, 2010, Haiti earthquake[END_REF][START_REF] Douilly | Simulation of broadband strong ground motion for a hypothetical M w 7.1 earthquake on the Enriquillo Fault in Haiti[END_REF]. Using the same source as these authors, which satisfies observations of static coseismic ground displacement, we find the largest horizontal peak ground acceleration at 0.18 g for the E-W component, while the mean horizontal peak ground acceleration is 0.17 g (Fig. 11). These two ground motion estimates are consistent with each other, though one must keep in mind that these calculations do not account for the more complex near-surface sediment structure likely to be present in the area of the submarine landslide.

Although there is no direct evidence that the submarine landslide described above was triggered by the 2010 Haiti earthquake, it appears that seafloor acceleration in that area was within the range that permits slope failure in marine sediments. It is therefore plausible that a submarine landslide off the southern coast of Haiti is the secondary source necessary to explain, at the same time, the far-field tide gauge and DART buoy observations and the exceptionally large run-up heights at a few locations along the southern coast of Haiti given the limited coseismic displacement there.

A number of similar examples have been documented, often in subduction contexts. For instance, the 1998, M w 7.0 Papua New Guinea earthquake was followed by a tsunami with run-up heights of over 7 m on a 20 km long coastal segment, too large to be accounted for by the coseismic ground displacement alone [START_REF] Okal | Comment on "Origin of the 17 July 1998 Papua New Guinea Tsunami: earthquake or landslide?[END_REF]Tappin et al., 2008). Numerical simulations showed that the tsunami could be attributed to a submarine landslide with a volume of 4 km 3 located 20 km offshore [START_REF] Heinrich | Numerical modeling of tsunami generation and propagation from submarine slumps: the 1998 Papua New Guinea event[END_REF]. Similarly, the 2006 Java tsunami produced run-up heights up to 20 m at Permisan, indicative of a submarine landslide off the Nusa Kambangan Island [START_REF] Fritz | Extreme runup from the 17 July 2006 Java tsunami[END_REF][START_REF] Hébert | The 2006 July 17 Java (Indonesia) tsunami from satellite imagery and numerical modelling: a single or complex source?[END_REF]. In the northern Caribbean, [START_REF] López-Venegas | Submarine landslide as the source for the October 11, 1918 Mona Passage tsunami: Observations and modeling[END_REF] argue that the tsunami associated with the 11 October 1918 M w 7.5 normal faulting earthquake in the Mona Passage between Hispaniola and Puerto Rico was the result of a submarine landslide at a ∼ 2000 m depth. Run-up heights in western Puerto Rico reached up to 6 m about 20 km from the landslide, whose volume was estimated at 10 km 3 , 5 times that of the landslide described here. In a tectonic context similar to the 2010 Haiti earthquake, the 1989, M w 6.9 Loma Prieta earthquake excited tsunamis in the nearby Monterey Bay that required a secondary offshore landslide, assumed to have occurred in shallow waters [START_REF] Ma | The origin of the tsunami excited by the 1989 Loma Prieta earthquake faulting or slumping?[END_REF].

In northern Haiti, the M w 7.6 to 8.0, 7 May 1842 earthquake was followed by a tsunami that killed ∼ 300 people, with wave heights of 2 m reported at a few locations along the northeastern Haitian coast and a maximum run-up of 4.6 m in the city of Port-de-Paix [START_REF] Scherer | Great earthquakes in the island of Haiti[END_REF]. Some far-field effects have been reported, but their reliability remains questionable [START_REF] O'loughlin | Caribbean Tsunamis: A 500-Year History from 1498-1998[END_REF]. Because it is the fastest-slipping active fault closest to the affected area, the Septentrional Fault is commonly assumed to be the source of that earthquake. However, a recent study shows that an 1842-like strike-slip earthquake on that fault cannot explain the tsunami run-up heights observed [START_REF] Gailler | Tsunami scenarios and hazard assessment along the northern coast of Haiti[END_REF]. As a consequence, these authors suggest that event may have occurred on the North Hispaniola Fault, a slow-slipping reverse fault ∼ 50 km offshore. The MMI IX intensities estimated along the northern Haiti coastline for the 1842 earthquake suggest accelerations on the order of 0.3 g [START_REF] Murphy | The correlation of peak ground acceleration amplitude with seismic intensity and other physical parameters[END_REF]. It is therefore plausible that a strike-slip event on the Septentrional Fault would have triggered submarine slope failures -possibly multiple ones -along the steep and near-shore slopes that characterize the bathymetry of the northern Haiti margin [START_REF] Leroy | Segmentation and kinematics of the North America-Caribbean plate boundary offshore Hispaniola[END_REF].

Underwater mass movements have been suggested for more than 100 years to trigger tsunamis in the near field of large earthquakes (e.g., [START_REF] Verbeek | Kort verslag over de aarden zeebeving op Ceram, den 30sten September 1899[END_REF][START_REF] Gutenberg | Tsunamis and earthquakes[END_REF]. They are now considered to be significant contributors to tsunami hazard (e.g., [START_REF] Synolakis | The slump origin of the 1998 Papua New Guinea tsunami[END_REF][START_REF] Brink | Size distribution of submarine landslides and its implication to tsunami hazard in Puerto Rico[END_REF][START_REF] Okal | The 1956 earthquake and tsunami in Amorgos, Greece[END_REF]. It has been noted that tsunamis in strike-slip tectonic regimes were more frequent than expected. Such events may be accompanied by submarine or subaerial landslides (e.g., [START_REF] Ma | The origin of the tsunami excited by the 1989 Loma Prieta earthquake faulting or slumping?[END_REF][START_REF] Imamura | Numerical Simulation of the 1992 Flores Tsunami: Interpretation of Tsunami Phenomena in Northeastern Flores Island and Damage at Babi Island[END_REF][START_REF] Yalçıner | Tsunamis in the Sea of Marmara: Historical documents for the past, models for the future[END_REF][START_REF] Hornbach | High tsunami frequency as a result of combined strike-slip faulting and coastal landslides[END_REF]Hoffman et al., 2014) or associated with a large coseismic deformation [START_REF] Frucht | A Fresh View of the Tsunami Generated by the Dead Sea Transform, 1995 M w 7.2 Nuweiba Earthquake, along the Gulf of Elat-Aqaba, Seismol[END_REF]. That local earthquakes can trigger shallow submarine slope failures is well known and was documented following the 2010 Haiti earthquake, in its epicentral area [START_REF] Hornbach | High tsunami frequency as a result of combined strike-slip faulting and coastal landslides[END_REF][START_REF] Fritz | Twin tsunamis triggered by the 12 January 2010 Haiti earthquake[END_REF]. That deep submarine landslides can also produce significant tsunamis is less clear, though it was for instance proposed for the Puerto Rico, M l 7.5, 1918 earthquake and tsunami [START_REF] López-Venegas | Submarine landslide as the source for the October 11, 1918 Mona Passage tsunami: Observations and modeling[END_REF].

Ten [START_REF] Brink | Mysterious tsunami in the Caribbean Sea following the 2010 Haiti earthquake possibly generated by dynamically triggered early aftershocks[END_REF] recently proposed that the tsunami observations described here were caused by a doublet of dynamically triggered aftershocks with magnitudes of M6.8 and M6.5 located 85 km southwest of the epicenter, with reverse faulting on south-dipping faults. We first note that these magnitudes are significantly larger than the M5.4 and 5.1 quoted in the original publication describing these aftershocks [START_REF] Fan | Local near instantaneously dynamically triggered aftershocks of large earthquakes[END_REF]. We computed the coseismic displacements that such large-magnitude events close to the coast would have caused at GPS sites onshore and found that they would have values (up to 5 cm in their preferred model) and directions (northward) that are not observed in the GPS data [START_REF] Calais | Transpressional rupture of an unmapped fault during the 2010 Haiti earthquake[END_REF]. Finally, our analysis of the detailed bathymetry and recently acquired seismic reflection profiles at the location of the alleged south-dipping reverse fault [START_REF] Leroy | Segmentation and kinematics of the North America-Caribbean plate boundary offshore Hispaniola[END_REF] indicates that they are actually north-dipping. We therefore argue that the explanation proposed by ten [START_REF] Brink | Mysterious tsunami in the Caribbean Sea following the 2010 Haiti earthquake possibly generated by dynamically triggered early aftershocks[END_REF] does not hold in the face of available data.

Conclusion

We have shown that observations of run-up heights up to 2 m along the southern coast of the Haiti at distances up to 100 km from the epicenter, as well as tide gauge and DART buoy records at distances up to 600 km from the 12 January 2010 epicenter, require a secondary source and can be explained by a submarine landslide. We have identified a landslide scar 30 km from the epicenter off the southern coast of Haiti at a depth of 3500 m, where ground acceleration would have been sufficient to trigger slope failure in soft sediments. This candidate source, 2 km 3 in volume, matches observations remarkably well assuming that the sediment collapse obeys a viscous flow with an initial apparent viscosity of 2 × 10 5 Pa s. Although that particular source cannot be proven to have been activated in 2010, our results add to a line of evidence that earthquake-triggered submarine landslides can cause significant tsunamis in areas of strike-slip tectonic regime. This result contributes to explaining why many tsunamis occur along strike-slip fault systems, especially when steep submarine slopes are within reach of significant coseismic ground acceleration, as is the case along most strike-slip fault systems in the Caribbean.

Data availability. The bathymetric data are available on demand [START_REF] Leroy | HAITI-SIS cruise, RV L'Atalante[END_REF].

Author contributions. AP, EC, PH and HH designed the study. AP and PH preformed the tsunami model calculations. SL acquired the bathymetric data, and SL and MR analyzed and interpreted the bathymetric data. HA and RD computed seismic ground motions. AP and EC led the manuscript writing. All co-authors contributed to the interpretation of the results and to the manuscript writing.

Competing interests. The authors declare that they have no conflict of interest.

Figure 1 .

 1 Figure 1. Regional tectonic context. The yellow star indicates the epicenter of the 2010, M w 7.1 Haiti earthquake. Other yellow stars indicate the approximate epicenter of large historical tsunamigenic earthquakes, in particular the 1842 event of northern Haiti, which is discussed in the text. White stars show the approximate location of other large historical earthquakes. Tsunami run-up heights observed after the 2010 event are indicated by squares: < 1 m (yellow), 1 to 2 m (orange), > 2 m (red), draw down (black). The locations of DART buoy 42407 and the Santo Domingo tide gauge are indicated by white triangles. The inset shows the large-scale tectonic framework with current plate velocities relative to the Caribbean shown in centimeters per year.

Figure 2 .

 2 Figure 2. Initial coseismic seafloor displacement shown as dashed lines (red is uplift, blue is subsidence, and contours are shown every 0.05 m) with maximum wave height shown as a colored map in the background in the near field of the 2010 Haiti earthquake. H: Hayes et al. (2010); M: Meng et al. (2012); S: Symithe et al. (2013); SF: Saint Fleur et al. (2015); F: Fritz et al. (2013). Calculations are done at a 30 bathymetric grid resolution. (f) Enlargement of the Port-au-Prince area for the F source model calculated using a 25 m resolution grid.

Figure 3 .

 3 Figure 3. Comparison between observed and simulated water heights at the Santo Domingo tide gauge (a) and at DART buoy 42407 (b).Their locations are provided in Fig.1. The high-amplitude signal on the DART record before ∼ 0.7 h is seismic surface waves.

Figure 4 .

 4 Figure 4. Water height comparison at a synthetic gauge located in the center of the Bay of Jacmel at 30 m water depth. Calculation on a bathymetric grid uses a 30 s resolution.

Figure 5 .

 5 Figure 5. Detailed bathymetry 20 km due south of Jacmel (Fig. 1) where the landslide scar discussed in the text has been identified. (a) Raw observations. (b) Geological interpretation. Swath bathymetric data are from Haiti-sis cruise (https://doi.org/10.17600/12010070,[START_REF] Leroy | HAITI-SIS cruise, RV L'Atalante[END_REF]. The white area corresponds to missing bathymetry of the boat campaign[START_REF] Leroy | Segmentation and kinematics of the North America-Caribbean plate boundary offshore Hispaniola[END_REF].

Figure 6 .

 6 Figure 6. Landslide velocity as a function of time for the dynamic viscosities v tested here.

Figure 7 .

 7 Figure 7. Time series calculated at the synthetic gauge G 1 as a function of time for the dynamic viscosities v tested here.

Figure 8 .

 8 Figure 8. Maximum modeled wave height using a viscosity v = 2 × 10 5 Pa s for the sliding material.

Figure 9 .

 9 Figure9. Time series calculated at synthetic gauge G 2 located on the coast of the Bay of Jacmel (green circle in Fig.8). This calculation uses a bathymetric grid with 25 m resolution.

Figure 10 .

 10 Figure 10. Comparison between observed and simulated water heights at the Santo Domingo tide gauge (a) and at DART buoy 42407 (b) for a dynamic viscosity v = 2 × 10 5 Pa s. High amplitudes observed before 0.6 h are due to the earthquake detected by the DART buoy.

Figure 11 .

 11 Figure 11. Estimated horizontal peak ground accelerations due to an earthquake of M w 7.1 from four ground motion prediction equations: Abrahamson and Silva (2008; blue line), Boore and Atkinson (2008; solid red line with ± 1 standard deviation in dashed lines), Campbell and Bozorgnia (2008; yellow line), and Chiou and Youngs (2008; green line). The error bar symbol shows more accurate estimates derived from hybrid broadband ground acceleration calculation at the site of interest, as proposed in Douilly et al. (2017). This combines the deterministic low-frequency part from the dynamic rupture simulation of[START_REF] Douilly | Three-dimensional dynamic rupture simulations across interacting faults: The M w 7.0, 2010, Haiti earthquake[END_REF] with the stochastic high-frequency one calculated through the specific barrier model[START_REF] Papageorgiou | A specific barrier model for the quantitative description of inhomogeneous faulting and the prediction of strong ground motion. Part I. Description of the model[END_REF] for a crossover frequency of 1 Hz.

Nat. Hazards EarthSyst. Sci., 20, 2055-2065[START_REF] Boore | NGA08_GM_TMR, Next generation attenuation ground motions for specified period (T ), magnitude, and distance (R) code, revised version[END_REF] https://doi.org/10.5194/nhess-20-2055-2020

https://doi.org/10.5194/nhess-20-2055-2020 Nat. Hazards Earth Syst. Sci., 20, 2055-2065, 2020

Acknowledgements. The authors are grateful to colleagues from the Institut de Physique du Globe for sharing their source model for the 2010 Haiti earthquake.

Financial support. This research has been supported by the "Yves Rocard" joint laboratory between the Ecole Normale Supérieure, the Commissariat à l'Energie Atomique and the CNRS as well as by the Interreg Caraibes project PREST (European Union and FEDER) and the Institut Universitaire de France.

Review statement. This paper was edited by Maria Ana Baptista and reviewed by Frederic Dias and two anonymous referees.


Shubham Pande

Shon Yadav

Suresh Balanethiram
email: sureshbalanethi-ram@gmail.com

Deleep R Nair

Sebastien Fregonese
email: sebastien.fregonese@ims-bordeaux.fr

Senior Member, IEEE Thomas Zimmer
email: thomas.zimmer@ims-bordeaux.fr.

Member, IEEE Anjan Chakravorty

S Yadav Is With Globalfoundries

Bangalore S Fregonese

An Efficient Thermal Model for Multifinger SiGe HBTs under Real Operating Condition

Keywords: SiGe HBTs, multifinger transistor, electrothermal effect, thermal modeling, self-heating, thermal resistance

In this work, we present a simple analytical model for electrothermal heating in multifinger bipolar transistors under realistic operating condition where all fingers are heating simultaneously. The proposed model intuitively incorporates the effect of thermal coupling among the neighboring fingers in the framework of self-heating bringing down the overall model complexity. Compared to the traditional thermal modeling approach for an n-finger transistor where the number of circuit nodes increases as n 2 , our model requires only n-number of nodes. The proposed model is scalable for any number of fingers and with different emitter geometries. The model is validated with 3D thermal simulations and measured data from STMicroelectronics B4T technology. The Verilog-A implemented model simulates 40% faster than the conventional model in a transient simulation of a five-finger transistor.

I. INTRODUCTION

S ILICON germanium heterojunction bipolar transistors (SiGe HBTs) are popularly used as power amplifiers in the RF front end modules. Power amplifiers are often expected to have large emitter area in order to allow large amount of currents. However, large emitter widths lead to a higher base resistance resulting in a lower maximum oscillation frequency (f max ) of the device. Therefore, it has been a common practice to partition a large emitter into smaller fingers each having small enough emitter width (W E ) leading to a multi-finger transistor. Although in such a structure, each emitter finger is electrically isolated by shallow trenches (ST), they are thermally coupled through the common Silicon substrate.

Self heating is a serious problem in modern bipolar transistors where lateral dimensions are significantly scaled down and additional trench isolations are used. In case of multi-finger transistors, additional thermal coupling from nearby fingers further increases the device temperature. Conventionally, this additional increment in temperature is captured by considering thermal coupling effects from other nearby fingers. For the calculation of the overall temperature at a finger, thermal effect of each finger is considered at a time and finally all effects are added up assuming the validity of superposition. The temperature dependence of thermal conductivity makes the temperature-power relationship nonlinear. This makes the straightforward application of superposition to include both the self-heating and thermal coupling effects in calculating the overall finger temperature questionable. The real operating condition within a power amplifier circuit exciting all the fingers in a multi-finger transistor together instead of exciting one finger at a time was elaborated in [START_REF] Dwivedi | Effects of BEOL on self-heating and thermal coupling in SiGe multi-finger HBTs under real operating condition[END_REF]. A state-of-theart static thermal model to cater the self-heating as well as thermal coupling effects in an n-finger transistor requires n 2 number of nodes as reported in [START_REF] Walkey | Equivalent circuit modeling of static substrate thermal coupling using VCVS representation[END_REF]. Besides, modeling thermal coupling using voltage controlled-voltage-sources (VCVS) in series with the self-heating resistances degrades the speed performance of a thermal network as illustrated in [START_REF] Balanethiram | Efficient modeling of static self-heating and thermalcoupling in multi-finger SiGe HBTs[END_REF].

In this work, we present an intuitive thermal model to predict the overall temperature at each finger in a multifinger transistor when all the fingers are excited together. Unlike the conventional methods, our model uses no VCVS as the thermal coupling effects are considered within the framework of selfheating. Eventually, the use of superposition can be avoided altogether resulting in just 'n' number of nodes for an n-finger transistor. The proposed model also considers the temperaturedependent thermal conductivity of the semiconductor material and does not need to use any superposition theorem to calculate the overall device temperature at any finger. The paper is organized as follows. Section II presents the elaborate model formulation with a hint towards model implementation. Section III presents a detailed model validation against 3D TCAD simulation and experimental data. The speed performance of the proposed model is also compared with the conventional thermal model. Finally, we present our conclusions in section IV.

II. MODEL FORMULATION

Electrothermal effect in HBTs causes heat generation at the base-collector junction. Most of the generated heat flows down towards the substrate contact because of the high thermal resistance offered by the upward path due to the interlayer dielectric of the back-end-of-line (BEOL) structures as reported in [START_REF] Dwivedi | Effects of BEOL on self-heating and thermal coupling in SiGe multi-finger HBTs under real operating condition[END_REF], [START_REF] Sahoo | Impact of back-end-of-line on thermal impedance in SiGe HBTs[END_REF], [START_REF] Esposito | Dedicated test-structures for investigation of the thermal impact of the BEOL in advanced SiGe HBTs in time and frequency domain[END_REF]. Therefore, the upward heat-flow is neglected in the initial analysis and is added later by adding an effective BEOL thermal resistance in parallel. In the front-endof-line (FEOL) portion, the base and emitter regions are also neglected because of their negligible thickness. This allows us to model the emitter finger as a rectangular heat source on a semi-infinite substrate as shown in Fig 1 . Note that the effect of base and emitter region can always be clubbed with the BEOL thermal resistance. For modeling purpose, it is assumed that the bottom of the substrate is maintained at an ambient temperature (T amb ) and the substrate extends to infinity in the lateral directions. In case of multifinger transistor, multiple heat sources are to be considered simultaneously under real operating condition. The modern application circuits such as power amplifiers are equipped with temperature insensitive bias techniques to ensure a near constant operating current [START_REF] Hua | High-linearity and temperature-insensitive 2.4 GHz SiGe power amplifier with dynamic-bias control[END_REF]- [START_REF] Sub | PCS/W-CDMA dual-band MMIC power amplifier with a newly proposed linearizing bias circuit[END_REF]. Adopting such biasing techniques in a multifinger transistor yields similar amount of collector current (I C ) through all the fingers at a given collector-emitter voltage (V CE ). This leads to an identical amount of power dissipation (P diss ) in each finger (since P diss ≈ I C × V CE ). Accordingly the modeling framework presented in this work assumed that identical amount of P diss is generated at each finger. Fig. 1 depicts the heat spread from a single heat source in a semi-infinite environment. Vertical position (z) dependent temperature variation, T (z), inside the system from the heat source (at a temperature T j ) to the heat sink (at T amb ) is indicated by the isothermal contours. In order to obtain a simplified model for T (z) in such a system, a single heatspreading angle (θ) under conical approximation is used defining the imaginary thermal boundary [START_REF] Holway | Approximate formulas for the thermal resistance of IMPATT diodes compared with computer calculations[END_REF]. If P diss is the power dissipated by the heat source, using average thermal conductivity formulation [START_REF] Balanethiram | Accurate modeling of thermal resistance for onwafer SiGe HBTs using average thermal conductivity[END_REF], the temperature variation along the z-direction (T (z)) can be written as

T (z) = T 1-α ref + P diss f G (z) 1 -α β 1/(1-α) (1)
where T ref = T amb at z = H (signifying a heat-sink at the substrate contact). α, β are the parameters of temperature dependent thermal conductivity in Silicon, κ si (T ) = βT -α [START_REF] Nghiêm | Electro-thermal simulation based on coupled Boltzmann transport equations for electrons and phonons[END_REF]. f G (z) signifies the position dependent geometry factor of the heat spread [START_REF] Yadav | A pragmatic approach to modeling selfheating effects in SiGe HBTs[END_REF]. First, T (z) is evaluated at the heat sink (z = H) with T ref = T amb in [START_REF] Dwivedi | Effects of BEOL on self-heating and thermal coupling in SiGe multi-finger HBTs under real operating condition[END_REF]. Eventually, for calculating T (z -∆z), T (z) from the previous calculation is taken as T ref in [START_REF] Dwivedi | Effects of BEOL on self-heating and thermal coupling in SiGe multi-finger HBTs under real operating condition[END_REF]. This is repeated to obtain the temperature variation for all z in the system from heat sink to the heat source.

The above formulation can be generalized for an asymmetric heat spread of a planar heat source with W as width in the x-direction and L as length in the y-direction. In that case, the position-dependent cross-sectional area of the heat spread can be written as [START_REF] Gao | Emitter ballasting resistor design for, and current handling capability of Al-GaAs/GaAs power heterojunction bipolar transistors[END_REF] A(x, y, z)

= (W + (H -z)Θ x ).(L + (H -z)Θ y ) (2)
with

Θ x = tan(θ x,l ) + tan(θ x,r ) and Θ y = tan(θ y,l ) + tan(θ y,r ) (3
) where H is the total height of the heat spread. Here θ x,l (θ y,l ) and θ x,r (θ y,r ) are the heat spreading angles on the left and right sides along W E (L E ) respectively. Resultant geometry factor can be obtained as

f G (z) = z 0 dz A(x, y, z ) = ln W (L+zΘx) L(W +zΘy) W Θ y -LΘ x . (4)
This equation presents a generic formulation and can be used to estimate the geometry factor of any asymmetric heat spread. Note that (4) need not be used in case of vertical heat spread, where a simple ratio between the height and the cross sectional area yields f G (z). Accordingly, the depth-dependent temperatures (including that at the heat source, z = 0) can be obtained using [START_REF] Dwivedi | Effects of BEOL on self-heating and thermal coupling in SiGe multi-finger HBTs under real operating condition[END_REF]. Note that our aim is to find out the finger temperature i.e., at z = 0 and as such there is no need to compute any other T (z). However, this calculation is exercised in order to showcase the model's capability to capture the underlying physics and accordingly gain a confidence on the proposed model. Estimation of the temperatures becomes difficult for a system with two or more heat sources. Consider a case in which two heat sources (having same area) each dissipating a power of P diss are kept close to each other (with centreto-centre spacing s in between them) along the x-direction as shown in Fig. 2(a). Similar to the single-heat-source structure, total heat spread of the system can be shown by dashed lines with spreading angle θ (which may be equal to 45 • ) [START_REF] Holway | Approximate formulas for the thermal resistance of IMPATT diodes compared with computer calculations[END_REF], [START_REF] David | Computerized thermal analysis of hybrid circuits[END_REF]. For such a system, one can obtain the rise in temperature by considering this overall heat spread and total power dissipation of 2P diss . Alternatively, one can also obtain temperature rise at one finger by defining individual heat spread. As the structure is symmetric, one can intuitively say that T 1 (z) = T 2 (z). Hence, the total heat spread in the x-z plane must be divided into identical halves as shown by the dashed line in Fig. 2(a). The vertical thermal boundary or zero spreading angle from the z-point where the spread of the two heat sources intersect essentially ensures equal heat flow volume in this case. This divided heat-flow volume with power dissipation P diss yields the same finger temperatures. In case of a system with three heat sources as shown in Fig. 2(b), the thermal boundaries of heat source in the middle are governed by two adjacent heat sources. On the other hand, for the first and third heat sources, only one thermal boundary depends on the adjacent (middle) heat source and the other boundary is defined by a heat spreading line of θ. Depending on the spreading angles for individual heat sources, this framework can allow a higher temperature rise at the second heat source compared to the other two. However the first and third heat sources have an equal rise in temperature (T 1 = T 3 ) due to identical boundary conditions. This is also evident from the identical isothermal lines (obtained from TCAD simulation) under the first and third heat sources as shown in Fig. 2(b). Under the second heat source, these lines appear more crowded indicating a higher temperature than the rest.

In order to accurately predict the temperature at each finger, an effective heat spreading angle (θ 1 ) has to be defined between the adjacent heat sources as shown by the dashed lines in Fig. 2(b). Prior to the intersection of the heat sources, the spreading angles of all the heat sources were same θ. If θ 1 is estimated accurately, separate geometry factors corresponding to all heat sources can be obtained. Since all the fingers are dissipating the same amount of power and since the heat flux is continuous, the same amount of heat flux lines must be crossing the x-y plane at the depth of intersection of the flux lines from different fingers. Therefore, estimation of θ 1 can be guided by a principle based on the symmetric total spreading as elaborated below. In case of two-finger device, the first heat spread in Fig. 2(a) is characterized by θ on one side and 0 • on other. A similar approach follows for the second heat spread as well. Therefore, both the spreads have effective spreading angle of θ + 0 • . In case of three-finger device (Fig. 2(b)), heat spread for the first finger has an angle of θ in one side and -θ 1 in the other. Similarly for the second finger, the angle is θ 1 in both the sides. Equating the effective total angles of the two fingers results into θθ 1 = 2 × θ 1 yielding θ 1 = θ/3. Therefore, the essential principle is to find out the total heat spreading angle associated with each finger and to equate them in order to obtain the effective heat spreading between the adjacent heat sources. Extending this technique to a system with n number of fingers, one can obtain a single equation for heat spreading angles as

θ i = θ n (n -2i) i = 1, 2, 3, ... n 2 for even n, i = 1, 2, 3, ... n -1 2 for odd n. ( 5
)
Here θ is the angle of outer heat spreading for the corner fingers. Note that the constant spreading angle framework does not always ensure θ = 45 • as have been reported in several cases [START_REF] David | Computerized thermal analysis of hybrid circuits[END_REF]- [START_REF] Masana | A closed form solution of junction to substrate thermal resistance in semiconductor chips[END_REF]. While testing against TCAD simulation (as presented in the next section) we have found that in most of the geometries θ = 46 • yields excellent accuracy. Once the heat spreading angles are obtained, they are used to compute the cross-sectional area of each heat spread and the geometry factors corresponding to each heating finger using (2) to (4) and subsequently the temperature at each finger using [START_REF] Dwivedi | Effects of BEOL on self-heating and thermal coupling in SiGe multi-finger HBTs under real operating condition[END_REF].

In practice, each transistor finger is to be modeled using separate electrical model where a thermal sub-circuit is available in order to capture the self-heating effect. Conventionally the thermal coupling effect is captured by using a more complicated thermal network that uses voltage-controlled voltage sources [START_REF] Walkey | Equivalent circuit modeling of static substrate thermal coupling using VCVS representation[END_REF]. In the present work, since we have computed the geometry factor (f G ) for each heating finger, the corresponding thermal resistance is easily obtained and can be used within the already existing self-heating network. This way one can get rid of the thermal coupling network altogether and reduce the overall node-count of the thermal network from n 2 (required for conventional model) to just n for an n-finger transistor. The resulting improvement in the simulation speed will be discussed in the next section.

III. MODEL EVALUATION A. Comparison with 3D TCAD simulation

First, we test our proposed model against 3D TCAD thermal simulation results of multifinger SiGe HBT structures having no trench isolation. Using Synopsys Sentaurus [START_REF] Synopsys | Sentaurus: Sentaurus device user guide, release H-2013[END_REF], we simulated a single finger structure and four multifinger structures with two, three, four and five fingers, respectively, having an emitter area (A E ) of 0.2×5 µm 2 for each finger. Following the modeling framework, each heat source of area 0.2×5 µm 2 corresponding to each transistor finger is located at the top of the Silicon substrate. A constant power of 30 mW is dissipated at each heating finger. This simulation scenario is illustrated through a cross-sectional view (x-z plane) of a 5-finger structure in Fig. 3. Each finger is isolated from the neighboring finger with a centre-to-centre spacing of s = 2.5 µm. The total simulated area of the substrate is 300×300 µm 2 with a Silicon thickness of 100 µm. In the TCAD simulation, values for the temperature coefficients of thermal conductivity for Silicon, α = 1.263 and β = 2099 are used. For all the TCAD simulations, the bottom surface is maintained at ambient temperature (T amb =300 K) and the heat source injects a constant uniform heat flux into the silicon substrate. All other surfaces, including the part of the top surface not covered by the heat source, are considered adiabatic.

Fig. 4(a) compares our model results for the z-dependent temperature variation along the middle of the heat source (T (z)) with TCAD simulation data corresponding to a single finger structure and corner (first) fingers in case of multifinger structures, without any trench isolation. Similarly, Fig. 4(b) compares our model results for the T (z) variation with TCAD simulation data corresponding to the internal fingers in multifinger structures without any trench isolation. Note that for optimum model fitting, we have used θ = 46 • for all the structures. Accordingly the spreading angles of internal fingers are calculated using [START_REF] Esposito | Dedicated test-structures for investigation of the thermal impact of the BEOL in advanced SiGe HBTs in time and frequency domain[END_REF]. Excellent model agreement with the TCAD simulated T (z) data for all the structures builds up our confidence in the calculation of all the heat spreading angles corresponding to each finger and subsequently the imaginary boundaries. This motivates us to extend our modeling framework for shallow trench isolated (STI) isolated multifinger transistor structures as detailed below.

In case of multifinger structure with STI, each finger is housed within a STI with a constant finger-edge to trenchedge distance of 0.14 µm along both xand y-directions. Cross sectional view of an ST-isolated 5-finger HBT structure simulated in TCAD is shown in Fig. 5. Corresponding to STMicroelectronics B4T technology, the width (and height) of STI is chosen as 0.45 µm (and 0.36 µm) [START_REF] Chantre | Pushing conventional SiGe HBT technology towards "Dotfive" terahertz[END_REF]. The material chosen for the trench has very low thermal conductivity (SiO 2 with κ = 0.014 W/cm-K) compared to that of the substrate material (Silicon). The dimensions of the substrate, emitter finger(s) and finger spacing remain identical with the previous study carried out for no-trench devices. A constant power of 30 mW is dissipated at each finger for the TCAD simulation as well as for the model. In this case, the heat flow volume of each finger is divided into three regions as shown in Fig. 5. The symmetric trapezoidal volume within STI, followed by the cuboidal volume within STI, and the rest below STI. Note that in order to evaluate the temperature profile in the section below STI, one can apply the proposed model from the bottom of STI to the heat sink in the same fashion as explained in the no-trench case. Temperature profile within STI can be obtained following the technique generally applied for calculating the self-heating temperature since no effect of thermal coupling from other fingers need to be considered within this heat flow volume. A simple depth/area ratio is used to calculate f G (z) inside the cuboidal section and ( 4) is used for the remaining sections. Fig. 6(a) compares our modeling results for the T (z) variation with the 3D TCAD simulation data corresponding to the corner fingers in case of multifinger structures with STI. Thermal boundaries corresponding to the internal fingers from the bottom of the STI are estimated using [START_REF] Esposito | Dedicated test-structures for investigation of the thermal impact of the BEOL in advanced SiGe HBTs in time and frequency domain[END_REF]. Eventually using θ and θ i in (4) and then using (4) in (1) the overall temperature profile is obtained. In Fig. 6(b), we compare our modeling results for T (z) variation with the TCAD simulation data corresponding to the internal fingers in multifinger structures with STI. Excellent model agreement is observed both in Fig. 6(a) and (b). Although our model assumes the STI to be a prefect thermal insulator, TCAD simulation considers a slight heat leakage through the STI due to the non-zero κ of SiO 2 . However, model correlation with TCAD data justifies our approximation of the perfect insulator used for SiO 2 . This

T j2 T j3 T j4 T j5 T j1 θ θ 1 θ 2 θ 2 θ 1 θ STI STI z [µm] x [µm]
T j2 T j3 T j4 T j5 s s/2 θ θ z [µm]
x [µm] observation allows us to consider perfect thermal insulator for all the trenches in subsequent simulation and model. As a next step, we test our model for multifinger transistor system where each finger is individually surrounded by STI and the whole transistor is housed within a deep trench isolation (DTI). Like in previous case, the emitter finger and STI are separated with a constant edge-to-edge distance of 0.14 µm along both xand y-directions. The distance between the inside edge of DTI and the midpoint of the corner emitter finger in the x-direction (along the emitter width) is exactly half of the finger-to-finger spacing (s). However, STI and DTI are located side-by-side in the y-direction (along the emitter length). Cross sectional view of the DTI five-finger structure simulated in TCAD is shown in Fig. 7. The width and height of STI remain unaltered from the previous simulation. The width (and height) of DTI is taken as 0.5 µm (4.26 µm). The dimensions of the substrate, emitter fingers, spacing between the fingers and the amount of power dissipated at each finger are not altered from the previous study.

Here we divide the heat flow volume of each finger into six sections: two inside STI, three in-between STI and DTI (due to different relative distance between STI and DTI in xand y-directions), and one below DTI. Inside DTI, the temperature profile of each finger will be the same as the heat flow volume of individual finger is identical. The corresponding geometry factor f G (z) is evaluated with a symmetric lateral spread of θ (=46 • ) or by a simple depth/area ratio (as applicable in different sections) and eventually the T (z) profile is obtained using (1). For the section below DTI, it can be assumed that the bottom of DTI acts as a single heat source dissipating n × P diss power, where n is the number of fingers. Since the heat flow volume inside the DTI region is equally shared among the fingers, this assumption leads to the prediction of the same temperature at the bottom of DTI and identical T (z) for all fingers. Black line in Fig. 8(a) shows a prediction for our T (z) model under this assumption where TCAD results for different fingers of a five-finger structure with DTI are also plotted. Even though the model predicts the same temperature for all the fingers unlike TCAD simulation, the prediction is quite close. If we closely examine the isothermal lines just below DTI in Fig. 7, we find that the lines are not strictly horizontal, hence the T (x) cannot be the same just below the DTI till certain depth (say, δ). Therefore, we can apply our proposed angle theory conceptualized in [START_REF] Esposito | Dedicated test-structures for investigation of the thermal impact of the BEOL in advanced SiGe HBTs in time and frequency domain[END_REF] to take care of the asymmetric heat spreading within this vertical region (from the bottom of DTI to a depth of δ below DTI). The value of δ is taken as 5 µm and is found to remain unaltered for different finger dimensions and STI-to-DTI distances. A single section (with nP diss total power) is assumed below this region as shown in Fig. 9 in order to evaluate the overall temperature at each finger. Fig. 8(b) compares this refined model for T (z) variation with TCAD simulation data. It is observed that the model is in excellent agreement with the TCAD simulated results.

In order to investigate the model scalability, TCAD simulations of ST and DT-isolated transistor structures with different emitter dimensions and number of fingers were carried out.

B. Validation with Measurement

On-wafer measurement was carried out on a DTI five-finger SiGe HBT. Measurements were done at 300 K using a SUSS MicroTec probing station equipped with a thermal chuck. The emitter fingers of the transistor (each having W E ×L E =0.23× 5 µm 2 ) are electrically isolated but thermally coupled within the DTI. Each one of the five emitters has independent probing access (although all fingers are excited together in this study), the collector is common for all the fingers and the bases are all grounded as detailed in [START_REF] Weiss | Mutual thermal coupling in SiGe: C HBTs[END_REF].

The proposed thermal model is implemented in Verilog-A and is connected as the self-heating sub-circuit to each of the five transistors. HICUM/L2 parameters are used for the electrical network [START_REF] Schröter | Compact Hierarchical Bipolar Transistor Modeling with HICUM[END_REF]. The parameters for the thermal network including the BEOL thermal resistance are extracted following the methodology presented in [START_REF] Balanethiram | Extraction of BEOL contributions for thermal resistance in SiGe HBTs[END_REF]. Circuit simulations are carried out with circuit simulator QucsStudio [START_REF] Margraf | Qucsstudio -a free and powerful circuit simulator[END_REF].

Fig. 11 shows the measured (symbols) output characteristics of a 5-finger SiGe HBT structure fabricated in STMicroelectronics B4T technology [START_REF] Chantre | Pushing conventional SiGe HBT technology towards "Dotfive" terahertz[END_REF]. Solid lines show the circuit simulation results with proposed thermal sub-circuit and dotted lines show the same with self-heating turned-off. Excellent model correlation with the experimental data is observed. In the given constant voltage bias condition, different fingers dissipate different amounts of power because slightly different temperatures across the fingers affect the corresponding collector current differently. In this situation, the present modeling framework tends to overestimate the temperature difference between the fingers. However, the circuit simulation results reveal only 5% difference between the maximum and minimum finger temperatures. The actual difference can be even lower.

Note that for this five-finger transistor, only five self-heating thermal networks were needed, thus avoiding the need of any complicated thermal coupling circuit altogether. In order to quantify the speed improvement of our model over the stateof-the-art thermal model for multifinger transistor [START_REF] Walkey | Equivalent circuit modeling of static substrate thermal coupling using VCVS representation[END_REF], quasistationary and transient simulations of a 5-finger SiGe HBT are carried out for both the models using QucsStudio. HICUM/L2 parameters are used for electrical network of the transistor while only static part of the thermal subcircuit is implemented for both the models. In case of quasistationary sweep of V BE from 0 to 0.9 V with 10000 steps, our model simulates 10% faster than the state-of-the-art. For a transient simulation with V BE = 0.85 V V BE = 0.90 V V BE = 0.95 V Fig. 11. Output characteristics of a 5-finger SiGe HBT structure fabricated in STMicroelectronics B4T technology [START_REF] Chantre | Pushing conventional SiGe HBT technology towards "Dotfive" terahertz[END_REF]: comparison of circuit simulation results with proposed thermal model (solid line) with measured data (symbols). Dotted lines show circuit simulation results with no thermal model. 0.9 V pulse with ns width and 40 ns rise time, our model simulates faster than the state-of-the-art.

IV. CONCLUSION

We have presented a simple, analytical, thermal model for multifinger SiGe HBTs. The proposed model is highly accurate as it considers the temperature dependence of thermal conductivity of silicon and at the same time requires no extra circuit node to account for the thermal coupling effects between nearby fingers. Other than the dissipated power, the input for the model are the dimensions and relative locations of emitter fingers and different trenches in order to compute the temperature at each finger. The model is implemented within the framework of existing self-heating sub-circuit of the main electrical model of bipolar transistor. A high level of model accuracy is observed when compared against 3D TCAD simulation as well as actual electrical measurements. The model is found to simulate 40% faster than the stateof-the-art thermal coupling model while tested for transient simulation.

Fig. 1 .

 1 Fig. 1. Cross sectional view of a single finger bipolar transistor structure with no trench isolation showing heat source, heat sink, isothermal lines and the imaginary boundaries.

Fig. 2 .

 2 Fig. 2. Cross sectional view of (a) two-finger and (b) three-finger transistor structure with no trench isolation showing heat source, heat sink, isothermal lines and the imaginary boundaries.

Fig. 3 .Fig. 4 .

 34 Fig. 3. The imaginary boundary used in the model overlaid on the TCAD simulated temperature profile in the x -z plane of a 5-finger device without any trenches. θ 1 and θ 2 can be calculated using (5). For θ = 46 • , the resulting values are θ 1 = 27.6 • and θ 2 = 9.2 • .

Fig.Fig. 6 .

 6 Fig. The imaginary boundary used in the model overlaid on the TCAD simulated temperature profile in the x -z plane of a 5-finger device with STI. Inset shows the region near STI.

Fig. 7 .

 7 Fig. 7. TCAD simulated temperature contours in the x -z plane obtained by taking a cut-plane across the middle of the fingers for a 5-finger device with both STI and DTI.

10 - 4 Fig. 8 .Fig. 9 .

 489 Fig. 8. T (z) variation at the middle of the fingers in 5-finger structure with STI and DTI: a comparison between TCAD simulation (symbols) and (a) the proposed model with δ = 0 (line); (b) the proposed model with δ = 5 µm.

Figs. 10

 10 (a), (b) and (c) show the variation of finger temperatures with emitter length (for a fixed W E = 0.2 µm), emitter width (for a fixed L E = 5 µm) and number of emitter fingers (for fixed W E = 0.2 µm and W E = 5 µm). Figs.10 (a) and (b) present the results of three different fingers in 5-finger structures, whereas Fig.10 (c) depicts the temperatures of middle and corner fingers for different multifinger devices. The proposed model (solid lines) demonstrates an excellent agreement with the TCAD results (symbols). Maximum modeling error in Fig.10(c) appears to be less than 1.5%.

Fig. 10 .

 10 Fig. 10. Variation of finger temperature (a) with emitter length (for a fivefinger device with fixed W E = 0.2 µm), (b) with emitter width (for a fivefinger device with fixed L E = 5 µm), (c) with number of emitter fingers (for multifinger devices with W E = 0.2 µm and L E = 5 µm): comparison of modeling results (lines) with TCAD simulations (symbols).

This work was supported in part by the EU under Project Taranto, in part by ISRO project ELE/17-18/176/ISRO/ANJA and in part by Department of Science and Technology, India, under Project EMR/2016/004726.The authors would like to thank STMicroelectronics for providing the B9MW wafer.