Laurence Piche

Jean-Christophe Daigle

Rinaldo Poli

Jerome P Claverie
email: claverie.jerome@uqam.ca[

Investigation of Steric and Electronic Factors of (Arylsulfonyl)phosphane-Palladium Catalysts in Ethene Polymerization

Keywords: Phosphane ligands, Palladium, Polymerization

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

Introduction

The evolution of olefin polymerization catalysis since Ziegler's discovery in 1953 has involved a prolific coupling of polymer science with organometallic chemistry. However, there are still no commercially viable catalysts for the controlled copolymerization of simple olefins with polar functional monomers.

Currently, commercial processes for the copolymerization of ethene with polar functional monomers such as acrylates employ free radical processes which require extreme pressures and afford little or no control over polymer architecture (tacticity or crystallinity, blockiness, molecular weight distribution), and thus limit the range of material performances. A need exists for new molecular catalysts capable of polymerizing polar monomers with controlled microstructure under mild conditions. [1] A significant advance was reported by Johnson et al [2] who discovered that cationic palladium diimine complexes can copolymerize ethene and acrylates to afford branched copolymers where the acrylate is placed in a terminal position. In 2002, Drent et al [3] disclosed that an ill-defined catalytic system containing a phosphine sulfonate and a palladium complex, either tris(dibenzylideneacetone)dipalladium(0) or palladium(II) acetate, permits the preparation of ethene-acrylate copolymers where the acrylates are incorporated in main chain positions. Well-defined palladium catalysts containing a phosphine aryl sulfonate ligand were then disclosed by Hearley et al., [4] Goodall et al., [5,6] Kochi et al., [7,8] Liu et al., [9] Skupov et al., [10] Luo et al., [11] Vela et al. [12] and most recently Guironnet et al. [13] Among those reports, acrylate copolymerization with ethene was mentioned by Goodall, [5,6] Skupov, [10] and Guironnet. [13] These studies employ the catalyst (o-Ar 2 PC 6 H 4 SO 3 )PdMe(L) with Ar = o-OMePh, which corresponds to the ligand originally presented by Drent. [3] The role of the ancillary ligand L (L = pyridine, [14] lutidine, [15] DMSO, [13] allyl group [9] ) on the catalytic activity has been studied in detail. However, at this time, little is known on the influence of the aryl phosphine sulfonate structure. We recently reported that the introduction of the bulky and electron-rich aryl groups (Ar = -[o- (2',6'-(OMe) 2 C 6 H 3 )-C 6 H 4 ]) resulted in a very active catalyst which affords polyethylene of high molecular weight [10] but with a modest propensity to incorporate any other monomer than ethene. We infer that this behavior stems from the steric hindrance which precludes the facile coordination of any olefin larger than ethene. Thus, it appears that there might be a trade-off between, on one side, the high activity and high molecular weights favored by bulky and electron-rich P^O sulfonated aryl ligands and, on the other side, the propensity to incorporate polar comonomers, which is observed with less bulky phosphines. To clarify this issue, we have turned our attention toward catalysts based on polyaromatic sulfonated phosphines o-Ar 2 PC 6 H 4 SO 3 H where Ar is phenyl, naphtyl (Np), phenantryl (Pa) or anthracenyl (An). Non-sulfonated polyaromatic phosphine analogs were initially developed by Müller et al, who demonstrated that their properties are changed by altering the number of aromatic rings associated with the phosphine. [16] These phosphines become better donors as the number of aromatic rings increases and their Tolman cone angle increases from 145 o for PPh 3 to 177 o for PNp 2 Ph and 186 o for PAn 2 Ph. Thus, the larger phosphines are the better donors, and we should expect that catalysts based on the larger phosphines would be more active and generate polyethylene of high molecular weight.

Results and Discussion

The synthesis of the sulfonated arylphosphine is a one-pot procedure (Scheme 1). For the phenyl substituted phosphines 1 and 2, the dilithiated salt of benzene or toluene sulfonic acid is reacted with commercial diphenylchlorophosphine. For the other phosphines (3-5), the sulfonated benzyl group is introduced first upon reaction of trichlorphosphine with the lithiated salt. The resulting dichlorophosphine salt is not isolated, but it is reacted directly with two equivalents of the lithium salt of the desired aryl group. This procedure was found to be very rapid and reproducible as long as the benzyl sulfonic acid was sufficiently anhydrous. The n-BuLi concentration needs also to be carefully adjusted, as an excess of n-BuLi leads to the formation of n-butyl phosphines (as shown by MS), and a defect of n-BuLi leads to the isolation of phosphine oxides (R 2 POH, also shown by MS). Although the MS indicates the presence of a single compound for ligands 3 and 4, the 31 P and 13 C NMR spectra clearly show the presence of two distinct species [START_REF]None of the ligands show 31 P-1 H coupling (expected JPH ~ 500 Hz) by either 1 H or proton coupled 31 P NMR, probably because the intermolecular proton exchange is fast[END_REF] characterized by very similar spectroscopic properties. Calculations by DFT indicate that these species correspond to two rotational isomers (Figure 1) with syn and anti conformations of the aryl groups across the P atom (Figure 1), with energies differing by less than 2 kcal/mol, which is in good agreement with the 80:20 proportion found by 31 P NMR at room temperature. The activation barrier is above 10 kcal/mol, which is high enough for the structures to appear as distinct species in the NMR timescale, even at higher temperature (no coalescence was observed at T = 120 o C). Bis(phenantryl)phenyl phosphine [START_REF] Müller | [END_REF] shows only one resonance in the 31 P NMR spectrum indicating that the presence of the ortho sulfonic acid group contributes to the slow conversion between both rotamers. From the phosphine structures (optimized by DFT), we have calculated Tolman cone angles, [19] that is to say the apex angle of a cylindrical cone with origin 2.28 Å from the center of the phosphorous atom whose sides just touch the Van der Waals surfaces of the outermost atoms of the organic substituents.

The Tolman angles for the syn conformers of 3 and 4 are respectively 192 and 190 o whereas they are 206 and 207 o for the anti conformers. This is significantly higher than tert-butyl phosphine (182 o ), but slightly smaller than the highly hindered tris(2,4,6-trimethylphenyl)phosphine (212 o ). Therefore, these sulfonated aryl phosphines exhibit considerable bulk. Based on the pioneering work of Mingos, [20] the availability of the P lone pair increases when the size of the aryl group increases, thus the larger phosphines are better electron donors. Thus, the order of basicity of these phosphines is expected to be 1 ~ 2 < 3 ~ 4 < 5. It also corresponds to the observed ranking for the 31 P chemical shifts, which decreases from 4 ppm (1 and 2) to -30 ppm (5). The catalyst synthesis proceeds smoothly following the procedure highlighted in reference. [10] The yields are in the following order: 1Pd ~ 2Pd > 3Pd > 4Pd > 5Pd, which does not follow the expected basicity of those ligands, indicating that steric factors are the dominant influence in determining the reactivity of these sulfonated phosphines towards Pd centers. Catalysts 3Pd, 4Pd and 5Pd are sparingly soluble in most common solvents except DMSO. The overall structure observed for 1Pd resembles those of other (P^O)PdMe(L) complexes (Figure 2 12) in pseudo equatorial positions. This half-boat conformation has been reported for the majority of aryl sulfonate catalysts [9], [12], [15], [21] except for bulkyl aryl groups or when pyridine is replaced by DMSO. [10], [13] The complete characterization of catalysts 3Pd and 4Pd is complex because each rotamer reacts to give a separate catalyst, resulting in a doubling of all phosphine resonances (Figure 3). The analysis is further complicated by the presence of two distinct exchange processes. [14] The first one is the exchange between bound pyridine and free pyridine (if added in excess, Figure 3). The second one is the inversion of the six member ring -Pd-O-S-C=C-P-. The non-sulfonated phosphine aryl substituents are respectively occupying a pseudo equatorial and pseudo axial position. For catalyst 1Pd and 2Pd, the exchange between axial and equatorial positions is fast in the NMR time scale, even at -90 o C in CD 2 Cl 2 . For catalysts 3Pd to 5Pd, the exchange is slow. This ring inversion process is very sensitive to steric bulk : for example, for Ar = Ph All catalysts are able to polymerize C 2 H 4 at 85 o C (P = 300 psi), with the activity decreasing from 1Pd to 5Pd. The resulting polymers are highly linear, as shown by 13 C NMR and by examination of the Mark-Houwink plot in triple-detection GPC.

Unexpectedly, the drastic increase in steric hindrance from 1Pd to 4Pd results in a decrease of the average molecular weight. Thus, the least bulky and more acidic phosphine yields a catalyst with the highest activity and producing polymers with the highest molecular weights. Contrarily to what was reported by us, [10], [START_REF]a by-product during ligand workup in our[END_REF] we found that 1Pd and 2Pd are also able to copolymerize acrylates with ethylene with activities and molecular weights which are nearly identical to those obtained with MePd(pyridine)P(3-Me-6-SO 3 -C 6 H 3 )(o-OMe-Ph) 2 . For example, at P = 100 psi, T = 100ºC and for a concentration of tert-butyl acrylate of 1.70 mol/L, an insertion of 6% was obtained with catalyst 1Pd, while for similar conditions (P = 100 psi, T = 100ºC and monomer concentration = 0.85 mol/L), an insertion of 15% of tert-butyl acrylate was observed with catalyst 2Pd. However, catalysts 3Pd, 4Pd and 5Pd do not yield any copolymer under comparable conditions.

Conclusions

Phosphine sulfonate palladium complexes were prepared and used as catalysts for ethene polymerization without the need of activation. Linear polyethylenes were obtained with these catalysts, but acrylate -ethene copolymers could only be obtained with 1Pd and 2Pd. Surprisingly, the introduction of steric hindrance in the catalyst scaffold results in lower molecular weights and lower activities. The origin of these phenomena is not totally clear at this moment. We believe that several other ligand structures will need to be prepared and characterized before being able to derive structure-property relationships.

Experimental Section

General Considerations. All manipulations were done under inert atmosphere using standard Schlenk techniques. Solvents were degassed and dried over activated molecular sieves. Benzene and toluene sulfonic acid were dried by azeotropic distillation with benzene.

Dimethyl(N,N,N',N'-tetramethylethylenediamine)palladium(II), PdMe 2 (TMEDA), was prepared according to de Graaf. [START_REF] De Graaf | [END_REF] All acrylic monomers were purified by sparging them with argon and passing them over a bed of inhibitor-remover resin (Aldrich) as acrylic monomers are usually inhibited with quinones which interfere with the catalyst. The monomers were then spiked with tert-butyl catechol (0.25% wt:wt) in order to prevent spontaneous radical polymerization of the acrylate during the polymerization process. 1 H, 13 C and 31 P NMR spectra were recorded on a Varian Inova 600 MHz spectrometer at ambient temperature except for the polymers which were analyzed in deuterated tetrachlorethane at 115 o C. The molecular weight distributions were determined by gel permeation chromatography (GPC) using a Viscotek HT GPC equipped with triple detection operating at 160 o C. The eluent was 1,2,4-trichlorobenzene, and separation was performed on three PolymerLabs Mixed B(-LS) columns. The dn/dc of pure linear polyethylene was found to be 0.106 mL/g at this temperature. Electrospray mass spectra (ESI-MS) of organic compounds were recorded on a Agilent 6210 LC-MSD TOF mass spectrometer. Standard numbering of polyaromatic C and H was used below.

Preparation of ligand 1, 2-Diphenylphosphanyl-4-methylbenzenesulfonic acid. To a solution of dry toluenesulfonic acid (1.03 g, 6 mmol) in THF (30 mL) was added nBuLi 2.5 M in hexanes (4.8 mL, 12 mmol) at 0°C. After stirring for 1 h at room temperature, the solution was added dropwise to a solution of bis(phenyl)chlorophosphine (1.32 g, 6 mmol) in THF (20 mL) at 0°C. After stirring for 4 h at room temperature, the solvent was removed in vacuo leaving a white solid. The solid was dissolved in dichloromethane (50 mL) and extracted with acidic water (2 mL of concentrated HCl in 30 mL of water) and then twice with water (30 mL). The organic solvent was removed in vacuo. The product was then recrystallized from dichloromethane/diethylether at -32°C. The resulting white crystals were dried in vacuo. Yield = 0.9 g (42 %). 1 H NMR (CDCl 3 ): δ: 8.15 (s, 1H, C(P)-C(SO 3 )=CH-), 7.70-7.45 (m, 11H, H 4 -Ph, H 4 -ArSO 3 , H 2 -Ph, H 3 -Ph), 6.96 (d, J PH = 14 Hz, 1H, C(P)-CH=C(Me)), 2.28 (s, 3H, CH 3 -ArSO 3 ). 13 2-Diphenylphosphanylbenzenesulfonic acid. To a solution of dry benzenesulfonic acid (0.80 g, 5 mmol) in THF (25 mL) was added nBuLi 2.5 M in hexanes (4.0 mL, 10 mmol) at 0°C. After stirring for 1 h at room temperature, the solution was added dropwise to a solution of bis(phenyl)chlorophosphine (1.10 g, 5 mmol) in THF (15 mL) at 0°C. After stirring for 4 h at room temperature, the solvent was removed in vacuo leaving a white solid. The solid was dissolved in dichloromethane (40 mL) and extracted with acidic water (2 mL of concentrated HCl in 30 mL of water) and then twice with degased water (30 mL). The organic solvent was removed in vacuo. The product was then recrystallized from dichloromethane/diethylether at -32°C. The resulting white crystals were dried in vacuo. Yield = 0.9 g (53 %). 1 H NMR (CDCl 3 ): δ: 8.31 (s, 1H, C(P)-C(SO 3 )=CH-), 7.75-7.42 (m, 12H, H 4 -Ph, H 4 -ArSO 3 , H 5 -ArSO 3 H 2 -Ph, H 3 -Ph), 7.22 (m, 1H, C(P)-CH=CH-in ArSO 3 ). 13 C NMR (CDCl 3 ): 151.8 (CSO 3 , J PC = 12.0 Hz), 134.6 (Cipso in phenyl, J PC = 9.2 Hz), 133.5 (-C(P)-CH-in phenyl, J CP = 12.8 Hz), 133.3 (C(P)-CH=CHin C 6 H 4 SO 3 ), 132.8 (C(P)-C(SO 3 ), 129.8 (-CH-C(SO 3 )=C(P), J PC = 10.1 Hz), 129.3 (C(P)-CH=CH-CH-in phenyl, J PC = 11.9 Hz), 129.3 (C(P)-CH=CH-CH-in phenyl), 128.7 (C(P)-CH=CH-in ArSO 3 ), 128.6 (-CH=CH-C(SO 3 )=C(P)), J PC = 8.2 Hz). 31 P NMR (CDCl 3 ): δ: 4.3 (s). MS: found 342.0488, calc 342.0480.

Preparation of ligand 3, 2-(Di-naphthalen-1-yl-phosphanyl)benzenesulfonic acid.

To a solution of benzenesulfonic acid (0.8 g, 5 mmol) in THF (20 mL) was added nBuLi 2.5 M in hexanes (4.2 mL, 10.5 mmol) at 0°C. The 0.5 mmol in excess were used to quench 0.5 mmol of residual water in benzenesulfonic acid. After stirring for 2 h at room temperature, this solution was added dropwise to a mixture of PCl 3 (0.69 g, 5 mmol) in THF (20 mL) maintained at -78°C. The resulting whitish suspension was stirred for 1 h. In a separate Schlenk flask, nBuLi 2.5 M in hexanes (4 mL, 10 mmol) was added to 9-bromo-naphthalene (2.07 g, 10 mmol) in THF (30 mL) at 0°C. This mixture was left for one hour at room temperature and then introduced dropwise to the whitish suspension. After stirring for 2 h at room temperature, the solvent was removed in vacuo, leaving a purple solid. After dissolution in dichloromethane (40 mL), acidic ion-exchange resin (Amberlite IRC-50 (H) 16-50 mesh, 10 g) was added and stirred for 3 hours. The supernatant was dried in vacuo. The resulting solid, dissolved in acetonitrile, was stirred for 3 hours. After filtration, the solvent was removed. The resulting white crystals were dried in vacuo. Yield = 1.4 g (63 %). 1 To a solution of benzenesulfonic acid (0.8 g, 5 mmol) in THF (20 mL) was added nBuLi 2.5 M in hexanes (4.2 mL, 10.5 mmol) at 0°C. After stirring for 2 h at room temperature, this solution was added dropwise to a mixture of PCl 3 (0.69 g, 5 mmol) in THF (20 mL) maintained at -78°C. The resulting whitish suspension was stirred for 1 h. In a separate Schlenk flask, nBuLi 2.5 M in hexanes (4 mL, 10 mmol) was added to 9bromophenanthrene (2.57 g, 10 mmol) in THF (30 mL) at 0°C. This mixture was left for 1 h at room temperature and then introduced dropwise to the whitish suspension. After stirring for 2 h at room temperature, the solvent was removed in vacuo, leaving a purple solid. After dissolution in dichloromethane (40 mL), acidic ion-exchange resin (Amberlite IRC-50 (H) 16-50 mesh, 10 g) was added and the mixture was stirred for 3 h. The supernatant was dried in vacuo. The resulting solid, dissolved in acetonitrile, was stirred for 3 h. After filtration, the solvent was removed. The resulting pale yellow crystals were dried in vacuo. Yield = 1.0 g (37 %). 1 H NMR (DMSO-d 6 ): δ: 8.80-8.72 (m, 4H, H 6-5 -Pa), 8.65 (m, 2H, H 7 -Pa), 8.12 (m, 1H, H 3 -ArSO 3 ), 8.02 (d, J = 7.7 Hz, 1H, H 5 -ArSO 3 ), 7.92 (dd, 3 J = 7.9 Hz, 3 J = 1.2 Hz, 1H, H 4 -ArSO 3 ), 7.67 (m, 2H, H 10 To a solution of dry benzenesulfonic acid (0.92 g, 5.8 mmol) in THF (20 mL) was added nBuLi 2.5 M in hexanes (5.2 mL, 13 mmol) at 0°C. After stirring for 2 h at room temperature, the solution was added dropwise to a solution of PCl 3 (0.787 g, 5.8 mmol) in THF (20 mL) at -78°C and stirred for 1 h. In a separate Schlenk flask, nBuLi 2.5 M in hexanes (4.64 mL, 11.6 mmol) was added to 9-bromoanthracene (3.00 g, 11.6 mmol) in THF (30 mL) at 0°C. This mixture was left for 1 h at room temperature and then introduced dropwise to the whitish suspension. After stirring for 2 h at room temperature, the solvent was removed in vacuo, leaving a purple solid. After dissolution in dichloromethane (40 mL), acidic ion-exchange resin (Amberlite IRC-50 (H) 16-50 mesh, 12 g) was added and the mixture was stirred for 3 h. The supernatant was dried in vacuo. The resulting solid, dissolved in acetonitrile, was stirred 3 h. After filtration, the solvent was removed. The resulting dark yellow crystals were dried in vacuo. Yield = 1.9 g (49 %). 1 H NMR (DMSO-d 6 ): δ:8.45 (s, 2H, H 10 -An), 8.40 (dd, 3 J = 3.5 Hz, 3 J = 9.1 Hz, 4H, H 1,8 -An), 7.89 (d, J = 8.4 Hz, 4H, H 4,5 -An), 7.44 (m, 4H, H 3,6 -An), 7.31 (t, J = 7.0 Hz, 1H, C(P)-C(SO 3 )= CH-CH-), 7.22 (m, 4H, H 2,7 -An), 6.95 (t, J = 7.0 Hz, 1H, C(P)-C(SO 3 )= CH-), 6.90 (m, 2H, C(P)-CH=CH-). 13 C NMR (CDCl 3 ): δ: 153.9 (CSO 3 ), 135.3 (C(P)-C(SO 3 )=CH-), 134.6 (C(P)-C(SO 3 )=CH-CH=), 134.2 (C(P)-CSO 3 , J PC = 19.9 Hz), 134.2 (C(P)-CH=CH-in ArSO 3 , J PC = 14.9 Hz), 131.0 (C 4,5 -An), 129.1 (C 1,8 -An), 128.6 (C 8a,9a -An), 127.9 (C(P)-CH=CH-in ArSO 3 ), 127.7 (C 4a,10a -An), 126.7 (Cipso in An, J PC = 22.9 Hz), 125.4 (C 10 -An), 125.1 (C 2,7 -An), 124.6 (C 3,6 -An). 31 . PdMe 2 (TMEDA) (0.063 g, 0.25 mmol) and ligand 1 (0.089 g, 0.25 mmol) were dissolved in 10 mL of dry THF under an inert atmosphere and stirred for 30 min. Pyridine (0.0965 g, 1.25 mmol) was then added followed by stirring for another 30 min. During the stirring, a white precipitate was formed. After adding 25 mL of Et 2 O, the precipitate was collected, washed with Et 2 O and dried under vacuum. Yield = 0.100 g (72 %). 1 H NMR (CDCl 3 ): δ: 8.81 (d, J = 4.9 Hz, 2H, Hortho pyridine), 8.17 (dd, 3 J = 4.5 Hz, 3 J = 8.0 Hz, 1H, C(P)-C(SO 3 )=CH-), 7.86 (t, J = 7.7 Hz, 1H, Hpara pyridine), 6.63 (m, 4H, Hortho phenyl) 7.51 (m, 2H, Hmeta pyridine), 7.46 (m, 6H, Hmeta + Hpara phenyl), 7.33 (d, J = 7.9 Hz, 1H, C(P)-CH=), 6.80 (d, J = 9.5 Hz, 1H, C(P)-CH=C(Me)-CH-), 2.25 (s, 3H, ArCH 3 ), 0.49 (d, J PH = 2.63 Hz, 3H, Pd-Me). 13 C NMR (CDCl 3 ): δ: 150.5 (N-C=C), 146.9 (CSO 3 , J PC = 13.7 Hz), 140. PdMe 2 (TMEDA) (0.113 g, 0.44 mmol) and ligand 2 (0.152 g, 0.44 mmol) were dissolved in 10 mL of dry THF under inert atmosphere and the resulting solution was stirred for 30 min. Pyridine (0.04 g, 0.50 mmol) was then added followed by stirring for another 60 min. During the stirring, a white precipite formed. After adding 25 mL of Et 2 O, the white precipitate was collected, washed with Et 2 O and dried under vacuum. Yield = 0.110 g (81%). 1 . PdMe 2 (TMEDA) (0.063 g, 0.25 mmol) and ligand 3 (0.111 g, 0.25 mmol) were dissolved in 10 mL of dry THF under an inert atmosphere and the resulting solution was stirred for 30 min. Pyridine (0.02 g, 0.30 mmol) was then added followed by stirring for another 60 min. After adding 10 mL of Et 2 O, the purple precipitate was collected, washed with Et 2 O and dried under vacuum. Yield = 0.081 g (50%). 1 . PdMe 2 (TMEDA) (0.063 g, 0.25 mmol) and ligand 4 (0.136 g, 0.25 mmol) were dissolved in 10 mL of dry THF under an inert atmosphere and the resulting solution was stirred for 30 min. Pyridine (0.02 g, 0.30 mmol) was then added followed by stirring for another 60 min. After adding 10 mL of Et 2 O, the light brown precipitate was collected, washed with Et 2 O and dried under vacuum. Yield = 0.076 g (41%). 1 H NMR (DMSO-d 6 ): δ: 8.52 (d, J = 8.3Hz, 2H, H ortho pyridine), 8.40-8.15 (m, 4H, C(SO 3 )-CH=CH-CH=CH), 8.05-7.65 (m, 18H, HPa), 7.57 (t, J = 7.5 Hz, 1H, H para pyridine), 7.30 (dd, 3 J = 7.7 Hz, 3 J = 6.9 Hz, 2H, H meta pyridine), 0.63 (d, J PH = 3.4 Hz, 3H, Pd-Me). 13 . PdMe 2 (TMEDA) (0.063 g, 0.25 mmol) and ligand 5 (0.136 g, 0.25 mmol) were dissolved in 10 mL of dry THF under an inert atmosphere and the resulting solution was stirred for 30 min. Pyridine (0.02 g, 0.30 mmol) was then added followed by stirring for another 60 min. After adding 10 mL of Et 2 O, the yellow precipitate was collected, washed with Et 2 O and dried under vacuum. Yield = 0.069 g (37%). 1 H NMR (DMSO-d 6 ): δ: 9.26 (d, J = 9.3 Hz, 2H, Hortho pyridine), 8.96 -8.65 (m, 4H, C(SO 3 )-CH=CH-CH=CH), 8.82 (dd, 3 J = 9.0 Hz, 3 J = 3.5 Hz, 2H, H meta pyridine), 8.52-7.50 (m, 18H, HAn), 7.39 (t, J = 7.5 Hz, 1H, H para pyridine), 0.66 (d, J PH = 3.0 Hz, 3H, Pd-Me). 13 Polymerizations. Polymerizations were carried out in a stainless steel reactor (100 or 450 mL, Parr). Catalyst, toluene and eventually comonomer were added to a Schlenk flask in a nitrogenfilled glove box. The reactor, which was first dried and kept under nitrogen, was loaded with the toluene solution by cannula transfer from the Schlenk flask under nitrogen. The reactor was then sealed, pressurized with ethene, stirred and heated. The polymerizations were performed at constant pressure in the feed reactor and the activities were calculated from the rate of ethene consumption which was monitored by the decrease of the ethene pressure in the feed tank. Once the reaction was over, the reactor was cooled down to room temperature and slowly depressurized. The polymers were precipitated in four volumes of methanol, collected by centrifugation or filtered, washed with methanol and dried under vacuum.

Computational Details. All geometry optimizations were performed with the Gaussian03 suite of programs [24] using the B3LYP functional, which includes the three-parameter gradientcorrected exchange functional of Becke [25] and the correlation functional of Lee, Yang, and Parr, which includes both local and nonlocal terms. [26] The basis set chosen was the standard 6-31+G**, which includes both polarization and diffuse functions.

  Scheme 1

Figure 1 .

 1 Figure 1. Enthalpic changes between the syn and anti conformations of phosphines 3 and 4 vs the dihedral angle C1 P C2 C3 (D, indicated with stars). Only the lowest transition state TS is shown : the other TS (located at D ~ 0 o ) is at least several kcal/mol higher in energy.

  ), with the Pd atom in a square planar environment and the Me group trans to the sulfonate group. The six-member cycle Pd(1)-P(1)-C(131)=C(132)-S(1)-O(11)-adopts a half-boat conformation, with C(111) and O(13) in pseudo axial positions and C(121) and O(

  (o-OMe), the inversion barrier (measured by NMR) of PdMe(py)(PAr 2 PhSO 3 ) is 5.7 kcal/mol in CD 2 Cl 2 , whereas for Ar = Ph(o-C 6 H 3 (2,6-OMe) 2 ) the barrier is 8.3 kCal/mol in CD 2 Cl 2 .

Figure 2 .

 2 Figure 2. ORTEP view of 1Pd. Thermal ellipsoids are shown at the 50% probability level. Selected bond lengths (Å) and angles ( o ): Pd(1)-C(1) 2.057(10) ; Pd(1)-O(11) 2.164(7); Pd(1)-P(1) 2.229(3) ; Pd(1)-N(1) 2.110(8); C(1)-Pd(1)-N(1) 90.6(4); C(1)-Pd(1)-P(1) 89.1(3); C(1)-Pd(1)-O(11) 174.7(4); N(1)-Pd(1)-P(1) 172.4(3).

Figure 3 .

 3 Figure 3. Superposition of 13 C NMR spectra (downfield region). For 1Pd, bound pyridine (BP) is in rapid exchange with free pyridine (FP) whereas for 5Pd, the exchange is intermediate in the NMR time scale. For 3Pd, the exchange is slow, and the two rotamers are observed, as indicated by the C-SO3 resonances.

  C NMR (CDCl 3 ): δ: 150.0 (C-CH 3 ), 140.6 (CSO 3 , J PC = 12.4 Hz), 135.7 (C(P)-CH=C(Me)), 134.8 (Cipso in phenyl, J PC = 11.4 Hz), 134.7 (C(P)-C(SO 3 ), J PC = 11.0 Hz), 134.0 (C(P)-CH-in phenyl, J PC = 11.4 Hz), 130.1 (C(P)-CH=CH-CH-in phenyl, J PC = 13.1 Hz), 130.1 (C(P)-CH=CH-CH-in phenyl, J PC = 13.1 Hz), 129.4 (-CH-C(SO 3 )=C(P)), 129.3 (-CH=CH-C(SO 3 )=C(P)), 21.5 (ArCH 3 ). 31 P NMR (CDCl 3 ): δ: 3.6 (s). MS: found 356.0626, calc 356.0636.

  H NMR (DMSO-d 6 ): δ: 8.14 (dd, 3 J = 7.36 Hz, 3 J = 4.21 Hz, 1H, H 3 -Ar-SO 3 ), 8.10 (d, J = 7.85 Hz, 1H, H 4 -Ar-SO 3 ), 7.88-7.83 (m, 2H, H 3 -Np), 7.82 (d, J = 8.09 Hz, 2H, H 2 -Np), 7.46 (dd, 3 J = 6.28 Hz, 3 J = 3.11 Hz, 2H, H 4 -Np), 7.43-7.40 (m, 2H, H 5 -Np), 7.33 (m, 4H, H 7,6 -Np), 7.25 (m, 2H, H 8 -Np), 7.13 (m, 1H, H 5 -Ar-SO 3 ), 7.01 (m, 1H, H 6 -Ar-SO 3 ). 13 C NMR (DMSOd 6 ): δ: 154.0 (C(SO 3 ), J PC = 28.8 Hz), 136.8 (C(P)-C(SO 3 ), J PC = 20.9 Hz), 136.7 (C(P)-CH=CH-in ArSO 3 ), 135.7 & 135.5 (C 5 -Np), 133.8 & 133.7 (C 4a -Np), 133.7 (C(P)-CH=CH-in phenyl), 132.5 ((C(P)-CH=CH-CH-in phenyl), 129.2 (C 8a -Np), 129.2 (Cipso in naphtyl, J PC = 15 Hz), 128.4 (C 7 -Np), 127.8 (C(P)-C(SO 3 )-CH=, J CP = 5.0 Hz), 127.2 (C 2 -Np, J PC = 26 Hz), 126.6 (C 4 -Np), 126.4 (C 3 -Np), 126.19 & 126.21 (C 6 -Np), 134.4 (C 8 -Np). 31 P NMR (DMSO-d 6 ): δ: -23.0 (s), -26.8 (s). MS: found 442.0796, calc 442.0793. Preparation of ligand 4, 2-(Di-phenanthren-9-yl-phosphanyl)benzenesulfonic acid.

  -Pa), 7.62 (m, 2H, H 4 -Pa), 7.45 (dd, 3 J = 16.2 Hz, 3 J = 7.8 Hz, 2H, H 3 -Pa), 7.38 (d, J = 7.5 Hz, 2H, H 2 -Pa), 7.35 (d, J = 7.7 Hz, 1H, H 6 -ArSO 3 ), 7.23 (d, J = 3.3 Hz, 2H, H 1 -Pa), 7.15 (m, 2H, H 8 -Pa). 13 C NMR (DMSO-d 6 ): δ: 152.5, 134.8, 132.3, 132.2, 132.1, 132.0, 130.4, 129.8, 128.9, 128.5, 127.7, 127.6, 127.3, 127.2, 126.1, 125.7, 125.5, 125.4, 124.4, 122.0, 121.7, 121.6. 31 P NMR (DMSO-d 6 ): δ: -22.4 (s), -24,8 (s). MS: found 542.1113, calc 542.1106. Preparation of ligand 5, 2-(Di-anthracen-9-yl-phosphanyl)benzenesulfonic acid.

  2 (C-CH 3 , J PC = 6.6 Hz), 138.5 (Cipso in phenyl, br.), 135.1 (Cpara in pyridine), 134.4 (C(P)-CH=CH-in ArSO 3 , J PC = 12.1 Hz), 131.8 (C(P)-CH=C(Me)-), 131.1 (C(P)-C(SO 3 )=CH-), 130.4 (C(P)-C(SO 3 )=CH-CH-), 130.0 (C(P)-CH=CH-CH in phenyl), 128.8 (C(P)-CH=CH-CH-in phenyl, J PC = 11.1 Hz), 125.2 (Cmeta in pyridine), 21.6 (ArCH 3 ), 0.9 (CH 3 -Pd). 31 P NMR (CDCl 3 ): δ: 28.9 (s). Preparation of 2Pd, [MePd(pyridine)P(-6-SO 3 -C 6 H 3 )(Ph) 2 ].

  H NMR (CDCl 3 ): δ: 8.80 (d, J = 4.2 Hz, 2H, Hortho pyridine), 8.28 (m, 1H, -C(SO 3 )-CH-), 7.86 (t, J = 7.4 Hz, 1H, Hpara pyridine), 7.63 (m, 4H, Hortho phenyl), 7.49 (d, J = 7.1 Hz, 2H, Hmeta pyridine), 7.45 (m, 6H, Hmeta + Hpara phenyl), 7.36 (t, J = 7.4 Hz, 2H, C(SO 3 )-CP-CH, C(SO 3 )-CH=CH), 7.05 (t, J = 8.6 Hz, 1H, C(SO 3 )-CP-CH=CH), 0.50 (d, J PH = 2.4 Hz, 3H, Pd-Me). 13 C NMR (CDCl 3 ): δ: 150.2 (N-C=C), 149.2 (CSO 3 , J PC = 13.0 Hz), 138.1 (Cipso in phenyl, br.), 134.5 (Cpara in pyridine), 134.2 (CP-CH=CH-in ArSO 3 , J PC = 12.2 Hz), 130.9 (CP-CH=CH-CH in phenyl), 129.9 (CP-CH=CH in ArSO 3 ), 129.8 (PC=CH in ArSO 3 , J PC = 6.9 Hz), 129.6 (-CH-C(SO 3 )=CP), 128.7 (C(P)-CH=CH-CH-in phenyl, J PC = 11.2 Hz), 128.6 (-C(P)-C(SO 3 )-, J PC = 7.8 Hz), 125.0 (Cmeta in pyridine), 0.6 (CH 3 -Pd). 31 P NMR (CDCl 3 ): δ: 29.2 (s).

  online: ((will be filled in by the editorial staff))Entry for theTable of Contents ((Please choose one layout.)) for Table of Contents -max. 350 characters; not the same text as the Abstract))

Table 1 .

 1 Ethene polymerization data (T = 85 o C, P = 300 psi)

	Cat	[Cat]	TON	Pol. wt	Mn [a]	PDI [b]
		(µmol/L)	(molE / molPd)	(g)	(g/mol)
	1Pd	47	43 10 3	11.3	9600	1.8
	2Pd	54	24 10 3	7.4	9300	1.7
	3Pd	76	1.7 10 3	0.71	5000	1.4
	4Pd	43	14 10 3	3.4	3100	1.5
	5Pd	85	4000	0.38	3000 b) 35000	1.2 4

[a] Determined by GPC analysis at 160°C in 1,2,4 trichlorobenzene. [b] bimodal distribution.

Table 2

 2

	Cat	[Cat] (µmol/L)	[TBA] (mol/L)	TON (molE /molPd)	Mn [a] (g/mol)	PDI [a]	TBA mol% [b]
	1Pd	94	1.7	704	5170	1.4	6
	2Pd	185	0.85	0	3000	1.2	15

Ethene -tert-butyl acrylayte (TBA) copolymerization data (T = 100 o C, P = 100 psi).

[a] Determined by GPC analysis at 160 °C in 1,2,4 trichlorobenzene. [b] Determined by NMR analysis at 110 °C in tetrachlorethane d2

Preparation of 4Pd [MePd(pyridine)P(-6-SO 3 -C 6 H 3 )(phenan- threne) 2

  H NMR (DMSO-d 6 ): δ: 8.58 (d, J = 7.1 Hz, 2H, Hortho pyridine), 8.51 (d, J = 7.6 Hz, 1H, C(SO 3 )-CH-), 8.38-8.21 (m, 3H, C(SO 3 )-CH=CH-CH=CH-), 8.08-7.71 (m,14H, H-Np), 7.39 (t, J = 8.9 Hz, 1H, Hpara pyridine), 7.34 (dd, 3 J = 7.2 Hz, 3 J = 6.50 Hz, 2H, H meta pyridine), 0.62 (d, J PH = 2.8 Hz, 3H, Pd-Me).13 C NMR (DMSO-d 6 ): δ: 154.0, 153.8, 151.8, 151.6, 150.5, 150.3, 136.6, 135.6, 134.3, 133.7, 133.4, 133.1, 132.4, 130.6, 129.7, 129.3, 128.8, 128.3, 127.9, 127.5, 127.0, 126.5, 126.2, 126.0, 125.7, 0.5. 31 P NMR (DMSO-d 6 ): δ: -18.0 (s), -22.2 (s).

MePd(pyridine)P(-6-SO 3 -C 6 H 3 )(anthra- cene) 2

	C NMR
	(DMSO-d 6 ): δ: 152.3, 149.5, 135.3, 134.8, 132.3, 132.1, 130.4, 129.8, 128.9, 128.6, 127.8, 127.6, 127.3, 127.2, 126.4, 126.2, 125.7, 125.4, 124.4, 123.2, 122.5, 122.0, 121.7, 0.6. 31 P NMR (DMSO-
	d 6 ): δ: -8.88 (s), -11.45 (s).
	Preparation of 5Pd [

  C NMR (DMSOd 6 ): δ: 152.9, 152.8, 152.2, 149.9, 147.3, 139.7, 136.7, 135.9, 135.7, 135.0, 134.9, 134.6, 131.8, 131.7, 130.4, 129.9, 129.8, 129.7, 129.4, 128.7, 127.6, 127.4, 127.3, 126.8, 126.7, 126.4, 126.2, 125.9, 125.7, 125.4, 0.5. 31 P NMR (DMSO-d 6 ): δ: -19.1 (s).

Acknowledgments

We are grateful to C. Tessier at Université Laval for the crystallographic analysis. We also thank the Centre de Calcul Midi Pyrénées (CALMIP) for granting free CPU time for the theoretical calculation. This work was supported by NSERC. We thank the FQRNT for a travel fellowship to L.P. ____________


Dimitri Neaux

Barbara Blanc

Katia Ortiz

Yann Locatelli

Renate Schafberg

Anthony Herrel

Vincent Debat

Thomas Cucchi

Constraints associated with captivity alter craniomandibular integration in wild boar

Keywords: domestication, morphological integration, modularity, geometric morphometrics, cranium, skull

The domestication process is associated with substantial phenotypic changes through time.

However, although morphological integration between biological structures is purported to have a major influence on the evolution of new morphologies, little attention has been paid to the influence of domestication on the magnitude of integration. Here, we assessed the influence of constraints associated with captivity, considered as one of the crucial first steps in the domestication process, on the integration of cranial and mandibular structures. We investigated the craniomandibular integration in Western European Sus scrofa, using three-dimensional (3D) landmark-based geometric morphometrics. Our results suggest that captivity is associated with a lower level of integration between the cranium and the mandible. Plastic responses to captivity can thus affect the magnitude of integration of key functional structures. These findings underline the critical need to develop integration studies in the context of animal domestication to better understand the processes accountable for the setup of domestic phenotypes through time.

Introduction

Domestication is the ongoing process of the intensification of interactions between humans and other animals [START_REF] Vigne | The origins of animal domestication and husbandry: a major change in the history of humanity and the biosphere[END_REF][START_REF] Zeder | The domestication of animals[END_REF] associated with substantial phenotypic changes through time [START_REF] Zeder | Core questions in domestication research[END_REF][START_REF] Sánchez-Villagra | The taming of the neural crest: a developmental perspective on the origins of morphological covariation in domesticated mammals[END_REF][START_REF] Lord | The history of farm foxes undermines the animal domestication syndrome[END_REF]. Identifying the mechanisms responsible for the emergence of domestication is crucial to understand its role in the trajectories of human societies over the last 10,000 years [START_REF] Zeder | Why evolutionary biology needs anthropology: Evaluating core assumptions of the extended evolutionary synthesis[END_REF] and the emergence of humans as a new evolutionary selective force (Erlandson & Braje, 2013;[START_REF] Smith | The onset of the Anthropocene[END_REF]. The initial morphological changes associated with the first responsive steps of animal populations to anthropogenic environments, prior to the emergence of selected breeds, are largely unknown and remain unidentified. Controlling the behaviour of wild animals, where they are removed from their natural habitat and moved into an anthropogenic environment, is generally considered as a first step and a catalyst of the domestication process [START_REF] Vigne | Early domestication and farming: what should we know or do for a better understanding?[END_REF][START_REF] Zeder | Core questions in domestication research[END_REF]. Previous studies have shown that a lifetime in captivity can induce changes in the functional demands of wild animals (e.g. locomotor, foraging, or feeding behaviours), modifying the shape of craniomandibular [START_REF] Hartstone-Rose | The three-dimensional morphological effects of captivity[END_REF][START_REF] Selvey | Cranial responses to captivity in Lemur catta and Propithecus verreauxi in natural history museum collections[END_REF]Neaux et al., 2020) and postcranial bony structures [START_REF] Morimoto | Exploring femoral diaphyseal shape variation in wild and captive chimpanzees by means of morphometric mapping: a yest of Wolff's law[END_REF][START_REF] Panagiotopoulou | Foot pressure distribution in White Rhinoceroses (Ceratotherium simum) during walking[END_REF][START_REF] Harbers | The mark of captivity: plastic responses in the ankle bone of a wild ungulate (Sus scrofa)[END_REF] and that captivity can leave an anatomical print on the musculoskeletal system, beyond the phenotypic variation range observed in animals in their natural habitat.

For a comprehensive understanding of these processes, it is crucial to take into account that morphological structures, such as the cranium and the mandible, may respond to constraints in a coordinated fashion as they are morphologically integrated [START_REF] Olson | Morphological Integration[END_REF]Cheverud, 1982;[START_REF] Klingenberg | Morphological integration and developmental modularity[END_REF]. This coinheritance of character complexes (Cheverud, 1995) has been described as the consequence of shared genetic processes, developmental pathways, functional selective pressures, and/or phylogenetic constraints [START_REF] Marcucio | Mechanisms that underlie co-variation of the brain and face[END_REF][START_REF] Parsons | Mind the gap: genetic manipulation of basicranial growth within synchondroses modulates calvarial and facial shape in mice through epigenetic interactions[END_REF][START_REF] Martínez-Abadías | Chondrocranial growth, developmental integration and evolvability in the human skull[END_REF]. Morphological integration, defined as the tendency of different traits to vary jointly in a coordinated manner [START_REF] Klingenberg | Morphological integration and developmental modularity[END_REF], has been suggested as having a major influence on morphological evolution [START_REF] Wagner | Perspective: complex adaptations and the evolution of evolvability[END_REF][START_REF] Schlosser | Modularity in Development and Evolution[END_REF][START_REF] Klingenberg | Developmental constraints, modules, and evolvability[END_REF]. Indeed, a high degree of covariation between structures (i.e. strong level of integration) channels morphological variation along specific trajectories of shape space, reducing the range of potential phenotypic diversity by constraining evolution along "lines of evolutionary least resistance" (Fig. 1.a; [START_REF] Schluter | Adaptive radiation along genetic lines of least resistance[END_REF][START_REF] Marroig | Cranial evolution in sakis (Pithecia, Platyrrhini). II: Evolutionary processes and morphological integration[END_REF][START_REF] Wagner | The road to modularity[END_REF][START_REF] Goswami | The macroevolutionary consequences of phenotypic integration: from development to deep time[END_REF]. Conversely, a low covariation reduces the constraints on morphological variation. In this case, the evolution of traits is equally possible in all directions of the shape space as the different structures can respond independently to selective forces, increasing the extent of potential phenotypic diversity. A link between environmental factors, acting during the life of an individual, and modifications in the magnitude of integration has been suggested (Cheverud, 1995;[START_REF] Klingenberg | Studying morphological integration and modularity at multiple levels: concepts and analysis[END_REF]. It has been hypothesized that morphological integration is labile in response to changes in environmental conditions and that the correlations between phenotypic traits can be altered by the environment [START_REF] Pigliucci | Phenotypic Evolution: A Reaction Norm Perspective[END_REF]. Indeed, as integration results from the coordinated plastic responses of several traits to variation in environmental factors, it is likely that changes in these factors may cause variation in the magnitude of integration. In this sense, the need for studies disentangling the relationship between morphological changes, due to environmental factors, and the level of integration has already been underlined [START_REF] Klingenberg | Studying morphological integration and modularity at multiple levels: concepts and analysis[END_REF].

To assess the impact of a lifetime of growth in a captive artificial environment on morphological integration in an ungulate, we used an experimental approach focusing on the skull of wild boar.

We collected weaned wild boar piglets from a genetically homogenous population and raised them in a captive anthropogenic environment close to their initial habitat (100 km away). In this experimental farm, the piglets were separated into two groups where their natural foraging behaviour was suppressed (100 m 2 stall with no possibility of foraging) or drastically limited (3,000 m 2 pen with limited possibility of foraging due to the lack of space); they were fed primarily on processed dry food pellets, developed for pig farming. We compared the level of morphological integration in the captive wild boar specimens with wild-caught wild boar populations. The captive wild boar had little possibility to forage and were fed on a diet requiring little mechanical demands. We hypothesized that the constraints of captivity during their growth, by reducing the range of functions performed and relaxing the need for functional integration, may be linked to a significant reduction in the magnitude of integration. integration is concentrated within certain parts of a structure (the modules) but is relatively weak between these modules. Modularity therefore means that integration in a structure is compartmentalised, with strong integration within modules and weak integration between modules. Modified after [START_REF] Klingenberg | Morphological integration and developmental modularity[END_REF][START_REF] Klingenberg | Evolution and development of shape: integrating quantitative approaches[END_REF].

Material and methods

Material

The dataset was composed of 46 adult European wild boar and pig skulls belonging to four different groups (see Supporting Information Data S1). We chose specimens from a limited region (i.e. Western Europe) to reduce the confounding effects of geographic and climateinduced morphological variation known to exist in Sus scrofa (Albarella et al., 2009). The first two groups consisted of wild boar from the DOMEXP project: a multidisciplinary experiment aiming to assess the effect of captivity on the musculoskeletal system (http://anrdomexp.cnrs.fr/). To test the plastic response of mobility reduction on the shape of a wild ungulate skull, we relied on a control population of wild boar living in a 100,000 m 2 fenced forest in Urciers (France). From this population, we sampled 24 piglets that were divided into two groups of 12 specimens of equal sex ratio (6 males and 6 females). These groups were raised from 6 to 24 months at the Zoological Reserve of La Haute Touche (France) in two different contexts of mobility reduction: an indoor stall of 100 m² ('stallcaptive' group) offering no possibility of natural foraging, and a 3,000 m² wooded pen ('enclosurecaptive' group) with only limited natural foraging possible. We supplied both groups with the same processed dry food pellets, including 15.5% of raw protein adapted for pig diets. This experiment received ethics approval from the French Ministère de l'Enseignement Supérieur et de la Recherche (APAFIS#5353-201605111133847). The relatively small sample sizes for the 'stallcaptive' and the 'enclosurecaptive' groups are inherent to the experimental nature of the study. As sample size can affect the results of integration studies [START_REF] Rohlf | Use of two-block partial least-squares to study covariation in shape[END_REF]Bookstein et al., 2003), we choose to use similar sample sizes for the other studied groups. In addition to the two captive groups, we also sampled adult free-ranging specimens ('wildcaught' group). This group included four individuals from the initial free-ranging herd of Urciers, (i.e. the same population as the captive ones). These specimens came from a wild boar farm, where human interactions are intentionally kept to a minimum in order to ensure that the behaviour of the wild boar remains as natural as possible. They are free to forage for food in the woods. In addition to the specimens of the DOMEXP project, the 'wild-caught' group also included seven free-ranging wild boar from the same geographic and climatic environment (i.e. temperate central France) as the DOMEXP specimens. Like most wild boar in Western Europe, these free ranging specimens had an omnivorous diet consisting mostly of vegetable foods, e.g. acorns, roots and crops [START_REF] Schley | Diet of wild boar Sus scrofa in Western Europe, with particular reference to consumption of agricultural crops[END_REF]. All these specimens were wild-caught between one and two years of age. We included a fourth group of long-term domesticated populations of German, Polish, and French Landrace pigs ('Landrace' group), i.e. locally adapted traditional breeds [START_REF] Negri | European landrace conservation: An introduction[END_REF]. They were raised in stalls, with a strong mobility reduction, and were between one and nine years of age.

Data acquisition and analyses

We used 94 homologous landmarks and 67 semilandmarks placed on three-dimensional (3D) surfaces to describe the cranial and mandibular shape (Supporting Information Data S2). We digitised the anatomical landmarks and semilandmarks using IDAV Landmark v3.0 software [START_REF] Wiley | Evolutionary Morphing[END_REF]. We performed all the analyses in the R environment (R Core Team, 2019).

To remove variation related to their initial arbitrary position along the curves, the semilandmarks were slid along the tangent of the curves minimising bending energy [START_REF] Gunz | Semilandmarks: a method for quantifying curves and surfaces[END_REF]. These were then superimposed with the fixed landmarks using a generalised Procrustes superimposition [START_REF] Rohlf | Extensions of the Procrustes method for the optimal superimposition of landmarks[END_REF], implemented in the gpagen function of the package 'geomorph' (Adams et al., 2019) to obtain a new set of shape variables (Procrustes coordinates) and the centroid size (CS). The cranial and mandibular landmarks were subject to separate Procrustes superimpositions in order to avoid the increase of covariance and spurious results (Cardini, 2018).

Allometry is known to significantly affect the level of morphological integration as sizedependent shape changes contribute to produce integration between structures [START_REF] Klingenberg | Evolutionary covariation in geometric morphometric data: analyzing integration, modularity, and allometry in a phylogenetic context[END_REF]; therefore, we performed Procrustes ANOVAs [START_REF] Klingenberg | Geometric morphometrics of developmental instability: Analyzing patterns of fluctuating asymmetry with Procrustes methods[END_REF] with permutation procedures to quantify the allometry, with size computed as the decimal logarithm of CS (log CS;Collyer et al., 2015). This test was performed with the procD.lm function of the package 'geomorph' (Adams et al., 2019). We also tested the difference between the allometric slopes of the studied groups. Assuming that these differences were not significant, all the following analyses were computed on both the raw shape data and on size-corrected shape data, which are the residuals from the global multivariate regression of the shape against log CS, to account for the effect of allometry [START_REF] Monteiro | Multivariate regression models and geometric morphometrics: the search for causal factors in the analysis of shape[END_REF].

We performed a principal component analysis (PCA) using gm.prcomp ('geomorph') on all groups to assess the overall morphological variation and the distribution of individuals in the shape space. We evaluated the significance of shape differences among groups by performing a Procrustes ANOVA on aligned Procrustes coordinates using procD.lm.

To quantify the shape covariation, partial least squares (PLS) analyses [START_REF] Rohlf | Extensions of the Procrustes method for the optimal superimposition of landmarks[END_REF]Bookstein, 1991) and covariance ratios (CR; Adams, 2016;Adams & Collyer, 2016) were used jointly, as recommended by Adams (2016). We quantified the covariation as a proxy for the integration of cranium and mandible for each pair of axes by a correlation coefficient rPLS using integration.test ('geomorph'). This coefficient is supported by a permutation test for the null hypothesis that the distribution of specimens on one axis has no bearing on the distribution of the other axis. We computed the heatmap of shape deformations along the PLS axes to assess the location and the intensity of covariations using meshDist from the 'Morpho' package [START_REF] Schlager | Morpho: calculations and visualisations related to geometric morphometrics[END_REF]. In addition, differences in integration patterns were assessed by examining the general orientation of each group's distribution on the PLS scores [START_REF] Mitteroecker | The evolutionary role of modularity and integration in the hominoid cranium[END_REF][START_REF] Singh | Morphological evolution through integration: a quantitative study of cranial integration in Homo, Pan, Gorilla and Pongo[END_REF][START_REF] Neaux | Morphological integration of the cranium in Homo, Pan, and Hylobates and the evolution of hominoid facial structures[END_REF]. For this purpose, we tested for differences in the regression slopes between the studied groups on the between-group PLS. We assessed the overall modularity between cranium and mandible modules using the CR from modularity.test ('geomorph'). Modularity exists if integration is compartmentalised, i.e. concentrated within certain parts of a structure (the modules) but relatively weak between modules (Fig 1 .b; see Supporting Information Data S3). The value of CR provides a measure for characterising and evaluating the degree of modularity in biological data sets (Adams, 2016;Adams & Collyer, 2016). Morphological integration and modularity were assessed including all groups (between-group covariation) and within groups (within-group covariation).

Results

Variation analyses

Allometry explains nearly 20% of the shape variation in the cranium (p < 0.01; 19.74% of the total variance) and the mandible (p < 0.01; 17.39% of the total variance). The allometric slopes did not differ between the studied groups for the cranium (p = 0.44) or the mandible (p = 0.16).

In addition to raw shape, we computed the size corrected shape variables for further analyses.

On the PCA, PC1 accounted for 56.41% and 31.24% of the total variance for the cranium and mandible respectively (Supporting Information Data S4). For both structures, PC1 was driven by the strong divergence between the wild boar phenotype towards the negative side of the axis and the Landrace pigs towards the positive side. For the cranium, PC2 mainly separates the wild-caught from the captive wild boar. It is noteworthy that the plastic effect displayed on PC2 is different from the shape divergence between the wild boar and pigs, displayed on PC1, as the two shape changes are located on different PCs. We found significant (p < 0.05) pairwise differences of raw cranial shapes between all groups and the 'Landrace' but not between the wild boar groups (Supporting Information Data S4). We found the same results for the allometry-free cranial shapes. We found significant (p < 0.05) pairwise differences of mandibular raw shapes between all groups except between the 'stallcaptive' and 'enclosure captive' groups. For the allometry-free mandibular shapes, the difference between the 'stall captive' and the 'wild-caught' groups was also not significant.

Between-group covariation analyses

The correlation coefficient of the first pair of PLS axes (PLS1) between the cranium and the mandible for all the studied specimens is strong and significant for raw (rPLS = 0.89; p < 0.01; Table 1, Fig. 2.a) and allometry-free shapes (rPLS = 0.88; p < 0.01). The PLS1 pairs of axes account respectively for 86.65% and 92.67% of the total covariation. The main deformation associated with PLS1 is located in the anterior part of the nasal, in the nuchal crest region, in the zygomatic process of the frontal and in the tip of paroccipital processes for the cranium (Fig. 2.b). For the mandible, they are visible in the maximum of curvature between the mandibular ramus and corpus, in the inner part of the gonial angle region, on the insertion of the lower canines, and on the ventral part of the symphysis. The regression slopes between the studied groups were not different between the studied groups for the cranium (p = 0.44) or the mandible (p = 0.16) between-group PLS 1. The CR for all the studied specimens indicates a significant modularity between the cranium and mandible for raw (CR = 0.81; p < 0.01; Table 1) and allometry-free shapes (CR = 0.70; p < 0.01).

Within-group covariation analyses

The PLS computed for each studied group showed a significant level of integration for the 'enclosurecaptive' and 'wild-caught' groups for raw shapes (rPLS = 0.89; p = 0.04; Table 1) and (rPLS = 0.97; p < 0.01), and allometry-free shapes (rPLS = 0.88; p = 0.04) and (rPLS = 0.97; p < 0.01). The correlation coefficients of PLS1 are not significant for the 'stallcaptive' and 'Landrace' groups. The main deformation associated with PLS1 includes important changes in the anterior extremity of the rostrum, the occipital region, the lateral side of the ramus and the symphysis region for both the 'enclosurecaptive' (Fig. 3.a) and 'wild-caught' (Fig. 3.b) groups. Deformations include changes in the ventral edge of the zygomatic arch and in the pterygoid fossa region for the 'wild-caught' group. The CR values for all the studied groups indicates a significant modularity between the cranium and the mandible (Table 1).

Discussion

Our analyses confirm that captivity imposed on wild boar during their growth is linked to a reduction in the magnitude of integration. The results obtained from allometry-free shape data reveal similar tendencies, indicating a relatively low impact of allometry on patterns of covariation. This result underlines that changes in environmental factors can affect the magnitude of integration. Previous results on the same experimental sample (DOMEXP project) found that the shape of cranium and mandible are affected by changes in the functional demands associated with captivity [START_REF] Neaux | How changes in functional demands associated with captivity affect the skull shape of a wild boar (Sus scrofa)[END_REF]. Indeed, modifications in foraging and feeding behaviours have been identified as potential factors able to modify skull shape.

Furthermore, morphological integration between the cranium and the mandible is considered as a classic example of functional integration, where two structures interact in the same functional context [START_REF] Klingenberg | Studying morphological integration and modularity at multiple levels: concepts and analysis[END_REF]. Indeed, the upper and lower jaws need to be coordinated to achieve proper occlusion and perform functions, such as biting and chewing [START_REF] Hautier | Patterns of covariation in the masticatory apparatus of hystricognathous rodents: implications for evolution and diversification[END_REF]Figueirido et al., 2013). Therefore, our results show that captivity, inducing changes in foraging and feeding behaviour which likely reduce the need for functional integration [START_REF] Neaux | How changes in functional demands associated with captivity affect the skull shape of a wild boar (Sus scrofa)[END_REF], also diminishes the magnitude of integration between the cranium and the mandible, i.e. the structures performing these functions. In this sense, several studies have empirically shown that morphological integration can be highly variable over short timescales in response to environmental changes acting on shared developmental and functional processes (Beldade et al., 2002;[START_REF] Young | Serial homology and the evolution of mammalian limb covariation structure[END_REF][START_REF] Monteiro | Adaptive radiations, ecological Specialization, and the evolutionary integration of complex morphological structures[END_REF].

In our study, most of the wild-caught and captive wild boar groups did not display significant differences in terms of shape disparity. This similarity underlines that though captivity modifies functional demands in wild animals [START_REF] Neaux | How changes in functional demands associated with captivity affect the skull shape of a wild boar (Sus scrofa)[END_REF][START_REF] Hartstone-Rose | The three-dimensional morphological effects of captivity[END_REF][START_REF] Harbers | The mark of captivity: plastic responses in the ankle bone of a wild ungulate (Sus scrofa)[END_REF], it does not affect their potential range of morphological variation.

Our analyses also showed that, as for the group of captive wild boar raised in a stall, integration is also not significant for the group composed of Landrace pigs. These traditional breeds of pigs share several features with the captive wild boar from our experiment. They were given daily rations, mainly composed of agricultural products and food waste, allowing the relaxation of environmental constraints associated with the necessity to find and process food. This relaxation in one of the main functions performed together by the cranium and mandible (i.e. mastication), may result in a weaker morphological integration between these structures in Landrace pigs, as well as in captive wild boar raised in a stall. Furthermore, these two groups share the impossibility to perform foraging and rooting as they were both raised in stalls, i.e.

on artificial solid grounds. When possible, foraging and rooting are activities that both wild boar (Blasetti et al., 1988) and pigs (Buckner et al., 1998) spend a lot of time doing.

Modifications in rooting frequency, impacting the development of the muscles in the neck regions, may be associated with changes in cranial shape [START_REF] Owen | The zooarchaeological application of quantifying cranial shape differences in wild boar and domestic pigs (Sus scrofa) using 3D geometric morphometrics[END_REF]. Therefore, the impossibility for both captive wild boar raised in a stall and Landrace pigs to perform such functions may also explain the non-significant integration between the cranium and the mandible observed in these two groups. This confirms that a reduction in the range of functions available is linked to a significant reduction in the magnitude of integration. Although we found differences in the integration level between the studied groups, we did not find differences in integration patterns, suggesting that changes in constraints due to captivity affect the level of covariation between structures but not the way they covary. This result was expected, as previous studies have shown that integration patterns are fairly conservative, even at high taxonomic levels [START_REF] Goswami | Morphological integration in the carnivoran skull[END_REF][START_REF] Porto | The evolution of modularity in the mammalian skull I: morphological integration patterns and magnitudes[END_REF][START_REF] Neaux | Basicranium and face : assessing the impact of morphological integration on primate evolution[END_REF].

For both the between-group and within-group analyses, we also found significant modularity between the cranium and the mandible, corroborating the presence of two basic independent phenotypic modules in the skull (one cranial and one mandibular). The modularity between the cranium and mandible is likely explained by their respective functional roles. Indeed, whereas the morphology of the mandible is closely related to feeding behaviour [START_REF] Taylor | Feeding behavior, diet, and the functional consequences of jaw form in orangutans, with implications for the evolution of Pongo[END_REF]Daegling & McGraw, 2007;Anderson et al., 2014), the shape of the cranium is also affected by a multiplicity of other functions unrelated to food consumption (e.g. vision, respiration, mastication, brain protection; [START_REF] Lieberman | The Evolution of the Human Head[END_REF]. Our results confirm that even if the cranium and the mandible can be considered as two distinct modules (i.e. integration is stronger within these structures than between them), there is still a significant relationship between them, at least in our between-group analysis, which can be defined as intermodule integration [START_REF] Klingenberg | Cranial integration and modularity: insights into evolution and development from morphometric data[END_REF].

Conclusion

Our results support the hypothesis that behavioural changes associated with captivity, considered a catalyst of the animal domestication process [START_REF] Vigne | Early domestication and farming: what should we know or do for a better understanding?[END_REF][START_REF] Zeder | Core questions in domestication research[END_REF], do result in a reduction of the integration between the cranium and the mandible. However, this work will need to be expanded further using a greater dataset, as the relatively small sample size, inherent to the experimental nature of our study, could have partly biased these results [START_REF] Thiese | P value interpretations and considerations[END_REF]. Further studies would also help confirm our results that a weak integration could be the morphological response to anthropogenic changes in the functional demands associated with captivity, constituting possible future new markers for the domestication process that could be explored in the archaeological record.

Rabier, Patrick Roux and Colin Vion) for their help during the setup of the experimental structures, the care they provided to the experimental specimens and during the data acquisition. 1 Italicized sexes and ages were estimated based on osteological observations, using respectively the morphology of canine cross section (Mayer & Brisbin, 1988) and the mandibular tooth eruption and wear stages in occlusal view (Grant, 1982).

Daegling

Figure 1 .

 1 Figure 1. (a) A high integration level channels morphological variation, reducing the range of

Figure 2 .

 2 Figure 2. (a) First pair of partial least squares analysis axes (PLS1) between cranial and

Figure 3 .

 3 Figure 3. Heatmap of the shape covariation intensity of partial least squares analysis axes (PLS

  Figueirido B, Tseng ZJ, Martín-Serra A (2013) Skull shape evolution in durophagous carnivorans. Evolution 67,1975-1993.

  2 Abbreviations: MNHN = Muséum national d'Histoire naturelle (Paris, France); ZNS = Zentralmagazin Naturwissenschaftlicher Sammlungen (Halle/Saale, Germany). Grant A. 1982. The use of tooth wear as a guide to the domestic ungulates. In: Wilson B, Grigson C, Payne S, eds. Ageing and Sexing Animal Bones from Archaeological Sites, UK: BAR British Series, 991-108 Mayer JM & Brisbin IL. 1988. Sex identification of Sus scrofa based on canine morphology. Journal of Mammalogy 69:408-412 Data S2: Digitisation and definitions of landmarks a. Digitisation protocol All specimens were scanned using a Computed Tomography (CT) scanner with a spatial resolution of between 100 and 500 μm. The five wild boar from Urciers were scanned as live specimens at the Chirurgie et Imagerie pour la Recherche et l'Enseignement (CIRE) platform of the Institut National de Recherche pour l'Agriculture, l'Alimentation et l'Environnement (INRAE). Other individuals were scanned as dry specimens using a CT scanner close to the collections they were housed in. We segmented the bones using the segmentation tools of the Avizo v8.0 software, and then converted the volumes into a three-dimensional PLY format. b. Definitions of cranial (1-70) and mandibular (71-94) homologous landmarks. of supraorbital process of the frontal bone 27, 28 Meeting point of the parietal-frontal suture and temporal line 29, 30 Most anterior, dorsal point of the zygomatic process of the squamosal Most anterior point of the cheek-tooth row (excluding P1) 43, 44 Most posterior point of the cheek-tooth row 45 Most posterior point of the posterior nasal spine on the palatine bone 46, Most ventral, lateral point of the pterygoid process of the sphenoid 48, Most posterior point of the pterygoid hamulus 50, Meeting point of the pterygoid process with the ridge of the lateral pterygoid plate 52, Meeting point of the pterygoid hamulus with the ridge of the medial pterygoid plate 54 Most posterior point of the vomer in contact with the sphenoid 55, Most ventral, lateral, posterior point of the sphenoid-squamosal suture 57, Most ventral, medial, posterior point of the sphenoid-squamosal suture 59, Most posterior, medial point of the petro-occipital fissure 61, Most lateral point of the occipital condyle 63 Most anterior, ventral midline point of the premaxilla 64 Most posterior midline point of the nuchal crest 65, Most posterior, lateral point of the nuchal crest 67, Most lateral point of the foramen magnum 69 Most posterior, dorsal point of the foramen magnum 70 Most anterior point, ventral of the foramen magnum 71, Most anterior, lateral point of the lower canine alveolus 73, Most anterior point of the cheek-tooth row (excluding P1) 75, Most lateral point at the maximum of curvature between the mandibular ramus and corpus 77, Most lateral point at the maximum of curvature between the coronoid process and boar (Sus scrofa) cranium showing the homologous landmarks (red dots) and semilandmarks (blue dots) used in the study. e. Wild boar (Sus scrofa) mandible showing the homologous landmarks (red dots) and semilandmarks (blue dots) used in the study.Data S3: Summary of the covariance ratio (CR) testThe CR coefficient is compared to a distribution of values obtained by randomly assigning landmarks into subsets. A significant modular signal is found when the observed CR coefficient is small relative to this random distribution. A CR value between zero and one (e.g. CR 1) indicates that the degree of covariation is higher within than between modules, characterising a modular structure. A CR value close to one (e.g. CR 2), within a distribution of values obtained by randomly assigning landmarks into subsets (histogram), describes a structure where the covariations within and between modules are similar, characterising a random set of variables. A CR value larger than one (e.g. CR 3) defines a greater covariation between than within modules, i.e. an integrated structure. Modified after Adams (2016).Adams, D.C. 2016. Evaluating modularity in morphometric data: challenges with the RV coefficient and a new test measure. Methods in Ecology and Evolution. 7: 565-572. mandible, the divergence from wild to domestic animals was characterised by three main traits: (1) a taller and more upright ramus, (2) an anteroposteriorly shorter and taller corpus and (3) a reduced mandibular angle. The cranial shape changes along PC2, from wild-caught to captive wild boar involves three main shifts: (1) an anteroposteriorly longer cranium, (2) more robust zygomatic arches and (3) a more concave cranium. The mandible shape change along PC2 was characterised by (1) a decrease of the corpus length, (2) a taller ramus and (3) a wider mandible in the superior view. b. Pairwise ANOVA distance and p-values of Procrustes coordinates computed for the cranium and mandible.

Table 1 .

 1 Values of PLS, covariance ratios and coefficients for raw shapes and allometry-free shapes. rPLS: PLS coefficient of the first pair of PLS axes, %EC: percentage of covariation explained by the first pair of PLS axes, CR: Covariance Ratio.

		rPLS	p-value	%EC	CR	p-value
	raw shapes
	all groups	0.89	< 0.01	85.65	0.81	< 0.01
	stall -captive	0.82	0.51	59.33	0.71	< 0.01
	enclosure -captive	0.89	0.04	76.74	0.84	< 0.01
	wild-caught	0.97	< 0.01	72.86	0.95	< 0.01
	Landrace	0.88	0.06	69.31	0.81	< 0.01
	allometry free-shapes
	all groups	0.88	< 0.01	92.67	0.70	< 0.01
	stall -captive	0.89	0.29	56.90	0.88	< 0.01
	enclosure -captive	0.88	0.04	76.90	0.84	< 0.01
	wild-caught	0.97	< 0.01	56.14	0.97	< 0.01
	Landrace	0.84	0.52	32.16	0.88	< 0.01
	Significant values (p < 0.05) are in bold.

  b. List of specimens. M: male, F: female.

	2017-559 2017-579 S_bay_lds_4	F F M	18 33	enclosure -captive wild-caught Landrace	Réserve de la Haute Touche Chambord Bayerisches Landschwein (German)	MNHN MNHN ZNS
	Catalogue number 2017-561 2017-580	Sex 1 Age M F 18	Status enclosure -captive wild-caught	Location/Breeds Réserve de la Haute Touche Chambord	Location 2 MNHN MNHN
	2017-563 2017-581	M F	(months) 1 19	enclosure -captive wild-caught	Réserve de la Haute Touche Chambord	MNHN MNHN
	2017-557	F		stall -captive	Réserve de la Haute Touche	MNHN
	2017-564 COMP_2013-1262 F M	16-18	enclosure -captive wild-caught	Réserve de la Haute Touche Compiègne	MNHN MNHN
	H285	M		stall -captive	Réserve de la Haute Touche	MNHN
	2017-565 COMP_2013-1265 M F	21	enclosure -captive wild-caught	Réserve de la Haute Touche Compiègne	MNHN MNHN
	2017-560	M		stall -captive	Réserve de la Haute Touche	MNHN
	2017-566 COMP_2013-1273 M F	36-60	enclosure -captive wild-caught	Réserve de la Haute Touche Compiègne	MNHN MNHN
	2017-562	M		stall -captive	Réserve de la Haute Touche	MNHN
	2017-567 S_bay_lds_1	F F	13	enclosure -captive Landrace	Réserve de la Haute Touche Bayerisches Landschwein (German)	MNHN ZNS
	2017-555	F		stall -captive	Réserve de la Haute Touche	MNHN
	2017-568 S_bay_lds_3	F M	33	enclosure -captive Landrace	Réserve de la Haute Touche Bayerisches Landschwein (German)	MNHN ZNS
	2017-556	F		stall -captive	Réserve de la Haute Touche	MNHN
	2017-570 S_hv_br_6	F F	22	enclosure -captive Landrace	Réserve de la Haute Touche Hannover-Braunschweig Landschwein (German)	MNHN ZNS
	2017-569	F		stall -captive	Réserve de la Haute Touche	MNHN
	2017-572 S_hv_br_9	M F	51	enclosure -captive Landrace	Réserve de la Haute Touche Hannover-Braunschweig Landschwein (German)	MNHN ZNS
	H319	F		stall -captive	Réserve de la Haute Touche	MNHN
	2017-573 S_kr1	M F	18-20	enclosure -captive Landrace	Réserve de la Haute Touche French (Corsican)	MNHN ZNS
	2017-554	F		stall -captive	Réserve de la Haute Touche	MNHN
	PRA_172 S_kr2	F M	18-20	wild-caught Landrace	Urciers French (Corsican)	MNHN ZNS
	2017-571	M		stall -captive	Réserve de la Haute Touche	MNHN
	2017-583 S_pol_2	M F	36-60	wild-caught Landrace	Urciers Polnisches Landschwein (Polish)	MNHN ZNS
	2017-574	M		stall -captive	Réserve de la Haute Touche	MNHN
	2017-585 1850-435	F F	16-18	wild-caught Landrace	Urciers French	MNHN MNHN
	2017-575	M		stall -captive	Réserve de la Haute Touche	MNHN
	PRA_188 1860-43	F M	16-18	wild-caught Landrace	Urciers French	MNHN MNHN
	2017-558	M		enclosure -captive	Réserve de la Haute Touche	MNHN
	2017-577 DUP_C	M M	16-18	wild-caught Landrace	Chambord French	MNHN MNHN

Table 1 .

 1 Values of PLS, covariance ratios and coefficients for raw shapes and allometry-free shapes. rPLS: PLS coefficient of the first pair of PLS axes, %EC: percentage of covariation explained by the first pair of PLS axes, CR: Covariance Ratio.

	enclosure -captive	wild-caught	Landrace

Acknowledgements

We are most grateful to the staff of the Réserve Zoologique de la Haute-Touche (Christophe Audureau, Jérémy Bernard, Christophe Jubert, Sandrine Laloux, Emmanuel Maréchal, Régis

We thank the CIRE platform at INRAE (Hans Adriensen, Frédéric Elbout, Christian Moussu and Luc Perrigouard), CT teams from the Leibniz-IZW in Berlin (Guido Fritsch and Juliane Kühne) and the University Hospital Halle/Saale (Silvio Brandt). We would like to thank Jill Cucchi for copy editing this manuscript. We are grateful to two anonymous reviewers for their valuable comments and advices on an earlier version of the manuscript. This research was funded by ANR, through the DOMEXP project (ANR-13-JSH3-0003-01), the LabEx ANR-10-LABX-0003-BCDiv, in the programme 'Investissements d'avenir' ANR-11-IDEX-0004-02, programme Emergence SU-19-3-EMRG-02, Muséum national d'Histoire naturelle (Paris), and the CNRS INEE institute. We have no conflict of interest to declare.

Adams DC (2016) Evaluating modularity in morphometric data: challenges with the RV coefficient and a new test measure. Methods Ecol Evol 7, 565-572. Adams DC, Collyer M, Kaliontzopoulou A (2019) geomorph: geometric morphometric analyses of 2D/3D landmark data. Adams DC, Collyer ML (2016) On the comparison of the strength of morphological integration across morphometric datasets. Evolution 70, 2623-2631. Albarella U, Dobney K, Rowley-Conwy P (2009) Size and shape of the Eurasian wild boar (Sus scrofa), with a view to the reconstruction of its Holocene history. Environmental Archaeology 14, 103-136. Anderson PS, Renaud S, Rayfield EJ (2014) Adaptive plasticity in the mouse mandible. BMC Evolutionary Biology 14, 85. Beldade P, Koops K, Brakefield PM (2002) Developmental constraints versus flexibility in morphological evolution. Nature 416, 844-847. Blasetti A, Boitani L, Riviello MC, et al. (1988) Activity budgets and use of enclosed space by wild boars (Sus scrofa) in captivity. Zoo Biology 7, 69-79. Bookstein FL (1991) Morphometric Tools for Landmark Data: Geometry and Biology, Cambridge. UK: Cambridge University Press. Bookstein FL, Gunz P, Mitteroecker P, et al. (2003) Cranial integration in Homo: singular warps analysis of the midsagittal plane in ontogeny and evolution. J Hum Evol 44, 167-187. Buckner LJ, Edwards SA, Bruce JM (1998) Behaviour and shelter use by outdoor sows. Applied Animal Behaviour Science 57, 69-80. Cardini A (2018) Integration and Modularity in Procrustes Shape Data: Is There a Risk of Spurious Results? Evol Biol .46, 90-105. Cheverud JM (1982) Phenotypic, genetic, and environmental morphological integration in the cranium. Evolution 36, 499-516. Cheverud JM (1995) Morphological integration in the saddle-back tamarin (Saguinus fuscicollis) cranium. The American Naturalist 145, 63-89. Collyer ML, Sekora DJ, Adams DC (2015) A method for analysis of phenotypic change for phenotypes described by high-dimensional data. Heredity 115, 357-365.

Data Availability

The data that support the findings of this study are available from the corresponding authors upon reasonable request.


Ph Blanchard

D Gandolfo

L Laanait

J Ruiz

H Satz

On the Kertész line: Thermodynamic versus Geometric Criticality

Keywords: pacs: 05.50.+q, 64.60.C, 75.10.H, 05.70.Fh, 05.10Ln

The critical behaviour of the Ising model in the absence of an external magnetic field can be specified either through spontaneous symmetry breaking (thermal criticality) or through cluster percolation (geometric criticality). We extend this to finite external fields for the case of the Potts' model, showing that a geometric analysis leads to the same first order/second order structure as found in thermodynamic studies. We calculate the Kertész line, separating percolating and nonpercolating regimes, both analytically and numerically for the Potts model in presence of an external magnetic field.

Introduction

The critical behaviour in certain spin systems, such as the Ising model, can be specified in two equivalent, though conceptually quite different ways. In the absence of an external magnetic field, decreasing the temperature leads eventually to the onset of spontaneous symmetry breaking and hence to the singular behaviour of derivatives of the partition function. On the other hand, the average size of clusters of like-sign spins also diverges at a certain temperature, i.e., there is an onset of percolation. The relation between these two distinct forms of singular behaviour has been studied extensively over the years, and it was shown that for the Ising model on the lattice d , with d ≥ 2, implemented with a suitable cluster definition using temperature dependent bond weights, the two forms lead to the same criticality: the critical temperatures T c as well as the corresponding critical exponents coincide in the two formulations [7,4].

In the presence of an external field H, the Z 2 symmetry of the Ising model is explicitly broken and hence there is no more thermodynamic critical behaviour. Geometric critical behaviour persists, however; for T ≤ T p (H), there is percolation, while for T > T p (H), the average cluster size remains finite. In the T -H plane, there thus exists a line T p (H), the so-called Kertész line, separating a percolating from a non-percolating "phase" [11]. Given the mentioned correct cluster definition, it starts at T p (0) = T c , i.e., at the thermodynamic critical point.

We want to show here that the equivalence of thermodynamic and geometric critical behaviour can be extended to the case H = 0. Since in the case of continuous thermodynamic transitions, such as those of the Ising model, the introduction of an external field excludes singular behaviour (for the case d = 2 this is shown analytically [9,6]), our problem makes sense only for first order transitions, for which the discontinuity remains over a certain range of H, even though for H = 0, the symmetry is broken. The ideal tool for such a study is the q-state Potts' model on a lattice d , with q ≥ 3 and d ≥ 3. In this case, we have a thermodynamic phase diagram of the type shown in Fig. 1a, with a line of first order transitions starting at T c (0) and ending at a second order point T c (H c ) [10]; the transition at this endpoint is found to be in the universality class of the 3-d Ising model. In terms of the energy density ǫ(H) of the system (the energy per lattice volume), the phase diagram has the form shown in Fig. 1b; for H = 0, the coexistence range ǫ 2 ≤ ǫ(0) ≤ ǫ 1 corresponds to the critical temperature T c (0). The average spin m(ǫ) as order parameter vanishes for ǫ ≥ ǫ 1 and becomes finite for smaller ǫ. We want to show that in the temperature range T c (0) ≤ T (H) ≤ T c (H c ), the corresponding Kertész line T p (H) (see Fig. 1c) coincides with that of the thermal discontinuity and that it also leads to the same first order/second order phase structure. Let us begin with a conceptual discussion of the situation.

The q-state Potts' model in the absence of an external field provides q + 1 phases: the disordered phase at high temperature and q degenerate

H T T (0) c c c H ε ε ε 1 2 T (H ) H c H c (a) (b) non-percolating T H percolating T (0) c Kertesz Line T (H) p (c)
Figure 1: Thermodynamic and geometric phase structure for a first order transition ordered low-temperature phases. Spontaneous symmetry breaking has the system fall into one of these as the temperature is decreased. Turning on a small external field H aligns the spins in its direction and thus effectively removes the q -1 "orthogonal" low-temperature phases. Hence now only two phases remain: the ordered low-temperature state of spins aligned in the direction of H, and the disordered high-temperature phase. The two are for T < T c (H c ) separated by a mixed-phase coexistence regime. At the endpoint T = T c (H c ), there is a continuous transition from a system in one (symmetry broken) ordered phase to the corresponding (symmetric) disordered phase. The behaviour at H = H c in Fig. 1b is thus just that of the Ising model, and hence the endpoint transition is in its universality class.

In the geometric formulation for H = 0, with decreasing temperature or energy density there is formation of finite clusters of q different orientations; the clusters here are defined using the temperature-dependent F-K bond weights. At ǫ(0) = ǫ 1 , the Z q symmetry is spontaneously broken: for one of the q directions, there now are percolating clusters, and the percolation strength P (ǫ) becomes finite for ǫ < ǫ 1 . However, the disordered phase also still forms a percolating medium (for d ≥ 3). A further decrease of the energy density reduces the fraction of space in disordered state, and for ǫ(0) ≤ ǫ 2 , there is no more disordered percolation. Embedded in the disordered phase are at all times finite clusters of a spin orientation "orthogonal" to the one chosen by spontaneous symmetry breaking. In our treatment, we will therefore divide the set of clusters into three classes: disordered, ordered in the direction of symmetry breaking, and ordered orthogonal to the latter.

While for H = 0, any of the q directions could be the given orientation, for H = 0, the external field specifies the alignment direction, making the q -1 sets of "orthogonal" clusters essentially irrelevant. It is for this reason that at the endpoint of a line of first order transitions one generally encounters the universality class of the Ising model. Whatever the original symmetry of the system was, at the endpoint there remains only the aligned and the disordered ground states. The plan of the paper is as follows. In the next section, we recall the cluster treatment of the Potts' model and specify our method to identify the different cluster types. This will be followed by an analytic study valid for small external fields and by numerical calculations for different q up to asymptotic values of H. Formal details of the analytic calculation are given in the appendix.

The model

We consider a finite-volume q-state Potts model on the lattice d (d ≥ 2), at inverse temperature β = 1/T and subject to an external ordering field h. It is defined by the Boltzmann weight

ω Potts (σ) = i,j e β(δσ i ,σ j -1) i e hδ σ i ,1 , (1)
where the spins σ i take on the values of the set {1, . . . , q}, and where the first product is over nearest neighbour pairs (n,n). If we want to study the behaviour of clusters, in the sense of F-K clusters, we turn to the corresponding Edwards-Sokal formulation [5], given by the Boltzmann weight

ω ES (σ, η) = i,j [e -β δ η ij ,0 + (1 -e -β )δ η ij ,1 δ σ i ,σ j ] i e hδ σ i ,1 , (2)
where the edge variables η ij belong to {0, 1}. This "site-bond" model can be thought of as follows. Given a certain spin configuration, one puts between two neighbouring sites σ i = σ j an edge or bond with the probability 1e -β , and no edge with the probability e -β ; for σ i = σ j , no bond is present. When the field is infinite, all σ i = 1, and we are left with a classical bond percolation problem, while for finite field, one has a random bond percolation model in the random media given by the spin configuration.

In the presence of an external field, we find it convenient for the study of the Kertész's line to consider a modified version of the Edwards-Sokal formulation. We have three different types of spin combination: (0) two adjacent spins i, j are not equal, (1) two adjacent spins i, j are equal and parallel to h (we denote this direction as 1), or (2) two adjacent spins i, j are equal but not parallel to h. Correspondingly, we "color" the edge between i and j in three different colors n i,j , where the edge variables n ij belong to {0, 1, 2}. The resulting Boltzmann weight becomes

ω CES (σ, n) = i,j [e -β δ n ij ,0 + (1 -e -β )δ n ij ,1 χ (σ i =σ j =1) + (1 -e -β )δ n ij ,2 χ (σ i =σ j =1) ] i e hδ σ i ,1 , (3)
where the characteristic function χ(σ i = σ j = 1) is unity for σ i = σ j = 1 (parallel spins in the direction of h) and zero otherwise, while χ(σ i = σ j = 1) is unity for parallel spins not in the direction of h and zero otherwise. The summation over the spin variables then leads to the following Tricolor-Edge-Representation

ω TER (n) = i,j e -βδ n ij ,0 (1 -e -β ) (δ n ij ,1 +δ n ij ,2 ) × e hS 1 (n) (q -1) C 2 (n) (q -1 + e h ) |Λ|-S 1 (n)-S 2 (n) . (4)
Here, S 1 (n) and S 2 (n) denote the number of sites that belong to edges of color 1 and color 2, respectively, while C 2 (n) denotes the number of connected components of the set of edges of color 2, and |Λ| is the number of sites of the lattice under consideration.

Let us mention that such kind of graphical representation has already been considered for various spin models in presence of an external field [3].

Let p Λ (i ↔ j) be the probability that the site i is connected to j by a path of edges of color 1. As (geometric) order parameter we will consider the following mass-gap (inverse correlation length)

m(β, h) = -lim |i-j|→∞ 1 |i -j| ln lim Λ↑ d p Λ (i ↔ j) (5)
where i and j belong to some line parallel to an axis of the lattice. As (thermodynamic) order parameter, we shall consider the mean energy

E(β, h) = -1 β ∂ ∂β f (β, h)
, where f (β, h) is the free energy of the model4 .

Analytic results

Let us first have a look at the diagram of ground state configurations of the TER representation which are the translation invariant configurations maximizing the Boltzmann weight (4).

For the color a = 0, 1, 2, let b a be the value of the Boltzmann weight of the ground state configuration of color a per unit site. One finds b 0 = e -βd (q

-1 + e h ), b 1 = (1 -e -β ) d e h , b 2 = (1 -e -β ) d .
Notice that b 0 = b 1 on the line When q is large enough and h not too large, the TER representation (4) can be analyzed rigorously by a perturbative approach. Namely, by using the standard machinery of Pirogov-Sinai theory, we will show that the model undergoes a thermodynamic first order phase transition in the sense that the mean energy (as well as the magnetization) is discontinuous at some β c (h) ∼ β 0 (h). We also find for these values of the parameters, that the phase

β 0 (h) = ln[1 + (1 + (q -1)e -h ) 1/d ] (6
h β β 0 (h) 0 1 2
Figure 2: Diagram of ground state configurations: All the ground state configurations coexist at (0, β 0 (0)). Below β 0 (h), the 0-state dominates. Above β 0 (h), the 1-state dominates; it coexists with the 0-state on the line β 0 (h), and with the 2-state on the line h = 0, β ≥ β 0 (0). diagram of this model reproduces the diagram of ground state configurations (Fig. 2), see Appendix for more details.

In addition, the model exhibits a geometric (first order) transition, in the sense that, on the same critical line, the mass gap is discontinuous.

Theorem 1. Assume d ≥ 2, q and h such that

c d (1 + (q -1)e -h ) -1/2d < 1 (7)
holds, where c d is a given number (depending only on the dimension), then there exists a unique

β c (h) = β 0 (h) + O(1 + (q -1)e -h ) -1/2d ) such that 1. ∆E(β c (h), h) = E(β - c (h), h) -E(β + c (h), h) > 0 2. m(β, h) > 0 for β ≤ β c (h) and m(β, h) = 0 for β > β c (h).
The proof is given in the appendix.

Let us recall that it has already been shown that the Potts model (1) undergoes, for q large and h small, a first order phase transition on a critical line [1], where both the mean energy and the magnetization are discontinuous. Since, as already mentioned, the free energies of models ( 1) and ( 4) are the same, this critical line coincides with the one mentioned in the theorem.

In the absence of an external magnetic field, the statements of the theorem have been shown previously [13,14,12].

Condition (7) restricts the range of values of parameters to which our rigorous analysis applies. Moreover, we do not expect thermodynamic first order transitions when h is sufficiently enhanced. In the next section, we turn to numerical study on a wider range of values.

Numerical simulations

We have implemented a generalization of the Swendsen-Wang algorithm for our colored Edwards-Sokal model (3).

First, given a spin configuration, we put between any two neighbouring spins of the same color, an edge colored 0 with probability (w.p.) e -β , and w.p. 1e -β , an edge colored 1 if these spins are of color 1, and colored 2 otherwise. When two neighbouring spins disagree, the corresponding edge is colored 0.

Then, starting from an edge configuration, a spin configuration is constructed as follows. Isolated sites (endpoints of 0-bonds only) are colored 1 w.p. e h /(q -1 + e h ) and colored c ∈ {2, ..., q} w.p. 1/(q -1 + e h ). Nonisolated sites are colored 1 (w.p. 1) if they are endpoints of 1-bonds and colored c ∈ {2, ..., q} w.p. 1/(q -1).

This algorithm allows us two compute both quantities associated to spins configurations and those associated to edges configurations.

The numerical results for d = 2 are presented in Fig. 3. For q ≤ 4, we found a whole geometric transition line β c (h) for which m(β, h) > 0 when β < β c (h), and m(β, h) = 0 when β ≥ β c (h). The mass gap is continuous at β c (h). For β < β c (h), the mean cluster sizes remain finite, while for β ≥ β c (h) the size of 1-edge clusters diverges. The energy density as well as the magnetization do not show any singular behavior.

For q ≥ 5, some critical h c appears for which the geometric transition line β c (h) becomes first order when h < h c , i.e. m(β, h) > 0 for β ≤ β c (h) and m(β, h) = 0 for β > β c (h). In addition, on this part of the line, we find that the mean energy is discontinuous 5 . When h ≥ h c , only a geometric transition occurs and the scenario is the same as for q ≤ 4. Thus our numerics show that the geometric and thermodynamic transitions coincide up to h c , similarly to what we got analytically but only at (very) small field (and large q), see Fig. 3.

Let us mention that the numerics are in accordance with the theory for vanishing and infinite fields: β c (0) = ln(1 + √ q) and β c (∞) = ln 2. Figure 3: β c (h) for several values of q, with "first order" behavior in red, "second order" in blue. The first order behavior is both thermodynamic and geometric. The second order behavior is only geometric.

The system size in these calculations was L = 50, d = 2. The "first order" part of the transition lines has been determined via Binder cumulants [2]. The Hoshen-Kopelman algorithm [8] was used to study cluster statistics. For each value of q, more than 2 × 10 5 iterations were performed. Data have been binned in order to control errors in measurements.

Concluding remarks

For the Potts model in the presence of an external magnetic field, we have shown that when the Kertész line is first order, it coincides with the usual thermodynamic critical line. This property holds up to some critical point (h c , β c (h c )), beyond which the thermodynamic transition disappears. Such behavior may well appear also for a broader class of models exhibiting first order transition in the presence of an external field. We believe that the behavior at the above critical point also belongs to the universality class of the Ising model, as it is the case in the 3-state Potts model in three dimensions [10].

Acknowledgements

It is a pleasure to thank János Kertész for interesting comments. The BIBOS research center (Bielefeld) and the Centre de Physique Théorique (CNRS Marseille) are gratefully acknowledged for warm hospitality and financial support. The authors are indebted to an anonymous referee for constructive remarks.

Appendix

We first introduce the partition function of the TER representation with boundary conditions a ∈ {0, 1, 2} in a box Λ6 :

Z a (Λ) = n i∈Λ ω i (n)q C 2 (n)-δ a,2 i∈∂Λ j∼i δ n ij ,a (8)
where the sum is over all configurations n = {n ij } ij∩Λ∅ , ∂Λ is the boundary of Λ (set of sites of Λ with a n.n. in d \ Λ), the notation i ∼ j means that i and j are n.n., and

ω i (n) = (1 -e -β ) (δ n ij ,1 +δ n ij ,2 )/2 e -βδ n ij ,0 /2 e hχ(i∈"1 ′′ ) (q -1 + e h ) j∼i δ n ij ,0 (9)
where χ(i ∈ "1 ′′ ) means that the site i belongs to some edge of color 1. Next, consider a configurationn on the envelope of Λ:

E(Λ) = { i, j ∩ Λ = ∅}.
A site i ∈ Λ is called correct if for all j ∼ i, n ij takes the same value, and called incorrect otherwise. Denote I(n) the set of incorrect sites of the configuration n. A couple Γ = {Supp Γ, n(Γ)} where the support of Γ (Supp Γ) is a maximal connected subset of I(n), and n(Γ) the restriction of n to the envelope of Λ is called contour of the configuration n (here, a set of sites is called connected if the graph that joins all the sites of this set at distance

d(i, j) = max k=1,...,d |i k -j k | ≤ 1 is connected). A couple Γ = {Supp Γ, n(Γ)
Z a (Λ) = θ={Γ 1 ,...,Γn}ext b | Ext Λ θ| a n k=1 ρ(Γ k ) m = 0, 1, 2 Z m (Int m Γ k ), (10)
where ρ(Γ) = i∈supp Γ ω i (n Γ )q C(n Γ )-δ a,2 . From (10), we get

Z a (Λ) = b |Λ| a {Γ 1 ,...,Γn}comp n k=1 z a (Γ k ), (11)
where the sum is now over families of compatible contours and the activities z a (Γ) of contours are given by z a (Γ) = ρ(Γ)b

-| supp Γ| a m =a Zm(Int m Γ)
Za(Int m Γ) . It is easy to prove the following Peierls' estimate,

ρ(Γ)( max a=0,1,2 b a ) -| Supp Γ| ≤ e -τ | Supp Γ| , (12)
where e -τ = (1 + (q -1)e -h ) -1/2d . Indeed, first notice that an incorrect site i is either of color 1 or of color 2. In the first case one has j∼i (δ n ij ,0 + δ n ij ,1 ) = 2d, so that ω i (n Γ )/b 1 = (e β -1) -( j∼i δ n ij ,0 )/2 , implying

ω i (n Γ )/ max a=0,1,2 b a ≤ (1 + (q -1)e -h ) -( j∼i δ n ij ,0 )/2d .
Thus since 1 ≤ j∼i δ n ij ,0 ≤ 2d -1, each incorrect site of color 1 gives at most a contribution e -τ to the L.H.S. of (12). In the second case, one has

j∼i (δ n ij ,0 +δ n ij ,2 ) = 2d, so that w i (n Γ )/b 2 = (e β -1) -( j∼i δ n ij ,0 )/2 , implying ω i (n Γ )/ max a=0,1,2 b a ≤ (q -1 + e h ) -( j∼i δ n ij ,0 )/2d .
We then use again that 1 [12]) to obtain that each incorrect site of color 2 gives at most a contribution (e h + q -1) -1/2+1/2d ≤ e -τ to the L.H.S. of (12).

≤ j∼i δ n ij ,0 ≤ 2d -1 and that C 2 (n Γ ) ≤ i∈Supp Γ χ(1 ≤ δ n ij ,2 )/2 j∼i δ n ij ,2 (see
When the assumptions of the theorem are satisfied, the Peierls' estimate (12) provides a good control of the system by using Pirogov-Sinai theory [START_REF] Ya | Theory of Phase Transitions: Rigorous Results[END_REF]. We introduce the truncated activity

z ′ a (Γ) = z a (Γ) if z a (Γ) ≤ (c 0 e -τ ) | Supp Γ| (c 0 e -τ ) | Supp Γ| otherwise,
where c 0 is a numerical constant, and we call a contour stable if z a (Γ) = z ′ a (Γ). Let Z ′ a (Λ) be the partition function obtained from (11) by leaving out unstable contours, i.e., by taking the activities z ′ a (Γ) in (11), and let us introduce the metastable free energies f met a (β, h) =lim Λ↑ d (1/|Λ|) ln Z ′ a (Λ). The leading term of these metastable free energies equalsln b a . The corrections can be expressed by free energies of contour models which can be controlled by convergent cluster expansions. As a standard result of Pirogov-Sinai theory, one gets that the phase diagram of the system is a small perturbation of the diagram of ground state configurations. Namely, there exits a unique point β c (0) given by the solution of

f met 0 (β, h) = f met 1 (β, h) = f met 2 (β, h
) for which all contours are stable and such that Z a (Λ) = Z ′ a (Λ) for a = 0, 1, 2. There exists a line β c (h) given by the solution of f met 0 (β, h) = f met 1 (β, h) when h > 0 and such that, Z a (Λ) = Z ′ a (Λ) for a = 0, 1. For β < β c (h) one has Z 0 (Λ) = Z ′ 0 (Λ), and for β > β c (h) one has Z 1 (Λ) = Z ′ 1 (Λ). For h = 0 and β ≥ β c (0), one has in addition Z 2 (Λ) = Z ′ 2 (Λ). For the color a = 0, 1, 2, denote by • a (β, h) the expectation value under the a-boundary condition. As a consequence of the above expansions and analysis, we obtain by standard Peierls' estimates that for h ≥ 0

δ n ij ,1
1 (β, h) ≥ 1 -O(e -τ ) for β ≥ β c (h) (13) δ n ij ,0 0 (β, h) ≥ 1 -O(e -τ ) for β ≤ β c (h) (14) while in addition we also get for h = 0:

δ n ij ,2 2 (β, 0) ≥ 1 -O(e -τ ) for β ≥ β c (0) (15)
By definition of the mean energy, one has that ∆E = E(β -, h) -E(β + , h) is proportional to the difference δ n ij ,0 0 (β -, h)δ n ij ,0 1 (β, h), and the first statement of the theorem follows immediately from these properties.

To prove the second statement, we remark that if one imposes that the site i is connected to j by a path made up of edges of color 1, then under the boundary condition 0, there exists necessarily an external contour that encloses both the sites i and j. As a consequence of the above analysis the probability of external contours Γ decays like (c 0 e) -τ | Supp Γ| when the 0-contours are stable, i.e. when Z 0 (Λ) = Z ′ 0 (Λ). One thus gets p Λ (i ↔ j) ≤ (Cte e -τ ) |i-j| when β ≤ β c (h) from which the first statement of the theorem follows. On the other hand under the boundary condition 1, the probability that the site i is not connected to j can be bounded from above by a small number O(e -τ ) when Z 1 (Λ) = Z ′ 1 (Λ). This follows also from a Peierls type arguments and implies that the probability that the site i is connected to j under the boundary condition 1 is greater than 1 -0(e -τ ) for β ≥ β c (h). It gives also that the probability p Λ (i ↔ j) for the site i to be connected with j under the boundary condition 0 is also greater than 1 -0(e -τ ) for β > β c (h), implying the second statement.

  ) and that b 0 = b 1 = b 2 at the point β 0 (0). The diagram of ground state configurations, inferred from the values of the weights b 0 , b 1 , b 2 is shown in Fig. 2 (in the (h, β) plane).

  }where Supp Γ is a connected set of sites is called contour if there exists a configuration n such that Γ is a contour of n. For a contour Γ, let n Γ denote the configuration having Γ as unique contour, Ext Γ denotes the unique infinite component of d \ (Supp Γ), Int Γ = d \ (Ext Γ ∪ Supp Γ), and Int m Γ denote the set of sites of Int Γ corresponding to the color m ∈ {0, 1, 2} for the configuration n Γ . Two contours Γ 1 and Γ 2 are said to be compatible if their union is not connected and are called external contours if furthermoreInt Γ 1 ⊂ Ext Λ Γ 2 and Int Γ 2 ⊂ Ext Λ Γ 1 . For a family θ = {Γ 1 , . . . , Γ n } ext of external contours, let Ext Λ θ denote the intersection Λ ∩ n k=1 Ext Λ Γ k .With these definitions and notations, one gets the following expansion of the partition functions over families of external contours,

Note that all partitions functions, and hence the free energies, of models (1 )-(4) coincide

The Swendsen-Wang algorithm allows to compute both associated order parameters (mass-gap and mean energy).

The reader should not be confused by the fact that (8), called diluted partition function in PS-theory, differs from the usual one by an unimportant boundary term which makes the expansions(10) and(11) easier to write.