Francois Jeanneret

Frank Kirchner

Alain Clappier

Hubert Van Den Bergh

Bertrand Calpini

Total VOC reactivity in the planetary boundary layer 1. Estimation by a pump and probe OH experiment

The reactivity of hydroxyl radical (OH) with total volatile organic compounds (VOC), Z kvoc[VOC ], is directly obtained in the planetary boundary layer (PBL). The method is based on a perturbation technique where high initial OH concentrations are created by flash photolysis of ozone and subsequent reaction of O(•D) with H20. Laserinduced fluorescence is used to measure the residual OH concentration at different time delays after the perturbation (<100 ms) and obtain a direct estimate of the OH lifetime rOH for different atmospheric conditions. For specific experimental conditions the transport by the wind may be neglected, and the chemical processes governing the OH decay may be expressed by a detailed box model. With a simple chemical equation derived in this paper using roll and complementary measurements of CO, 03, and NOx, an in situ estimate of Z kvoc[VOC ] in the PBL is retrieved with an uncertainty of less than 20% in comparison to the detailed box model calculations. This analysis is applied to laboratory measurements with three synthetic NOx/VOC/O 3 gas mixtures, and the retrieved OH lifetimes and total VOC reactivity are discussed against model predictions.

1.

Introduction

Tropospheric ozone is a criteria pollutant formed by photooxidation of VOC in the presence of nitrogen oxides (NOx -NO + NO2) and solar radiation. In the planetary boundary layer (PBL), ozone is a major constituent of each photochemical episode [START_REF] Staffelbach | Photochemical oxidant formation over southern Switzerland, 1, Results from summer 1994[END_REF] and can be harmful to both humans and vegetation [START_REF] Calpini | Air pollution: Measuring techniques and impact on our environment[END_REF]Neffel, 1999;[START_REF] Fuhrer | Effects of air pollutants on plants[END_REF]. The oxidation sequence leading to ozone production is initiated by VOC reaction with the hydroxyl radical (OH). The relative behavior of VOC and NOx in ozone formation can be understood in terms of competition for OH. The latter participates in both the NOx and VOC reaction cycles, thus playing a central role in the atmospheric chemistry.

Here we summarize the most important steps in this atmospheric chemistry, while the reader is referred to Kirchner et al. thus inducing the ozone production by (R4) and (R5).

At high ratios of NOx to VOC concentration, OH will react mainly with NOx, and at a low ratio, the VOC reactions can predominate. However, the VOC and NOx lifetimes against the loss by reaction with OH will depend on the VOC-OH and NOx-OH rate constants. The OH rate constants with different VOCs differ for each individual VOC species. Hundreds of different VOC species are present in the atmosphere. They may be measured independently but with VOC detection techniques that are slow and request in particular a long period of integration time.

This paper presents ideas on what we can gain with the pump and probe OH method, a method where on the contrary high-frequency measurements are possible, designed to measure in situ the OH lifetime and consequently its reactivity with the sum of all VOC present in an air parcel (5• koH,voc[VOC]).

The latter parameter is shown to be specific to each ambient condition. Total VOC reactivity will give more significant indications on the ozone production chemistry than the VOC respective concentrations [Kirchner et al., this issue]. In section 1 a summary of the experiment is presented. Section 2 underlines the effect of the transport of the pumped OH out of the probed volume by the wind. Section 3 shows how the complete atmospheric chemical mechanism can be simplified if one considers our special initial condition of a pumped OH concentration of up to 5 orders of magnitude higher than its natural concentration and analyze the OH lifetime on a timescale shorter than 20 ms. Section 4 compares the total VOC reactivity with OH predicted by box model with experimental results for three different NOx/VOC gas mixtures.

Pump and Probe OH Technique

Because of its high reactivity, the hydroxyl "steady state" concentration remains much below one part per trillion depending on the local chemical composition, making its detection particularly demanding. The probe beam is generated by a frequency doubled pulsed dye laser pumped by a XeC1 excimer laser. This tunable laser source is aligned on one of the most intense absorption bands of OH around 282 nm (Q 1(2) transition). The fine wavelength tuning on the selected OH transition was performed initially using an intense source of OH from a methane/oxygen gas burner. The LIF signal induced by the presence of OH is directly expressed in terms of OH concentration, with an ac-curacy of better than 30%. This overall uncertainty is the sum of the uncertainties due to the different optical parameters influencing the LIF signal (in particular, the OH absorption cross section and the quantum yield of OH fluorescence), the geometrical factors such as volume sampled for both pump and probe beams, the optical detection efficiency, and the statistical variability of the fluorescence signal. According to an experiment where a known amount of OH was produced and its absolute concentration was retrieved by this method with a relative difference of less than 6%, one may consider this 30% uncertainty as an upper limit [START_REF] Calpini | Direct measurement of the total reaction rate of OH in the atmosphere[END_REF].

Note that the probe wavelength itself with an energy of 0.6 mJ/pulse, a quantum yield of 0.95 [Finlayson-Pitts and Pitts, 1986], and an ozone absorption cross section 0-03(282 nm) = 325 x 10 -2ø cm 2 [START_REF] Molina | Absolute absorption cross section of ozone in the 185-to 350-nm wavelength range[END_REF]] generates a "secondary" OH formation. For typical conditions this effect reaches a maximum secondary OH production of less than 5 % with respect to [OH]o. This effect is constant in time and appeared not to perturb the OH relaxation analysis at the applied conditions. A pulse box generator triggers the two lasers independently and allows controlling the time delay between the pump and the probe laser pulses.

The OH LIF is collected perpendicularly to the plane given by the gas flow and the two coaxial laser beams. The OH LIF signal is focused by a combination of two UV fused silica lenses with f/3.5 numerical aperture at the entrance of a double grating monochromator, and detected by a side-on photomultiplier tube.

Wind Effects

The time evolution of OH after the pump step is controlled by both the physical and chemical properties of the atmosphere. The atmospheric dispersion induced by diffusion and transport due to the wind blows the pumped OH out of the probed volume. This directly affects our measure of the OH time decay. As the aim of this study is to measure the atmospheric OH relaxation time in a "purely chemical" regime, the geometrical configuration of the experiment must be adapted so to avoid any wind effects, at least for a given time delay after OH was formed.

Initially, the pump laser beam creates OH in the volume defined by its direction and its cross-sectional area. The OH concentration is measured in a smaller volume delimited by the probe beam geometry as seen in Figure 1 When the overlap is complete (At < Ate), the transport has no effect on the OH time decay, and term A in (1) can be neglected. The measured OH relaxation is in this case exclusively due to the chemical mechanisms and the diffusion. On the other hand, when At > At2, there is no overlap, and the perturbed volume is totally out of the LIF zone in which case the air mixture in the probed volume is renewed. This condition will determine the laser frequency f < 1/At of the pump and probe laser beam in order to avoid any memory effect in the pumped OH volume. Thus the beam radii and the experimental repetition rate are defined by At < Ate: the chemical OH decay is not biased by transport effect; and a laser frequencyf = 1/At where At > At 2 to avoid any memory effect.

For calm atmospheric condition (U = 2 m/s) a typical experiment will be performed with the pump and the probe radii at the experimental site set at 85 and 5 mm, respectively. This would lead to a time window At• -40 ms where the transport would not influence the OH relaxation. In this case, no memory effect would be observed for At 2 > 45 ms or an experimental repetition rate typically lower than 20 Hz.

With At < Ate, term A in (1) is neglected: with the time evolution of any atmospheric species depending only on the chemical reactions. With these assumptions the pump and probe method may be used to measure the OH lifetime in a "pure chemical regime" or without transport effect. However, this also underlines that the pump and probe method is limited to the analysis of the OH time decay for short time delay only, typically below 100 ms. Section 3 will demonstrate that even under such conditions, large variations in the OH relaxation time are obtained for different atmospheric conditions and as a result the total VOC reactivity with OH can be retrieved.

OCk/Ot = V(KVCk) + &•. ( 4

OH Chemical Reactions in the Atmosphere

The atmospheric photochemistry is a complex mechanism including homogeneous and heterogeneous chemical reactions. This study will be restricted to to obtain realistic mixing ratios for peroxy radicals and a real- 2 on a time window of 100 ms. In this simulation the maximum initial OH concentration is about 3 x 10 •ø molecules/cm 3 for both the rural and moderately polluted cases since the ozone concentration was in both cases around 80 ppb as input before pump (see Table 1). In the two other cases the initial OH concentration is lower due to respectively lower initial ozone concentrations. 1. One may consider as negligible any OH production term. This assumption is shown not to be valid any more in highly polluted conditions where OH in the presence of high nitric oxide may generate a secondary OH production term not constant in time. case, the respective concentration of isoprene, d-limonene, a-pinene, or formaldehydes and higher aldehydes (input before the pump, see Table 1) were predicted to change between t = 0 and t = 20 ms in the following way: isoprene from 8.96 ppb to 8.76 ppb, d-limonene from 5.59 ppb to 5.53 ppb, and a-pinene from 2.14 ppb to 2.06 ppb, while the other VOC species were even more stable. With the different reaction rate used in RACM, these relatively small changes in the VOC concentration did not affect the analysis of the OH lifetime, and the assumption of constant VOC concentrations was confirmed.

In the polluted condition a fast relaxation is predicted during

One may use a constant value

In Figure 4 on short timescales can be much simplified. A pseudo firstorder reaction for OH has then been confirmed as a good estimate of the complete chemical code under the particular condition of this pump and probe experiment. This formalism can be checked against experimental results obtained in a controlled environment using a flow tube reaction chamber and different types of gas mixtures with synthetic air conditions.

Comparison Between Predicted and Measured

E kvo c [VOC] in Different VOC/NO x Environments

In this pump and probe experiment a gas mixture is continuously flown through the reaction cell. This synthetic air is formed of NO x, VOCs, and 03 added to the mr:O2:H20 gas stream at atmospheric pressure and temperature. VOCs are composed of propene, n-butane, and ethene in a ratio of 1:8:3 for total VOC concentrations varying between 125 and 210 ppb in accordance with atmospheric conditions reported by Chameides et al. [1992]. NO x and 03 concentrations are controlled so that three stable gas mixtures are tested with values following Finlayson-Pitts and Pitts [1986] for moderately polluted, ruralN%, and ruralvo c conditions, respectively. The experimental concentrations are listed in Table 2. We used the term ruralN% for an air parcel in rural conditions characterized by the presence of high NO x concentrations or similar to an air pollution condition downwind from a large urban area, and ruralvo c for a rural condition with low NO x concentrations. Hence, for NO x = NO + NO 2 concentrations, the experimental NO2/NO does not exhibit a realistic ratio. This is due to the dilution bench where 03 and NO x are mixed together before the reaction chamber in order to assure a homogeneous gas mixture at the reaction site. As ozone reacts very fast with nitric oxide, a large fraction of NOx at the reaction site is NO 2. For this reason, 03 and NO x are measured directly at the reaction site by UV differential absorption and chemiluminescence technique, respectively. The three VOC species are obtained directly from a gas calibration tank. The uncertainty on the VOC concentrations is estimated below 10%, while NO x and 0 3 gas detectors have both a precision of _+ 1 ppb.

The measurements were performed by time steps of 1 ms for a time delay varying from 1 to 20 ms. 3.

I i i i i i i i i i i i i i i i i i i i i I i i ! i i i ! ! i
The agreement between experimental and simulated OH lifetimes is very good for the moderately polluted case. For the ruralNo x case the experimental value is 20% higher than the value calculated by the model and 12% higher than the value estimated by applying (17). Considering the error limits, there is an agreement between experiment and model. For the ruralvo c there is a large difference between Tex p on the one hand and Tre f and 7RACM on the other hand. However, one has to take into account also the large error limits especially of Tre f. This uncertainty in Tre f results only from uncertainties in the kinetic data. They affect rRACM in a similar manner so that even for the ruralvo c case the OH lifetimes agree inside their error limits. We may now use ( 17 One should also mention that the actual experiment itself gives better results in highly polluted conditions where very short OH lifetimes are expected. This is partially linked to the limit of detection of the apparatus giving a lowest OH concentration detection of some 108 OH molecules/cm 3. In order to have a sufficient dynamic range for each of the three gas stream conditions, the experiments were performed with [OH]o slightly above the proposed domain of validity of (18), namely, with initial OH concentrations up to 3.75 X 10 TM OH molecules/cm 3 for the ruralvo c case. In this case the linearization errors defined in section 3 are larger with the most important effect for the longer OH lifetimes. In addition, the time window of 20 ms was set so to fulfill the condition of a "pure chemical regime" (no transport effect due to the gas flow). However, this time window appeared to be short especially for the detection of OH lifetime above typically 30 ms as is the case for the ruralvo c case.

The uncertainties defined in

Finally, it is important to underline that the experiments were performed with different NO2/NO ratios higher than what is the case in the real atmospheric environment. Reaction (R3) (HO 2 + NO --> OH + NO2) being an important source of secondary OH formation and the experiments being performed with low NO concentration (due to its rapid reaction with ozone in the dilution bench and conversion to NO2) , no secondary OH was perturbing rexp in this case. Nevertheless, this secondary OH formation could eventually have had a strong impact on the OH lifetime if the experiments were performed under more realistic NO2/NO ratios. In order to address this problem, box model runs with RACM were repeated for the moderately polluted case keeping any conditions similar to the experiment, except the NO2/NO ratio changed to 350 ppb/150 ppb. However, even in this case, no significant secondary OH formation was predicted during the first 20 ms.

Conclusion

With an in situ "pumped" OH concentration obtained by 03 flash photolysis, the OH time relaxation can be used to obtain a first estimate of the total volatile organic compound reactivity with OH in the atmosphere. This OH time decay is approximated for short time delays by a simplified chemical system. On the one hand, the interference due to transport by the wind and diffusion is avoided adapting the repetition rate and beam geometry of the experiment. On the other hand, the only significant chemical processes governing the OH loss are the • kVoc[VOC]RACM obtained from the complete RACM mechanism. For more rural-like experimental conditions, while the different results were still in accordance (within the error limits), both the chemistry with secondary OH interference and the apparatus with its actual OH detection limit were limiting the final accuracy of the results.

The upgrade of the experimental layout will consist first of using a Kr excimer laser for the ozone photodissociation. A higher ozone absorption efficiency will be gained at 248 nm, together with a higher energy per pulse (600 mJ) and a wider beam size (typically 10 x 25 mm without beam expander). The latter point is of importance since it will assure a "pure chemical regime" (no transport effect) over longer timescales. Furthermore, the OH detection will be performed with a gated photodetector and a LIF laser source tuned around 308 nm. This will essentially prevent the experiment from any secondary OH formation due to the probe laser beam itself. With these changes an OH detection limit of the order of 106 to 107 OH molecules/cm 3 is expected thus giving the chance for pump and probe OH measurements performed at much lower initial OH concentration (typically below 10 •ø OH molecules/cm3).

The pump and probe technique will also be applied as a remote sensing technique (lidar configuration).

Z kvoc[VOC] is the reactivity-weighted sum of all the different volatile organic compounds and shall give an important new criterion for our understanding of the control of the ozone formation in a photochemical smog. A companion paper by with these new experimental parameters measured in atmospheric conditions.

Figure 1 .

 1 Figure 1. Pump and probe laser beams configuration: wind effects at different time delay to < t l < t2 between the pump and the probe laser beams. The pump volume is moved away from the probed LIF zone: (a) full overlap, (b) partial overlap, and (c) no overlap.

  See, for example, collection of papers on the Tropospheric OH Photochemistry Experiment, Fritz Peak/Idaho Hill, Colorado, fall 1993 in Journal of Geophysical Research, 102(D5), 6169-6510, 1997. [see also Brune et al., 1998; Ehhalt et al., 1998; Mauldin et al., 1999]. Here we use a pump and probe method that makes it possible to directly measure the OH kinetics in different atmospheric conditions from remote to polluted urban air masses. It consists of creating a high initial OH concentration by flash photolysis of ozone in the Hartley band to form the excited oxygen atom O(1D), the latter reacting with water vapor to form OH. In this way, high OH concentration is initially created in the volume defined by the pump laser beam geometry. A second laser source is used for laser-induced fluorescence (LIF) detection of OH (probe beam). Thus the OH concentration is monitored at different times after the initial ozone flash photolysis, and its relaxation kinetics is followed in different environments. The pump and probe experiment was described in detail by Calpini et al. [1999], where it has been applied first in a reaction cell containing synthetic air at atmospheric pressure for the study of the reaction rate k(OH+CH4) and then in atmospheric conditions. The pumped OH concentration is obtained using a 4th harmonic Nd:YAG laser pulse at 266 nm with an energy of up to 80 mJ/pulse and a repetition rate of 10 Hz. ozone forms O(1D). The latter reacts with water vapor to form OH, but is also deactivated to O(3p) due to collisions with other air molecules or water vapor itself. Depending on the experimental conditions and the ozone and water vapor concentrations, the pumped OH concentration after the ozone flash photolysis, denoted below as [OH]o at t = 0, reaches typical values of up to 10 TM molecules/cm 3.

  Figure 1, and to accurately control their alignment on the same axis as well as the beam divergence of the two beams. The wind transports the pumped OH volume out of the probe volume. The two beam radii rpump and rprob e must assure a time window At, during which the experiment is not biased by the wind. During At the wind also fluctuates because of turbulent effects. Thus we have to consider Umax and Umin to be the maximum and minimum wind speed during the OH relaxation. The experimental conditions can be defined in three steps: full overlap (Figure la) of the probe beam in the pumped area if At < At• with •tl = If'pump-rprobel/$ .... (2) no overlap (Figure lc) if At > At 2 with At2 = Irpump + rprobel/Smin, (3) partial overlap (Figure lb) if At• < At < At 2. When the overlap is complete (At < Ate), the transport has no effect on the OH time decay, and term A in (1) can be

  check for the respective influence of the chemistry and the diffusion effects on the OH lifetime, (6) is expressed on a cylindrical coordinate grid with the diffusion terms con-sidered as production and loss terms for the chemical solver. For time delays below 100 ms a sensitivity analysis,

  the gas phase chemistry only. The chemical mechanism used in this study is the Regional Atmospheric Chemistry Mechanism (RACM) [Stockwell et al., 1997], a mechanism used to simulate the atmospheric chemical composition in the troposphere from the Earth's surface through the upper troposphere for remote to polluted urban conditions. The RACM mechanism is a further development of the RADM2 [Stockwell et al., 1990]. This lumped mechanism which was intensively tested against smog chamber experiments contains 77 chemical species among which 28 are radicals, with 234 reactions. In our case, the reactions including OH + OH were added to the original mechanism. Their importance will be discussed later. Because of the fast OH chemistry, the numerical runs predicting the OH time decay are operated on time steps of 0.1 ms/step. The chemical solver is based on the method of Gong and Cho [1993]. It is a hybrid implicit/explicit scheme, with an explicit scheme used for species with relatively long lifetimes and the implicit Newton-Raphson iterative scheme used for the case was taken from Fehsenfeld et al. [1992]. In the remote case the VOC split resulting from Kuhn et al. [1998] was modified by increasing the percentages of biogenically emitted VOC, in particular, isoprene and terpenes. The mixing ratios are summarized in Table 1 as "initial" values for each trace gas or family and for the four chemical regimes. The simulation was performed in two steps. The initial values of

Figure 2 .Figure 3 .

 23 Figure 2. Natural logarithm of the predicted OH concentration decay versus time for the four atmospheric conditionsdefined in Table 1.

  Figure 3 and associated discussion about the linearization error), (16) would also remain valid as long as the implicit assumption of constant concentration [Ci]o versus time would be verified. This assumption will be particularly sensible in the remote case where high amount of very reactive VOC species are present and may be depleted during the course of a measurement. On the basis of the box model run of the remote

  Figure 4. Species reactivity with OH at t = 0 for the polluted, rural, and remote atmospheric conditions defined in Table1. Note the change in the vertical scale and the expanded scale by a factor of 1000, 500, and 500, respectively, for the less reactive species. The sum of these low-reactive species with OH contributes to less than 0.6% with respect to the overall reactivity with OH.

  For each time delay the LIF signal is averaged over typically 1000 laser shots acquired on an oscilloscope. Each data point is converted to an [OH] equivalent concentration value taking into account the overall LIF detection efficiency. The pump and the probe beam diameters were set at 6 and 3 mm, respectively, and the flow speed of the synthetic air stream was 5.2 cm/s at the reaction site. The conditions (2) and (3) are fulfilled for time delays At• = 29 ms and At 2 = 86 ms. Therefore the measured OH relaxation for time delays below At• is obtained without per-

Figure 6 .

 6 Figure 6. Natural logarithm of the pump and probe OH concentration versus time measured in the reaction cell for the three different NOx/VOC/O3 gas mixtures defined in Table 2: (a) experimental values and (b) RACM-predicted results. OH lifetimes r are given by the slope of the linear fits. The experimental error bars represent the overall uncertainty of the statistical variability of the fluorescence signal and variations in the trace gas concentrations, assumed to be random.

  ) and get directly the total VOC reactivity with OH: E kvoc[gOr]exp = 1/•'OH-E kNox[NOx]oof the error sources on • kvoc[gOC]exp are due to • kNox[NOx]o, while the pump and probe experimental determination of 1/roll and 2koH[OH]o contributes to less than 20% of the error. In other words, the largest contribution to the error reported on • kvoc[VOC]exp is in this case due to uncertainties that are not related to the pump and probe method. This figure changes in the ruralvo c case where the term 2koH[OH]o plays the major role on the • kvoc[gOC]exp uncertainty. The total VOC reactivity • kvoc[VOC]ref is the reference value obtained as the sum of the propene, ncan see that for all three cases the error resulting from the simplification (difference between • kvoc[VOC]ref and • kVoc[VOC]RACM) is considerably smaller than the uncertainty in the kinetic data (error limits of • kvoc[VOC]ref ). Estimating the uncertainty of • kVoc[VOC]RACM being equal to the uncertainty of • kvoc[VOC]ref , there is an agreement with the experimental values inside the error limits. • kvoc[VOC]exp is calculated by (18) with error limits that become larger if rexp has larger uncertainties and also if the proportion of the OH loss by reaction with VOC gets smaller compared to the overall OH loss. In particular, in the case of high NOx concentrations the contribution of uncertainties in the kinetic data needed in (18) to the error limits of • kvoc[VOC]exp is high, for example, more than 50% in the ruralNox case.

  radical reactions with VOCs, NOx, 03, CO, and OH itself. Therefore a first in situ estimate of the total VOC reactivity with OH Z kvoc[VOC] can be obtained measuring the pumped OH concentration [OH]0 at t = 0 and the inverse OH lifetime with the pump and probe experiment, with additional measurements of NOx, CO, and 03 concentrations by standard trace gas detectors. Under these conditions a simple scheme was defined that gives E kvoc[VOC] with an uncertainty smaller than 20% with respect to the RACM chemical mechanism for an OH decay time below 12 ms in a highly polluted case, while in less polluted conditions a better accuracy is obtained even for longer delays. On the basis of this analysis a series of pump and probe OH measurements were performed in a controlled environment with gas mixtures that were comparable to three different levels of pollution. In the case of a polluted environment the pump and probe OH experiment gave an estimate of the experimental Z kvoc[VOC]exp in good agreement with the referenced value Z kvoc[VOC]ref calculated directly as the sum of the VOC species in the gas stream times their respective rate constant with OH. The associated uncertainty was in large part due to the poorly known rate constant of NOx with OH rather than the method itself. This result was confirmed by

Table 1 were used to initiate the box model for 20 min in order

 1

Table 1 . Initial Conditions for Four Typical Atmospheric Cases a

 1

		Polluted	Moderately		Rural	Remote
			Input Before		Input Before		Input Before		Input Before
	Species	Initial	Pump	Initial	Pump	Initial	Pump	Initial	Pump
	Inorganics
	03	50	50.8	80	84.7	72	81.3	50	48.6
	H202	2	2.0	2	2.0	2	3.0	2	2.0
	CO	1000	1000	200	200	200	200	200	200
	H20	1 X 107	1 X 107	1 X 107	1 X 107	1 X 107	1 X 107	X 107	1 X 107
	NOx	500	493	5.0	4.1	6.0	4.7	0.10	0.036
	NO	150	142	1	0.7	1	0.7	'"	0.002
	NO 2	350	351	4	3.4	5	4.0	0.10	0.034
	Alkanes
	CH 4	1700	1700	1700	1700	1700	1700		1700
	NMVOC	1000	996.5	100.0	98.59	102.6	105.0	30.00	28.64
	ETH (ethane)	53	53	5.3	5.3	......		0.53	0.53
	HC3 (low-reactive alkanes)	429	428	42.9	42.2	41.7	41.6	4.29	4.28
	HC5 (middle-reactive alkanes)	88	88	8.8	8.5	2.2	2.2	0.88	0.88
	HC8 (high-reactive alkanes)	25	25	2.5	2.3	0.1	0.1	0.25	0.25
	Alkenes
	ETE (ethene)	30	30	3.0	2.8	......		0.3	0.3
	OLT (terminal olefines)	2	1.97	0.2	0.16	..-	2.0	0.02	0.49
	OLI (internal olefines)	1	0.67	0.1	0.03	......		0.01	0.03
	ISO (isoprene)	............				38.5	32.0	10.00	8.96
	API (a-pinene)	............				3.4	2.5	6.67	5.59
	LIM (d-limonene)							3.33	2.14
	Aromatics
	TOL (toluene)	56	56	5.6	5.4	0.5	0.5	0.56	0.56
	XYL (xylene)	8	8	0.8	0.7	0.2	0.2	0.08	0.08
	Carbonyls
	HCHO	60	56	6.0	5.9	1.0	5.0	0.6	0.9
	ALD (higher aldehydes)	50	51	5.0	5.1	6.0	6.8	0.5	1.2
							2.7	...	0.8
	KET (ketones)	198	199	19.8	20.2	9.0	9.6	1.98	2.45
	Radicals
	peroxy radicals	.-.	5.5 x 10 -4	"'	4.0 x 10 -2	-..	9.4 x 10 -2	...	1.2
	HO 2	-..	3.0 x 10 -4	"'	1.8 x 10 -2	...	2.7 x 10 -2	..-	1.9 x 10 -3
	RO 2	...	2.5 x 10 -4	...	2.2 x 10 -2	...	6.7 x 10 -2	...	1.2

.................. Methacrolein ...............

amll concentrations in ppb. istic ratio between NO and NO 2. The result (input before the pump) was then used as initial conditions for the chemistry at the time when OH is pumped in the atmosphere. This is the first time step of the pump and probe simulation with this model. The polluted atmosphere is characterized by high con- centrations of pollutant precursors of nonmethane volatile or- ganic compounds (NMVOC) and NOx. The remote case, a fairly clean atmosphere, is composed of low precursor concen- trations except for some typical OH-reactive remote hydrocar- bons such as isoprene. The moderately polluted and rural conditions are cases that are similar with moderate NOx and high biogenic VOC conditions but with different VOC split. The species acronyms of Table 1 are defined by Stockwell et al. [1997]. The natural logarithms of the absolute OH concentration time decays predicted by RACM are presented in Figure

Table 2 . Gas Concentration in the Experimental Gas Mixture a

 2

	Moderately
	Compounds	Polluted	RuralNox	Ruralvoc
	C3H 6	17	10	11.6
	C4H10	143.3	84.3	97.3
	C2H 4	51.2	30.1	34.8
	NO	5	4	1
	NO 2	148	86	2
	03	71	80	120
	amll values in ppb.

Table 3 did not include errors resulting from the simplifications made for obtaining (17). To evaluate them, Tre f is directly compared with rRACM using the modified RACM mechanism where (OH + OH) and the up- dated rate constant for reaction (HO + NO2) were added. The results Tref/TRACM are presented in Table 3. The OH loss rate calculated with the modified RACM mechanism is 3-7% smaller than the results applying (17). Besides the formation of

Table 3 . OH Lifetimes and VOC Reactivity for the Three Experimental Conditions

 3

		Moderately Polluted	RuralNox	Ruralvoc
	?exp, ms	12 _+ 1.8 (R = 0.97)	16.3 _+ 3.2 (R = 0.96)	25.3 _+ 6.2 (R = 0.97)
	q'ref, ms	11.5 ___ 3.1	18.3 _+ 4.9	36.9 _+ 9.9
	TRACM , ms	11.9	19.6	39.8
	q'ref/q'RACM , %	96.9	93.3	92.8
	• kvoc[VOC]exp , S -1 • kvoc[VOC]ref , s -1	26.9 _+ 6.9 30.3 +_ 8.2	24.5 _+ 12.5 17.8 _+ 4.8	33.1 _+ 13.5 20.6 _+ 5.6
	• kvoc[VOC]R•CM, s -1	27.6	14.2	18.7

OH from secondary reactions, the other reason for this differ- ence is the change in the concentration of the reactants during the experiment and in particular the OH concentration de- crease itself. Replacing the OH initial concentration in (17) by the "mean" OH concentration (i.e., after 10 ms) lowers this overestimation of the OH loss using (17) to 1.5-3%. Never- theless, in all three cases the errors caused by the simplification remained below 10%.

Acknowledgments. The authors are grateful to the Swiss National Science Foundation (grant 21-50861.97) and the CTI SATURN project (3260.1) for supporting this work. The authors would like to thank each of the referees for their helpful comments and suggestions. The authors are also grateful to Ed Browell and his team at NASA Langley for stimulating discussions.


Christophe Jacquet

Akash Sharma

Matthieu Ehlinger

Jean-Noël Argenson

Sebastien Parratte

Matthieu Ollivier
email: matthieu.ollivier@ap-hm.fr

Maxime Fabre

Patient-specific high-tibial osteotomy's 'cutting-guides' decrease operating time and the number of fluoroscopic images taken after a Brief Learning Curve

Keywords: High-tibial osteotomy, Patient-specific instruments, Accuracy, Correction, Fluoroscopy, Operative time

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

In patients with moderate tibiofemoral osteoarthritis and varus alignment of the leg, an opening-wedge high-tibial osteotomy (HTO) is an ideal joint sparing operative solution [START_REF] Asik | High tibial osteotomy with Puddu plate for the treatment of varus gonarthrosis[END_REF][START_REF] Hantes | Satisfactory functional and radiological outcomes can be expected in young patients under 45 years old after open wedge high tibial osteotomy in a long-term followup[END_REF][START_REF] Schallberger | High tibial valgus osteotomy in unicompartmental medial osteoarthritis of the knee: a retrospective follow-up study over 13-21 years[END_REF]. The angular correction must however be optimal in all three spatial planes for the osteotomy to be effective. In general we have seen Computer-assisted surgery being used intra-operatively to control the overall frontal alignment of the lower limbs axis, however, it has not been able to demonstrate an ability to control the degrees of correction required within the osteotomy itself. In fact, the pin trackers of the navigation system cannot be implanted in each segment of the tibia to control the exact opening and therefore produce the exact correction required in the frontal and sagittal planes [START_REF] Schröter | Surgical accuracy in high tibial osteotomy: coronal equivalence of computer navigation and gap measurement[END_REF]. Furthermore, the navigation system requires a significant capital investment and the system might be cumbersome to use in the operative room [START_REF] Hankemeier | Navigated open-wedge high tibial osteotomy: advantages and disadvantages compared to the conventional technique in a cadaver study[END_REF][START_REF] Saragaglia | Navigated osteotomies around the knee in 170 patients with osteoarthritis secondary to genu varum[END_REF]. The recent introduction of patient-specific cutting guides (PSCGs) based on pre-operative CT-scan templating now offers the possibility of providing instrumentation specific guides to each patient allowing the accurate control of the correction during the procedure for both the frontal and the sagittal planes [START_REF] Arnal-Burró | Three dimensional-printed patient-specific cutting guides for femoral varization osteotomy: do it yourself[END_REF][START_REF] Munier | Can three-dimensional patient-specific cutting guides be used to achieve optimal correction for high tibial osteotomy? Pilot study[END_REF]. The accuracy of these systems has been demonstrated in different in vivo and in vitro studies, with superior correction accuracy when compared to conventional techniques [START_REF] Donnez | Are three-dimensional patientspecific cutting guides for open wedge high tibial osteotomy accurate? An in vitro study[END_REF][START_REF] Jacquet | More accurate correction using "patientspecific" cutting guides in opening wedge distal femur varization osteotomies[END_REF][START_REF] Kroes | Numerical optimization of alignment reproducibility for customizable surgical guides[END_REF]. These studies using PSCGs have also shown a decrease in operating time and a reduction in the number of fluoroscopic images per intervention compared with conventional techniques [START_REF] Arnal-Burró | Three dimensional-printed patient-specific cutting guides for femoral varization osteotomy: do it yourself[END_REF][START_REF] Corona | Preliminary results using patient-specific 3d printed models to improve preoperative planning for correction of post-traumatic tibial deformities with circular frames[END_REF][START_REF] Munier | Can three-dimensional patient-specific cutting guides be used to achieve optimal correction for high tibial osteotomy? Pilot study[END_REF]. As for each and every new technological advance that emerges, this system may have a learning curve, but this has not been evaluated in the literature [START_REF] Kayani | Robotic-arm assisted total knee arthroplasty has a learning curve of seven cases for integration into the surgical workflow but no learning curve effect for accuracy of implant positioning[END_REF][START_REF] Sodhi | The learning curve associated with robotic total knee arthroplasty[END_REF]. From our early experience with systems, it was our hypothesis that PSCGs for opening-wedge high-tibial osteotomies (HTO) would have a short learning curve with regards to the operating time and the use of fluoroscopy without any influence on the accuracy of correction. The objectives of this study were to determine the surgeon's learning curve for PSCGs in HTO surgery in order to identify the point of the surgeon's transition from the learning phase to the proficiency phase based on the analysis of the following: (1) the intra-operative parameters measured as operating time, surgeon's comfort levels using the State-Trait Anxiety Inventory (STAI) scoring system, the number of fluoroscopic images and (2) the radiological results measuring the differences between the targeted correction and the achieved correction on pre and post-operative CT-scans, respectively and (3) the functional outcomes at 1 year.

Materials and methods

All patients undergoing medial opening-wedge HTO with PSCGs between February 2014 and November 2017 were enrolled in this single-centre, prospective, observational study. The inclusion criteria were: patients < 60 years old with medial knee osteoarthritis (Ahlbäck ≤ 3 [START_REF] Ahlbäck | X-ray classification and examination technics in gonarthrosis[END_REF]) with significant genu varum (HKA < 176°). The exclusion criteria were: advanced osteoarthritis (Ahlbäck > 3). Patient consent was obtained pre-operatively after they were informed of the procedure in accordance with the principles of the Declaration of Helsinki. All osteotomies were performed by three different surgeons. The first surgeon specialized in sports surgery of the knee and hip was referenced in this study as the "Sports" surgeon (SP). The second surgeon who is a specialist in knee and hip arthroplasty surgery and more specifically reconstructive surgery was referenced as the "Recon" surgeon (MFA). The last surgeon was trained in sports surgery and arthroplasty; we referenced him as the "Combi" surgeon (MO). Details regarding all surgeon's procedures are presented Table 1.

Prior to the start of the enrollment period, all three surgeons had significant experience performing osteotomies (more than 50 cases each) using the conventional technique with the same fixation plate, with no experience of using PSCGs for osteotomy previously. Preoperatively, all patients had conventional radiographs taken (weight-bearing longleg, A/P and lateral views) and underwent a CT-scan using a dedicated protocol including the hip, the knee and the ankle. The targeted correction was defined by the surgeon based on the physical examination and the radiological analysis of the hip-knee-ankle (HKA) angle, medial proximal tibia angle (MPTA) and the posterior proximal tibia angle (PPTA). An order form was then filled by the surgeon to convey the requested corrections to the engineers.

Virtual templating

The CT-scan protocol consisted of acquiring images centred on the femoral head, the knee (allowing the distal femur and 15 cm of the proximal tibia to be captured), and one centred over the ankle. The slice thickness was 0.625 mm for the knee and 2 mm for the hip and ankle (GE Light Speed VCT64). All measurements were constructed using algorithm-calculated constructions on the corresponding tibia which were then mapped. Previous accuracy and reproducibility analysis estimated that this system allows automated measurement of upper femur anatomy with a margin of error < 2 mm and < 1° [START_REF] Banerjee | Innovations in hip arthroplasty threedimensional modeling and analytical technology (SOMA)[END_REF]. The concept was then to create a 3D model of the tibia using the CT-scan images and to mimic the HTO correction following the request of the surgeons with regard to the degrees of correction required in both the frontal and sagittal planes. After simulating the osteotomy on the 3D tibial model, the plate (Activmotion HTO plate, Newclip Technics ® , Haute-Goulaine, France) was virtually positioned in the ideal anatomical location to stabilize the osteotomy. The placement of the screws and screw sizes were also calculated. Once the final construct was virtually created, the pre-osteotomy cutting jig was designed to guide the cut and to drill the final screws hole. Additionally, two k-wires were virtually positioned to protect the lateral hinge during the cut and the opening of the wedge. Intra-operatively, the surgeon positioned the PSCG on the tibia, drilled the holes, performed the cut, and carried out the opening of the wedge. The plate was positioned on the tibia and when the plate holes match the pre-drilled holes on the tibia, the correction was achieved both on the frontal and the sagittal plane using this method (Fig. 1a-f).

Surgical technique

After performing the standard surgical approach, the PSCG was positioned on the tibia and secured by using two k-wires [START_REF] Dessyn | Adding a protective K-wire during opening high tibial osteotomy increases lateral hinge resistance to fracture[END_REF]. The final position of the cutting guide was confirmed with fluoroscopy. The screw holes were then pre-drilled using the jig and dedicated pin were inserted to secure the position of the jig on the tibia. Using the dedicated slot of the jig, the saw blade was guided during the cut. The upper part of the PSCG was then removed to finalize the cut and another slot of the jig was used to guide the bi-planar anterior cut below the tibial tuberosity. The osteotomy cut was then gradually opened/distracted using a laminar-spreader until having the holes of the plate aligned with the predrilled holes on the tibia. The screws were then inserted in the plate to secure the construct. As the screw length was measured on the virtual pre-operative planning, sizing the screws intra-operatively was no longer required. The osteotomy gap was then filled with an injectable phosphocalcic cement (Quickset, Graftys ® , Aix-en-Provence, France).

Post-operative clinical and radiological follow-up

Postoperatively, the same rehabilitation protocol was applied for all patients. They were allowed free range of knee motion and weight-bearing as tolerated protected with the use of crutches. One month post-operatively patients underwent both clinical and radiological examination which included the same standardized CT-scan protocol used pre-operatively. Subsequently this was then used to evaluate the accuracy of the correction compared to the pre-operative planning. Further follow-up visits were scheduled at 3, 6, 9 and 12 months which included a clinical evaluation and standard radiographic analysis (A/P and lateral at 3, 6, 9 and 12 months and long-leg standing at 6 and 12 months).

Patients were allowed to return to their recreational activities at 6 months.

Outcome measures

1. Operating time Operating time was defined as the time from the surgical incision to final wound closure.

Surgeon anxiety levels The Spielberger State-Trait

Anxiety Inventory (STAI) questionnaire is a validated subjective assessment tool to quantify an individual's stress levels with individual traits arising from the clinical environment [START_REF] Marteau | The development of a six-item short-form of the state scale of the Spielberger State-Trait Anxiety Inventory (STAI)[END_REF]. The six-item questionnaire has a 4-point rating scale and total scores range from 6 to 24, with higher values indicating higher levels of stress. The STAI questionnaire was completed by each surgeon at the end of each intervention.

Number of fluoroscopic images taken intraoperatively at

the end of each surgery.

Accuracy of post-operative alignment correction

The accuracy of the postoperative alignment correction was defined by the difference between the targeted correction defined preoperatively and the correction obtained postoperatively measured on the CT-scan (ΔHKA, ΔMPTA, ΔPPTA). For this evaluation the intra-and inter-observer intraclass correlation coefficient were, respectively, 0.86 and 0.81. 5. Functional outcomes To evaluate the functional outcomes, the Knee Osteoarthritis Outcome Score was used [START_REF] Roos | The Knee injury and Osteoarthritis Outcome Score (KOOS): from joint injury to osteoarthritis[END_REF]. Patients were asked to complete this score before the surgery and post-operatively at

Statistical analysis

Sample size calculation was performed using operative time as the primary outcome measure and published data on operative times with similar surgical techniques [START_REF] Arnal-Burró | Three dimensional-printed patient-specific cutting guides for femoral varization osteotomy: do it yourself[END_REF][START_REF] Munier | Can three-dimensional patient-specific cutting guides be used to achieve optimal correction for high tibial osteotomy? Pilot study[END_REF]. The minimal clinical difference was set at 10 min and standard deviation at 10 min [START_REF] Arnal-Burró | Three dimensional-printed patient-specific cutting guides for femoral varization osteotomy: do it yourself[END_REF][START_REF] Munier | Can three-dimensional patient-specific cutting guides be used to achieve optimal correction for high tibial osteotomy? Pilot study[END_REF]. This study required 20 patients in each arm to detect this minimum difference between surgeons in terms of operative time using a two-tailed, twosample t tests with a power of 80% and significance level of 5%. The CUSUM sequential analysis tool was used to assess learning curves for operative time, surgeon stress levels as assessed using the STAI questionnaire, Number of fluoroscopic images taken intraoperatively, and accuracy of lower-limb alignment correction and functional outcomes. Standardized target values for the CUSUM analyses were set using the overall mean values for these outcome measures. CUSUM values represent a running total of the differences between the value of each data point and the standardized target values for each outcome. Learning curves for accuracy of lower-limb alignment correction was assessed using the calculation of the root mean square error values for radiological outcomes. Categorical data were compared using the Chi square test and Fisher's exact test where greater than 25% of cells had less than 5 cases. Normally distributed continuous variables were compared using independent t tests for unpaired variables, paired t test for paired (matched) variables, and one-way ANOVA for multiple variables, The Mann-Whitney test was used for non-parametric data. Statistical significance was set at P < 0.05 for all statistical tests.

Results

The mean operating time for all the patients was 26.3 min ± 8.8. No significant difference was observed between the three surgeons (Table 2). CUSUM analysis for operating time revealed a sharp inflexion point after the first 10 cases for every surgeon (10 cases for "Sport", 7 cases for "Combi" and 9 cases for "Reconstruction"). This analysis helped to identify two distinct phases in the learning curve. Phase 1 represents the initial learning segment and phase 2 represents the proficiency stage (Fig. 2).

The mean surgeon's anxiety level for all three surgeons was 8 ± 3. No significant difference was observed for this parameter between the three surgeons (Table 2). CUSUM analysis of preoperative stress levels as assessed using the STAI questionnaire revealed an inflexion point after 8 cases for every surgeon (5 cases for "Sport", 7 cases for "Combi" and 10 cases for "Reconstruction") (Fig. 3).

The mean number of fluoroscopic images for all the surgery was 4.3 ± 1.2. No significant difference was observed for this parameter between the three surgeons (Table 2). CUSUM analysis of number of fluoroscopic images revealed an inflexion point after 9 cases for every surgeon (10 cases for "Sport", 9 cases for "Combi" and 8 cases for "Reconstruction") (Fig. 4).

The mean ΔHKA was 1.0° ± 1.0°, the mean ΔMPTA 0.5° ± 0.6°, and the mean ΔPPTA was 0.4° ± 0.8°. No significant differences were observed between all three surgeons for these three parameters (Table 2) and the CUSUM analysis did not show any learning effect for these parameters (Fig. 5). The subgroup analysis of the first 10 cases compared with the next 10 (as well as the following 9 for "Sports") for each surgeon, the results did not show any significant difference in postoperative alignment accuracy (Table 3).

No significant difference was observed between the three surgeons for all the items of the KOOS at 1 year post-operative (Table 2). Moreover, the CUSUM analysis did not show any learning effect for this parameter. An average improvement of 29 ± 24 for the KOOS Pain, 29 ± 26 for the KOOS Symptom, 29 ± 27 for the KOOS ADL, 27 ± 34 for the KOOS Sport/Reco and 29 ± 31 for the KOOS QQL was observed (all P < 0.001).

Discussion

The main findings of this study were; using PSCGs during HTO was associated with a learning curve of 10 cases for a decrease in the operating time, 8 cases for an improvement in surgeon comfort levels and 9 cases for a reduction in the number of fluoroscopic images needed during HTO. There was no learning curve to achieve an accurate correction when compared to pre-operative planning and furthermore there was no learning curve to observe an improvement in functional outcome.

With the advance of new technologies in orthopaedic surgery aimed at improving surgical accuracy and improving the comfort of the operating surgeon, the learning curve and the impact on the outcome is often unknown. Currently there are no data on the learning curves associated with PSCG and HTO surgery. Recently some studies have explored the learning curve associated with robot-assisted surgery in knee arthroplasty [START_REF] Kayani | Robotic-arm assisted total knee arthroplasty has a learning curve of seven cases for integration into the surgical workflow but no learning curve effect for accuracy of implant positioning[END_REF][START_REF] Kayani | The learning curve associated with robotic-arm assisted unicompartmental knee arthroplasty: a prospective cohort study[END_REF][START_REF] Sodhi | The learning curve associated with robotic total knee arthroplasty[END_REF]. This present study recognizes the use of PSCG in HTO surgery has a similar short learning phase akin to robot-assisted surgery in knee arthroplasty, and as a whole takes the advantage of exploiting the benefits without adversely affecting the accuracy of the outcome.

Therefore, a comparison of the different parameters within the surgeons learning curve with conventional techniques is necessary. Few studies have evaluated operating time with a conventional HTO technique. Shi et al. [START_REF] Shi | Three dimensional patient-specific printed cutting guides for closing-wedge distal femoral osteotomy[END_REF] in their comparative study of PSCGs versus conventional technique, measured a mean operating time of 96.5 min for the conventional group and 77.7 min for the PSCGs group. In another comparative study of Pérez-Mañanes et al. [START_REF] Pérez-Mañanes | 3D surgical printing cutting guides for open-wedge high tibial osteotomy: do it yourself[END_REF] the mean tourniquet time for the conventional technique was 92 min versus 62 min for the PSCGs group. The analysis of the learning curves obtained in this study showed a benefit in terms of operating time even during the initial learning segment using PSCGs. During the second phase (proficiency stage), the comparison of the results of this study with those previously published showed that a potential shortening of 70% of the operating time can be expected compared to conventional technique. The operating time learning curves were also relatively similar between the three surgeons (10 cases for "Sport", 7 cases for "Combi" and 9 cases for "Reconstruction"), although their experiences and their fields of activity were different. It shows that the relative ease of use of PSCGs is not related to the surgeon's Fig. 3 CUSUM analysis charts demonstrating the mean learning curve time for the surgeon's anxiety level for the three surgeons; a sharp inflexion was observed after 8 cases Fig. 4 CUSUM analysis charts demonstrating the mean learning curve for the number of fluoroscopic images for the three surgeons; a sharp inflexion was observed after 9 cases Fig. 5 CUSUM analysis charts for the ΔHKA for the three surgeons. No learning curve effect was observed background. In this study, improvements in the surgeon's anxiety levels with PSCGs followed the same trend than the trend observed for the operating time with a baseline STAI scores reached after eight cases. Six-item State-Trait Anxiety Inventory is a short-form version of the state scale, consisting of six items chosen for reliability and validity. When compared with the full-form of the STAI, the sixitem version offers a briefer and just as acceptable scale for subjects while maintaining results that are comparable to those obtained using the full-form of the STAI [START_REF] Marteau | The development of a six-item short-form of the state scale of the Spielberger State-Trait Anxiety Inventory (STAI)[END_REF].

The use of fluoroscopy is necessary in the conventional technique as well as in the use of PSCGs. However, in a study by Pérez-Mañanes et al. [START_REF] Pérez-Mañanes | 3D surgical printing cutting guides for open-wedge high tibial osteotomy: do it yourself[END_REF] the mean number of fluoroscopic images required for each intervention was 55 for the conventional technique. In this present study, the mean number of images taken were 4 ± 1 with a learning curve of 9 cases. Therefore, the use of PSCGs drastically reduces radiographic exposure to the patients.

Another objective of this study was to analyse the reliability of this technique to reach the targeted correction. There was no learning curve for achieving the planned correction when compared to the post-operative surgical correction. Few studies have been published on the use of PSCGs for HTO. Victor et al. [START_REF] Victor | Virtual 3D planning and patient specific surgical guides for osteotomies around the knee: a feasibility and proof-of-concept study[END_REF] reported accuracy of 1° in the frontal plane and 2° in the sagittal plane; however, the analysis was performed with radiographs only. Perez-Mananes et al. [START_REF] Pérez-Mañanes | 3D surgical printing cutting guides for open-wedge high tibial osteotomy: do it yourself[END_REF] ]had similar results in the frontal plane with an accuracy of 2°, but the tibial slope was not measured nor was a CT-scan performed. Munier et al. [START_REF] Munier | Can three-dimensional patient-specific cutting guides be used to achieve optimal correction for high tibial osteotomy? Pilot study[END_REF] in their pilot study observed an accuracy measured on a CT-scan of less than 2° in the frontal and sagittal plane. A recent study [START_REF] Jud | Malpositioning of patient-specific instruments within the possible degrees of freedom in high-tibial osteotomy has no considerable influence on mechanical leg axis correction[END_REF] demonstrated that incorrect positioning of patient-specific instruments within the possible degrees of freedom in HTO has no considerable influence on mechanical leg axis correction which further reinforces the reproducibility of this procedure.

The use of navigation is another established technique that theoretically improves the accuracy of postoperative limb alignment correction. Saragaglia et al. [START_REF] Saragaglia | Navigated osteotomies around the knee in 170 patients with osteoarthritis secondary to genu varum[END_REF] reported a 96% accuracy for a planned HKA angle of 184° ± 2° in their study with computer-assisted surgery. The results of this study compares favourably to those of Saragaglia. However, navigation lengthens the operating time (average time of operative procedure: 82 ± 9.3 min for Hankemeier and 74.6 ± 11.4 min for Saragaglia) [START_REF] Hankemeier | Navigated open-wedge high tibial osteotomy: advantages and disadvantages compared to the conventional technique in a cadaver study[END_REF][START_REF] Saragaglia | Navigated osteotomies around the knee in 170 patients with osteoarthritis secondary to genu varum[END_REF]. Published radiological results often point to an under correction after HTO surgery with conventional techniques. The problem is particularly apparent with opening-wedge procedures as the under correction is about 2° more than with closing wedge [START_REF] Brouwer | Osteotomy for medial compartment arthritis of the knee using a closing wedge or an opening wedge controlled by a Puddu plate. A one-year randomised, controlled study[END_REF][START_REF] Nerhus | Radiological outcomes in a randomized trial comparing opening wedge and closing wedge techniques of high tibial osteotomy[END_REF]. Marti et al. [START_REF] Marti | Accuracy of frontal and sagittal plane correction in open-wedge high tibial osteotomy[END_REF] reported a 50% accuracy between the correction achieved and the one planned after openingwedge HTO. This was vastly different from this study where the use of PSCGs led to a reliable correction that was similar to the planned correction.

Regarding the functional results at 1 year evaluated by the different items of KOOS, a significant improvement of these scores without a learning curve effect and without significant difference between the three surgeons was observed. The average gains observed in this study for each KOOS item seem superior to those using conventional osteotomy techniques [START_REF] Duivenvoorden | Comparison of closing-wedge and opening-wedge high tibial osteotomy for medial compartment osteoarthritis of the knee: a randomized controlled trial with a six-year follow-up[END_REF][START_REF] Lind-Hansen | Open-wedge high tibial osteotomy: RCT 2 years RSA follow-up[END_REF][START_REF] Toksvig-Larsen | Ten-year results of physical activity after high tibial osteotomy in patients with knee osteoarthritis[END_REF]. However, the short decline and the absence of a control group do not currently allow this study to conclude on the superiority of PSCGs in terms of functional results. It will be necessary to continue to monitor these patients to determine whether an improved accurate correction affects the functional outcomes of patients in the future. There are several limitations of this study that must be considered when interpreting these findings. Firstly, the surgeons who practised these procedures were familiar with HTO surgery fixation and had already used custom cutting guides in knee arthroplasty surgery. Their learning curve may not be directly transferrable to other less experienced surgeons. Secondly, follow-up time was limited to 1 year following surgery and so long-term data on functional outcomes and revision rates were not available. Thirdly, it is a non-comparative study and although data exist on the accuracy of correction, the operative time and the use of fluoroscopy during conventional surgery, it would have been valuable to be able to make a direct comparison using a control group. Despite these limitations, we believe that this project is original and fairly designed to analyse the learning curve in using PSCGs in HTO surgery and the accuracy of the correction obtained. Moreover, a significant number of patients (71) compared to studies using PSCGs in HTO were included (10 for Munier et al. [START_REF] Munier | Can three-dimensional patient-specific cutting guides be used to achieve optimal correction for high tibial osteotomy? Pilot study[END_REF] 8 for Perez-Mananes et al. [START_REF] Pérez-Mañanes | 3D surgical printing cutting guides for open-wedge high tibial osteotomy: do it yourself[END_REF] and 14 for Victor et al. [START_REF] Victor | Virtual 3D planning and patient specific surgical guides for osteotomies around the knee: a feasibility and proof-of-concept study[END_REF]). The use of a CT-scan to measure the accuracy of the postoperative correction made it possible to limit the uncertainty of the measurement on standards radiographs [START_REF] Boonen | Inter-observer reliability of measurements performed on digital long-leg standing radiographs and assessment of validity compared to 3D CT-scan[END_REF]. The HTOs were performed by three different surgeons with different experiences of knee surgery, which schematically represents the different types of practitioners likely to use PSCGs in their daily practice.

The results of this study should reassure surgeons wanting to start using PSCGs during HTO in their daily practice. Indeed, the use of PSCGs requires a short learning curve to increase the efficiency of operating time and reduction in the use of fluoroscopic imaging. But, there is no impact of cumulative experience on accuracy of limb alignment, which is important for the safe implementation of this procedure into routine surgical practice.

Conclusion

In conclusion, the use of PSCGs requires a short learning curve to increase the efficiency of operating time and reduction in the use of fluoroscopic imaging. There is no learning curve associated with improving the accuracy of postsurgical HTO correction and the functional results at 1 year.

Fig. 1

 1 Fig. 1 Preoperative 3D templating illustrating the operative technique of custom cutting guides. a, b Optimal positioning of the cutting guide with two pins (one to protect the lateral hinge during the osteotomy and one to guide the saw). c Pre-drilling the eight holes of the

Fig. 2

 2 Fig.2CUSUM analysis charts demonstrating the mean learning curve for operative time for the three surgeons; a sharp inflexion was observed after 10 cases separating the learning phase to the proficiency phase

Table 1

 1

	Pre-operative data in PSCGs/all the P value between patients undergoing HTO with		(n = 21) Combi's surgeon	(n = 29) Sport's surgeon	surgeon (n = 21) Reconstruction's	Total (n = 71)
	the three surgeons were non-significant	Age (years) Male gender (%)	46.1 ± 4.2 44	42.1 ± 9.0 47	46.2 ± 4.3 47	44.3 ± 6.8 45
		HKA (°)	172.1 ± 2.6	171.7 ± 2.7	170.6 ± 2.6	171.6 ± 2.8
		MPTA (°)	80.8 ± 2.6	79.7 ± 3.2	77.9 ± 2.8	79.8 ± 3.1
		PPTA (°)	9.4 ± 2.5	8.7 ± 3.3	10.7 ± 2.6	9.3 ± 2.9
		KOOS Pain (/100)	52 ± 17	57 ± 16	55 ± 16	56 ± 15
		KOOS Symptoms (/100)	49 ± 15	55 ± 17	52 ± 16	53 ± 16
		KOOS ADL (/100)	50 ± 14	54 ± 16	51 ± 15	52 ± 16
		KOOS Sport/Recon (/100)	44 ± 17	38 ± 21	41 ± 19	40 ± 21
		KOOS QQL (/100)	43 ± 23	50 ± 19	46 ± 23	47 ± 22
		HKA hip-knee-ankle angle, MPTA medial proximal tibial angle, PPTA posterior plateau tibial angle

  1 year. The KOOS is a validated knee specific PROM. It holds 42 items in five separately scored subscales: KOOS Pain, KOOS Symptoms, Function in daily living (KOOS ADL), Function in Sport and Recreation (KOOS Sport/Rec), and kneerelated Quality of Life (KOOS QOL). The functional improvements were measured as the differences between the pre-operative values and values measured at 1 year postoperatively for each sub-category of the KOOS and defined as: ΔPain, ΔSymptoms, ΔADL, ΔSport/Rec and ΔQQL for each patient.During the enrolment period, 71 patients were included (21 performed by "Combi", 21 by "Reconstruction" and 29 by "Sport"). The mean patient age was 44.3 ± 6.8 No patient was lost to follow-up at 1 year. Pre-operatively the mean HKA angle was 171.6 ± 2.8; the mean MPTA was 79.8° ± 3.1°; the mean PPTA was 9.3° ± 2.9° (Table1). Local Ethical Committee approval was obtained prior to study's initiation (Comité Informatique et Liberté (CIL)/ Assistance Publique des Hopitaux de Marseille (AP-HM)/ Registration Number: 2014-36).

Table 2

 2 Post-operative data in patients undergoing HTO with PSCGs/all the P value between the three surgeons were nonsignificant (n.s.) STAI The Spielberger State-Trait Anxiety Inventory, Post-operative bleeding bleeding was defined by the quantity of blood in the drain on day 2, ΔKOOS Pain, ΔKOOS Symptoms, ΔKOOS ADL, ΔKOOS Sport/ Reco and ΔKOOS QQL difference between the value obtained in the pre-operative questionnaire and that obtained at 1 year postoperatively

		Combi's sur-	Sport's sur-	Reconstruction's	Total (n = 71)
		geon (n = 21)	geon (n = 29)	surgeon (n = 21)
	Operative time (min)	28.9 ± 9.7	25.2 ± 6.6	27.0 ± 10.1	26.3 ± 8.8
	Number of fluoroscopy	4.7 ± 1.2	4.1 ± 1.1	4.6 ± 1.5	4.3 ± 1.2
	Surgical surgeon anxiety's level	7 ± 2.7	7 ± 2	9 ± 3.5	8 ± 3
	ΔHKA	0.7 ± 0.8	1.2 ± 1.0	0.9 ± 1.0	1.0 ± 1.0
	ΔMPTA	0.4 ± 0.5	0.7 ± 0.7	0.4 ± 0.5	0.5 ± 0.6
	ΔPPTA	0.3 ± 0.5	0.4 ± 0.7	0.6 ± 0.9	0.4 ± 0.8
	ΔKOOS Pain	33 ± 27	28 ± 22	27 ± 25	29 ± 24
	ΔKOOS Symptoms	26 ± 28	29 ± 24	33 ± 23	29 ± 26
	ΔKOOS ADL	28 ± 25	28 ± 27	30 ± 29	29 ± 27
	ΔKOOS Sport/Reco	19 ± 31	32 ± 35	33 ± 29	27 ± 34
	ΔKOOS QQL	33 ± 31	20 ± 36	32 ± 29	29 ± 31

Table 3

 3 Accuracy of postoperative limb alignment for the three surgeons, subgroup analysis of 10 patients HKA hip-knee-ankle angle, MPTA medial proximal tibial angle, PPTA posterior plateau tibial angle/ the accuracy of the postoperative alignment was defined by the difference between the desired correction defined preoperatively and the correction obtained postoperatively measured on CT-scan (ΔHKA, ΔMPTA, ΔPPTA)

		Cases 1-10	Cases 11-20	Cases 21-29	P value
	Sport's surgeon
	ΔHKA (°)	1.3 ± 0.7	0.8 ± 1.0	1.2 ± 0.9	n.s
	ΔMPTA (°)	0.9 ± 1.1	0.7 ± 0.6	0.6 ± .0.8	n.s
	ΔPPTA (°)	0.4 ± 0.4	0.4 ± 0.6	0.2 ± 0.4	n.s
	Combi's surgeon
	ΔHKA (°)	0.7 ± 0.7	0.8 ± 0.9	-	n.s
	ΔMPTA (°)	0.4 ± 0.7	0.5 ± 0.5	-	n.s
	ΔPPTA (°)	0.3 ± 0.6	0.3 ± 0.4	-	n.s
	Reconstruction's surgeon
	ΔHKA (°)	0.7 ± 0.7	0.8 ± 0.9	-	n.s
	ΔMPTA (°)	0.4 ± 0.7	0.5 ± 0.5	-	n.s
	ΔPPTA (°)	0.3 ± 0.6	0.3 ± 0.4	-	n.s

International Knee and Joint Centre, Abu Dhabi, UAE

Funding No funding was required for this study.

Compliance with ethical standards

Conflict of interest Some of the authors disclosed potential conflict of interest.

Ethical approval

Local ethic committee approval was obtained prior to study Initiation (N°2014-36).


F S Gentile

Rérat M § , D'arco

Substitutional carbon defects in Silicon. A quantum mechanical characterization through the IR and Raman spectra

Keywords: Silicon, substitutional carbon defect, Ab-initio calculation, IR and Raman spectra

The Infrared and Raman spectra of eight substitutional carbon defects in silicon are computed at the quantum mechanical level by using a periodic supercell approach based on hybrid functionals, an all electron Gaussian type basis set and the Crystal code. The single substitutional C s case and its combination with a vacancy (C s V and C s SiV) are considered first. The progressive saturation of the four bonds of a Si atom with C is then examined. The last set of defects consists of a chain of adjacent carbon atoms C i s , with i = 1 to 3. The simple substitutional case, C s , is the common first member of the three sets. All these defects show important, very characteristic features in their IR spectrum. One or two C related peaks dominate the spectra: at 596 cm -1 for C s (and C s SiV, the second neighbor vacancy is not shifting the C s peak), at 705 and 716 cm -1 for C s V, at 537 cm -1 for C 2 s and C 3 s (with additional peaks at 522, 655 and 689 for the latter only), at 607 and 624 cm -1 , 601 and 643 cm -1 and 629 cm -1 for SiC 2 s , SiC 3 s , SiC 4 s , respectively. Comparison with experiment allows to attribute many observed peaks to one of the C substitutional defects. Observed peaks above 720 cm -1 must be attributed to interstitial C or more complicated defects.

The Infrared and Raman spectra of eight substitutional carbon defects in silicon are computed at the quantum mechanical level by using a periodic supercell approach based on hybrid functionals and all electron Gaussian type basis set. All these defects show important, very characteristic features in their vibrational spectrum. Comparison with experiment allows to attribute many observed peaks to one of the C substitutional defects.

INTRODUCTION

Electronic and mechanical properties of silicon are known to be significantly affected by the presence of impurities and defects. Carbon, along with oxygen, represents one of the most common impurity. Carbon is substitutionally dissolved as isovalent impurity and can interact with other point defects. It was found that implantation-generated excess of silicon self-interstitials after boron implantation is greatly reduced in the presence of carbon 1 . Carbon substitutional is considered immobile 2 and characterized by single vibrational mode at 608 cm -1 , while carbon interstitial is highly mobile 3 . Carbon substitutional reacts with silicon self-interstitials, producing a new carbon interstitial defect. EPR measurements suggest that this defect has C 2v symmetry and a doublet electronic ground state.

In the present study we investigate eight defects in silicon, limiting the analysis to substitutional carbon. The aim is to explore whether the IR and/or Raman spectra provide enough information for the clear identification of the most important of these defects. As these defects are EPR inactive (they have a closed shell ground state, or the most stable state is a singlet, as for the two here considered that couple C s and V, a vacancy), the vibrational analysis is the most general tool available for this characterization. The extension to defects involving interstitial carbon, possibly in combination with N or O will be the object of a forthcoming paper 4 .

The eight investigated defects can be grouped in three sets, the first one includes the simple substitutional C s defect, and its perturbed cases C s V and C s SiV, in which C s is a first or second neighbor of a vacancy. The second set includes 4 cases, to be indicated as SiC i s , with i=1 to 4, in which the four neighbors of a Si atom are progressively substituted by C atoms.

The third set includes three defects, C i s , with i=1 to 3, in which a chain of adjacent carbon atoms has been taken into account.

The supercell scheme (cubic supercells containing 64 and 216 atoms will be considered), the very effective hybrid B3LYP 5,6 functional, an all electron Gaussian type basis set and the Crystal code have been used. As outlined previously [7][8][9][10][11][12] , the B3LYP hybrid functional, with 20% non-local exact exchange, leads to a proper description of (electronic) spin polar-ization, and accurate vibrational spectra.

The paper is structured as follows. Section METHODOLOGY contains the details of the computational procedures and conditions. Section RESULTS describes the geometry, the charge and spin (when applies) density, and the vibrational features of the eight defects. In Section DISCUSSION the results are discussed and compared with experiments. A few conclusions are drawn in section CONCLUSIONS. A Supplementary Material Section contains some of the technical details concerning the vibrational part, and shows the Raman spectra of the eight defects.

METHODOLOGY

Calculations have been performed by use of the B3LYP global hybrid functional, 5,6 as implemented in the Crystal17 rogram. 13 An all electron Gaussian type basis set (derived from Pople's standard 6-31G 14 ) as described in D'Arco et al. 15 has been adopted for silicon: the contraction is 6-5111G* (22 atomic orbitals, AOs), where the first shell is of s type, followed by 4 sp and one d shell; the exponents of the two outermost sp shells have been set to 0.45 Bohr -2 and 0.15 Bohr -2 , respectively. For the carbon atom, the 6-31G* 16 contraction has been adopted without modifications (14 AOs both).

The truncation of the Coulomb and exchange infinite lattice series is controlled by five thresholds T i (see Crystal manual, Ref. 17, for more details), which have been set to 8 (T 1 -T 4 ) and 16 (T 5 ). The convergence threshold on energy for the self-consistent-field (SCF) procedure has been set to 10 -8 and 10 -11 E h for the structural and frequency calculations, respectively.

The DFT exchange-correlation contribution to the Fock matrix has been evaluated by numerical integration over the unit cell volume. Radial and angular points for the integration grid are generated through Gauss-Legendre radial quadrature and Lebedev two-dimensional angular point distributions. The default pruned grid with 74 radial and 974 angular points has been used.

A periodic supercell approach is used to simulate different defect concentrations. Cubic super-cells containing 64 and 216 atoms have been considered (to be referred to in the following as S n , with n= 64 or 216). They correspond to the expansion of the bulk silicon conventional cell, containing 8 atoms, by a factor 2 and 3 along the three lattice parameters.

The defects reduce the symmetry; the deformation of the cell from the cubic shape is however negligible. Overall, the size of the basis set (number of AOs) is (about) 1400 in S 64 and (about) 4744 in S 216 . Reciprocal space has been sampled by using a regular sub-lattice with a shrinking factor 4 for S 64 and 2 for S 216 ; this corresponds to a number of k-points ranging from 4 to 36, according to the point symmetry of the defective system.

The experimental Raman spectrum of pristine silicon 18 is characterized by a single, three fold degenerate, peak at 521 cm -1 . The calculated peak with the present functional and basis set is at 533 cm -1 , 12 cm -1 above. So any observed IR signal or any Raman peak other than the one at 521 cm -1 must be attributed to defects.

As the mass of carbon is less than half the one of silicon (12.01 vs 28.09 u), we can anticipate that many of the modes to which C contributes will appear well above the phonon band of pristine silicon, that has a calculated higher limit at 533 cm -1 . Besides that, the C-Si bonds are polar (typical electronegativity values are 1.9 for Si and 2.6 for C); as a consequence, the IR defect peaks are expected to be very intense, and to dominate the spectrum.

The identification of the nature of the modes (which atom is contributing? what kind of movement is described by the mode?) can easily be performed through two tools implemented in the CRYSTAL code: a) the isotopic shift (see the note following eq. 1 in the Supplementary Material Section ); at negligible cost the mass of a subset of atoms can be altered arbitrarily. Note that at variance with respect to experiment, i) only very particular atoms of the cell can be involved in the isotopic substitution; ii) it is not necessary that the attributed mass corresponds to an existing isotope. Very large masses can be attributed to atoms in order to magnify the mass effect. b) the graphical animation of the modes. This permits to have a direct and clear evidence of the atoms and groups providing important contributions to the eigenvectors; it is available at www.crystal.unito.it in the Defects in Silicon section.

RESULTS

In this section we consider eight defects containing substitutional carbon, organized in three series.

The first series contains the isolated substitutional carbon, C s , and the combination of C s with a vacancy V: C s SiV, when C is a second neighbor of V, and C s V, when C s is sitting on the border of V. The ground state of C s is a singlet closed shell, whereas in the two other cases the four first neighbors of the vacancy carry a uncoupled electron, so that a quintuplet q (S z =2), a triplet t (S z =1) and a singlet s (S z =0) spin states are possible. The singlet is more stable than the triplet by 0.10 (C s SiV) and 0.72 eV (C s V), respectively. The q state is less stable than the triplet (see Table 2).

In Figure 1, top, these three defects are schematically represented, with atom-atom distances R and Mulliken charges q and bond populations b. The electronegativity of C is larger than the Si one (2.6 vs 1.9), so that the Mulliken net charge q of C is -0.29 |e| in C s ; it increases The IR spectrum of C s is very simple (see Figure 3, left). It is dominated by a single peak at 596 cm -1 , whose intensity is 2015 km/mol. In C s SiV the peak splits in three, and blue shifts by a few cm -1 (600, 601 and 617 cm -1 ); the sum of the intensities of the three peaks is about the same as the one of C s . When carbon is on the border of the vacancy, the shift is larger (more than 100 cm -1 ). The lower symmetry splits the single peak of C s in two peaks at close wavenumbers, 705 and 716 cm -1 , with a much smaller intensity.

The second series includes C s , C 2 s and C 3 s , a chain of 1 to 3 adjacent substitutional carbon atoms. The ground state is always a closed shell; charges and distances are shown in Figure 1, second row. The carbon net charge is as high as -0.55 and -0.58 |e| in C 2 s , and in C 3 s (on the two terminal atoms; for the central one it is -0.31 |e|). The C-C distances are much larger than in diamond (1.81 and 1.97 Å in the double and triple substitutional case), due to the rigidity of the bulk Si lattice.

As regards the IR spectrum, in C 2 s and C 3 s the 596 cm -1 C s peak red shifts by 60 cm -1 to 537 cm -1 . In C 3 s , three additional peaks of lower intensity appear at 585 (11 cm -1 below the C s peak), 655 and 689 cm -1 .

The former corresponds to a C-C stretching, the latter to the C-Si stretching of the central carbon atom, as the graphical animation of teh modes confirms. The 585 cm The IR spectrum of the three SiC x s defects, with x =2 to 4, is shown in Figure 3, right column. The peak at 596 cm -1 of C s splits in two peaks at 607 and 624 cm -1 (SiC 2 s ), and at 601 and 643 cm -1 (SiC 3 s ), by effect of the reduced symmetry (C 2v and C 3v , respectively). In SiC 4 s a single peak (symmetry is T d ) at 629 cm -1 dominates the spectrum, intermediate between the two peaks of SiC 2 s and SiC 3 s . The isotopic substitution 12 C → 13 C red shifts these peaks by 18-20 cm -1 , confirming that they are related to carbon. Less intense peaks in the range 520-548 cm -1 appear in the 3 cases. Their isotopic shift is only 3-5 cm -1 large, indicating that they are due to the perturbation of the Si lattice around the defect, rather than being directly connected to C. The defect formation energies E f were calculated according to the following equation:

E f = (E D -n * E C ) -(E P -m * E Si ) (1)
where E D and E P are the total energies of the defective (D) and pristine (P ) silicon supercells (S 64 or S 216 ), n and m are the numbers of carbon atoms added to, or silicon atoms removed from, the defective supercells. Formation energies for all studied defects are reported in Table 2. E Si is the total energy of the silicon atom in the pristine silicon structure (that is, half the electronic total energy of the primitive cell, that with the present basis set and functional is -289.46807 E h per atom). As regards the definition of E C , we have followed two different approaches, denoted by superscripts 1 and 2. In the first one (E 1 f ), E C is obtained in a similar way from pristine diamond (-38.08851046 E h per atom). In the second definition,

E 2
f is the energy of an isolated carbon atom. The difference between the two is obviously nothing else than the binding energy of the diamond crystal. E 1 f for C s is 1.5 eV. The one of C 2 s is not twice the one of C s , but about three times, and three more eV are necessary to create C 3 s . This discontinuity is confirmed by a parallel discontinuity in the C s Si distances, and C s net charges.

If we compare the two systems with two or three C atoms, the table shows that localizing the perturbation around a single Si atom is strongly favoured: 3.11 (SiC 2 s ) vs 4.57 (C 2 s ) eV, and 4.84 (SiC 3 s ) vs 7.71 (C 3 s ) eV. Looking at the SiC i s series, the progressive insertion of one, two, three and four C atoms requires 1.5, 1.61, 1.73 and 1.68 eV. We considered also the formation energy of C-Si-Si-C, when the two carbon atoms are linked to two first neighbors silicon atoms. The E 1 f formation energy is 3.88 eV, much lower than the C 2 s one (4.57 eV) but much higher that the SiC 2 s one (3.11 eV). It is clear that if n carbon atoms are present in bulk silicon, they will tend to aggregate.

The formation energy of a single silicon vacancy obtained from our model is about 4.7 eV (the experimental estimated value is 4.0 eV 19 ). E 1 f for C s SiV is 6.52 eV, slightly larger (0.32 eV) than the sum of the formation energy of the two defects (4.7 +1.5 =6.2 eV). E 1 f of C s V is, on the contrary, exactly equal to the sum of E 1 f for V and C s .

DISCUSSION

Associating an experimental IR, Raman or EPR signal to a specific defect is challenging and, very often, controversial. In this section we try to connect the microscopic structure of defects to experiments comparing calculated and measured IR spectra.

In 1965 Newman 2 published the IR spectrum of carbon doped monocrystalline silicon. The spectra collected at 290 K were characterized by a single peak at 608 cm -1 . For a sample containing both 12 C and an increased percentage of 13 C, two peaks at 604 and 586 cm -1

were observed; the second one was attributed to 13 C, with an isotopic shift of 18 cm -1 .

The dominant peak in our simulated spectrum of C s is at 596 cm -1 , that downshifts by 18 cm -1 to 578 cm -1 when 13 C is used instead of 12 C (see Table 1). The present scheme (B3LYP plus the basis set described in the computational section) underestimates the dominant Raman peak of perfect bulk silicon by 12 cm -1 (from 533 to 521 cm -1 ). It is reasonable to suppose that a similar underestimation is affecting the spectra of the defective systems.

Taking into account this red shift, it can be concluded that the peaks observed by Newman can be attributed to a single carbon substitutional defect, C s .

In the same experiment two more modes were observed when collecting spectra at 63 K. For 12 C plus 13 C doped samples, absorption bands at 1217 and 1175 cm -1 with much smaller intensities were attributed to possible second harmonics of the local mode vibrations. We are collecting some evidences that interstitial carbon defects could be, alternatively, responsible for these peaks 4 .

In 1970 Bean et al. 20 reported a large list of modes that could be associated to oxygen and carbon defects or their complexes. For samples containing only carbon, four low intensity modes in the range from 528 cm -1 to 552 cm -1 were reported, without attribution to specific local structures.

Our calculated IR spectrum of the C 2 s defect (two vicinal substitutional carbon atoms, see figure 4), is dominated by a single peak at 537 cm -1 , well within the range of the peaks reported by Bean et al.. Other C substitutional defects show peaks in the 528-552 cm -1 range: C 3 s , with the two dominant peaks at 537 (as for C 2 s ) and 527 cm -1 ; SiC 3 s and SiC 4 s present low intensity peaks in the range 522-532 cm -1 . Note however that in the case of SiC 3 s and SiC 4 s also the dominant peaks at 601 and 643 cm -1 for the first, and at 629 cm -1 for the second, should appear in the experimental spectrum. So the C 2 s and C 3 s defects are the best candidates for the attribution of these experimental peaks.

Bean et al. 20 reports also about a broad band in the range 704-735 cm -1 . The authors suggest that this band could be related to C-V centers. Inspection of figures 3 and4 shows that all C substitutional defects present peaks below 690 cm -1 , the exception being C s V, a carbon first neighbor of a vacancy, that in its singlet state (the most stable), has two dominant peaks at 705 and 716 cm -1 . We can then confirm the attribution by Bean et al..

CONCLUSIONS

The IR and Raman spectra, as well as the local geometry, the charge and spin (when applies) distribution of eight substitutional carbon defects are produced. The Raman spectrum is in all cases dominated by a peak at 530 cm -1 , corresponding to the only peak of pristine bulk silicon. The landscape of the IR spectra is much more varied, and allows in most of the cases a clear identification of the specific defect. C s has a single very intense peak at 596 cm -1 , whereas C s V is the only one with peaks above 700 cm -1 (705 and 716 cm -1 ).

With respect to C s , the SiC i s defects blue shift to 607 and 624 cm -1 , and red shift to 548 cm -1 (SiC 2 s ), blue shift to 601 and 643 cm -1 (SiC 3 s ), or to 629 cm -1 (SiC 4 s ), whereas C 2 s ) and C 3 s ) red shift to 523-533 cm -1 . The IR spectrum is then able to characterize all these defects. The most relevant experimental peaks have been attributed to some of these defects.

Peaks above 720 cm -1 are most probably attributable to interstitial C defects that will be the object of a forthcoming investigation.

The present study confirms that modern, high quality quantum-mechanical computer codes are extremely useful tools for the analysis of the properties of defects in crystalline compounds. They can be considered as standard tools, that complement in many cases, and substitute in others, experimental investigations, when the latter are difficult, or very expensive, or unable to discriminate between different hypotheses. -0.549 -0.549

(

  in absolute value) to -0.44 and -0.55 |e| in the two defects including the vacancy. The R distances, that are 1.56 Å in pristine diamond and 2.57 Å in silicon, are in between these values for the various C-Si distances. The bond populations b are intermediate between the ones of diamond and silicon.

  -1 mode corresponds to C-Si stretching of two outer carbon atoms. The C-Si peak is absent in C 2 s for symmetry reasons. It appears on the contrary in the Raman spectrum, see the Supplementary Material, Figure S1). The mode at 651 cm -1 corresponds to the asymmetric C-Si stretching. The third series includes SiC 1 s ≡ C s , SiC 2 s , SiC 3 s and SiC 4 s (second row in Figure 1), corresponding to substitute one, two, three and four neighbors of a central silicon atom. The carbon net charge, that in C s is -0.290 |e|, and is more than compensated by the four Si atoms (+0.110 |e|), jumps to -0.824, -0.825 and -0.896 |e| in the three other cases, whereas the central Si atom has a net charge of +0.475, +0.687 and +0.896 |e|. The net charge of the "external" Si atoms is constant at about +0.240 |e| for the three cases. The Si-C bond population reduces from +0.34 |e| for C s to about +0.26 |e| for the three other cases, as a consequence of the increased distance between the central Si atom and its neighbors from 1.82 to 2.11, 2.16 and 2.19 Å along the series.

Figure 1 :

 1 Figure 1: The local cluster around the substitutional carbon and vacancy defects here considered. Bond lengths (in Å), Mulliken net charges (in |e| and bold) and bond populations (in |e| and italic) are reported for nearest non equivalent atoms. Superscript s indicates an open shell singlet. Carbon atoms are in grey, the vacancy position is identified by the white circle.

Figure 2 :Figure 3 : 3 Figure 4 :

 2334 Figure 2: The local cluster around the C x s defects here considered. Bond lengths (in Å), Mulliken net charges (in |e| and bold) and bond populations (in |e| and italic) are reported for nearest non equivalent atoms. Carbon atoms are in grey.

  ), computed with reference to pristine diamond and silicon bulk, and to the isolated carbon atom. The S 64 supercell has been used for the first four defects, S 216 for the more extended defects, as indicated. 18 Institute of Solid State Physics, University of Latvia as the Center of Excellence has received funding from the European Union's Horizon 2020 Framework Programme H2020-WIDESPREAD-01-2016-2017-TeamingPhase2 under grant agreement No. 739508, project CAMART²

Table 1 :

 1 Wavenumber shift δ (in cm -1 due to the 12 C → 13 C isotopic substitution) of the most intense IR modes. The absolute intensities (in km/mol) are also reported. A semicolon separates the information concerning one mode from the one of another mode.

	Defect	Sym			ν( 12 C) →ν ( 13 C), δ; (intensity)
	C s	T d			596→578, 18 (2015)
	(C s -Si-V) s	C s		617→598, 19 (791); 601→583, 18 (632); 600→582, 18 (615)
	(C s -V) s	C s		716→693, 23 (344); 705→683, 22 (355)
	C 2 s	C 3v		651→628, 23 (R); 544→543, 1 (R); 537→537 0(941)
	C 3 s	C 2v 689→664, 25 (144); 655→631, 24 (196); 612→592, 20 (R); 585→567, 8 (25); 537→536, 1 (299
	SiC 2 s	C 2v		624→605,19 (1522); 612→592, 20 (R); 1294, 548→540, 8 (943)
	SiC 3 s	C 3v		642→623,19 (2298); 601→582, 19 (1850); 540→535, 5 (653)
	SiC 4 s	T d		629→609,20 (5692); 540→537, 3 (559)
	Method		E D		E 1 f (eV)			E 2 f (eV)
			(Ha)	s	t	Q	s	t	Q
	C s		-18274.52383	1.50			-6.55
	C s -Si-V		-17984.83849	6.52 6.61 7.42	-1.54 -1.44 -0.64
	C s -V		-17984.85291	6.19 6.92 7.02	-1.86 -1.14 -1.03
	C 2 s		-18023.03141	4.57			-11.53
	C 3 s (S 216 )		-61770.68856	7.71			-16.45
	SiC 2 s (S 216 )	-62022.23700	3.11			-12.99
	SiC 3 s (S 216 )	-61770.79382	4.84			-19.31
	SiC 4 s (S 216 )	-61519.35278	6.52			-25.69

Table 2 :

 2 Total (E D , in E h ) and formation E 1 f (in eV) energy of the substitutional defects in silicon (see equation 1

ACKNOWLEDGMENTS

Access to the HPC resources of CINES/IDRIS/TGCC obtained thanks to the grant 2018-[A0050810537] (Ph. D'Arco) made by GENCI is warmly acknowledged.

((Additional Supporting Information may be found in the online version of this article.))