Ludovic Castro

Gabriel Therukauff

Aurelien Vantomme

Alexandre Welle

Luc Haspeslagh

Jean-Michel Brusson

Laurent Maron

Jean-François Carpentier

Evgueni Kirillov

A Theoretical Outlook on the Stereoselectivity Origins of Isoselective Zirconocene Propylene Polymerization Catalysts

For the Table of Contents entry

Origins of stereoselectivity in propylene polymerization for a series of isoselective metallocenes belonging to two different families, namely C1-symmetric {Cp/Flu}-type and C2-symmetric {SBI}-type systems, were apprehended by theoretical calculations using an improved DFT model, and meso/rac pentad distributions were predicted with good precision.

INTRODUCTION

Highly isotactic polypropylene (iPP) can be efficiently produced using two main families of zirconocene precatalysts: (i) C1-symmetric one-carbon-bridged cyclopentadienylfluorenyl ({R 1 R 2 C-(Flu)(Cp)}, hereafter referred to as {Cp/Flu}, 1-R and 2, Scheme 1) and (ii) C2-symmetric silicon-bridged ansa-bis(indenyl) ({R 1 2Si-(2-Me-4-R 2 -Ind)2} or {SBI}, 3). 1 Understanding the actual mechanism of stereocontrol operating with those systems, which allows engineering of more stereoselective catalysts, has been the subject of numerous studies. 2 For example, some general aspects of regio-and stereochemistry of propylene insertions for rac-{bridge-(1-Ind)2}ZrX2 (where bridge = R2Si, R2C or C2H4; X = Cl, alkyl), 3 in particular, rac-{Me2Si(2-Me-4-Ph-Ind)2}ZrCl2 (3), 4 have been already studied computationally using various levels of theory. Thus, enantiomorphic site control (ESC) via ubiquitous chain migratory insertion has been corroborated as the operative mechanism for this system, which resulted in the design of new, more stereoselective {SBI}-type metallocenes. 4c Scheme 1. Zirconocene precatalysts 1-R (Me, Ph), 2 and 3 used as models in this study. 10.1002/chem.201801438

Accepted Manuscript

Chemistry -A European Journal On the other hand, for the {Cp/Flu}-type metallocenes, the situation is less obvious.

While the chain migratory insertion (CMI) mechanism has been generally accepted for the syndioselective Cs-symmetric systems, 5,3i two mechanisms were hypothesized for their isoselective C1-symmetric analogues (Scheme 2). First, chain "stationary" insertion (or site epimerization), involving monomer insertion only on the more stereoselective site, followed by site epimerization (backskip), has been suggested by Razavi et al. 6 as the main operating mechanism of formation of isotactic sequences. In addition, an alternating mechanism has been proposed by Marks et al., 7 which involves both a lesser and a more stereoselective sites of the molecule for monomer insertion. Bercaw et al. have discussed a case where both above mechanisms could be functioning simultaneously. 8 In the latter study, metallocene precatalyst 2 that features superior isoselectivity was developed; it was proposed that this catalyst enchains propylene via solely the site epimerization (backskip) mechanism, completely suppressing the alternating one. 8b Possible stereocontrol mechanisms for C1-symmetric systems analogous to 1-R and 2 (Scheme 1) have been computationally probed by Ziegler, Razavi et al. 9 and others. 10 Yet, the obtained thermochemical data were found not sufficiently explicit to allow drawing straightforward conclusion on the operative stereocontrol mechanism.

Scheme 2. Proposed mechanisms of stereocontrol with C1-symmetric metallocene catalysts. 8 10.1002/chem.201801438

Accepted Manuscript

Chemistry -A European Journal

In this study, we aimed at rationalizing the mechanism of regio-and stereocontrol of the highly isoselective polymerization of propylene with complexes 1-R 11 and 2 8b (Scheme 1) by theoretical computations of the first, second and third insertion steps. For straight comparison and benchmarking, DFT studies on the same level of theory were employed for system 3, which is well known to operate via an enantiomorphic site control mechanism involving chain migratory insertion reactions (vide supra). 2 Besides, for the above catalytic systems, a probabilistic modeling based on the procedure described by Farina et al. 12 and Chien et al. 13 was carried out using the computed thermochemical data; 14 meso/rac pentad distributions for iPP samples obtained at various temperatures were predicted and found to be in good to acceptable consistency with experimental data. this is due to the unclear and controversial structure of MAO, its dynamic behavior 18 and obvious computational costs related to the large structure of this system. The nature of the counterion may play a significant role in stereochemistry of propylene insertion, namely, influencing the propagation, site epimerization (backskip) and misinsertion kinetics. 5b,19 At the same time, the direct influence of the counterion is neglected in many theoretical models applied to stereochemistry computations. 14,20 This can be rationalized in the view of the following two possible phenomena: (1) in case of tightly bound ion pairs incorporating strongly coordinating counterions (e.g., [MeB(C6F5)3] -), extraction and substitution of the counterion upon coordination of a monomer molecule is an endothermic process, which activation barrier magnitude will be barely independent of the site involved in coordination; 9c

RESULTS AND DISCUSSION

General

(2) in late polymerization steps, the growing polymer chain may expulse the counterion on some distance from the coordination site, thus making its further involvement minimal.

In a recent contribution, we have developed a DFT theoretical model competent for quantitatively describing propylene insertion into metal-alkyl bonds of neutral and cationic metallocenes including the "naked" cation [{Ph(H)C(3,6-tBu2-Flu)(3-tBu-5-Me-Cp)}ZrMe] + derived from 1-Me. 21 In this model, ΔH/ΔH # values are calculated, which were shown for a binary system involving bis(cyclopentadienyl)yttrium alkyl and propylene for which experimental data are available, to be very close to the corresponding ΔG/ΔG # values calculated using both dispersion correction and SMD solvation model. This DFT model (see SI) has been used throughout this study for all computations.

In order to assess the regio-and stereoselectivities of catalysts, all possible mechanisms of monomer insertion were considered. First, two main types of insertions are possible: the primary (1,2-; pr) and secondary (2,1-; sec) insertions, which lead, for the first polymerization step, to iso-butyl and sec-butyl products (R = Me), respectively (Scheme 3). Scheme 3. Primary and secondary insertions of propylene into Zr-R bond.

Accepted Manuscript

Chemistry -A European Journal Second, eight different insertion reactions have to be considered for C1-symmetric systems 1-R and 2, as depicted in Scheme 4. Indeed, these systems can adopt two different configurations, as the Zr-polymeryl chain can reside either in anti or syn orientation with respect to the tBu group of the Cp ligand (these two possibilities are hereafter denoted as A and S, respectively). Moreover, the propylene molecule can insert into the Zr-C(polymeryl) bond via either the re or si face. Hence, for example, primary insertion si into the anti isomer of 1-R and 2 will be noted as A-pr-si. It is noteworthy that for insertion into the anti isomer, propylene binds at the more crowded site of the catalyst (i.e., the one accommodating the tBu-Cp substituent) and that the resulting product will use in turn the syn coordination site for the forthcoming insertion step. Correspondingly, for the insertion into the syn isomer, propylene coordinates at the less crowded site and the resulting product will use the anti coordination site. For the C2-symmetric system 3, only four different insertions, namely pr-si, pr-re, sec-si and sec-re, have to be considered. were calculated (Table 1), and eight possible pathways via the formation of the corresponding propylene π-adducts were considered at each insertion step. For the 5-Ph-Cp-substituted analogue [{Ph(H)C(3,6-tBu2-Flu)(3-tBu-5-Ph-Cp)}ZrMe] + (1-Ph For each insertion step, the zero of energy corresponds to the anti isomer of the "naked" cation (more stable than the syn isomer by 2.5 kcal.mol -1 for the very first insertion step) and a propylene molecule on an infinite distance (Fig. 1).

On the basis of the results obtained for systems 1-R (R = Me, Ph) and 2, the following trends can be highlighted:

1)

Regardless the actual monomer insertion step, primary (1,2-) coordination and insertion of propylene into Zr-C(alkyl) bond are systematically highly selective for both systems, which is evidenced from the corresponding lower activation barriers (typically by 3-11 kcal•mol -1 ) than those calculated for secondary (2,1-) insertions. This appears to be driven in part by steric considerations. For example, the corresponding optimized structures of the transition states for the primary insertions in 1-Me (Fig. S2) revealed no particular steric repulsion between the methyl group of the coordinated propylene molecule and the ligand system, while those for the secondary insertions feature strong steric hindrance imposed by the ligand system. Also, the Zr-iBu and higher Zrpolymeryl products resulting from the 1,2-insertions are more thermodynamically stable; this apparently stems from the presence of b-C-H … Zr agostic interactions (Fig. 1) stabilizing the cationic center. 22 10.1002/chem.201801438

Accepted Manuscript

Chemistry -A European Journal

2)

As demonstrated for 1-Me, the very first propylene insertion into Zr-Me bond is not stereoselective. Indeed, all four primary insertion reactions feature very close barriers (7.1-8.5 kcal•mol -1 ) and afford thermodynamically stable (by 14.6 and 17.7 kcal•mol -1 ) products (Fig. 1).

3) Stereoselectivity control over monomer insertion for the C1-symmetric system 1-Me comes out from the second insertion step, that is propylene insertion in the Zr-C(iBu) bond. 23 First of all, in agreement with the mechanism of Razavi et al.,6,9d,8 A-pr-si insertion of propylene on the crowded site ("head down" into the free space of the central region of the fluorenyl ligand) is the preferred one in terms of kinetics. Also, the energy differences between the propylene insertion barriers for A-pr-re and A-pr-si computed for the second and third insertion steps (DDH ¹ = 4.6 and 3.9 kcal•mol -1 , respectively) suggest that an enantiofacial misinsertion is kinetically disfavored on the crowded site. As expected, the open site of the metallocene is significantly less selective: small energy differences (< 1 kcal•mol -1 ) between the barriers of the S-pr-re

Accepted Manuscript

Chemistry -A European Journal and S-pr-si insertions of propylene were calculated for the second and third insertion reactions.

Second, the iso-butyl and the higher products formed after the S-pr-si and S-pr-re second and third insertions, respectively, are more stable by 3-5 kcal•mol -1 than their counterparts formed after the A-pr-si and A-pr-re insertions. This difference is due to the relative positioning of the iBu (or polymeryl) groups with respect to the tBu substituent of the Cp moiety -anti for the S-pr-re and S-pr-si insertions products and syn for the A-pr-re and A-pr-si ones (Fig. 2).

Overall, these results suggest that, after the stereoselective A-pr-si insertion, the resulted constrained syn product with the propagating chain initially residing on the crowded site should rearrange via epimerization (Scheme 5) to a less sterically congested, more stable (by ca. 4 kcal•mol -1 ) anti isomer. 24 This sequence of insertionepimerization steps leads to the formation of iPP by a stationary propagation, i.e. a chain "stationary" insertion mechanism. 6 Scheme 5. Epimerization (backskip) process for C1-symmetric {Cp/Flu} metallocene catalysts.

10.1002/chem.201801438

Accepted Manuscript

Chemistry -A European Journal

4)

The data computed for the second insertion step with system 1-Ph (Table 1) are also in complete agreement with the epimerization mechanism. Yet, the 1-Ph system is computed to be less stereoselective than 1-Me, due to the smaller energy difference between the barriers of misinsertion A-pr-re and stereoselective A-pr-si insertion (DDH ¹ = 3.1 kcal•mol -1 ). Also, the higher propylene A-pr-si insertion barrier calculated for 1-Ph (DDH ¹ = 1.5 kcal•mol -1 ) calls for a lower polymerization activity than that observed More, epimerization of the A-pr-si insertion product to that of the S-pr-re insertion (Scheme 5) via presumably a very low-energy transition state 24 is much more favored by thermodynamics (by 12.7 kcal•mol -1 ) than that calculated for the 1-Me system (vide supra). The latter result and also the fact that the S-pr-re insertion features a much less stable propylene adduct (+13.3 kcal•mol -1 ) and a higher by energy transition state (17.6 kcal•mol -1 ) both suggest that the only stereocontrol mechanism operating for this system is via chain "stationary" insertion.

Using the three-parameter probabilistic model for C1-symmetric systems described by Farina et al. (see the SI), 12 we have predicted meso/rac distributions from the computed insertion enthalpy values (Table 2). This model appeared to be competent to estimate (within acceptable margins) pentad distributions at different temperatures. For example, for 1-Me, the model reproduces very well the results obtained at lower polymerization temperature (60-70 °C), while at higher temperature (100 °C) the deviation (quadratic mean, RMS) was found higher. Apparently, the latter discrepancy can be a result of a higher amount of regioerrors (issued from secondary insertions) formed at higher temperature and whose account cannot be included in the given probabilistic model. 25 A similar trend was obtained for system 1-Ph (Table 3), yet with more significant discrepancy between the experimental and simulated results at higher polymerization temperatures. Of note, a much lower regioregularity (primary

Accepted Manuscript

Chemistry -A European Journal insertions = 91.8%) was determined experimentally in PPs produced with this system at higher polymerization temperature (Tpol = 100 °C).

The DFT calculated probabilistic descriptors a, b and k are all consistent with general isoselective behavior of the studied metallocene systems. The near-one calculated values of descriptors b and k are consistent with the site epimerization mechanism of Razavi et al. 6 when, after isoselective insertion on the crowded site, the propagating chain immediately

migrates to the open site of the catalyst. On the other hand, the relatively high value of a (along with inequality of the (1-k) descriptor to 0) is diagnostic that the stereocontrol alternating mechanism, 8 in which propylene insertion occurs on the lesser stereoselective open site, may also be valid for that given system. Both results are in line with the general assumption by Bercaw et al. 8 that both alternating and site epimerization mechanisms of stereocontrol can operate simultaneously for one and the same catalyst. More, the predominance of the more selective chain "stationary" mechanism over the less selective alternating one for this particular metallocene system is essentially determined by the steric hindrance of the 3-cyclopentadienyl substituent and, to a lesser extent, of the substituents on the fluorenyl ligand. Thus, the bulkier the 3-cyclopentadienyl substituent (e.g., 2-methyl-2adamantyl >> tBu >> Me), the higher the stereocontrol via chain migratory mechanism and the lower the role of the alternating mechanism. The latter aspect is clearly demonstrated on system 2, which features the highest degree of stereocontrol in terms of [m 4 ] pentads content (>99%) among C1-symmetric {Cp/Flu} catalysts. 8b In this case, the "stationary" stereocontrol mechanism is consistent with instant site "epimerization" (k = 1) realized after a stereoselective (b ® 1) monomer insertion on the more stereoselective (crowded) site, while the stereoselectivity of the lesser stereoselective (open) site (not used for the coordination) has utterly no influence (0 < a < 1). a The values for the insertion barriers were calculated as DH ≠ ins = HTS -HReactant. b The corresponding TS geometries were not located (n.l.) or not calculated (n.c.).

Accepted Manuscript

Chemistry -A European Journal C2-Symmetric System (3). Similar energetic data for the four possible structures for propylene insertion at the first, second and third steps of polymerization were determined for the cationic system [rac-{Me2Si(2-Me-4-Ph-Ind)2}ZrMe] + (3); representative data are summarized in Table 5. By analogy with the {Cp/Flu} counterparts, system 3 revealed a strong tendency to afford regioselective primary insertions and also exhibited no stereopreference at the very first insertion step. For instance, the activation barrier for the sec-si and sec-re insertions are much higher (by 5-10 kcal.mol -1 ); this is because the methyl group of propylene points toward indenyl ligands in both structures, and in the sec-si structure it induces the rotation of the phenyl ring in order to avoid steric repulsion (for the second insertion step, see Fig. 3). Also, formation of the corresponding primary insertion products proved to be favored thermodynamically, apparently stabilized by agostic interactions.

Starting from the second insertion step, system 3 becomes highly stereoselective. The latter fact stems from a significant difference (DDH ¹ = 4.3-4.6 kcal•mol -1 ) between the barriers of the stereoselective (pr-si) insertion and misinsertion (pr-re). Indeed, in the structure of the transition state for the pr-re insertion (Fig. 3), the methyl group of propylene points towards the 4-indenyl phenyl substituent and lies in cis position with respect to the iPr substituent of the iso-butyl group (representing the growing polymeryl chain), which in turn generates major steric congestion. In the corresponding structure for the stereoregular pr-si insertion, these repulsions are minimized. These results are in agreement with an operative enantiomorphic site control mechanism.

Of note, system 3 features significantly lower propylene insertion barriers into Zrpolymeryl bond (DDH ¹ = 3-6 kcal•mol -1 ) than those calculated for the {Cp/Flu}-type systems 1-R and 2. This fact is diagnostic of a higher propensity of the former system to polymerize propylene, which is actually in line with the observed superior productivities of {SBI}-based catalysts as compared to those of {Cp/Flu}-based systems. 11 10.1002/chem.201801438

Accepted Manuscript

Chemistry -A European Journal We also probed modeling of the microstructure of iPP produced with system 3 using a one-parameter model (Table 6). The model describes accurately enough (in terms of RMS) the experimental pentad distributions for samples obtained at high polymerization temperature (80-95 °C). The near-one values of the probabilistic descriptor a, DFT computed for the second and third insertion steps at different temperatures, are all in agreement with very high isoselectivity of the catalytic system. However, the theoretical model predicts much higher stereoselectivity than that achieved experimentally. This discrepancy again may stem from insertion side-processes, like regioirregular secondary insertion, which are not included in the theoretical statistical model. a The values for the insertion barriers were calculated as DH ≠ ins = HTS -HReactants. The values in brackets for the insertion barriers were calculated as DH ≠ = HTS -HAdduct.

Table 6. Pentad distributions (%) and corresponding probability parameters determined experimentally for system 3, and those simulated using a one-parameter model and DFT computed. a

The validity of chain "stationary" insertion mechanism (i.e., site epimerization by backskip) for isoselective C1-symmetry {Cp/Flu}-type metallocenes was demonstrated using a recently developed theoretical model 21 for the investigation of stereocontrol mechanism of propylene polymerization. However, an alternative, yet less stereoselective alternating mechanism of stereocontrol can be operational at the same time. In line with previously reported experimental results, 8b steric bulkiness imposed by the 3-cyclopentadienyl substituent (i.e., 2-methyl-2-adamantyl vs tBu) is shown to be beneficial for the suppression of this alternative mechanism. Eventually, this makes the more stereoselective site epimerization mechanism to be unique for the heavily 3-cyclopentadienyl-substituted {Cp/Flu}-type metallocenes. On the other hand, introduction of a bulkier group (Ph vs Me) into the 5-cyclopentadienyl position facilitates the alternating mechanism, eventually affording a less stereoselective catalyst. For the C2-symmetric {SBI}-type system 3, 4 the well-accepted chain migratory stereocontrol mechanism was corroborated using the new theoretical model. 21 Implementation of the computed thermochemical data into the probabilistic analysis allowed reproducing and predicting statistical pentad distributions for iPPs obtained with both families of metallocene catalysts operating via dissimilar stereocontrol mechanisms. Further studies on the application of this theoretical model for modeling/understanding of other metallocene systems, as well as computer-aided design of new more performing/stereoselective polymerization catalytic systems, are underway in our laboratories. Farina, M.; Di Silvestro, G.; Terragni, A. Macromol. Chem. Phys. 1995, 196, 353-367.

(e) Di Silvestro, G.; Sozzani, P.; Terragni, A. Macromol. Chem. Phys. 1996, 197, 3209-3228. 13 (a) Yu, Z.; Chien, J. C. W. J. Polym. Sci. A: Polym. Chem. 1995, 33, 125-135. (b) It should be noted that these agostic interactions may be absent when a counterion (not considered by this theoretical model) is present and coordinated to the cationic metal center.

23 Noteworthy, starting from the second insertion step, monomer coordination onto the metal center of {Cp/Flu} systems became endothermic.

24

The corresponding state of the epimerization process can hardly be located because of a very low and smooth energetic profile.

25

As proposed in the previous theoretical study, 9a other possible factors, like entropic factors, back-side insertion and chain epimerization may also be responsible for the lower precision of statistical modeling of pentad distributions in polymers obtained at higher temperatures.

10.1002/chem.201801438

Accepted Manuscript

Chemistry -A European Journal

  Considerations used for Theoretical Models. It is generally accepted that ionic complexes of general structure [Cp2ZrMe] + [Me-MAO] -(Cp = cyclopentadienyl-type ligand), which form during the activation step in a reaction between parent Cp2ZrX2 (X = halogen or alkyl) and excess MAO, constitute a loosely associated form of a "true" catalysta "naked" cation [Cp2ZrMe] + and the corresponding counterion [Me-MAO] -. 15,16 Although different theoretical models treating complex systems incorporating molecular versions of anions such as [MeB(C6F5)3] -, [B(C6F5)4] -, and [Me{MeAl(µ3-O)}3] -were generally successful, 17 implementation of the MAO-based anion in the calculations is still problematic;

Scheme 4 .

 4 Scheme 4. Eight possible propylene π-adducts of cationic complexes for C1-symmetric

  Systems. Energetic data for the three consecutive propylene insertion steps for the cationic precursor [{Ph(H)C(3,6-tBu2-Flu)(3-tBu-5-Me-Cp)}ZrMe] + (1-Me)

  ) only the second insertion step was computed, and for the significantly bulkier [{Me2C(Oct)(3-(2-methyl-2adamanthyl)-Cp)}ZrMe] + (2), only the second and the third insertion steps were calculated.

Figure 1 .

 1 Figure 1. Optimized geometries of anti (left) and syn (right) isomers of the iso-butyl product

Figure 2 .

 2 Figure 2. Optimized transition state geometries for the second insertion of propylene.

a

  Experimental data taken from ref. 11; a and b -the probabilities of stereoselective insertion on the open and crowded metallocene faces, respectively; k -the probability of site epimerization. b RMS = ((Σ(Iobs -Icalc)²/9) 0.5 . c Calculated as:

a

  Experimental data taken from ref. 11; a and b -the probabilities of stereoselective insertion on the open and crowded metallocene faces, respectively; k -the probability of site "epimerization". b RMS = ((Σ(Iobs -Icalc)²/9) 0.5 . c Calculated as:

a

  Experimental data taken from ref. 8b; a and b -the probabilities of stereoselective insertion on the open and crowded metallocene faces, respectively; k -the probability of site "epimerization". b RMS = ((Σ(Iobs -Icalc)²/9) 0.5 . c has no influence as k = 1. d Calculated as:

Figure 3 .

 3 Figure 3. Optimized transition state structures for the four possible propylene insertions in

Table 1 .

 1 Energetic data calculated for different insertion steps of propylene for systems 1-R (R = Me, Ph) and 2. a

	10.1002/chem.201801438

Table 2 .

 2 Pentad distributions (%) and corresponding probability parameters 12 determined experimentally for 1-Me, and those simulated using a three-parameter model and DFT computation data. a

	Calculated from DFT data

Table 3 .

 3 Pentad distributions (%) and corresponding probability parameters 12 determined experimentally for 1-Ph, and those simulated using a three-parameter model and DFT computation data. a

	Calculated from DFT data

Table 4 .

 4 Pentad distributions (%) and corresponding probability parameters 12 determined experimentally 8b for 2, and those simulated using a three-parameter model and DFT computation data. a

					Calculated from DFT data
		Exp		Calc	2 nd insertion	3 d insertion
	T, [° C]	0	20		0	20	0	20
	Mn [×10 3 ]	370	425	-
	1,2 ins	n.d.		-	-	-	-	-
	[mmmm]	>99	>99	99.0	97.7	96.6	99.5	99.2
	[mmmr]	n.d.	n.d.	0.4	0.9	1.3	0.2	0.3
	[rmmr]	n.d.	n.d.	0.0	0.0	0.0	0.0	0.0
	[mmrr]	n.d.	n.d.	0.4	0.9	1.3	0.2	0.3
	[mrmm] +[rmrr]	n.d.	n.d.

Table 5 .

 5 Energetic data calculated for the first three propylene insertion steps in system 3. a

	Reaction		DH (1 st insertion)	DH (2 nd insertion)	DH (3 d insertion)
		Adduct	-5.3	-3.9	-3.8
	pr-si	TS	4.7 (10.0)	5.8 (9.8)	5.4 (9.2)
		Product	-15.3	-18.4	-15.5
	pr-re	Adduct	-5.0	-1.7	-1.7

  Yu, Z.; Chien, J. C. W. J. Polym. Sci. A: Polym. Chem. 1995, 33, 1085-1094. 14 Angermund, K.; Fink, G.; Hensen, V. R.; Kleinschmidt, R. Chem. Rev. 2000, 100,

	1457-1470.

15

Bochmann, M. Organometallics 2010, 29, 4711-4740. 16

For a recent report on the reactivity of dichloro metallocenes 1-Me and 3 with MAO and molecular activators, see: Theurkauff, G.; Bader, M.; Marquet, N.; Bondon, A.; Roisnel, T.; Guegan, J.-P.; Amar, A.; Boucekkine, A.; Carpentier, J.-F.;

Kirillov, E. Organometallics 2016, 35, 258-276. 22

10.1002/chem.201801438 Accepted ManuscriptChemistry -A European Journal

ACKNOWLEDGMENTS

This work was supported by TOTAL SA and TOTAL Research & Technology Feluy. L. Maron and L. Castro are grateful to the CALMIP for computation time.

Accepted Manuscript

Chemistry -A European Journal


Fabien Destic
email: fabien.destic@isae.fr

Christophe Bouvet

Impact damages detection on composite materials by THz imaging

This paper presents a Non-Destructive Testing (NDT) method based on the penetration properties of terahertz (THz) waves. A CW raster-scanning THz imaging setup, using a 3.8 THz Quantum Cascade Laser as a source, is used to perform NDT of polypropylene/polypropylene composite samples. Results from transmission and reflection THz imaging are compared to ultrasound C-scan. THz images in reflection give similar results to C-scan whereas THz transmission images provide more information about delaminations and cracks in the fiber fabrics.

Introduction

The field of THz waves lies in the electromagnetic spectrum in the frequency range between 300 GHz and 10 THz, or in wavelengths from 30 μm to 1 mm. On the one hand, the energy of THz photons is very low, below 40 meV. On the other hand, most non-metallic and non-polar media are "transparent" to THz waves. These two properties make THz waves a potential candidate for Non-Destructive Testing in a variety of fields, such as art conservation science [START_REF] Gallerano | THz-ARTE: non-invasive terahertz diagnostics for art conservation[END_REF], plastic weld joints [START_REF] Wietzke | Terahertz imaging: a new non-destructive technique for the quality control of plastic weld joints[END_REF], quality control on chocolate bars [START_REF] Jördens | Towards real-time terahertz quality assurance of chocolate products[END_REF] or the absorption of water into polyamide and wood plastic composite [START_REF] Jördens | Investigation of the water absorption in polyamide and wood plastic composite by terahertz time-domain spectroscopy[END_REF].

A composite material is a combination of two (or more) materials working together that give the composite its unique properties. This denomination is applied to many materials and many sectors: glass, aramide or carbon fibre reinforced composites used in the aeronautics and automotive industries or sports equipment, and wood composite in furniture manufacturing, among others.

NDT using THz waves can be an alternative when "traditional" methods are not very efficient (composite made from soft epoxy resin, for instance) or to improve spatial resolution. In a simplified manner, there are two kinds of THz radiation technology: pulsed and Continuous Wave (CW).

In pulsed systems, a femtosecond Ti:Sa or mode-locked fiber laser (duration ∼100 fs) produces an optical-pulse train sent on a THz emitter (photoconductive antenna or nonlinear crystal) to form electromagnetic pulses. After being transmitted through or reflected by the sample, the resulting waveform is coherently detected in time-domain using a delay line and Fourier Transformed to be analyzed in the frequency domain, providing a broad spectrum that covers typically the frequency range [0.1-3] THz [START_REF] Zhang | Introduction to THz wave photonics[END_REF]. The family of CW or quasi-CW systems can be divided into two types depending on the detection used: coherent detection, generally used with low power tunable sources, also uses an optical delay line and has the same characteristics as pulsed systems, and direct detection where phase information is not available but operations are faster and easier.

Over the past ten years, some research has been conducted on characterization or imaging, in the THz range, of various composite materials, most of which use pulsed THz-Time Domain Spectroscopy systems (THz-TDS). For instance, THz-NDT have been applied to titanium nanospheres in low density polyethylene (LDPE) [START_REF] Rutz | Quality control of polymeric compounds using terahertz imaging[END_REF] and glass-fibre reinforced polymers [START_REF] Rutz | Non-destructive testing of glass-fibre reinforced polymers using terahertz spectroscopy[END_REF][START_REF] Jördens | Terahertz spectroscopy to study the orientation of glass fibres in reinforced plastics[END_REF]. Fire damage testing on carbon or glass fibre was also explored [START_REF] Redo-Sanchez | Damage and defect inspection with terahertz waves[END_REF][START_REF] Karpowicz | Fire damage on carbon fiber materials characterized by THz waves[END_REF][START_REF] Stoik | Nondestructive evaluation of aircraft composites using transmissive terahertz time domain spectroscopy[END_REF]. The European FP7 project, DOTNAC, is dedicated to the development of a TDS setup to inspect aeronautics composite materials. Defects on glass fiber reinforced plastics (foreign material inserts, delaminations, or moisture contamination) can be visualized or coatings analyzed on carbon fiber samples [START_REF] Ospald | Aeronautics composite material inspection with a terahertz time-domain spectroscopy system[END_REF]. Amenabar et al. propose a review to THz Non-Destructive Testing of composite matter [START_REF] Amenabar | In introductory review to THz non-destructive testing of composite mater[END_REF].

Continuous wave THz inspection is also possible as demonstrated on sprayed-on foam insulation (SOFI) with a gas laser at 1.63 THz [START_REF] Redo-Sanchez | Damage and defect inspection with terahertz waves[END_REF] or carbon fibre at 0.6 THz with Gunn diodes [START_REF] Karpowicz | Fire damage on carbon fiber materials characterized by THz waves[END_REF]. Composite materials diagnostic using CW THz waves produced by Quantum Cascade Laser on thin samples of glass fibre or Kevlar composite have been studied [START_REF] Destic | THz QCL-based active imaging applied to composite materials diagnostic[END_REF]. This approach seems to be very promising as QCLs are THz sources emitting a few tens of mW, at frequencies in the [START_REF] Gallerano | THz-ARTE: non-invasive terahertz diagnostics for art conservation[END_REF][START_REF] Wietzke | Terahertz imaging: a new non-destructive technique for the quality control of plastic weld joints[END_REF][START_REF] Jördens | Towards real-time terahertz quality assurance of chocolate products[END_REF][START_REF] Jördens | Investigation of the water absorption in polyamide and wood plastic composite by terahertz time-domain spectroscopy[END_REF] THz range, with a constantly improving beam quality.

Experimental setup

The raster-scanning imaging experimental setup (Fig. 1) uses a 3.8 THz (λ = 79 μm) Quantum Cascade Laser (QCL) as a THz source. The laser, based on the Resonant-Phonon depopulation scheme, has been designed by R. Colombelli's group at Université Paris-Sud [START_REF] Chassagneux | Limiting factors to the temperature performance of THz quantum cascade lasers based on the resonant-phonon depopulation scheme[END_REF]. When cooled at ∼ 10 K by a closed cycle cryocooler, it emits 35 mW peak power in pulsed operation [START_REF] Destic | THz absolute power measurement: a simple and reliable method[END_REF]. An ILX Lightwave LDX-36010-35 pulsed current source drives the QCL at 1.4 A by 950 μs long pulses at a 21 Hz repetition rate. L 1 is a TPX lens that works both as a lens and a cryostat window. The collimated laser beam diameter is ∼13 mm in horizontal and ∼15 mm in vertical, divergence being ∼5 mrad. A second TPX lens, L 2 , is used to focalize the laser beam in a 800 μm diameter spot where the sample under test is placed. The average power is ∼130 μW.

The sample to be imaged is moved in X and Y by two Newport SMC100 motorized translation stages. The setup can operate either in transmission mode with the off-axis parabolic mirrors inside the red dotted rectangle or in reflection mode by means of a HRFZ-Si beam splitter and mirrors in the black dotted rectangle. At each point or pixel, the mean and standard deviation of 10 intensity values from the detector (Tydex Golay cell, NEP = 113 pW.Hz -1/2 ) is acquired via a Signal Recovery 7270DSP lock-in amplifier and, so, a transmission or reflection image can be constructed. Using a metallic resolution target, a spatial resolution of 0.5 mm in vertical and horizontal has been measured for the two arrangements of the imaging setup. The acquisition of a 30 mm × 30 mm image with 0.5 mm steps lasts about 1 h 25 min.

In order to characterize the THz images as objectively as possible, we will use two figures of merit: where S i,j and N i,j = σ (S i,j ) are the signal amplitude and the noise at pixel (i, j).

Impact damage detection

Samples

The samples used to demonstrate damage detection were 3 mm thick polypropylene composite plates, commercialized by Goodfellow under reference PP403300. This kind of composite has applications in automotive components, industrial cladding, audio products, personal protective equipment and sports goods.

The composite deals with 3 mm thickness plate manufactured from 20 woven plies oriented at 0 • of thermoplastic self-reinforcing composite produced by compacting polypropylene fabric. This type of composite possesses a unique combination of impact/crash resistance and stiffness/strength. The density is 0.92 g/cm 3 , the tensile modulus 5.0 GPa and the tensile stress is 180 MPa. In self-reinforced composite, the polymer matrix is reinforced with high-tenacity fibers of the same polymer family: in the present case, it deals with PP matrix and PP fibers. During the hot compacting manufacturing process, part of the fibers are melted and recrystallized upon cooling to bond the structure together. The global morphology of the obtained composite is not as clear as a classical composite with well separated plies, and, consequently, the resulting impact damage does not consist in delamination and matrix cracking (Fig. 2(c)) but rather in a damaged zone including matrix cracking and fiber/matrix debonding.

Damage creation

Fig. 2(a) shows the falling weight impact testing system used. The falling weight (mass = 2.03 kg) is in free fall guided in a tube. The impact energy is known from the velocity, measured by sensors [START_REF] Bouvet | Dynamic deformation, damage and fracture in composite materials and structures[END_REF]. Two 100 mm × 150 mm plates (Fig. 2(b)) have been impacted, in their center, with energies at levels of 10 J (sample S1) and 20 J (sample S2). The impact damage (Fig. 2(c)) consists of:

• dent depth due to the impactor contact • matrix cracking • delamination (failure of interface between consecutive plies) • failure of the ply located back side Of course, the higher the impact energy level, the greater the damage.

Results

Reflection THz imaging

For both samples, a reflection image is acquired with the setup and conditions described in Fig. 2. For each sample, S1 and S2, two images are taken: one with the sample illuminated on the impact side, so-called front side or recto, one illuminating the sample opposite to the impact i.e. back side or verso. The THz images are then compared to ultrasonic maps from a Krautkramer USD30 equipment. To perform these maps, called C-scan images, the plate is fully immersed in water. A single transducer of 3 MHz frequency (∼500 μm in wavelength) operating as a transmitter-receiver is scanned in a plane parallel to the composite plate surface and only the reflection from the back face is monitored. This technique allows to measure both the amplitude of the reflection and the time of flight, i.e. the time which is needed for an emitted pulse reflected from the back face of the sample to the receiver, and then to evaluate the depth of the first damage. As our setup does not give phase information, the amplitude C-scan ultrasounds images are compared to the THz reflection images (after converting intensity to amplitude). On all images, amplitude is normalized relative to its maximum and displayed on a color scale.

For each sample, paired images are shown and compared: THz reflection image and ultrasonic amplitude C-scan, sorted according to illumination side and impact energy.

In order to easily compare the images obtained with THz and C-scan techniques, the damage was highlighted with a dotted white form and the impact location was highlighted with an arrow (Figs. 34567). Of course the damage form was the same for every image of the same impact, and symmetry was done for verso observations (Figs. 567). The sole purpose of the damage form was to compare the damage obtained with the different observations but it cannot be considered as the exact damaged zone. In fact, the complementarity of the different observations was used in order to plot the damaged zone as accurately as possible.

Recto

Regarding the damage created on sample S1 by the 10 J impact (Fig. 3), the overall size of the observed damage zones is similar. In the THz image (Fig. 3(a)), the central spot, corresponding to the impact point, shows the maximum amplitude signal, while it is the most damaged zone. It is probably due to the dent depth created by the impactor (Fig. 2(c)), which may act as a concave reflector and as a light focalizer.

Fig. 4 summarizes the results obtained on the plate impacted with an energy of 20 J. The THz image (Fig. 4(a)) shows no significant difference with Fig. 3(a), corresponding to the 10 J impact. Indeed, in both cases, the damage area is expressed by a maximum signal located within a circle of 3-4 mm diameter which is probably due to the dent depth of the impactor in the "soft" material and, so, the creation of a concave reflector that may act as a light focalizer. There is not any significant size difference on the blue crescent-shaped areas (signal minimum). In the ultrasonic image Fig. 4(b), the damage is contained within a zone of 20 mm high and 15 mm wide. Within this zone, three maximum spots are seen. They are assumed to correspond to non-damaged areas, which seems unlikely regarding their locations and sizes. Rather they are probably measurement artifacts Finally, for both the THz image and for the C-scan image, the observed damage zones are difficult to interpret.

Verso

Results for the back side or verso observations are reported in Fig. 5 and Fig. 6 for the 10 J and 20 J impacts, respectively. In all cases, the verso observations are clearer than recto observations, for both THz and for ultrasonic investigations. This is most likely due to due to the shape of the back side which is less deformed than the front side. Indeed the dent depth left by the impactor at the front side probably disturbs the wave penetration in the plate. At the back side, in spite of the ply failure (Fig. 2(c) and Fig. 8), the plate shape is smoother.

At the two impact energies, the signal levels (eq. colors), shapes and sizes produced either by the THz setup or the ultrasound C-scan are quite similar: an elliptical blue-green zone corresponding to a low signal for sample S1 in the 10 J impact case (Fig. 5(a) and Fig. 5(b)), a "potato-like" zone for sample S2 in the 20 J case. But, both for 10 J impact and for 20 J impact, the damaged zone evaluated using THz investigation seems clearer than with the ultrasonic one.

Transmission THz imaging

Then, a transmission image is acquired for both samples. The sample is illuminated by the QCL on the recto. The average attenuation due to sample absorption and reflectance is about 30 dB. Fig. 7(a) shows the resulting image for sample S1, impacted at 10 J. It presents a transmission minimum (red and brown) in an area of about 2 mm wide and 7 mm long (dimensions are determined from the intensity profile, not shown here). It corresponds to a non-impacted side crack, caused by the traction, which follows the weaving pattern. This crack, in the thickness of the plate, can be likened to a slot that diffracts terahertz radiation thus the signal is reduced. In the case of

Considerations about THz images

Uncertainties are difficult to evaluate. In fact, the C-scan is a very classic method which is used here as the reference. Unfortunately, in this case, the C-scan of woven PP composite gives uncertain results. This is most probably due to the complex morphology of the impact damage due to the hot compacting manufacturing process which melts the fibers of consecutive plies together.

Table 1 summarizes the DR and SNR obtained for the different THz images. On the one hand, dynamic range is higher for reflection images (∼19-20 dB) than for transmission, probably due to the shape of samples. For recto illumination, the dent depth of the impactor focalize the reflected light and the signal is high. Recto illumination were also used for transmission images but in these cases, the light is absorbed and diffracted during propagation through the material. In the damaged zone, diffraction effects are much higher and so, the signal intensity is lower.

On the other hand, SNR is better for transmission than for reflection experiments. This behavior is a result of the lock-in amplifier response to signal dynamic range: given that the amplifier input range is fixed regarding the high level, then the low level signal measurement is noisier and the SNR decreases as the DR increases.

Microscope observations of the damage

In order to validate our results, the two samples were cut at the impact point, perpendicularly to the crack. After polishing, the damage was observed with a Scanning Electron Microscope (SEM). Corresponding images are reported in Fig. 8.

In both cases, a crack is visible on the back side of the impact (verso) and, around this crack, a delamination zone. In picture 8(b) of damage to 20 J, it is clearly visible that the rear face of the plate is slightly convex. This back side curvature causes the decrease in the reflected signal and forms observed both on the THz images and on C-scan images. Weave fibers are deformed, up to delaminations. As confirmed by the 8(a) photography, plate S1, impacted at 10 J, is less damaged than S2, impacted at 20 J. The curvature of the back side is also lighter. The SEM observation shows a delamination damage ∼2 mm for the sample impacted at 10 J and about 6 to 7 mm for the second sample. The observation of the two plate cuts are in good agreement with the THz transmission images where damage widths of the same order of magnitude are observed.

Conclusion

Impact damage detection on composite fiber polypropylene/polypropylene plates was performed by THz imaging in transmission and in reflection. The results were compared with ultrasonic C-scan NDT, followed by microscopic observation of induced damage. The C-scan reveals damage parallel to the plate, thus delamination. In the case of unidirectional laminates, this damage is sharp and flat and as such it is easily detected by C-scan. In the case of fabrics, where the damage is more complex with cracks between pleats but also in the fold, between the two directions of weave called meta-delamination [START_REF] Vieille | About the impact behavior of woven-ply carbon fiber-reinforced thermoplastic-and thermosetting-composites: a comparative study[END_REF], the C-scan gives an indication about the damaged area but without specifying the nature and location of the damage. The THz reflection reverse side appears to give the same indication as the C-scan while the THz transmission appears to highlight the fiber cracks. The two THz imaging techniques are complementary and can provide a NDT with information about the location and type of the damage, crucial for understanding damage in composite structures.

Fig. 1 .

 1 Fig. 1. Experimental setup.

•

  the dynamic range, DR = 10.log Max min where Max and min are the maximum and the minimum amplitude values in the image • the signal-to-noise ratio, SNR = 10.log [ S i,j N i,j ]

Fig. 2 .

 2 Fig. 2. Damage creation.

Fig. 3 .

 3 Fig. 3. Sample S1 recto.

Fig. 4 .

 4 Fig. 4. Sample S2 recto. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

Fig. 5 .

 5 Fig. 5. Sample S1 verso. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

Fig. 6 .

 6 Fig. 6. Sample S2 verso.

Fig. 7 .Fig. 8 .

 78 Fig. 7. THz transmission image. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

Table 1

 1 Figures of merit of the THz images.

	Fig.	Test parameters	DR	SNR
	3(a)	S1 refl. recto	19.2 dB	17.8 dB
	4(a)	S2 refl. recto	19.2 dB	17.8 dB
	5(a)	S1 refl. verso	19.4 dB	13.6 dB
	6(a)	S2 refl. verso	19 dB	10.8 dB
	7(a)	S1 trans.	19.8 dB	14.7 dB
	7(b)	S2 trans.	19 dB	10.8 dB




Abstract-Recent telecommunication systems exploit the degrees of freedom of the radio propagation channel more and more. Deeper investigations of system performance require realistic channel models. For this purpose, a new MIMO UWB channel sounder is introduced. The novelty is its capacity to estimate simultaneously direction of arrival and departure in 3D, including polarization state information. The second novelty is a new and faster algorithm to estimate the multipath properties of the propagation channel. An extensive measurement campaign is performed in different residential environments and statistical results are provided.

Index Terms-Channel models, direction of arrival estimation, radio propagation, ultra wideband technology.

I. INTRODUCTION

HE radio propagation channel has been investigated for many years. Early studies were restricted to statistical path loss models. Rapidly, wideband channel models were required to design digital air interfaces. The simplest one is the tap delay line (TDL) propagation channel model. It was introduced by COST action 207 [START_REF] Cost | Digital Land Mobile Radio Communications[END_REF]. The first application was the GSM standard, with a set of channel models using 6 taps: typical urban (TU), bad urban (BU), rural area (RA), etc. TDL propagation models were tuned for most telecommunication systems. The number of taps is proportional to the system bandwidth: 18 for Wi-Fi [START_REF] Medbo | Channel Models for HIPERLAN/2 in different Indoor Scenarios[END_REF] and up to several hundred for ultra wide band (UWB) systems [START_REF] Molisch | IEEE 802.15.4a channel model -final report[END_REF]. To support new multiple input multiple output (MIMO) telecommunication systems, this model was enhanced by including direction of arrival (DOA) and direction of departure (DOD) information on each tap [START_REF] Erceg | TGn channel models[END_REF].

In the past few years, massive MIMO antenna systems have become a hot topic in wireless research. This concept aims to drastically improve spectral efficiency and reduce transmitted power [START_REF] Rusek | Scaling Up MIMO: Opportunities and Challenges with Very Large Arrays[END_REF]. Experimental research platforms are already available for realistic investigations [START_REF] Shepard | ArgosV2: A Flexible Many-Antenna Research Platform[END_REF]. To reduce complexity and improve energy efficiency, spatial modulation could be This paragraph of the first footnote will contain the date on which you submitted your paper for review. This work was supported in part by Orange Labs.

P. Pajusco is with the Microwave Department, Telecom Bretagne Engineering School, Brest, France (e-mail: patrice.pajusco@telecombretagne.eu).

N. Malhouroux is with the Orange Labs, Networks and Carriers, Orange Labs, Belfort, France (e-mail: Nadine.malhouroux@orange.com).

G. El Zein is with Institute of Electronics and Telecommunications of Rennes, Rennes, France (e-mail: ghais.el-zein@insa-rennes.fr). also another promising candidate for these large array antennas [START_REF] Di Renzo | Spatial Modulation for MIMO Wireless Systems[END_REF]. Whatever the air interface principle, future telecommunication systems will use more and more antennas. Thus, they will be able to exploit simultaneously most of the degrees of freedom of the radio channel. A comprehensive characterization of the radio channel is therefore required to enhance existing propagation channel models. This paper investigates the different dimensions of the static channel simultaneously: complex amplitude, delay, frequency, 3D direction and polarization. In this context, an original UWB-MIMO channel sounder was realized. Its novelty consists in its ability to estimate simultaneously direction of arrival and departure in 3D, including polarization state information. The second novelty is a new and faster algorithm to estimate the multipath properties of the propagation channel. All these results were used experimentally to provide new statistical results of space-time parameters of the propagation channel

The paper is divided as follows. In Section II, a new UWB channel sounder is introduced. Section III proposes a new methodology to extract deterministic and diffuse multipath components. An extensive measurement campaign in three residential environments is detailed in Section IV, in which different examples are provided. The last section is dedicated to the statistical analysis, followed by the conclusion.

II. THE UWB CHANNEL SOUNDER

Experimental investigation of UWB indoor propagation channel requires a channel sounder. An overview of measurement setups is provided in [START_REF] Pagani | Ultra-Wideband Radio Propagation Channels[END_REF] and [START_REF] Ahmadi-Shokouh | Ultra-wideband (UWB) communications channel measurementsa tutorial review[END_REF]. In most cases, a vector network analyzer (VNA) is used due to ease of implementation. However, this solution is limited to static configurations. Dynamic channel measurements require dedicated equipment. MIMO measurements are a bit more complex because several links have to be measured at the same time. One solution involves several receivers and orthogonal signals at the transmitter. Switching antennas at the receiver and/or at the transmitter is an alternative to the purely parallel structure. This solution requires a perfect synchronization between the transmitter (TX) and the receiver (RX). However, it also significantly reduces the cost and complexity of the channel sounder and makes it possible to have many antennas. If the overall duration is short enough, the measurements can be considered real-time measurements. Otherwise, the channel environment has to be static. A practical implementation, combining parallel receivers and switching transmitters, is provided in [START_REF] Kmec | A novel Ultra-Wideband Real-Time MIMO Channel Sounder Architecture[END_REF].

In this study, the UWB MIMO channel sounder investigates spatial and polarization domains simultaneously.

Comprehensive Characterization of the Double Directional UWB Residential Indoor Channel P. Pajusco, N. Malhouroux, and G. El Zein T The frequency bandwidth was set between 2.5 and 12.5 GHz. For the sake of clarity, the main features of the equipment, already introduced in [START_REF] Malhouroux | Capacity of UWB propagation channel in outdoor to indoor configuration[END_REF], are reported here. An overview of the equipment and picture of the system are depicted in Fig. 1 and Fig. 2, respectively. The channel transfer function is estimated using a VNA, which imposes a static environment during measurements. The antennas are identical at the transmitter and the receiver. They combine a virtual circular array and a 5-element array. Such an architecture enables us to electronically beamform in any direction with an arbitrary polarization [START_REF] Pajusco | Novel Array Structure for Space-Time Characterization of the UWB Channel[END_REF]. All the equipment was piloted by a laptop using Matlab™.

III. CHANNEL ESTIMATION PARAMETERS

A. High-Resolution Algorithm

The simplest way to estimate the multipath radio channel is to use linear methods. For instance, the inverse Fourier transform of the transfer function provides the channel impulse response. Similar transformations are possible to estimate physical paths in the angular domain. These methods are robust whatever the channel, but resolution is limited.

Accuracy can be improved by considering a priori knowledge about the channel. The aim of high-resolution methods is to estimate the parameters of the model accurately. The most common channel model is a set of discrete paths. In practice, the radio channel also includes a diffuse component which can reduce accuracy. In indoor environments, the ratio of the diffuse component ranges from 10 to 30 % in line of sight (LOS) [START_REF] Haneda | Double Directional UltraWideband Channel Characterization in a Line-of-Sight Home Environment[END_REF][START_REF] Tsuchiya | UWB Indoor Double-Directional Channel Sounding for Understanding the Microscopic Propagation Mechanisms[END_REF] and can reach 90 % in non-line of sight (NLOS) [START_REF] Salmi | Incorporating diffuse scattering in geometrybased stochastic MIMO channel models[END_REF]. Even if the accuracy is not as high as expected, high-resolution techniques are often used because the representation of the channel with a limited number of parameters is very convenient for channel modeling.

The MUSIC, ESPRIT, CLEAN and SAGE algorithms are the best-known high-resolution methods. The MUSIC algorithm [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF] relies on a subspace decomposition method. Model parameters are obtained by minimizing a numerical function. The resolution is costly in terms of computation and limits the number of problem dimensions to two. When the antenna array has translation invariance, ESPRIT algorithm [START_REF] Roy | ESPRIT-estimation of signal parameters via rotational invariance techniques[END_REF] is very efficient because it provides an analytical solution without requiring any iterative searches. Moreover, the complexity is linear with respect to the number of problem dimensions. However, the set of solutions in each dimension has to be paired, which is not so trivial. The CLEAN method [START_REF] Cramer | On the analysis of UWB communication channels[END_REF], introduced in astronomy, is often used in UWB. With this iterative method, the path parameters of the strongest path are estimated and removed from the original data. The CLEAN method is close to the successive interference cancelation (SIC) method. The last method of this brief overview is the SAGE algorithm [START_REF] Fleury | Channel parameter estimation in mobile radio environments using the SAGE algorithm[END_REF]. This maximum likelihood method iteratively re-estimates each path parameter, considering the other paths to be perfectly known.

The solution is improved after each step. However, this method takes time and requires an initial solution.

B. The CLEAN SVD Algorithm

There is no constraint on the shape of the antenna array when using the CLEAN algorithm. However, above a certain number of dimensions, the algorithm is no longer applicable. A new extension of the CLEAN algorithm is proposed here to deal with the large amount of data provided by this channel sounder. The basic idea is to exploit the strong relationship between the physical properties of paths and the subspace channel radiation pattern introduced in [START_REF] Pajusco | Relationship between physical and subspace propagation channels, Simulation and experimental study[END_REF]. Actually, there is no one-to-one relationship between a single physical path and a singular value of the MIMO channel matrix. For instance, the radiation pattern of a singular value can beam in different directions corresponding to different paths. However, it always beams approximately in the physical direction of one path. The proposed algorithm is detailed in Fig. 3.

The input data are the MIMO channel impulse responses. The first stage is to find the excess delay of the main path using the power delay profile (PDP). Then, the singular value decomposition (SVD) is applied to the MIMO channel at the delay. It provides the far field associated with the first singular value and enables a rough 3D detection of the direction of departure and arrival of a potential path. Due to the limited beam former resolution, a local optimization is performed to refine the direction estimates. Then, the two arrays are beamed toward these optimal directions to compute the four SISO channels corresponding to each polarization:

, , ,

VV HV VH HH

h h h h . These four channel impulse responses are interpolated along the delay to get an accurate 22 matrix gain of the path. If the amplitude is not high enough, the path is dropped and this excess delay will no longer be scanned. Otherwise, the wideband path contribution is computed and removed from MIMO input data. The procedure is repeated until there is no longer excess delay available.

C. Diffuse Multipath Component Modeling

The CLEAN-SVD algorithm provides a discrete representation of the radio channel. However, the residual between the measurement and the model is not equal to zero. This error relies on different factors. The most important factor involves the limitations of the propagation model, which is not simply a set of discrete components. However, a part of the residual also comes from the algorithm itself. For instance, the input radiation patterns are not perfectly known. The whole residual, generally called the diffuse multipath component, has to be modeled. Indeed, this component can have a major impact on channel capacity. Reference [START_REF] Czink | Including Diffuse Multipath Parameters in MIMO Channel Models[END_REF] proposes a global modeling mixing specular components and an omnidirectional diffused component. In [START_REF] Salmi | Incorporating diffuse scattering in geometrybased stochastic MIMO channel models[END_REF], authors introduce another model taking into account polarization and different clusters.

The 802.11n MIMO channel model [START_REF] Erceg | TGn channel models[END_REF] represents the diffuse multipath component by a stochastic process including correlation properties between antenna signals. The principle is briefly recapped here: The TX R and RX R matrices can be computed directly from MIMO measurements, but the result depends on the antenna array characteristics. It is suitable to derive the correlation coefficients from the power angular spectrum (PAS). Formulas are generally provided in 2D and single polarization [START_REF] Erceg | TGn channel models[END_REF], but the 3D extension is possible, as follows. G  the 3D radiation pattern of the k th element.

      1/2 1/2 T RX TX IID H R H R  (1)
              * 2 , * 00 , , sin , , TX TX m n TX mn TX TX PAS PAS G G dd G G                             (2)
In this formula, each antenna pattern has the same phase reference location, i.e. the array factor is included in the antenna pattern. This principle simplifies the coefficient computation. The same computation can be performed at the receiver. To compute the experimental value of TX  and RX  matrices, the power angular spectrum at the TX and RX were estimated using beamforming. An average was performed for all excess delays. As a result, the correlation matrices are supposed the same whatever the excess delay.

D. Model Validation

The two previous sections detailed how to extract the channel parameters of the model, the specular components and the diffuse multipath components. Then, the channel model can be generated by the sum of these two basic propagation channel models. The whole procedure is recapped in Fig. 4. Different simulations were performed to validate this approach. According to the experimental setup, the MIMO configuration was set to 180180 with 500 MHz bandwidth.

The first example considered a channel with 20 paths discrete channel model. The path delays were arbitrary distributed within 10 ns delay. The direction of departure, the direction of arrival and the polarization of each ray were chosen arbitrary. Noise was applied to reach a 20 dB signal to noise ratio. The CLEAN-SVD algorithm was applied. Contrary to SAGE algorithm, the CLEAN-SVD does not require any iteration process. The results are provided in Table I. In this simple multipath configuration, a good agreement between model (in black) and estimate (in red) is found. Direction of arrival and departure resolution is better than a degree and amplitude accuracy is better than 0.5 dB. The channel capacity was computed and the result is the same with the channel model and its estimate with CLEAN-SVD.

      search arg max M pdp            
Subspace radiation pattern G and G of the first singular value of :,:,

TX RX M h        , , arg max H-H TX RX TX RX TX RX opt opt rebuilt                    0 arg max VV VH SISO SISO HV HH SISO SISO h h h h                              max max
arg max G and arg max G The second example is more realistic. Two components with identical power level are mixed: a discrete part and a diffuse multipath component. The discrete path is composed of 3 paths. Two of them are located at the same delay. The diffuse part is spread in delay and in angular domain. According to experimental results [START_REF] Czink | Including Diffuse Multipath Parameters in MIMO Channel Models[END_REF], an exponential and a von Mises-Fisher distribution were used in delay and angular domain respectively. This DOA-DOD configuration is recapped in Fig. 5. The three crosses represent the locations of the deterministic paths. The color enables pairing between transmitter and receiver. The estimate of the discrete paths (CLEAN-SVD) and diffuse multipath components is depicted in Fig. 6. The deterministic part is estimated with the same accuracy as the previous example. The two diffuse components are found at the right place. The spread are similar but the shape is not exactly the same as the model. Due to the random process, each simulation leads to a different results but the general shape is the same.

TX RX TX RX               New Path h h h      Start Yes       2 , , , TX RX k k TX RX pdp E h k k      search Removing from M     search   Significant
To quantify the accuracy of this model, the channel capacity was computed. The water filling capacity was performed considering the whole 180180 antenna array. The result is depicted in Fig. 7. It shows that the lack of diffuse multipath component leads to underestimate the channel capacity. Contrariwise, simple uniform distribution channel model overestimates the channel capacity. Taking into account a non-uniform diffuse multipath component provides approximately the same capacity. This channel decomposition will be used in next section to analyze MIMO channel measurements.

IV. PROPAGATION MEASUREMENTS

A. Measurement Setup

Propagation measurements were carried out with the channel sounder presented in Section II. Thus, the obtained data aim to provide new statistical results to improve existing channel models. Due to the massive deployment of personal telecommunication systems, we focused on residential environments. The frequency bandwidth was set between 2.5 and 12.5 GHz. This includes the whole allocated bandwidth of UWB systems. The TX and RX antenna system of the virtual array was placed in 60 locations along a circle. Only three of the five antennas were used to reduce the measurement duration. This configuration is equivalent to a 180180 MIMO transmission and required about 4 hours to be captured. To avoid variability, measurements were performed without people in the buildings, i.e. during the night or the weekend.

B. Measurement Description

For the current study, 32 measurement locations were investigated, representing more than 1 million channel impulse responses, 30 GB of raw data and 130 hours of actual measurements. To facilitate the reuse of such data, each environment was modeled in 3D. A panoramic picture was shot at each measurement location, which provides an accurate view of the environment during measurement. This makes it possible to superimpose the direction of the arrival/departure results on the panoramic pictures, which is very convenient in analyzing spatial spreading as explained in [START_REF] Conrat | A versatile propagation channel simulator for MIMO link level simulation[END_REF]. Such an approach was first used in [START_REF] Fleury | High-resolution channel parameter estimation for MIMO applications using the SAGE algorithm[END_REF] to analyze double directional outdoor measurements. The following sections detail the three different environments.

1) Meeting room

Because of its simplicity, the meeting room environment was used to validate the channel sounder. This room (10.6 m  6.8 m) is furnished mainly with a long table surrounded by chairs. The TX and RX face each other, with one at each end of the table.

3) Residential House

It is also important to characterize the residential house environment. Moreover, statistical data collected in this type of configuration are not often reported in the literature. Measurements were performed in a house built in the 1970s. The inner walls are made of brick and weight-bearing walls are made of concrete. The receiver was set on the first floor and the transmitter was placed in various locations on each floor. This made it possible to investigate radio links between floors, which is not so frequent in the literature. The description of the environment and the 14 measurement locations are reported in Fig. 8. The E1 location was measured twice to check the consistency of the measuring equipment.

C. Measurement Examples

Such measurements contain a large amount of information. To present some examples, the most important parameters were gathered on a single page. Some explanations are useful to benefit of the full information. The configuration of the measurement is given on in the top right of the figure.

Significant parameters are provided such as TX-RX distance, link condition, excess delay spread and angular spread at TX and RX. These angular spread parameters are computed in 3D according to [START_REF] Fleury | First-and second-order characterization of direction dispersion and space selectivity in the radio channel[END_REF]. This explains that these parameters have no unit.

The

1) LOS scenario

This measurement was performed in LOS situation in the meeting room. The whole results are provided in Fig. 10. As expected, the LOS path is the strongest one and comes from the opposite antenna. Several reflections are visible especially one from a metallic whiteboard (the other coming from walls). Even in LOS situation, the average power spatial distribution highlights diffuse phenomena around main paths.

2) NLOS scenario with TX and RX on the same floor

This second example was performed in NLOS situation in house environment. The whole results are provided in Fig. 11. The LOS path does not appear on the power delay profile. The strongest path comes from diffraction on the metallic door frame. In fact, most of paths come from frame of the opened door. This clearly appears on the average power spatial distribution where the shape of the power is similar with the shape of the door. The second contribution comes from wall reflection.

3) NLOS scenario with TX and RX on different floor

This last example was performed between floors. The whole results are provided in Fig. 12. The analysis at the first delay confirms that this path corresponds to the direct path. Indeed, the LOS direction fit to the TX-RX direction. The other delay examples present more paths and spatial spreading, even if they are at a local maximum of the power delay profile. In fact, electromagnetic waves experience much more complex mechanism which explains the increase of the scattering component. Nevertheless, main power comes approximately from the LOS direction.

V. STATISTICAL RESULTS

A. Path Loss

Statistical path loss channel models are very convenient to evaluate the system link budget without any geographical data base. A wide set of models are available in the literature. Most of them are dedicated to the single frequency band. This study focuses on frequency dependence, especially when the transmission link involves several floors.

The path loss can be modeled statistically by [START_REF] Molisch | IEEE 802.15.4a channel model -final report[END_REF]. The distance and frequency are each taken into account with a power exponent parameter.

   

22 00 0 , , . The frequency dependence of the path loss depends on the physical properties of the environment. It also depends on the type of interaction. For instance, the analytical expression of the diffraction attenuation involves the frequency parameter [START_REF] Qiu | Wideband wireless multipath channel modeling with path frequency dependence[END_REF]. The antenna pattern can also modify the path loss when the frequency changes [START_REF] Pajusco | Frequency Dependence of The UWB Indoor Propagation Channel[END_REF].

n o o df L f d L f d df                (3)
To estimate the  parameter in the different indoor environments, the path loss was estimated for each frequency by averaging over the different MIMO links. The antenna effects were minimized by removing the antenna gain of the LOS at each frequency. 13. The meeting room configuration, which is not reported, is similar to the apartment configuration.

In the apartment environment, the  parameter value is 1.1, which is consistent with the IEEE 802.15.4a channel model [START_REF] Molisch | IEEE 802.15.4a channel model -final report[END_REF]. In the house environment, the value is higher. It increases from 4.3 when the TX and RX are on the same floor, to 11.9 when they are on different floors. This increase relies on the properties of the building materials, whose losses increase with frequency, especially the concrete slabs [START_REF] Tesserault | Determination of Material Characteristics for Optimizing WLAN Radio[END_REF].

B. Delay Spreading

Delay dispersion is important for wideband channel models. The median delay spread was computed in the different environments at different frequencies. Measurement results in an anechoic chamber were also performed to get reference values. In the house environment, a single curve is reported, because results are similar on the different floors. The results are depicted in Fig. 14. A linear decrease of the delay spread can be observed. The same trend with frequency bands was reported in [START_REF] Nobles | Delay spread and received power measurements within a building at 2GHz, 5 GHz and 17 GHz[END_REF][START_REF] Janssen | Propagation measurements in an indoor radio environment at 2[END_REF]. In our experiment, the slope of the delay spread parameter was estimated at -2% per GHz. For a single location, the frequency correlation is very high. In the apartment and house environments, with 4 GHz spacing, the values are 93% and 87%, respectively. Other values of frequency spacing are available in Table II. From one location to another, the delay spread can change a great deal. Its value depends on the amplitude of the direct path versus the power of the multipath components. The path loss of the direct path is mainly due to transmission through walls and concrete slabs. If the multipath component is assumed constant whatever the location, the delay spread should be correlated with the path loss. The behavior is confirmed in Fig. 15, especially in the apartment environment. On the same floor in the house environment, the trend seems to be the same, but there are not enough measurements to conclude. In the case of measurement on different floors, the direct path is highly attenuated by concrete slab transmission. Thus, just the multipath component remains and the correlation no longer exists.

C. Spatial Spreading

Spatial spreading is the key parameter in the doubledirectional channel model. Many statistical results are available in 2D in the azimuth plane [START_REF] Cost | European Co-operation in Pervasive Mobile & Ambient Wireless Communications[END_REF]. However, results in the elevation plane or at different frequencies are rare. This experiment makes it possible to provide such information. Due to the symmetry of our configuration and the small number of measurement points, results at RX and TX are mixed. Azimuth and elevation spread are depicted in Fig. 16 and Fig. 17, respectively. The higher the frequency, the lower the azimuth spread and elevation spread. If a linear relation is considered, the slope is -0.76° per GHz and -1.76° per GHz for the elevation and azimuth spread, respectively. Due to these low values, it is possible to express the decrease in terms of percentage. In this case, the slope is about -3%/GHz for both parameters.

The frequency correlation parameter is very high especially in the house environment. Values for different frequency spacing are presented in Table II.

VI. CONCLUSION

A UWB MIMO channel sounder and the CLEAN-SVD algorithm were introduced. They make it possible to extract all the characteristics of the UWB static channel: path gain (including polarization state information), path excess delay, and 3D path directions at both the transmitter and the receiver. Measurement campaigns were performed in different residential areas. Examples were detailed to highlight propagation mechanisms. This statistical analysis in residential configuration provided new results. The higher the frequency, the lower the channel spreading parameter values: -2% per GHz for delay spread and -3% for azimuth and elevation spread. A significant correlation between path loss and delay spread was highlighted when transmitter and receiver are located on the same floor. A path loss model, including frequency dependence, was also proposed. These results will help to refine existing UWB MIMO channel models in residential environments. Future work will focus on outdoorto-indoor link configurations.

1GHz f   2 GHz f   4 GHz f   6 GHz f  

Fig. 1 .Fig. 2 .

 12 Fig. 1. Overview of the UWB channel sounder

Fig. 3 .

 3 Fig. 3. Principle of the CLEAN-SVD algorithm

  the RX correlation matrix and IID H a matrix of independent zero mean, unit variance, complex Gaussian random variables.

Fig. 4 .

 4 Fig. 4. Procedure to generate the channel model

Fig. 5 .Fig. 6 .

 56 Fig. 5. Configuration of the model

Fig. 8 .

 8 Fig. 8. Measurement location in residential house environment

Fig. 9 .

 9 Fig. 9. Measurement location in residential apartment environment

  figure a) represents the power delay profile, which is computed by averaging the 180180 instantaneous channel impulse responses. The vertical red line indicates the specific excess delay values used to analyze the direction of departure and arrival. The figure b) represents four examples of spatial distribution of power at the transmitter and the receiver. The results are superimposed on the panoramic picture at the transmitter and receiver. The geometrical direction of the transmitter/receiver is also indicated, which is very convenient in NLOS situation. The figure c) presents the average over all the excess delay values. At last, the figure d) shows the main discrete components obtained from the CLEAN SVD algorithm. The color circle enables pairing TX and RX paths.

Fig. 10 .Fig. 11 .Fig. 12 .

 101112 Fig. 10. Measurement example in meeting room environment

with n the distance dependence coefficient 0 d the reference distance  the frequency dependence coefficient 0 f

 00 the center frequency

Fig. 13 .

 13 Fig. 13. Path loss versus frequency The basic frequency dependence

Fig. 14 .Fig. 15 .

 1415 Fig. 14. Delay spread parameter versus frequency

Fig. 16 .

 16 Fig. 16. Azimuth spread parameter versus frequency

TABLE I RESULTS

 I OF CLEAN-SVD WITH NOISY DISCRETE CHANNEL MODEL

	Delay	DOA	DOD
	ns					dB dB dB dB
	10.57 76.1° -63.5° 138.0° -134.2° -6.5 -4.3 -1.1 -6.3
	10.58 76.0° -63.0° 138.0° -135.0° -6.4 -4.3 -1.1 -6.3
	10.58 86.4°	76.2° 110.1° -134.7° 0.0 -2.6 -1.3 -7.9
	10.58 86.0°	75.0° 110.0° -135.0° -0.1 -2.6 -1.1 -7.6
	11.20 133.9° -176.0° 88.1° -168.5° -0.7 1.3 -2.1 0.4
	11.18 134.0° -177.0° 88.0° -168.0° -0.4 1.2 -2.6 0.3
	11.26 50.4° -26.3° 111.4° -69.2° -6.3 -8.8 -4.5 -6.3
	11,26 50.0° -27.0° 111.0° -69.0° -6.3 -8.9 -4.5 -6.3
	16.01 88.6° -43.8° 169.6° 77.0°	1.7 -10.8 -0.2 -5.0
	16.03 89.0° -45.0° 170.0° 78.0°	1.3 -10.6 -0.1 -5.0

ACKNOWLEDGMENT

The authors wish to thank Y. Chartois and A. Bazi for their support in implementing part of the acquisition and post processing software respectively, C. Sabatier for designing the UWB monopole antenna and A. Johansen for her English review.