Pascal Pigeon

Siden Top
email: siden-top@chimie-paristech.fr

Anne Vessières

Michel Huché

Meral Goermen

El Mehdi

Marie-Aude Arbi

Michael J Plamont

Gérard Mcglinchey

Jaouen

Anne Vessie

` Res

Michel Huche´

Meral Go¨rmen

Mehdi El Arbi

Marie-Aude Plamont

Michael J Mcglinchey

Ge´rard Jaouen
email: gerard-jaouen@chimie-paristech.fr

A new series of ferrocifen derivatives, bearing two aminoalkyl chains, with strong antiproliferative effects on breast cancer cells

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The unusual reactivity of bioorganometallic species (complexes of biological interest incorporating a direct metal-carbon bond) is fuelling the emergence of a new organometallic-based approach to medicinal chemistry. [1][2][3][4][START_REF] Jaouen | [END_REF] We have recently illustrated the potential of such a topic by changing the tamoxifen series for ferrocifens. [6][7][8][9] OH-Tam 1, via its pro-drug tamoxifen, is the medication most widely used to cure hormone-dependent breast cancers (Chart 1). [10][11][12][13][14] The main mechanism of action of this molecule is an antiestrogenic effect which occurs at low concentration (in the sub-micromolar range) via a specific interaction with the alpha form of the estrogen receptor (ERa). However non-genomic interactions, i.e. an effect not mediated by ER, have also been described, 15 but in this case the IC 50 value found for OH-Tam on an ER-negative cell line such as MDA-MB-231 is high (29 mM) and not useful for therapeutic purposes. 16 We have prepared a series of ferrocenyl analogues of OH-Tam by substituting the phenyl ring of OH-Tam, and also varying the length of the amino side chain (n = 2, 3-5, 8), and found that the ferrocenyl derivative 2, with n =3 ,i s highly cytotoxic to both hormone-dependent and hormoneindependent breast cancer cells (IC 50 values in the range of 0.5 mM), while ferrocene itself has an IC 50 of 160 mM. [17][18][19] It has been hypothesized that the novel mechanism of action of this Fc-OH-Tam 2 complex could involve the generation of quinone methide owing to a conjugated redox ferrocenyl antenna. [20][21][22] Very recently we have fully characterized these quinone species and shown that they could be key metabolites explaining the different behavior between the organic and organometallic series. 23 In cancer therapy the discovery of a new type of mechanism is generally of interest. 24 Access to quinone methides requires the presence of a phenolic group in the original structure. 20 However, it is also possible to generate these antiproliferative species when the phenol is replaced by an aniline since quinone imines are accessible via the same redox process triggered by the ferrocenyl antenna. 25,26 This poses the question as to the unequivocal character of such a mechanism. What happens if one blocks the generation of quinone methides by protecting the phenolic functionality of Fc-OH-Tam, by incorporating a second -O(CH 2 ) 3 N(CH 3 ) 2 aminoalkyl chain identical to the one already present? This is the question we wished to examine.

The first product tested in this new series was, for chronological reasons, the organometallic complex 3. 7 However, when we found that the ferrocenophanes are often more active than their open chain counterparts, 27,28 we included the cyclic compound 4. Finally, since the ferrocene in these structures can only operate conjugatively with difficulty to activate the redox activation process seen in Fc-OH-Tam, 2, we chose to add the purely organic molecule 5, where the aryl group is lipophilic but less bulky than the ferrocenyl entity. It is the totality of our results in this domain, leading to unexpected biological effects, that we present herein.

Results and discussion

Synthesis of the compounds

The compounds 3, 4 and 5 were obtained by dialkylation of the corresponding phenols 6, 7 and 8, respectively, with dimethylaminopropyl chloride hydrochloride (Scheme 1). The phenols themselves were prepared by the McMurry coupling procedure that we have used profitably in previous syntheses. 7,9,29 The dialkylation of the diphenols 6 and 7 was carried out in DMF using sodium hydride as the base to form the diphenolate which, when heated at 90 1C with dimethylaminopropyl chloride hydrochloride, furnished 3 (20%), starting from 6, and 4 (53%) from 7. Dialkylation of the diphenol 8 was carried out in acetone at reflux using K 2 CO 3 as the base; the isolated yield of 5 was 64%.

Biochemical studies

The study of the effect of the three bis-dimethylaminoalkyl compounds 3, 4, 5 on the growth of hormone-independent (MDA-MB-231) and -dependent (MCF-7) breast cancer cells was performed together with the determination of their relative binding affinity (RBA) for ERa and of their lipophilicity. Results are reported in Table 1 and Fig. 1.

Study of the antiproliferative effect of 3, 4 and 5 on MDA-MB-231 cells. The compounds show a strong antiproliferative effect on MDA-MB-231 which can be attributed to their cytotoxicity. Quite surprisingly these IC 50 values are in the same range, around 0.40 mM, for the three compounds. This is absolutely not the case for their corresponding dimethylamino/ phenol molecules, respectively, 1, 2 (Chart 1) and 9 (Chart 2) which have IC 50 on MDA-MB-231 of 29, 0.5 and 0.015 mM. 16,18,31 For the organic compound 5, replacement of the OH of 1 by a dimethylamino chain increases dramatically the cytotoxicity of the compound (ratio = 85), it has almost no effect on the 3/Fc-OH-Tam couple (ratio = 1.1), and it noticeably decreases in the ferrocenophane series (ratio = 0.04).

Determination of the RBA values of the compounds for the alpha form of the estrogen receptor (ERa) and study of their effect on the growth of MCF-7 cells. All three molecules are recognized by the alpha form of the ER. This is a rather surprising result for these molecules lacking an OH group that is supposed to be one of the essential functional groups allowing the proper anchoring of a molecule in the ER binding site (Table 1). The high RBA value found for 5 was quite unpredictable and the decrease between 5 and 3 and 4 can be attributed to the bulkiness of the ferrocenyl unit. However, the molecular modeling studies (vide infra) provide an explanation for this result.

On MCF-7 cells, the three compounds show a slight, but reproducible, proliferative estrogenic effect at low concentrations (1 Â 10 À8 and 1 Â 10 À7 M; Fig. 1) after 72 h of incubation. In accordance to what is observed on MDA-MB-231 cells, they become highly cytotoxic at higher concentrations (between 1 Â 10 À7 M and 1 Â 10 À6 M). We have previously noted this consecutive dual effect for other series of related complexes albeit possessing diphenolic functions. 28

Molecular modeling

Having seen that the three molecules 3, 4, 5 are recognized by ERa, we decided to use molecular modeling to establish the details of the molecule-receptor interaction. The structure used was that of the ligand binding domain (LBD) of ERa occupied by OH-Tam. 32 Only the amino acids that make up the wall of the cavity were retained. OH-Tam was then removed and replaced successively with 3, 4 or 5.A n energy minimisation was then carried out, with the heavy atoms of the receptor immobilized, in order to establish the optimal position for the molecule within the active site. The Merck Molecular Force Field (MMFF) was employed for this purpose. Owing to the large number of atoms, in excess of the software limit of 600 for quantum mechanical computations, the calculations were carried out using molecular mechanics rather than quantum mechanics. Energy variation values DE corresponding to the binding of the molecules within the active site of the receptor are shown in Table 1. It is interesting to note that these energy values are all very negative, indicating that binding of these molecules in the active site is possible. The value obtained for the organic compound 5 is very close to that of OH-Tam, although less for compounds 3 and 4, and is easily explained by the steric crowding of the ferrocenyl group. We note also that there is good correlation between these values and the RBA values.

A visualization of the docking of 5 and of OH-Tam in the binding site of the estrogen receptor, in the antiestrogen conformation that corresponds to the interaction with OH-Tam, is shown in Fig. 2. This shows clearly that 5 and OH-Tam dock in a similar way inside the active site. The first anchor point of the ligand is an interaction between the nitrogen of the dimethylamino chain and the aspartic acid Asp 351. This is a strong hydrogen bond between the acid function and the terminal nitrogen of the chain, which, thanks to the difference in their pKa, exist in the form of zwitterions, as shown in the following equation:

In the case of OH-Tam the second anchor point is composed of two hydrogen bonds with Arg 394 and Glu 353 which, owing to their proximity, are also found as zwitterions. The first is a bond between the hydrogen atom of the phenol group and Glu 353 (in carboxylate form) and the second is a bond between the oxygen of the phenol group and Arg 394 (in the form of argininium) cf. below. Note that in Shiau's structure, a molecule of water is also present at this level. 32 In the case of 3, 4,or5, the phenol function no longer exists, which eliminates the first hydrogen bond, but the second can continue to exist, between the oxygen of the ether function of the amino chain and the argininium. Owing to a displacement of the molecule within the active site, this is however less strong than with OH-Tam. Nevertheless, the whole structure is stabilized by the positioning of the aliphatic chain -(CH 2 ) 3in a hydrophobic channel situated in this area. Overall, the replacement of the phenol function by the chain costs around 15 kcal mol À1 to the DE of the bond, but does not rule it out. These calculations provide a molecular basis for explaining both the RBA values of these compounds and their estrogenic effect as observed on hormone-dependent MCF-7 cells.

Discussion

We prepared a number of compounds selected from the ferrocifen and tamoxifen series, bearing two aminoalkyl chains of the type O-(CH 2 ) 3 N(CH 3 ) 2 , and found that they showed a strong antiproliferative effect, probably linked to their cytoxicity, on both hormone-dependent and hormone-independent breast cancer cells. In this series the ferrocenyl unit does not seem to play a specific redox activating role but rather behaves as a compact aromatic unit similar to a lipophilic aryl group. Even though a ferrocene is a compact, lipophilic aromatic metallocene (see the log P o/w values of the compounds in Table 1), it is bulkier than a simple organic aromatic. This concept of bioisosterism has recently come to the fore to explain in part the behaviour of ferrocenyl groups in biology. [START_REF] Metzler-Nolte | Ferrocenes: ligands, materials and biomolecules[END_REF] In addition, we found recently that 3 and 5 show the same significant bactericidal and fungicidal effects (on P. aeruginosa, S. aureus and C. albicans). [START_REF] Arbi | [END_REF] Therefore, the driving force of the toxicity of these molecules on cancer cells, bacteria and fungi seems to be connected exclusively to the presence of two aminoalkyl chains. This result differs from that reported for the ferrocenyl-diamine 10. In that case the ferrocenyl complex is active against the bacterium Mycobacterium tuberculosis while its organic equivalent, 11, is 8 times less effective. 35 The toxicity of these compounds seems closer to that found for polyamines derived from putrescine and spermine than to that of our compounds. 36 We note also that cobaltifen 12 (Chart 2) has been described and an IC 50 value of 2.5 mM was found on this same cell line, MDA-MB-231. 37 This value, more than six times higher than that of the compounds described here, can be explained by the presence of a very bulky organometallic group which dilutes the relative influence of the aminoalkyl chains.

Finally, Japanese researchers have synthesized and studied ridaifen-B 13 (Chart 2). 38,39 This compound shows a potent antiproliferative effect on both MDA-MB-231 and MCF-7 cells (IC 50 value of 1.3 mM after 48 h). In addition a COMPARE analysis performed on different cell lines revealed that the mechanism of action of this molecule differs from that of tamoxifen. They also showed that the toxicity of other symmetrical diamines was very similar. 38 This leads us to consider the origin of the significant toxicity demonstrated by this family of molecules.

Since molecules 3, 4, 5 and 12 can be excellent complexing agents, we used molecular modelling to estimate their affinities with Zn 2+ , starting from the idea that this metal is strongly present in the cells and that its displacement, thanks to the presence of basic pincer groups, could lead to their malfunction. AZn 2+ cation was placed between the two amino chains of 3, 4, and 5, and an energy minimisation was then carried out by the use of the Merck Molecular Force Field (MMFF) followed by semi-empirical quantum mechanical calculations (PM3). The D r H1 enthalpy of formation of the resulting complex was obtained, and then in each case, the complex was dissociated without geometric modification into two entities: a Zn 2+ cation and a molecule of 3, 4 or 5. The variations in D r H1 enthalpy of the complexation were calculated from the equation: D r H1 = D r H1(complex) À D r H1(molecule) À D r H1(Zn 2+ ). A D r H1 value around À275 kcal mol À1 was obtained for all these compounds. This result confirms that these molecules are likely to be good complexing agents for divalent cations such as Zn 2+ as well as Ca 2+ . This hypothesis seems to be verified by the formation of a complex between 3 and ZnCl 2 . In fact, addition of one equivalent of ZnCl 2 dissolved in THF to a THF solution of 3 led immediately to the formation of an orange oil at the bottom of the flask. Addition of water to this orange oil produced a yellow solution proving that the complex is partially soluble in water. The 1 H NMR spectrum of the complex in DMSO-d 6 shows a slight deshielding of the protons of the CH 2 CH 2 NMe 2 group (1.83-1.98, 2.46 and 2.48, 2.23 and 2.24 ppm, respectively) compared to that of 3 (1.78-1.94, 2.36 and 2.38, 2.16 and 2.17 ppm, respectively); the remaining protons are unchanged. Other complexing agents for dicationic metals such as hydroxamic acid, present at the end of the SAHA chain, have also been identified. [40][41][42]

Conclusion

We have shown that in the series of the OH-ferrocifens, where the formation of quinone methides can give rise to a strong cytotoxic effect on cancer cells, there may be another possibility of action if the phenol function is blocked by an aminoalkyl chain that prevents the generation of quinone methides. This possibility is not specific to organometallic species since it is conserved when the metallocene is replaced by an organic arene. This means of access to new cytotoxic compounds will require new targets to be identified. This will be the subject of further study, although several possible leads have been suggested here, such as the complexing role of two aminoalkyl chains on metallic acid cations like Zn 2+ or Ca 2+ . We note that these entities combine an antitumoral effect with an antibiotic effect. [START_REF] Arbi | [END_REF]

Experimental General remarks

All reactions took place under argon using standard Schlenk techniques. Anhydrous THF was obtained by distillation from sodium/benzophenone. Thin layer chromatography was performed on silica gel 60 GF254. Infrared spectra were obtained on an IRFT BOMEM Michelson-100 spectrometer equipped with a DTGS detector as a KBr plate. 1 H and 13 C NMR spectra were recorded on a 300 MHz Bruker spectrometer. Mass spectrometry was performed with a Nermag R 10-10C spectrometer. Elemental analyses were performed by the microanalysis service of CNRS at Gif sur Yvette. HRMS measurements were performed by a Thermo Fischer LTQ-Orbitrap XL apparatus with an electrospray source.

1,1-Bis[4-(3-dimethylaminopropoxy)phenyl]-2-ferrocenyl-but-1-ene, 3
Compound 6 (0.424 g, 1 mmol) was dissolved in 10 mL of DMF and sodium hydride (0.32 g, 8 mmol) was added. The mixture was stirred for 10 min, then 3-dimethylamino-1propyl chloride hydrochloride (0.375 g, 2.5 mmol) was added. The mixture was heated at 90 1C overnight. The mixture was cooled and 5 mL of ethanol were slowly added in order to destroy the remaining sodium hydride. Then, the mixture was concentrated under reduced pressure. The residue was dissolved in dichloromethane, washed twice with a diluted aqueous solution of sodium hydroxide followed with water. After drying on magnesium sulfate, the solution was concentrated under reduced pressure and the residue was chromatographed on silicagel with a 4/1 solution of chloroform/ triethylamine as eluent to yield product 3 as an oil in 20% yield. An alternative route of synthesis of 3 was published recently. [START_REF] Arbi | [END_REF] Compound 7 (0.422 g, 1 mmol) was dissolved in 10 mL of DMF and sodium hydride (0.32 g, 8 mmol) was added. The mixture was stirred for 10 min, then 3-dimethylamino-1-propyl chloride hydrochloride (0.375 g, 2.5 mmol) was added. The mixture was heated at 90 1C overnight, then more 3-dimethylamino-1-propyl chloride hydrochloride (0.375 g, 2.5 mmol) was added. The heating was continued for 4 h. The mixture was cooled and 5 mL of ethanol were slowly added in order to destroy the remaining sodium hydride. Then, the mixture was cooled and concentrated under reduced pressure. The residue was dissolved in dichloromethane and was washed twice with a diluted aqueous solution of sodium hydroxide followed with water. After drying over magnesium sulfate, the solution was concentrated under reduced pressure and the residue was chromatographed on silicagel with a 4/1 solution of chloroform/triethylamine as eluent. The residue was recrystallized from ethanol-water solution to yield product 4 in 53% yield. Mp: 119 1C. Compound 8 (1.58 g, 5 mmol) and potassium carbonate (2.764 g, 20 mmol) were stirred in 100 mL of acetone, then 3-dimethylamino-1-propyl chloride hydrochloride (1.739 g, 11 mmol) was added. The mixture was refluxed overnight, cooled and concentrated under reduced pressure. The residue was dissolved in dichloromethane then was washed twice with a diluted aqueous solution of sodium hydroxide followed with water.

After drying over magnesium sulfate, the solution was concentrated under reduced pressure and the residue was recrystallized from ethanol-water solution to yield product 5 in 64% yield. Mp: 91 1C. 1 Compound 3 (0.92 g, 1.55 mmol) was dissolved in 4 mL of THF. A solution of ZnCl 2 (0.21 g, 1.55 mmol) in THF (4 mL) was slowly added. An orange oil was immediately formed at the bottom of the flask. The solvent was removed to leave the Zn complex as an orange oil that was dried under vacuum.

Biochemical experiments Materials

Stock solutions (1 Â 10 À3 M) and serial dilutions of the compounds to be tested were prepared in DMSO just prior to use. Dulbecco's Modified Eagle Medium (DMEM), fetal calf serum, glutamine and kanamycin were obtained from Invitrogen. MCF-7 and MDA-MB-231 cells were from the Human Tumor Cell Bank. Glutamine, 17b-estradiol and protamine sulfate were from Sigma.

Determination of the Relative Binding Affinity (RBA) of the compounds for ERa RBA values were measured on ERa purchased from Pan Vera (Madison, WI, USA). A volume of 10 mL of the solution containing 3500 pmol mL À1 were added to 16 mL of buffer (10% glycerol, 50 mM bis-tris-propane pH 9, 400 mM KCl, 2 mM DTT, 1 mM EDTA, 0.1% BSA) in a silanized flask. Aliquots (200 mL) of this solution in polypropylene tubes were incubated for 3 h at 0 1C with [6,7-3 H]-estradiol (2 Â 10 À9 M, specific activity 1.62 TBq mmol À1 , NEN Life Science, Boston MA) in the presence of nine concentrations of the compounds to be tested. At the end of the incubation period, the free and bound fractions of the tracer were separated by protamine sulfate precipitation. The percentage reduction in binding of [ 3 H]-estradiol (Y) was calculated using the logit transformation of Y (logit Y:l n [ y/1 À Y] versus the log of the mass of the competing steroid. The concentration of unlabeled steroid required to displace 50% of the bound [ 3 H]-estradiol was calculated for each steroid tested, and the results were expressed as RBA. The RBA value of estradiol is by definition equal to 100%.

Culture conditions

Cells were maintained in monolayer culture in DMEM with phenol red/Glutamax I, supplemented with 9% of decomplemented fetal calf serum and 0.9% kanamycin, at 37 1Cina5% CO 2 air humidified incubator. For proliferation assays, cells were plated in 24-well sterile plates with 1.5 Â 10 4 cells for MDA-MB-231 and with 3 Â 10 4 cells for MCF-7 in 1 mL of DMEM without phenol red, supplemented with 9% of fetal calf serum desteroided on dextran charcoal, 0.9% Glutamax I and 0.9% kanamycin, and were incubated for 24 h. The following day (D0), 1 mL of the same medium containing the compounds to be tested diluted in DMSO, was added to the plates (final volumes of DMSO: 0.1%; 4 wells for each condition). After three days (D3), the incubation medium was removed and 2 mL of fresh medium containing the compounds was added. At different days (D3, D4, D5 and D6), the protein content of each well was quantified by methylene blue staining as follows. Cell monolayers were fixed and stained for 1 h in methanol with methylene blue (2 mg mL À1 ), and then washed thoroughly with water. Two millilitres of HCl (0.1 M) were then added, and the plate was incubated for 1 h at 37 1C. Then the absorbance of each well was measured at 655 nm with a Biorad microplate reader. The results are expressed as the percentage of proteins versus the control. Experiments were performed at least in duplicate.

Molecular modeling

Molecular modeling studies were carried out using the programs Spartan, Trident and Odyssey. [START_REF] Spartan | Spartan, Trident and Odyssey[END_REF]

Chart 1 Scheme 1

 11 Scheme 1 Synthesis of the bis-dimethylaminoalkyl compounds 3, 4 and 5.

Fig. 1 2 Fig. 2

 122 Fig. 1 Effect of various concentrations of the compounds on the growth of hormone dependent breast cancer cells MCF-7 after 72 hours in a medium without phenol red. Representative data of one experiment performed 3 times with similar results (six measurements AE limits of confidence).

1 -

 1 [Bis(4-(3-dimethylaminopropoxy)phenyl)methylidene]-[3]ferrocenophane, 4

Table 1

 1 Biochemical data of the compounds: IC 50 values on hormone-independent breast cancer cells MDA-MB-231, relative binding affinity values (RBA) for ERa, enthalpy variation (DE) of the association of the molecules in the ERa binding site and lipophilicity (log P o/w )

	Compound	IC 50 (mM) MDA-MB-231 a	RBA (%) on ERa b	DE (kcal mol À1 )o nE R a	log P o/w	d
	1	29	38.5 c	À151.2	3.2 (Z), 3.4 (E)
	3	0.45 AE 0	2.8 AE 0.1	À89.5	3.56
	4	0.40 AE 0.02	2.05 AE 0.08	À92.5	-
	5	0.34 AE 0.05	54 AE 12	À136.3	3.86

a After 5 days of culture. b Mean of two experiments performed on purified ERa, except for 1 (lamb uterus cytosol); the RBA value of estradiol, the compound of reference, is by definition equal to 100%. c RBA, value from ref. 9. d Measured as described in ref.

30.

  .3 (2 C). IR (KBr, n/cm À1 ): 2948, 2868, 2816, 2764 (CH 2 ,C H 3 ). MS (EI, 70 eV) m/z: 594 [M] + , 121 [CpFe] + , 86 [CH 2 CH 2 CH 2 NMe 2 ] + , 58 [CH 2 NMe 2 ] + . HRMS (ESI, C 36 H 47 FeN 2 O 2 :[ M+H ] + ) calcd: 595.29815, found: 595.29681. The monoalkylated Fc-OH-TAM 2 was also isolated in 36% yield. Anal. Calcd for C 36 H 46 FeN 2 O 2 (H 2 O): C, 70.54; H, 7.80; N, 4.67%. Found: C, 70.58; H, 7.90; N, 4.57%.

	1 H NMR (300 MHz, CDCl 3 ): d 0.94 (t, J = 7.4 Hz, 3H, CH 3 ), 1.79-1.95 (m, 4H, CH 2 ), 2.18 (s, 6H, NMe 2 ), 2.19 (s, 6H, NMe 2 ), 2.37 (t, J = 7.2 Hz, 2H, CH 2 N), 2.40 (t, J = 7.2 Hz, 2H, CH 2 N), 2.50 (q, J = 7.4 Hz, 2H, CH 2 ), 3.83 (t, J = 1.9 Hz, 2H, C 5 H 4 ), 3.89 (t, J = 6.6 Hz, 2H, CH 2 O), 3.91 (t, J = 6.6 Hz, 2H, CH 2 O), 3.98 (t, J = 1.9 Hz, 2H, C 5 H 4 ), 4.02 (s, 5H, Cp), 6.66 (d, J = 8.7 Hz, 2H, C 6 H 4 ), 6.77 (d, J = 8.7 Hz, 2H, C 6 H 4 ), 6.86 (d, J = 8.7 Hz, 2H, C 6 H 4 ), 7.02 (d, J = 8.7 Hz, 2H, C 6 H 4 ). 13 C NMR (75 MHz, CDCl 3 ): d 15.5 (CH 3 ), 27.5 (CH 2 ), 27.6 (CH 2 ), 27.9 (CH 2 ), 45.5 (2NMe 2 ), 56.5 (2CH 2 N), 66.1 (2CH 2 O), 67.9 (2CH C 5 H 4 ), 69.1 (5CH Cp), 69.3 (2CH C 5 H 4 ), 87.2 (C C 5 H 4 ), 114.1 (2CH C 6 H 4 ), 114.2 (2CH C 6 H 4 ), 130.4 (2CH C 6 H 4 ), 130.9 (2CH C 6 H 157

4 ), 136.5 (C), 137.2 (C), 137.3 (C), 137.4 (C),

  1 H NMR (300 MHz, CDCl 3 ):

	d 1.75-1.87 (m, 2H, CH 2 ), 1.87-2.13 (m, 2H, CH 2 ), 2.17 (s, 6H,
	NMe 2 ), 2.21 (s, 6H, NMe 2 ), 2.24-2.31 (m, 2H, CH 2 cycle),
	2.35 (t, J = 7.2 Hz, 2H, CH 2 N), 2.42 (t, J = 7.2 Hz, 2H,
	CH 2 N), 2.54-2.68 (m, 2H, CH 2 cycle), 3.83 (t, J = 6.3 Hz, 2H,
	CH 2 O), 3.89 (s, 4H, C 5 H 4 ), 3.94 (s, 2H, C 5 H 4 ), 3.96 (t, J =
	6.3 Hz, 2H, CH 2 O), 4.13 (s, 2H, C 5 H 4 ), 6.53 (d, J = 8.3 Hz,
	2H, C 6 H 4 ), 6.80 (d, J = 8.3 Hz, 2H, C 6 H 4 ), 6.86 (d, J =
	8.3 Hz, 2H, C 6 H 4 ), 7.05 (d, J = 8.3 Hz, 2H, C 6 H 4 ). 13 C NMR
	(75 MHz, CDCl 3 ): d 27.5 (CH 2 ), 27.6 (CH 2 ), 28.7 (CH 2 cycle),
	41.1 (CH 2 cycle), 45.4 (NMe 2 ), 45.5 (NMe 2 ), 56.4 (2CH 2 N),
	65.9 (CH 2 O), 66.1 (CH 2 O), 68.2 (2CH C 5 H 4 ), 68.6 (2CH
	C 5 H 4 ), 70.1 (2CH C 5 H 4 ), 70.2 (2CH C 5 H 4 ), 84.0 (C C 5 H 4 ),
	86.8 (C C 5 H 4 ), 113.1 (2CH C 6 H 4 ), 114.0 (2CH C 6 H 4 ),
	130.4 (2CH C 6 H 4 ), 131.7 (2CH C 6 H 4 ), 133.1 (C), 135.8 (C),
	136.2 (C), 140.2 (C), 157.1 (C), 157.7 (C). IR (KBr, n/cm À1 ):
	2943, 2856, 2814, 2763 (CH 2 ,C H 3 ). MS (EI, 70 eV) m/z: 592
	[M] + , 86 [CH 2 CH 2 CH 2 NMe 2 ] + , 58 [CH 2 NMe 2 ] + . HRMS
	(ESI, C 36 H 45 FeN 2 O 2 :[ M+H ] + ) calcd: 593.28250, found:
	593.28145. Anal. Calcd for C 36 H 44 FeN 2 O 2 (H 2 O) 0.25 : C, 72.42;
	H, 7.51; N, 4.69%. Found: C, 72.47; H, 7.54; N, 4.56%.
	1,1-Bis[4-(3-dimethylaminopropoxy)phenyl]-2-phenyl-but-1-ene, 5

  H NMR (300 MHz, CDCl 3 ): d 0.84 (t, J = 7.4 Hz, 3H, CH 3 ), 1.71-1.84 (m, 2H, CH 2 ), 1.84-1.98 (m, 2H, CH 2 ), 2.13 (s, 6H, NMe 2 ), 2.18 (s, 6H, NMe 2 ), 2.30 (t, J = 7.3 Hz, 2H, CH 2 N), 2.34-2.47 (m, 4H, CH 2 +CH 2 N), KBr, n/cm À1 ): 2953, 2871, 2813, 2762 (CH 2 , CH 3 ). HRMS (ESI, C 32 H 43 N 2 O 2 :[M+H] + ) calcd: 487.33191, found: 487.33075. Anal. Calcd for C 32 H 42 N 2 O 2 (H 2 O) 0.25 :C , 78.25; H, 8.72; N, 5.70%. Found: C, 78.26; H, 8.53; N, 5.46%.

	3.78 (t, J = 6.4 Hz, 2H, CH 2 O), 3.94 (t, J = 6.4 Hz, 2H,
	CH 2 O), 6.45 (d, J = 8.8 Hz, 2H, C 6 H 4 ), 6.67 (d, J = 8.8 Hz,
	2H, C 6 H 4 ), 6.79 (d, J = 8.7 Hz, 2H, C 6 H 4 ), 6.93-7.13 (m, 7H,
	C 6 H 5 +C 6 H 4 ). 13 C NMR (75 MHz, CDCl 3 ): d 13.6 (CH 3 ),
	27.5 (CH 2 ), 27.6 (CH 2 ), 29.0 (CH 2 ), 45.2 (NMe 2 ), 45.5 (NMe 2 ),
	56.4 (CH 2 N), 56.5 (CH 2 N), 65.9 (CH 2 O), 66.1 (CH 2 O), 113.3
	(2CH C 6 H 4 ), 114.0 (2CH C 6 H 4 ), 125.9 (CH C 6 H 5 ), 127.8
	(2CH arom ), 129.7 (2CH arom ), 130.6 (2CH arom ), 131.9 (2CH arom ),
	135.7 (C), 136.2 (C), 137.9 (C), 140.9 (C), 142.7 (C), 156.9 (C),
	157.7 (C). IR (Complexation of 1,1-bis[4-(3-dimethylaminopropoxy)phenyl]-2-
	ferrocenyl-but-1-ene, 3, with ZnCl 2

Acknowledgements

We thank Barbara McGlinchey for helpful discussions, and the Agence Nationale de la Recherche (ANR 2010BLAN706 1 blanc, ) Mecaferol *) for financial support. M. G. acknowledges a fellowship from the Ministe`re des Affaires Etrange`res (France).


Markus Muller

Philippe Dollfus

Michael Schroter

Markus Müller

Michael Schröter

-D Drift-Diffusion Simulation of Two-Valley Semiconductors and Devices

Keywords: III-V semiconductors, drift-diffusion, charge carrier transport, negative differential mobility, TCAD, HBT. 1D Drift

published or not. The documents may come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

Transistors based on III-V semiconductors provide high transit frequencies simultaneously with high breakdown voltages, owing to their higher low-field mobility and larger bandgap compared to silicon based transistors. The development of both process technology along with device design optimization and adequate physics-based compact models for technology deployment may benefit tremendously from numerical device simulation (a.k.a. technology computer aided design or TCAD) under the following conditions. (i) The associated TCAD tools provide a reasonably accurate prediction of the impact of structural changes on the electrical device characteristics. (ii) The computational effort enables obtaining useful results within reasonable time. Meeting these requirements is becoming increasingly difficult with the reduction of critical device dimensions and the use of materials that exhibit additional physical effects compared to silicon. Both of these attributes apply to III-V semiconductors and are the main reason for the lack of suitable TCAD tools in that field.

In contrast to silicon based devices, III-V based devices employ materials in which electron transport in more than one conduction band is relevant [START_REF] Butcher | Intervalley transfer of hot electrons in gallium arsenide[END_REF][START_REF] Ishibashi | Nonequilibrium electron transport in HBTs[END_REF][START_REF] Palankovski | Simulation of heterojunction bipolar transistors on Gallium-Arsenide[END_REF]. Transport in (typically) two valleys with quite different mobility is usually accounted for using the so-called single electron gas (SEG) approximation with a modified electron velocity-field equation that accounts for the negative-differential-mobility (NDM) effect [START_REF] Palankovski | Simulation of heterojunction bipolar transistors on Gallium-Arsenide[END_REF][START_REF] Gomila | Stationary states and phase diagram for a model of the gunn effect under realistic boundary conditions[END_REF][START_REF] Bonilla | The Gunn effect: Instability of the steady-state and stability of the solitary wave in long extrinsic semiconductors[END_REF]. This approach has been used in moment-based transport equations like the drift-diffusion (DD) and hydrodynamic transport formulations. Due to their computational efficiency, those formulations are being widely used by device simulators (e.g. [START_REF]Synopsis Sentaurus Device User Guide[END_REF][START_REF]Minimos-NT[END_REF][START_REF][END_REF]) in industry and academia. It has been shown though that the SEG approach does not provide satisfactory results for III-V devices for several reasons [START_REF] Shao | What silicon modelers should know about III-Vs for TCAD[END_REF][START_REF] Wedel | On the use of Drift-Diffusion and Hydrodynamic Transport Models for Simulating the Negative Differential Mobility Effect[END_REF][START_REF] Li | Physical modeling of degenerately doped compound semiconductors for high-performance HBT design[END_REF][START_REF] Zampardi | Industrial view of III-V devices compact modeling for circuit design[END_REF]. The NDM formulation of the mobility leads to severe convergence problems, prohibiting simulations in the operating region of interest for circuit applications [START_REF] Tatum | Simulation of Hot-Electron Effects with Multi-band Semiconductor Devices[END_REF][START_REF] Wedel | On the use of Drift-Diffusion and Hydrodynamic Transport Models for Simulating the Negative Differential Mobility Effect[END_REF][START_REF] Schroter | A Closed Form Solution for the Low-Current Collector Transit Time in Group IV and Group III-V HBTs[END_REF]. Moreover, the SEG approach is questionable due to its physical basis, which has extensively been discussed in [START_REF] Blotekjaer | Transport equations for electrons in two valley semiconductors[END_REF]. Fundamentally, such issues can be overcome by using the Boltzmann transport equation (BTE) (e.g., [START_REF] Tomizawa | Numerical Simulation of Submicron Semiconductor Devices[END_REF][START_REF] Wedel | A Box-Integration/WENO solver for the Boltzmann Transport Equation its Application to HighSpeed Heterojunction Bipolar Transistors[END_REF][START_REF] Hong | Deterministic solvers for the Boltzmann transport equation[END_REF]) or more complicated quantum transport formalisms (e.g., [START_REF] Wen | A multiscale TCAD approach for the simulation of InP DHBTs and the extraction of their transit times[END_REF]).

However, these approaches are not suitable for practical applications due to their far too large computational effort and resulting simulation times. Therefore, a self-consistent solution of moment-based transport equations is highly desirable due to their far lower computational cost and parameters that are easier to adjust to measured device characteristics. Despite the fact that their physical basis at small length scales becomes questionable, it has been shown that even the DD approach can still be employed if the mobility is properly adjusted [START_REF] Lundstrom | Carrier transport in BJTs: From ballistic to diffusive and off-equilibrium[END_REF]. This insight is very beneficial for maintaining a basic analytical understanding of the major transport effects in a device. In particular, the simplicity of the DD transport formalism preserves the bridge between BTE and physics-based compact models for circuit design, since all mainstream compact models for FETs and HBTs are based on 1D DD transport (e.g. [START_REF] Schröter | Compact Hierarchical Bipolar Transistor Modeling With HiCUM[END_REF][START_REF] Chauhan | BSIM Industry standard compact MOSFET models[END_REF]).

It is thus desirable to overcome the limitations of the SEG approximation within the DD transport formalism. In [START_REF] Tatum | Simulation of Hot-Electron Effects with Multi-band Semiconductor Devices[END_REF], a two-valley DD approach was proposed that uses a carrier temperature calculated from electrothermal considerations but still assumes a SEG. It has already been shown that two-valley hydrodynamic transport with two separate electron gas systems may be used to accurately reproduce BTE results [START_REF] Tomizawa | Transport equation approach for heterojunction bipolar transistors[END_REF]. However, the system in [START_REF] Tomizawa | Transport equation approach for heterojunction bipolar transistors[END_REF] requires a large number of parameters and is, in our opinion, not practical for the analysis and optimization of transistors that have both a spatially dependent material composition and doping density, such as heterojunction bipolar transistors (HBTs). In this work a two-valley (2v) formulation of 1D DD transport with two separate electron gas systems is proposed as a compromise between physical rigor and suitability for practical use of TCAD for device design, optimization and modeling. Section II introduces the rationale and physical basis for the 2v-DD equation system, which is then formulated in section III. The application of the new 2v-DD model is demonstrated in section IV.

The DD simulations shown herein have been obtained with the in-house simulator COOS [START_REF] Claus | COOS: A wave-function based Schrödinger-Poisson solver for ballistic nanotube transistors[END_REF], which has been extended for the simulation of III-V compound semiconductor heterojunction bipolar transistors with more than a single valley. All simulations herein have been conducted at 300 K.

II. FUNDAMENTALS

The BTE can be used to model carrier transport on a microscopic scale in bulk semiconductors and devices [START_REF] Tomizawa | Numerical Simulation of Submicron Semiconductor Devices[END_REF][START_REF] Wedel | A Box-Integration/WENO solver for the Boltzmann Transport Equation its Application to HighSpeed Heterojunction Bipolar Transistors[END_REF]. The bulk velocity field characteristics of GaAs, InP, GaSb and InAs obtained from the BTE are shown in Fig. 1 along with measured data. For a semiconductor with several valleys one may use the method of moments to derive a system of equations from the BTE [START_REF] Blotekjaer | Transport equations for electrons in two valley semiconductors[END_REF]. Conservation of particles is expressed as

, ( 1
)
where x is the transport direction, t is the time, n j is the carrier density and v j is the carrier velocity in valley j. The collision term in (1), denoted with subscript c, is of crucial importance for the simulation of many-valley semiconductors. It describes the coupling between the valleys due to scattering, including thermal carrier generation and recombination.

For the two valley (j = , L) equation system, numerical solutions based on the hydrodynamic and energy-transport approximations outlined in [START_REF] Blotekjaer | Transport equations for electrons in two valley semiconductors[END_REF] have been obtained in [START_REF] Tomizawa | Transport equation approach for heterojunction bipolar transistors[END_REF][START_REF] Grasser | A review of hydrodynamic and energy-transport models for semiconductor device simulation[END_REF]. However, to the best of our knowledge, such formulations have

E(kV cm -1 ) v(10 7 cm s -1 ) (c) (a) (b) (d) dn j dt ------- d v j n j   dx ----------------- + dn j dt -------     c =
not been implemented in any open-source or commercial TCAD simulator. This may be related to practical issues. First, the implementation and numerical effort is quite large. Second, the accuracy compared to measurements is not necessarily improved when more moments of the BTE are included, since more moments require more parameters that also need to be accurately determined [START_REF] Schröter | SiGe HBT technology: Future trends and TCAD-based roadmap[END_REF][START_REF] Korn | Device simulation of high-performance SiGe heterojunction bipolar transistors[END_REF]. If predictive device simulations are required, the BTE should be solved [START_REF] Schröter | SiGe HBT technology: Future trends and TCAD-based roadmap[END_REF].

III. MODEL FORMULATION

A. Two-Valley Drift Diffusion Model

To reduce the computational effort, the analysis is restricted to just the most important valleys. For the practically relevant III-V semiconductors, transport in the  and the L valley dominates for electric fields typically encountered in device and circuit applications. To be on the safe side, the BTE simulations used as references in this work also include the X valleys. The L valley results shown herein are the averages over the simulated L and X valleys.

From (1) and j = , L, one gets two continuity equations ,

Here J nj is the current density, q is the elementary charge, µ nj is the carrier mobility,  nj is the quasi-Fermi potential (QFP) and w j is the carrier energy in valley j. This formulation of the DD current density is better suited for the analysis of heavily doped devices compared to Blotekjaer's formulation and can also be derived directly from the BTE [START_REF] Marshak | Electrical current in solids with position-dependent band structure[END_REF]. For a 2v-DD equation system the collision terms can be written as

where j =  when k = L and vice-versa. The functions f are the intervalley transition frequencies that describe the intervalley scattering processes as a function of w j and w k . The collision terms cause the energy dependence of the carrier populations among the valleys. The intervalley transition frequencies can be obtained from BTE simulations as explained in [START_REF] Tomizawa | Numerical Simulation of Submicron Semiconductor Devices[END_REF].

Since the carrier energies are not directly available in a DD formulation,  nj and f jk are modeled as functions of the absolute value of a driving force F j . This force needs to be calculated from local quantities and is used to estimate the carrier energy in each valley. For a homogeneous bulk semiconductor the carrier energy is a function of the absolute value of the electrical field E [START_REF] Lundstrom | Fundamentals of Carrier Transport[END_REF]. This allows to find empirical model equations for  nj and f jk as functions of E from bulk BTE simulations. In non-homogeneous structures other fields are more suitable as driving forces, as will be discussed in Sec. IV.

B. Intervalley Transfer Rate Model

According to Fig. 2(a), f L shows a distinct dependence on the electric field, which is modeled by

dn j dt ------- 1 q --- dJ nj dx ---------- dn j dt -------     c + = J nj q nj w j  n j d nj dx ----------- - = n j d t d -------     c n j f jk w j   - n k f kj w k   + = , ( 5
)
where a f is a model parameter and E fL is a normalization constant (= 1 kV/cm). The critical field for the onset of intervalley transfer E ivt is modeled doping dependent as

, ( 6
)
where E ivt,0 , N ivt and  ivt are model parameters. In equilibrium, (4) yields the relationship , [START_REF]Synopsis Sentaurus Device User Guide[END_REF] with n L),equ as the equilibrium carrier concentration in the (L) valley. The parameter f L,equ in ( 5) is implicitly determined by the band structure, [START_REF]Synopsis Sentaurus Device User Guide[END_REF] and by the model for f L , which is introduced next. The results in Fig. 2(b) indicate that f L is almost independent of doping but a linear function of the field, hence .

(

The field dependence of f L has only minor influence on the distribution of carriers amongst the valleys, since the dynamic range of f L is much larger. The noise observed in the simulated f L results from the very small amount of carriers in the L valley at low fields.

C. Velocity Formulation

The carrier velocity in each band as a function of the electric field is shown in Fig. 3 for GaAs. The velocity field curve for carriers in the  valley shows a slight NDM effect for doping concentrations below about 5x10 16 cm -3 [START_REF] Arabshahi | Temperature and doping dependencies of electron mobility in InAs, AlAs and AlGaAs at high electric field application[END_REF][START_REF] Maloney | Transient and steady-state electron transport properties of GaAs and InP[END_REF], which has not been included in this work for the following reasons. (i) The low impurity concentrations where this effect occurs are of little relevance in advanced high-speed HBTs. (ii) For field values where this region would be relevant, the amount of carriers in the  valley is small. Neglecting this effect causes minor deviations between BTE and the 2v-DD model at high fields. 1.4x10 16 , 2.2x10 16 , 3.3x10 16 , 5.2x10 16 , 8.0x10 16 , 1.1x10 18 , 1.7x10 18 , 2.7x10 18 , 4.2x10 18 , 6.4x10 18 ]cm -3 . Comparison between BTE (symbols) and analytical approximations (lines) for (a) f L with ( 5) and (b) f L with [START_REF]Minimos-NT[END_REF]. The mobility model proposed herein has separate parameters for each valley j. It consists of a low-field mobility due to lattice and impurity scattering [START_REF] Caughey | Carrier mobilities in silicon empirically related to doping and field[END_REF] ,

f L E   f L,equ a f 1 E E ivt - E fL ---------------------     exp +     ln + = E ivt E ivt,0 1 N I N ivt ---------      ivt + = f L E 0 =   n L,equ n ,equ --------------f L E 0 =   = f L E   f L0  f E + = 0 5 10 15
where  j,L is the lattice mobility,  j,min is the mobility at high impurity concentrations and α j and N j,ref are model parameters.

Finally, the mobility field dependence is given as [START_REF] Shao | What silicon modelers should know about III-Vs for TCAD[END_REF] with v j,sat as the saturation velocity and β j as a model parameter. This mobility model is available in this or a similar form in most TCAD simulators.

IV. MODEL APPLICATION

A. Steady-State Bulk Velocity-Field Characteristics

The average carrier velocity in a homogeneous bulk semiconductor with two valleys is , [START_REF] Wedel | On the use of Drift-Diffusion and Hydrodynamic Transport Models for Simulating the Negative Differential Mobility Effect[END_REF] where n  +n L = N I . Under steady-state conditions the l.h.s. of (4) equals zero, resulting in

, ( 12
)
which can be evaluated using ( 5) and ( 8) [START_REF] Mccumber | Theory of negative-conductance and of Gunn instabilities in two-valley semiconductors[END_REF]. The parameters for the model equations of f L , f L ,  nL and  n have been obtained by fitting ( 5), ( 8) and [START_REF] Shao | What silicon modelers should know about III-Vs for TCAD[END_REF] to BTE simulation results. The bulk velocity field characteristics are then obtained using [START_REF] Wedel | On the use of Drift-Diffusion and Hydrodynamic Transport Models for Simulating the Negative Differential Mobility Effect[END_REF] and [START_REF] Li | Physical modeling of degenerately doped compound semiconductors for high-performance HBT design[END_REF]. The comparison of the 2v-DD model with BTE data is shown in Fig. 4. Quite good agreement, including the NDM effect, is observed over a wide range of doping levels 1 .

The model parameters are summarized in Appendix B. ----------------------------- ------------------

E(kV cm -1 ) v L (10 7 cm s -1 ) (a) (b) N I N I  j LI   j min   j L   j min  - 1 N I N j ref  -------------      j + -------------------------------- + =  nj  j LI  1  j LI  E v j sat  ------------------      j + 1  j  ---------------------------------------------------- = v n L v L n  v  + n L n  + -
= n L n  ------ f L E   f L E   -

B. Driving Force and Application to Test Structures

In this section different choices for the driving force F j , labeled (i)-(iii), are discussed using two n+-n-n+ test structures as examples. Common choices for driving forces are the gradients of the QFPs and the electric field. Approach (i) and (ii) are defined as:

(i) (ii) (i) is often used in DD simulations and has been proven useful in practice [START_REF] Gokhale | Numerical solutions for a one -dimensional silicon n-p-n transistor[END_REF]. Using (ii), the local energy in valley j is estimated from d(W j0 -W j )/dx, where W j is the band edge of valley j and the index 0 indicates its equilibrium value. Note that for a homogeneous bulk semiconductor, d(W j0 -W j )/dx is equal to the electric field 2 .

BTE results of the mean carrier energy, along with the gradients for (i) and (ii) are shown in Fig. 5. The QFPs in Fig. 5 have been extracted from the BTE assuming ,

1. Due to the lack of published experimental data over a sufficiently wide doping concentration range (GaAs, InP) or at all (GaSb, InAs) BTE data have been used as reference for calibrating the DD model. 2. In the general 3D case the d/dx operators in (i) and (ii) need to be replaced with . )

J nj J  nj  
N I N I N I N I F j d nj dx ----------- = F j d dx ------W j0 W j -   = n j n ij   nj - V T ------------------     exp =
as would be valid for DD transport; n ij is the effective intrinsic carrier density of valley j,  is the electrostatic potential and V T is the thermal voltage. All fields show a similar behavior and a significant deviation from the shape of the mean carrier energy. Thus, at least for this structure, the local fields are unsuitable for estimating the non-local energy acquisition. In Fig. 6 simulations with (i) and (ii) are compared to measurement data and BTE simulations. In both cases, the velocity saturation in the  valley had to be turned off to achieve convergence with (i) and (ii). As expected, both choices are unable to reproduce the BTE results satisfactorily and, thus, to describe non-local effects within the DD transport framework.

To overcome this limitation, another driving force needs to be found. This driving force is used in approach (iii) and will be introduced, next.

n Γ / (n L + n Γ ) (a) (b) (c)
Before (iii) is derived, the principal properties between energy and electric field are discussed based on the BTE results in Fig. 7. Regardless of the bias level, the energy peaks on the right side of the n-layer. Even though the field is almost constant for low bias, the energy shows a spatial dependence. At high bias, the mean energy tends to saturate around 0.12 eV, while the field is still significantly increasing.

To account for the non-local energy-field dependence, the local field needs to be transformed into an effective local field such that the local energy is estimated correctly. As shown in Appendix A, a simplified energy balance equation for the  valley reads .

(

) 14
Here,  w is the energy relaxation time, T 0 is the lattice temperature and T  is the average carrier temperature in the  valley. Near steady-state, ( 14) allows defining an effective local field , where   is a parameter related to the energy relaxation time and the saturation velocity (see Appendix A). In regions where carrier heating is relevant, but intervalley transfer is not, a rough solution of ( 14) may be obtained by assuming that the field is spatially constant and that carriers travel at saturation velocity in the transport direction x. Then, the solution of ( 14) is [START_REF] Slotboom | Non-local impact ionization in silicon devices[END_REF] , ----------------- x ( m)

d dx ------ 5 2 --- k B T  q ------------J n -     J -n d n dx ------------ 3 2 ---k B n  T  T 0 -  w ------------------ - = d n,eff dx - (16)
- 5k B 2q  -------------T  x d n dx ------------      T 0 -     = 0 0.2 0.
w (eV) (a) (b) V V T  x d  n dx ---------------      T 0 2q 5k B --------- d n dx ------------ exp  x -   -----------     d x b x  + =
where   is used as a fitting parameter and x b is the first point in the carrier heating region.

In n + regions, the gradient of the quasi-Fermi potential is negligible and hence T  equals T 0 . In the n region ( 16) can be used to estimate T  as long as intervalley transfer is negligible. When T  becomes sufficiently large, intervalley transfer becomes relevant and an analytical solution for T  is difficult to obtain. From the BTE solutions in Fig. 7(b) one can deduce that the energy, and thus also T  , is saturating beyond a certain critical field. Hence, a suitable driving force is , (iii-a)

where E crit and  are used as fitting parameters and ( 15)-( 16) are used to calculate |d neff /dx|. For the L valley it is assumed that , (iii-b) which overestimates the L valley energy. This error is not critical, since carriers transferred to the L valley travel near saturation velocity anyway. Furthermore, f L is only weakly dependent on F L . Equations (iii-a) and (iii-b) define approach (iii). Approach (iii) leads to quite good agreement for the electrostatic potential and carrier densities as shown in Fig. 6(b), (c). Also, the terminal current density exhibits qualitatively better agreement with BTE and experimental data at high fields (see Fig. 8(a)), but is now somewhat too large. With a proper adjustment of the L valley saturation velocity, the agreement can be improved further. The proposed use of the effective field in (iii-a) provides a practical approach towards capturing negative different mobility effects in a DD transport framework. The simulation results for an InP n+-n-n+ structure, shown in Fig. 8, further underline the suitability of this approach to capture non-local effects.

C. Application to GaAsSb/InP HBT Fig. 9 depicts simulation results for a GaAsSb/InP HBT [START_REF] Lövblom | Developement of Sub-Millimeter-Wave Inp/ GaAsSb Double Heterojunction Bipolar Transistors[END_REF], using the 2v-DD equation system. In addition to the two electron continuity equations (3) and the Poisson equation, a hole continuity and transport equation needs to be solved. Herein, the mobility model and corresponding parameters from [START_REF] Palanovski | Simulation of Heterojunction Bipolar Transistors[END_REF] are used. At the abrupt heterojunction, a special boundary condition is required [START_REF] Palanovski | Simulation of Heterojunction Bipolar Transistors[END_REF]. In the middle of the neutral base the hole quasi-Fermi potential is forced to equal the applied base voltage [START_REF] Schröter | A compact physics-based large-signal modelfor high-speed bipolar transistors with special consideration of high current densities and two-dimensional effects[END_REF].

To estimate the non-local energy, approach (iii) is used, where the collector is treated in the same way as the n region before. The material properties of the GaAsSb alloy are calculated using the parameters of its constituents and assuming bowing parameters equal to zero. The solution shows similar trends as the results in [START_REF] Wedel | A Box-Integration/WENO solver for the Boltzmann Transport Equation its Application to HighSpeed Heterojunction Bipolar Transistors[END_REF][START_REF] Tomizawa | Transport equation approach for heterojunction bipolar transistors[END_REF]. In particular, the L valley carrier density in the BC space-charge-region increases (s. Fig. 9(b)) and exceeds the valley carrier density with increasing bias. The impact of intervalley transfer on the transit frequency ------------------------------ ---------= can be seen in Fig. 9(c), where simulation results without intervalley transfer are shown for comparison. As the current density increases, more carriers get transferred to the L valley, decreasing the transit time due to the lower velocity in that valley. A calibration of the model to measured data for this transistor is left for future work. V. CONCLUSION A two-valley DD formulation for the simulation of III-V semiconductors has been developed. The transition of charge carriers between the two bands is described by a model for intervalley transfer that employs an effective electric field in order to capture the non-local energy acquisition of carriers. As opposed to standard single-valley transport models employed in most TCAD simulators, the new 2v-DD transport formulation captures the negative differential mobility effect in a natural way while at the same time maintaining computational efficiency and limiting the number of model parameters.

F  d n,eff dx  1 d n,eff dx  E crit -
-      + 1   -------------------------------------------------------------- = F L d L dx -
The parameters for widely used III-V materials have been obtained from bulk BTE simulations at room temperature. Examples for the application of the new 2v-DD model have been presented. Future investigations will focus on the application of the two-valley DD model to III-V HBTs and compact modeling support. Also, it would be useful to determine bowing parameters for III-V alloys.

Appendix A: Carrier temperature The energy-balance equation is given by ( 6) in [START_REF] Wachutka | Unified framework for thermal, electrical, magnetic, and optical semiconductor device modeling[END_REF]. It is also valid also for III-V semiconductors as long as intervalley transfer, and thus the associated energy flux, can be neglected. With

n Γ / (n Γ + n L ) 0 1 2 3 0 0.2 0.4 0.6 V (V) J mA μm -2 0 1 2 0 1 2 3 x (μm) Ψ (V) (b) (a) (c)
the assumptions given in [START_REF] Chang | Simplified energy-balance model for pragmatic multi-dimensional device simulation[END_REF] and using (3), one can derive [START_REF] Schroter | A Closed Form Solution for the Low-Current Collector Transit Time in Group IV and Group III-V HBTs[END_REF]. Assuming dJ n /dx = 0 and carriers traveling at saturation velocity allows to rewrite (14) as , [START_REF] Wedel | A Box-Integration/WENO solver for the Boltzmann Transport Equation its Application to HighSpeed Heterojunction Bipolar Transistors[END_REF] where v sat is the  valley saturation velocity and   = 5v sat  w /3 is the energy relaxation length as defined in [START_REF] Slotboom | Non-local impact ionization in silicon devices[END_REF]. This equation allows to derive ( 15) and ( 16) [START_REF] Slotboom | Non-local impact ionization in silicon devices[END_REF]. The parameters of the two-valley DD system are listed in the tables below. x/ m E/ (eV)

E Γ E L E V E FnΓ E FnL E Fp (a) (b)
(c)

Fig. 1 .

 1 Fig.1. Measured (symbols) and simulated (solid lines) bulk velocity field characteristics of intrinsic (a) GaAs[START_REF] Blakemore | Semiconducting and other major properties of gallium arsenide[END_REF], (b) InP[START_REF] Fawcett | High-field transport in indium phosphide[END_REF], (c) InAs and (d) GaSb at 300 K from the BTE. Note that no experimental data for (c) and (d) have been found for 300 K.

Fig. 2 .

 2 Fig.2. Intervalley transition frequency in GaAs forN I = [1.0x10 15 ,1.4x1016 , 2.2x1016 , 3.3x1016 , 5.2x1016 , 8.0x1016 , 1.1x1018 , 1.7x10 18 , 2.7x10 18 , 4.2x10 18 , 6.4x10 18 ]cm -3 . Comparison between BTE (symbols) and analytical approximations (lines) for (a) f L with (5) and (b) f L with[START_REF]Minimos-NT[END_REF].

Fig. 3 .

 3 Fig.3. Average velocity in (a) the  and (b) the L valley of GaAs from BTE (symbols) and from (10) (lines) for the same impurity concentrations N I as listed in Fig. 2.

=Fig. 4 .

 4 Fig.4. Bulk velocity-field characteristic of (a) GaAs, (b) InP, (c) InAs and (d) GaSb from BTE simulation (symbols) and the 2v-DD equation system (lines) for the same impurity concentrations N I as listed in Fig. 2.

Fig. 5 .Fig. 6 .

 56 Fig.5. Comparison of different fields (left axis) based on local quantities and the mean carrier energy (right axis, black solid line)for the GaAs n+-n-n+ structure of Fig.6. The applied voltage is 3 V.

Fig. 7 .

 7 Fig.7. (a) Electric field and (b) mean carrier energy from BTE simulations of the n + nn + GaAs structure for different applied voltages 0.25-3 V.

Fig. 8 .

 8 Fig.8. Comparison of (a) current density, (b) electrostatic potential and (c) carrier distribution between the valleys of an InP structure with the same doping profile as depicted in Fig. 6(c) from BTE (dots) and DD simulations with the energy estimator approach (iii) (solid lines) and (i) (dashed lines). The bias in (b) and (c) is V = 3 V.

Fig. 9 .

 9 Fig.9. Simulated GaAs 1-z Sb z /InP HBT with z = 0.5 in the base and z = 0 in all other regions. (a) Band diagram for V BE = 0 V. (b) Electron density in the (solid lines) and L valleys (symbols) and impurity concentration (dashed line) in the base and collector for operating points J C = [4x10 -3 , 0.18, 4.6, 7.4] mA/ m 2 . (c) Transit frequency vs. collector current density with intervalley valley transfer turned on and off. V BC = -0.2 V in all cases. The dashed line in (a) indicates the net doping concentration in log scale. Appendix B: Extracted Model Parameters

Table 1 :

 1 Parameters for the velocity field model in the  valley.

	parameter	GaAs	InP	InAs	GaSb
	 ,L / (cm 2 /Vs)	7060	3410	23300 20300

Table 2 :

 2 Parameters for the velocity field model in the L valley.

		10 19
	L /(cm -3	10 16
	Γ , n
	n	10 13
		0.15	0.2	0 .25	0.3
				x/μm
		400	ivt on ivt off
	f T /GHz	200
		10 -1 0		10 0	10 1
				J C /mA μm -2

)

Table 3 :

 3 Parameters for the intervalley transfer rate model.

	parameter	GaAs	InP	InAs	GaSb
	 ,min / (cm 2 /Vs) 457	36	324	81
	N ,ref /(10 17 cm -3 ) 0.52	1.42	0.383	0.735
	 	0.60	0.48	0.50	0.56
	 	4.17	6.60	4.17	1.38
	v ,sat /(10 7 cm/s)	2.30	2.35	2.60	3.76
	parameter	GaAs	InP	InAs	GaSb
	 L,L /(cm 2 /Vs)	377	185	2240	2230
	 L,min /(cm 2 /Vs)	62	48	254	103
	N L,ref /(10 17 cm -3 ) 4.09	8.34	4.87	1.43
	 L	0.62	0.65	0.51	0.62
	 L	1.73	1.98	2.27	1.98
	v L,sat /(10 7 cm/s)	0.92	0.80	1.18	0.84
	parameter	GaAs	InP	InAs	GaSb
	f L0 / THz	1.70	7.99	1.59	0.11
	δ f /(GHz/(kVcm)) 5.79	8.96	25.0	5.33
	a f / THz	0.36	0.19	0.32	0.64
	E ivt,0 /(kV/cm)	3.00	11.4	1.63	0.019
	N ivt /(10 16 cm -3 )	78.6	1.22x10 3 178	4.31x10 -3
	 ivt	0.46	0.58	0.54	0.48

This work has been partially supported by the EU project H2020-ECSEL project TARANTO and the German National Science Foundation (DFG SCHR695/14, /21).