Antoine Aerts

Philippe Carbonnière

Falk Richter

Alex Brown
email: alex.brown@ualberta.ca

HCOOD Antoine Dcooh

Aerts

P Carbonniére

Vibrational states of deuterated trans-and cis-formic acid

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

I. INTRODUCTION

Formic acid (HCOOH) is the simplest carboxylic acid and, hence, serves as a benchmark system for this class of molecules. Moreover, as the smallest carboxylic acid, it provides a test case for ab initio studies, including the determination of vibrational state energies. Formic acid is the most abundant organic acid in the atmosphere 1 and, also, is of great astronomical importance especially for its potential role as an organic building block. 2 Formic acid, which was first detected in the Sg B2 star-forming region in 1971 by Zuckerman et al., 3 has also been observed towards hot molecular cores. 4 Formic acid features two conformational isomers, i.e., trans and cis. Based on a recent full-dimensional anharmonic computational study by Richter and Carbonnière, 5 the cisisomer is about 1415 cm -1 higher in energy than the trans-isomer when comparing zeropoint-energies. Another computational study on formic acid by Tew and Mizukami 6 at a similar level of theory predicted this value to be 1412 cm -1 . The isomerism enriches the spectra and could manifest as interesting dynamics in the cis-trans isomerization process.

However, the trans form is about 800 times more abundant than the other rotameric form at room temperature. Therefore, the observation of the cis form in the gas phase requires a long absorption path (and no band-overlap) or more elaborate techniques such as demonstrated recently by Meyer and Suhm. 7 Another approach to detecting the cis-form involves using rare-gas matrix isolation, and this has been the primary method used for the experimental study of the cis form. However, this method induces environmental effects on the vibrational spectrum in going from the gas phase to the matrix. Moreover, the observed shifts in peak positions are dependent on the vibrational state and, thus, do not permit a systematic comparison to gas-phase theoretical results as shown by Meyer and Suhm. 7 Therefore, if one wishes to benchmark theoretical approaches for obtaining vibrational energies, gasphase measurements must be conducted directly for both isomers and their isotopologues.

Experimental data on the cis form are scarce [7][8][9][10] , especially for the deuterated species, and, to the best of our knowledge, results for the deuterated forms are available only from raregas matrices isolation measurements; 11,12 however, new gas phase measurements are being undertaken. 13 In the present work, we focus on the computation of fundamental vibrational states as well as combination bands and overtones of the deuterated forms of trans-formic acid, i.e., 2 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI:10.1063/1.5135571 HCOOD, DCOOH, and DCOOD. We also report the fundamental transition frequencies of the deuterated forms of cis-formic acid. Although the interest in the deuterated forms of formic acid results from the astrophysical interest in the non-deuterated species, spectroscopic evidence of HCOOD and DCOOH in the interstellar medium has been reported. 14,15 Although the doubly-deuterated form, i.e., DCOOD, is less-likely to be observed, doubly deuterated molecules have been observed in the interstellar medium, e.g., water (D 2 O), 16 ammonia (ND 2 H), 17,18 D 2 S, 19 and formaldehyde (D 2 CO). 20 Moreover, even triply deuterated species have been detected, e.g, ammonia (ND 3 ). 21 The vibrational states of formic acid (HCOOH) were determined in a number of theoretical studies. [5][6][7]11,12,[22][23][24][25][26][27][28] Tew and Mizukami 6 computed vibrational eigenstates using vibrational configuration interaction on a global potential energy surface (PES) constructed with 17076 CCSD(T)(F12*)/cc-pVTZ-F12 energies valid up to 15 000 cm -1 . Their theoretical trans fundamental energies have a root mean square deviation of 3.0 cm -1 with respect to the gas phase measurements. Maçôas et al. 11 assigned the overtone spectra of HCOOH and DCOOH using a correlation corrected vibrational self-consistent field technique on a PES constructed with 30816 points at the MP2/6-311++G(2d,2p) level of theory. Demaison et al. 23 computed a semidiagonal quartic force field for both isomers at the CCSD(T)/cc-pVTZ level of theory, derived semiexperimental equilibrium structures of HCOOH, and analysed some of its resonance polyads. Other studies provided equilibrium structures and the anharmonic ν OH stretching frequency using the MP2/6-311++G (3df,3pd) level of theory 24 as well as ν OD using the density functional tight binding formalism (DFTB). 25 The ν OH overtones were studied using a 1D potential computed at the MP2/cc-pVTZ level 26 and fundamentals, combination bands, and overtones were computed with vibrational self-consistent field (CC-VSCF) using second-order perturbation correction at the MP2/6-311++G(2d,2p) level for HCOOD. 12 The computational benchmark of Meyer and Suhm 7 with systematic comparison of band positions with the gas-phase measurements showed the limitation of density functional theory based methods while suggesting the double hybrid B2PLYP-D3(BJ)/augcc-pVQZ and MP2/aug-cc-pVQZ methods using second order vibrational perturbation theory (VPT2) to obtain the vibrational energies as a good compromise between cost and accuracy. Going beyond VPT2, most recently, Richter and Carbonnière 5 computed the vibrational eigenstates of formic acid using the block improved relaxation method [29][30][31][START_REF]Multidimensional Quantum Dynamics: MCTDH Theory and Applications[END_REF][START_REF] Meyer | [END_REF] in the Heidelberg Multi Configuration Time Dependent Hartree (MCTDH) method 29,[START_REF]Multidimensional Quantum Dynamics: MCTDH Theory and Applications[END_REF][START_REF] Worth | The MCTDH Package[END_REF][START_REF] Meyer | [END_REF][36][37] on a PES constructed based on CCSD(T)-F12a/aug-cc-pVTZ-F12 energies.

In the present work, methods of quantum dynamics coupled to a PES fit to quantum chemistry calculations are used to produce vibrational eigenstates for deuterated forms of formic acid using the Heidelberg implementation of the MCTDH method. 29,[START_REF]Multidimensional Quantum Dynamics: MCTDH Theory and Applications[END_REF][START_REF] Worth | The MCTDH Package[END_REF][START_REF] Meyer | [END_REF][36][37] We use the same method as Richter and Carbonnière 5 , their PES, and the (block-) improved relaxation method [29][30][31][START_REF]Multidimensional Quantum Dynamics: MCTDH Theory and Applications[END_REF][START_REF] Meyer | [END_REF] implemented in the MCTDH Heidelberg package which proved to be successful in assigning vibrational states for formic acid (HCOOH), 5 including the recently confirmed re-assignment of the Fermi-coupled ν 5 and 2ν 9 states. 38

II. METHODS

We compute the vibrational eigenstate energies of the deuterated forms of formic acid, HCOOD, DCOOH, and DCOOD, using the PES of Richter and Carbonnière 5 (referred to as the RC-PES) with the (block-) improved relaxation method [29][30][31][START_REF]Multidimensional Quantum Dynamics: MCTDH Theory and Applications[END_REF][START_REF] Meyer | [END_REF] in the Heidelberg MCTDH software package. 29,[START_REF]Multidimensional Quantum Dynamics: MCTDH Theory and Applications[END_REF][START_REF] Worth | The MCTDH Package[END_REF][START_REF] Meyer | [END_REF][36][37] The PES, valid up to about 6000 cm -1 above the zero point energy of the trans isomer, was generated in internal valence coordinates using the AGAPES method 39,40 and was based upon fitting of CCSD(T)-F12a/aug-cc-pVTZ-F12 ab initio energies. The PES has been converged to yield approximately a.5 cm -1 threshold for the average vibrational error for levels up to 6000 cm -1 above the trans-isomer zero point energy for HCOOH. The kinetic energy operator was generated using the TANA program. 41,42 The exploitation of the efficiency of MCTDH requires the Hamiltonian, i.e., the kinetic energy operator (KEO) and PES, to be in sum of product (SOP) form. Both the KEO and PES are constructed in SOP form by design. General details of the methods have been provided previously, [39][40][41][42] while specific details for formic acid, are given in the vibrational study of HCOOH by Richter and Carbonnière. 5 The valence coordinates, used for the PES and vibrational state computations are depicted in Figure 1, where the hydrogen isotopes Deuterium (D) or Hydrogen (H) are represented by x H. The grids, mode combinations, number and type of primitive functions, as well as the number of single particle functions (SPFs) used in the block improved relaxation computations within MCTDH are given in Table I. In comparison with the previous study of HCOOH, 5 we used the full torsional τ 2 grid and an increased number of SPFs. Assignments of the vibrational eigenstate energies to the fundamental, overtone, or combinational bands 4 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. Mode combinations (r 1 , -cos(π -θ 3 )) (r 2 ,r 3 ,r 4 ) (cos(θ 1 ),cos(θ 2 ))

PLEASE CITE THIS ARTICLE AS

( were determined by visualisation of the reduced densities.

In addition to the MCTDH computations, we carried out density functional theory (DFT) B3LYP 43,44 /aug-cc-pVTZ 45,46 VPT2 47 computations as implemented in Gaussian 16 48 for fundamental transition frequencies of both isomers of all considered species; a smaller number of overtones and combination bands are also reported. Although showing limitations in the predictions, especially in the case of resonant states and at higher transition energies, the functional proved to be in reasonable agreement with the gas-phase measurements with a maximal deviation of 49 cm -1 but a RMSD of 25 cm -1 for the (anharmonic) fundamental frequencies of trans-HCOOH. Since the goal of the present work is not to test functionals for accurate prediction of vibrational frequencies of formic acid, as has been done previously, 7 but rather to provide accurate prediction based on a full dimensional PES, a diversity of functionals has not been considered.

5

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI:10.1063/1.5135571

III. RESULTS AND DISCUSSION

The accuracy of the RC-PES 5 for the computation of fundamental, overtone, and combination vibrational bands of trans-HCOOH was assessed in the original work with the computed fundamental transition frequencies deviating by less than 5 cm -1 from known gasphase experimental transition frequencies. Further gas-phase experimental measurements for the cis-isomer fundamentals reported by Meyer and Suhm 7 showed a deviation of 8 cm -1 at most with the reported theoretical values from the RC-PES. In comparison, the environmental effect from the gas phase to the matrix is inconsistent, i.e., the frequency shift depends on the vibrational state, and is observed to be as large as 27 cm -1 , which is of similar order of magnitude as the difference between cis and trans fundamental transition frequencies stressing the need for accurate gas-phase measurements, including of the deuterated species.

Our theoretical results include fundamental, overtone, and combination vibrational band energies along with experimental results when available for all the deuterated trans-formic acid species: DCOOD, HCOOD, and DCOOH. For comparison, anharmonic vibrational energies computed at the B3LYP/aug-cc-pVTZ level of theory are presented. We also provide the harmonic frequencies and optimized geometries in the Supplemental Material in Tables S1 andS2. Meyer and Suhm 7 tested a number of functionals and basis sets for HCOOH and their best results, as compared to experiment, used a higher level of theory B2PLYP-D3(BJ)/aug-cc-PVQZ 45,46,49 which involves a double hybrid functional with Grimme's dispersion correction (D3) in combination with Becke-Johnson (BJ) damping 50 . The method of our choice deviates at most by 5 cm -1 from the theoretical results of Meyer and Suhm 7 for trans-HCOOH. In general, the vibrational energies from the anharmonic vibrational perturbation theory computations (VPT2) are mostly in good agreement with experimental data although the maximal deviation is of 46 cm -1 for the ν 2 (ν CD ) stretching of DCOOD; the good agreement for the relatively low-energy states could be attributed to the small size of the formic acid molecule and its relatively high rigidity. In the following sections, we discuss the vibrational states for each isotopologue, where cis-formic acid fundamental transition frequencies are given at the end; we do not discuss overtones and combination bands for the cis-isomer as experimental measurements are not yet available, while the fundamental frequencies should be available in the near future. 13 6 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. ) with respect to the ground state energy as determined using block-improved relaxation within MCTDH on the RC-PES (MCTDH) as well as those computed using VPT2 at the B3LYP/aug-cc-pVTZ level of theory (VPT2). Mode assignment is given as ν=stretching, δ=bending, γ=rocking, def.=deformation.

PLEASE CITE THIS ARTICLE AS

A. Trans-DCOOD

All experimental fundamental band centers of trans-DCOOD measured in the gas phase have been reported, [51][52][53][54][55][56][57] and the most recent results are provided in Table II. Experimental analysis including rotation has been completed for the ν 1 and ν 4 fundamentals, 55,56 i.e., the OD and C-O stretches, respectively. Two high-energy overtone and combination bands were also reported. 58 The higher energy stretching overtone 2ν 1 (2ν OD ) energy was obtained using the improved relaxation procedure; the theoretical value of 5173.0 cm -1 is in good agreement with the experimental value 58 of 5181 cm -1 and supports the validity of the PES at relatively high energy. The ν 2 + ν 3 combination band (ν CD +ν C=O ) 58 observed at 3981 cm -1 could not be converged in the present work due to the high density of states in that energy region and the mixed nature of the state (which is a function of the choice of coordinate system).

Furthermore, molecular geometrical changes observed 59 in HCOOH following the multiple excitation of the ν 1 mode (ν OD ) could not be studied in this work due to the limited energy range validity of the PES used.

Theoretical results compared to experimental measurement for the fundamental bands are gathered in older nature of the concerned experimental results (from 1969) which perhaps could benefit from updated measurements. As expected, the anharmonic vibrational perturbation theory results behave more poorly at higher energy due to anharmonicity, but show an RMSD of 21 cm -1 with measured fundamental transition frequencies and a stunningly good agreement for the five lowest energy fundamentals. The relatively good agreement for the OD stretching (ν 1 ) with a deviation of 29 cm -1 in comparison to the maximal deviation of 46 cm -1 of the CD stretching (ν 2 ) could be attributed to the "local" character of the ν 1 mode. 26 All other assigned overtones and combination bands are given in the Supplementary Material in Tables S3 andS4.

B. Trans-HCOOD

All fundamental bands were recorded in the gas phase for trans-HCOOD 51,53,54,[60][61][62] and in rare-gas matrices; data for overtone or combination bands is very limited in the gas phase 58,[65][66][67] , see Table V, but more significant in rare-gas matrices thanks to Marushkevich et al., 12 see Table S5 in the Supplementary Material. From gas-phase measurements, the rotational structure of the ν 3 , ν 4 , ν 5 , ν 6 , ν 7 , and ν 9 fundamentals have been studied 57,[65][66][67][68][69][70][71] as well as in the 2ν 9 overtone band. 65 Theoretical results compared to gas phase experimental and B3LYP/aug-cc-pVTZ VPT2 anharmonic computations are gathered in Table III for the fundamental transition frequencies. For the fundamentals, the RMSD is 7 cm -1 and the maximum deviation is 19 cm -1 for ν 8 but as has already been pointed out in the DCOOD case, the concerned experimental results were obtained a relatively long time ago. Therefore, given the experimental advances in high resolution spectroscopy, these could benefit from being remeasured. Again, the anharmonic vibrational perturbation results behave more poorly at higher energy but show an RMSD of 18 cm -1 with measured fundamental transition frequencies.

In Table IV, we gather the gas-phase fundamental transition frequencies alongside the rare-gas, i.e., Ne, Ar, Kr, and Xe, matrices measurements of Domanskaya et al. 63 in terms of spectral shift from the gas phase transition frequencies. The corresponding absolute values are repeated in the Supplementary Material in Table S6. In general, the shifts between the gas phase to the matrix seem to be consistent, i.e., the same direction and similar magnitude, for the various rare gases given a single vibrational mode; however there are some notable exceptions, e.g., ν C-O where Ne and Xe lead to small negative shifts while 9 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. Ar and Kr have larger positive shifts. This analysis emphasizes the need for benchmark gas-phase measurements to compare with theoretical predictions, e.g., for the cis-isomer discussed in Section III D. Most rare gas matrix measurements overestimate the gas phase frequencies for the modes ν 1 (ν CH ), ν 4 (γ CH ), and ν 8 (τ HCO ) and underestimate for modes ν 2 (ν OD ), ν 3 (ν C=O ), ν 6 (δ COD ), and ν 9 (τ COD ). The effect of the environment in matrix isolation for the ν 7 (OCO-COD def.) mode seems to be the weakest with a spectral shift with the gas-phase of +1.4 and -1.2 cm -1 at maximum for Ar and Xe, respectively.

PLEASE CITE THIS ARTICLE AS

The reported gas-phase combination and overtone band frequencies within the scope of this study are given in Table V. Our theoretical results agree with the measurements at lower frequencies of 2ν 9 (2τ COD ), ν 5 + ν 9 (ν C-O +τ COD ), and

ν 5 + ν 7 (ν C-O +OCO-COD def.
) with a maximum deviation of about 8 cm -1 for the ν 5 + ν 7 combination band. At higher frequencies, theoretical results for ν 2 + ν 3 (ν OD +ν C=O ), ν 1 + ν 3 (ν CH +ν C=O ), and 2ν 2 (2ν OD ) behave more poorly although the agreement with experiment is still fairly good with a maximum deviation of about 11 cm -1 for 2ν 2 . However, the reported measurements 58 at higher transition (>4400 cm TABLE VI: For half-deuterated trans-formic acid (trans-DCOOH), fundamental vibrational energies (in cm -1 ) with respect to the ground state energy as determined using block-improved relaxation within MCTDH on the RC-PES (MCTDH) as well as those computing using VPT2 at the B3LYP/aug-cc-pVTZ level of theory (VPT2). Gas phase and matrix (Ar) isolation measurements are given for comparison. Mode assignment is given as ν=stretching, δ=bending, γ=rocking, and def.=deformation.

C. Trans-DCOOH

All fundamental bands were recorded in the gas phase for DCOOH 51,53,54,60,73,74 and some have also been reported in rare-gas matrices; 11,64 the gas phase measurements and matrix shifts with respect to the gas phase measurements are given in Table VI. Two overtones were also reported: 2ν 2 (2ν CD ) at ≈ 4360 cm -1 , 58 where the MCTDH theoretical value is 4385.0 cm -1 and 2ν 5 (2ν C-O ) at 2254.244±0.011 cm -1 , 75 with a computed value of 2272.1 cm -1 .

From gas-phase measurements, the rotational structures of the ν 2 , ν 3 , ν 5 , ν 6 , ν 7 , ν 8 , and ν 9 fundamental bands as well as the 2ν 5 overtone have been studied. 73,[75][76][77][78][79] Theoretical results compared to gas phase experimental measurements and B3LYP/aug-cc-pVTZ anharmonic computations are gathered in Table VI. For the fundamentals, the MCTDH results have an RMSD of 3 cm -1 with a maximum deviation of about 4 cm -1 for mode ν 9 . As already pointed out for other isotopologues, the anharmonic vibrational perturbation theory results

behave more poorly at higher energy and, show an RMSD of 23 cm -1 when compared with measured fundamental transition frequencies.

In Table VI, we gather the gas-phase fundamental transition frequencies alongside the rare-gas (Ar) matrix measurements of Maçôas et al. Table S10. Shifts between the gas phase to the matrix follow the same trends as for the modes in HCOOD only for the γ CH (or γ CD ) rocking and ν OD (or ν OH ) stretching modes. These observations stress again 7 the need for gas phase measurements to compare to theoretical values which are typically computed in the gas phase. A Fermi resonance between the ν 3 stretch (ν C=O ) and the 2ν 8 (2τ DCO ) overtone in rare-gas matrices was reported by Maçôas et al. 11 for trans-DCOOH around 1723 and 1760 cm -1 , all values are given in Table S10. The Fermi resonance was captured in our computations as it can clearly be seen from the reduced density functions given in Figure 2, where the corresponding energies of the vibrational bands are 1722.8 and 1760.9 cm -1 , respectively. We assigned the fundamental ν 3 (ν C=O ) at 1722.8

cm -1 given the experimental ν C=O stretch value reported 73 at 1725.874974(34) cm -1 . All other assigned overtones and combination bands are given in the Supplementary Material in Tables S11 andS12.

D. Cis isomer

Formic acid is one of the simplest systems that exhibits cis-trans isomerism. The first observation of the thermally less abundant cis form was achieved with rotational spectroscopy by Hocking. 8 Gas-phase infrared spectra of the cis form of formic acid are scarce, [7][8][9][10] and to the best of our knowledge, not published for its deuterated forms. However, the cis isomer was isolated in rare-gas matrices in experimental studies on HCOOD 12,63 and DCOOH. 11,64 The block improved relaxation method implemented in MCTDH, coupled to the RC-PES proved again 5 to be successful in the estimation of vibrational frequencies of formic acid in 12

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. the case of its trans deuterated forms in this work. For this reason, we provide our theoretical values for the cis isomer using MCTDH on the RC-PES as well as anharmonic vibrational energies at the B3LYP/aug-cc-pVTZ level of theory in Tables VII, VIII, and IX for the fundamental vibrational frequencies of DCOOD, HCOOD, and DCOOH, respectively. We note that the vibrational frequencies reported by Tew and Mizukami 6 for the undeuterated HCOOH cis-isomer differed quite significantly from those determined using MCTDH on the RC-PES. 5 However, careful analysis by Richter ) with respect to the (cis) ground state energy as determined using block-improved relaxation within MCTDH on the RC-PES (MCTDH) as well as those computed using VPT2 at the B3LYP/aug-cc-pVTZ level of theory (VPT2). Mode assignment is given as ν=stretching, δ=bending, γ=rocking, and def.=deformation.

PLEASE CITE THIS ARTICLE AS

Available rare-gas matrices measurements for the cis isomer agree with the theoretical results within the maximum spectral shifts in comparison with the gas-phase measurements observed in the case of trans-formic acid; the matrix-induced shifts for trans-formic acid are given in Tables IV andVI. The present MCTDH prediction for the ν 9 fundamental of HCOOD (374.8 cm -1 ) is in very good agreement with the value (378.9 cm -1 ) predicted recently using a 1-dimensional torsional potential computed at the CCSD(T) level of theory extrapolated to the complete basis set limit. 38 A Fermi resonance between the ν C=O stretch and the 2τ DCO overtone in rare-gas matrices was reported by Maçôas et al. ) with respect to the (cis) ground state energy as determined using block-improved relaxation within MCTDH on the RC-PES (MCTDH) as well as those computed using VPT2 at the B3LYP/aug-cc-pVTZ level of theory (VPT2). Experimental values obtained from rare-gas matrix isolation are distinguished by the rare-gas used, Ne, Ar, Kr, and Xe; some bands split into several peaks so more than 1 value is provided. Mode assignment is given as ν=stretching, δ=bending, γ=rocking, and def.=deformation.

around 1720 and 1777 cm -1 . That Fermi resonance was not observed in our gas-phase computations, nor was it detected in the MP2 computations of Maçôas et al.. 11 However, as stressed in the preceding sections and by Meyer and Suhm, 7 theoretical estimations of the band positions need to be compared to gas-phase measurements due to inconsistency of the spectral shifts induced by the rare-gas matrix. For this reason, we provide the theoretical results in the hope that it will help the experimentalists in the location and assignment of fundamental vibrational modes belonging to half and fully deuterated cisisomers. Since gas-phase measurements on the cis isomer are challenging 7 and as such no direct comparison between theory and experiment can yet be made, we have not provided extensive lists of overtone and combination bands for the cis-isomer.

14 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. ) with respect to the respective (cis) ground state energy as determined using block-improved relaxation within MCTDH on the RC-PES (MCTDH) as well as those computed using VPT2 at the B3LYP/aug-cc-pVTZ level of theory (VPT2). Experimental values were obtained from raregas matrix isolation using Ar. Mode assignment is given as ν=stretching, δ=bending, γ=rocking, and def.=deformation.

PLEASE CITE THIS

IV. CONCLUSIONS

The PES of Richter and Carbonnière 5 coupled to the block improved relaxation method implemented in the MCTDH software package was used to determine vibrational eigenenergies of both trans and cis deuterated formic acid. All fundamental vibrational energies as well as overtones and combination bands were computed up to approximately 3000 cm -1

(relative to their respective ZPE) for trans-DCOOD, HCOOD, and DCOOH. We report an RMSD of 8, 7, and 3 cm -1 from estimations using MCTDH with the fundamental transition frequencies measured in gas-phase experiments for DCOOD, HCOOD, and DCOOH, respectively. Perhaps then, not surprisingly, excellent agreement is obtained between theory and experiment for the limited number of overtone and combination bands for the trans-isomers.

For some states with the largest discrepancies, we suggest that new experimental measurements may be warranted. Moreover, the present comparisons reinforce the conclusion 7 that for benchmarking theoretically determined results must be compared to gas-phase measurements rather than those from rare-gas matrices as shifts are strongly vibrational state dependent. Fundamental vibrational eigenenergies for all deuterated forms of the cis isomer are provided in the hope that they will spur the first experimental gas phase measurements

  FIG.1: Valence coordinates used to describe formic acid; the torsional angles τ 1 and τ 2 are respectively the x H1 and x H2 (or D1/D2 depending on deuterium substitution) out of plane motions relative to the O2CO1 body.

  11 in terms of spectral shift with the gas phase measurements. The absolute values are repeated in the Supplementary Material in 11 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE ASDOI:10.1063/1.5135571

FIG. 2 :

 2 FIG. 2: Contour plot of the reduced density functions in the {r 2 , r 3 } and {τ 1 , τ 2 } planes of the ν 3 (ν C=O ), 2ν8 (2τ HCO ) Fermi resonance for trans-DCOOH.

  DOI:10.1063/1.5135571

TABLE I :

 I Specifications of the parameters used for the MCTDH (block-)improved relaxation computations. HO refers to the harmonic oscillator basis. Bond distances are given in Bohr and the dihedral angles are given in radians.

DOI:10.1063/1.5135571

	Assignment Mode	Gas phase a	MCTDH	VPT2
	τ COD	ν 9 493.2254252(35) 57 489.4	492.3
	τ DCO	ν 8	873.0 54	871.5	872.8
	OCO-COD def. ν 7 554.4394726(50) 57 552.0	549.8
	δ COD	ν 6	945.0 54	943.7	939.2
	γ CD	ν 5	1042 52	1035.1	1041.4
	ν C-O	ν 4	1170.79980(2) 56 1170.4	1150.8
	ν C=O	ν 3	1742.0 54	1759.0	1762.9
	ν CD	ν 2	2213.0 54	2228.1	2167.1
	ν OD	ν 1 2631.87379(43) 55 2629.4 2602.7, 2569.9 b
		RMSD		8	21

a Error on the last digit(s) is indicated in brackets when published. b DFTB

25

TABLE II :

 II For fully-deuterated trans-formic acid (trans-DCOOD), fundamental vibrational energies (in cm -1

PLEASE CITE THIS ARTICLE AS DOI:10.1063/1.5135571

  This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. Error on the last digit(s) is indicated in brackets when published.b Value reported in60 , original paper unavailable(1973)

	Assignment Mode	Gas phase a	MCTDH VPT2
	τ COD	ν 9 508.1320569(45) 57 505.2 508.7
	τ HCO	ν 8	1010.8 72b	1029.4 1029.9
	OCO-COD def. ν 7 558.2722502(62) 57 555.8 553.7
	δ COD	ν 6	972.8520(1) 65,66	968.6 968.1
	ν C-O	ν 5	1177.09378(2) 69 1179.1 1156.6
	γ CH	ν 4 1366.48430(39) 68 1361.9 1360.5
	ν C=O	ν 3	1772.119(24) 67	1770.5 1775.2
	ν OD	ν 2 2631.63835(17) 62 2629.4 2602.8
	ν CH	ν 1	2938-2942 60	2934.4 2901.3
		RMSD		7	18

Table

II

with a root mean square deviation (RMSD) of 8 cm -1 and largest deviation of 17 cm -1 for the ν C=O and 15 cm -1 for the ν CD modes. However, we point out the 7 a

TABLE III :

 III For half-deuterated trans-formic acid (trans-HCOOD), fundamental vibrational energies (in cm -1 ) with respect to the respective ground state energy as determined using block-improved relaxation within MCTDH on the RC-PES (MCTDH) as well as those computed using VPT2 at the B3LYP/aug-cc-pVTZ level of theory (VPT2). Mode assignment is given as ν=stretching, δ=bending, γ=rocking, and def.=deformation.

  63,64 they are provided in TableIIIand matrix shifts with respect to the gas phase measurements are given in TableIV. We stress in this case that the ν 1 and ν 2 fundamental modes are "switched" compared to the HCOOH case, i.e., they represent the

	Assignment Mode	Gas phase a	Ne 63	Ar 63	Kr 63	Xe 63
	τ COD	ν 9 508.1320569(45) 57	-1.2	-1.4/-2.9 -1.1/-3.4/-4.5 -3.0/-5.6
	τ HCO	ν 8	1010.8 72	+26.6 64	+25.6/+24.3	+22.8	+19.6
	OCO-COD def. ν 7 558.2722502(62) 57	-0.4	+1.4/+0.4	+1.1/-0.1 +0.0/-1.2
	δ COD	ν 6	972.8520(1) 65,66	-1.7	-1.8	-1.2/-4.7/-5.4 -4.1/-8.1
	ν C-O	ν 5	1177.09378(2) 69	-0.2	+4.8	+9.1/+2.7 -0.8/-1.2
	γ CH	ν 4 1366.48430(39) 68	+1.6 64	+7.8/+4.6	+1.9	-1.7
	ν C=O	ν 3	1772.119(24) 67	-1.1	-4.9/-6.5	-9.9	-16.4/-15.7
	ν OD	ν 2 2631.63835(17) 62	-0.1	-12.1/-13.6	-21.6/-26.8 -31.9/-38.3
	ν CH	ν 1	2938∼2942 60	+9.5∼+13.5 +19.2∼+23.2 +11.8∼+15.8 +3.6∼-0.4

CH and OD stretching in HCOOD, respectively. For this reason, the equivalent movement is assigned alongside the mode number throughout this paper. The available experimental 8 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI:10.1063/1.5135571

a Error on the last digit(s) is indicated in brackets when published.

TABLE IV

 IV

	: Half-deuterated trans-formic acid (trans-HCOOD) fundamental vibrational energies
	(in cm -1 ) compared to matrix isolation values in different rare gas (RG) matrices; valus given as
	δν i = ν gas-phase i	-ν RG i . Mode assignment is given as ν=stretching, δ=bending, γ=rocking, and
	def.=deformation.

DOI:10.1063/1.5135571

  +τ COD ν 5 + ν 9 1680.962939(64) 67 1681.4 1662.7 ν C-O +OCO-COD def. ν 5 + ν 7 1732.084(24) 67

	Assignment	Mode	Gas phase a	MCTDH	VPT2
	2τ COD	2ν 9 1011.6766(1) 65,66 1006.2	1006.2
	ν C-O 1729.9	1705.7
	ν OD +ν C=O	ν 2 + ν 3	4404 58	4397.3 4376.0, 4345.8 b
	ν CH +ν C=O	ν 1 + ν 3	4702-4736 58	4715.3 4670.5, 4704.3 b
	2ν OD	2ν 2	5181 58	5169.7 5118.1, 5163.6 b
	a Error on the last digit(s) is indicated in brackets when published.
	b CC-VSCF 11 based on MP2/6-311++G (2d,2dp)

TABLE V :

 V For half-deuterated trans-formic acid (trans-HCOOD), overtone and combination bands energies (in cm -1 ) with respect to the respective ground state energy as determined using blockimproved relation within MCTDH on the RC-PES (MCTDH) as well as those computed using VPT2 at the B3LYP/aug-cc-pVTZ level of theory (VPT2) and CC-VSCF results of Maçôas et al..11 Mode assignment is given as ν=stretching, δ=bending, γ=rocking, and def.=deformation.

PLEASE CITE THIS ARTICLE AS DOI:10.1063/1.5135571

  -1 ) frequencies might have a large uncertainty and could benefit from updated measurements. Reported transitions in matrices are given in the Supplementary Material in TableS5and all other assigned overtones and combination bands are given in the Supplementary Material in Tables S7 to S9. 10 This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. Error on the last digit(s) is indicated in brackets when published.

	Assignment Mode	Gas phase a	Ar 11	MCTDH VPT2
	τ COH	ν 9 631.5437158(56) 75	-5.1/-5.6	627.7 629.0
	τ DCO	ν 8	870 74	+4.8/+5.0	871.9 872.9
	OCO-COH def. ν 7 620.5683857(57) 75	+1.7/+2.9	617.3 616.2
	δ COH	ν 6	-	-	1294.1 1283.5
	ν C-O	ν 5	1142.31075(2) 77	-0.2/+2.1/-0.6	1139.3 1125.4
	γ CD	ν 4	970 60	+6.2/+4.9/+4.1	970.9 965.7
	ν C=O	ν 3 1725.874974(34) 73	-2/-3	1722.8 1763.8
	ν CD	ν 2	2219.6896(2) 79	+5.5	2216.5 2183.0
	ν OH	ν 1	3566 60	-3.5/-16/-18/-5.1/-14 3567.8 3531.4
			RMSD		3	23

a

  and Carbonniére, see Section III.B.2 of Ref.5, suggested that the source of this discrepancy was not found in the underlying PESs nor the level of ab initio theory upon which they were based, but rather due to shortcomings in the VCI procedure at the high energies needed to obtain the cis vibrational states.

	Assignment Mode MCTDH	VPT2
		ν 0 1429.518 a
	τ COD	ν 9	369.8	376.0
	τ DCO	ν 8	859.7	863.7
	OCO-COD def. ν 7	615.1	616.9
	δ COD	ν 6	883.8	872.4
	γ CD	ν 5	1029.4	1023.1
	ν C-O	ν 4	1159.5	1141.8
	ν C=O	ν 3	1782.0	1792.5
	ν CD	ν 2	2142.6	2111.1
	ν OD	ν 1	2682.7 2657.6, 2605.9 b
	a with respect to the trans ground state energy
	b DFTB 25

TABLE VII :

 VII For fully-deuterated cis-formic acid (cis-DCOOD), fundamental vibrational energies (in cm -1

PLEASE CITE THIS ARTICLE AS DOI:10.1063/1.5135571

  This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.Assignment Mode Ne 63 Ar 63 Kr 63 Xe 63 MCTDH VPT2

		ν 0	-	-	-	-	1412.023 a
	τ COD	ν 9	-	-	-	-	374.8	381.7
	τ HCO	ν 8 1039.3 1041	-	-	1013.3 1017.3
	OCO-COD def. ν 7 629.4 633.0 633.0 632.2	622.5	624.4
						631.6
	δ COD	ν 6 906.8 911.2 911.7 911.6	905.3	892.1
				910.5 909.0 909.2
	ν C-O	ν 5 1159.0 1164.0 1163.3 1159.8 1161.8 1137.3
	γ CH	ν 4 1388.0 1395.4 1388.6 1385.1 1382.7 1381.2
			1394.5 1391.0
	ν C=O	ν 3 1818.2 1821.0 1821.1 1819.6 1814.8 1823.5
			1798.9 1799.3 1796.4 1791.1
			1781.2 1788.9 1785.8 1781.3
	ν OD	ν 2 2872.7 2895.3 2888.4 2876.2 2869.5 2849.3
	ν CH	ν 1 2685.5 2668.7 2656.7 2641.3 2683.4 2657.6

11

for cis-DCOOH 13 a with respect to the trans ground state energy

TABLE VIII :

 VIII For half-deuterated cis-formic acid (cis-HCOOD), fundamental vibrational energies (in cm -1

ARTICLE AS DOI:10.1063/1.5135571

  with respect to the trans ground state energy. b In Fermi resonance with a DCO wagging mode.

	Assignment Mode	Ar 11	MCTDH VPT2
			ν 0	-	1415.469 a
	τ COH	ν 9	500.7-498.1	490.9	493.0
	τ DCO	ν 8	-	861.0	864.6
	OCO-COH def. ν 7	655.2	652.0	647.3
	δ COH	ν 6 1248.9/1234.0/1229.2/1239.4 1234.7 1119.9
	ν C-O		ν 5	1143.4	1130.0 1216.8
	γ CD		ν 4	961.9/962.9	985.6	965.1
	ν C=O	b	ν 3	1777.2/1777.4/1778.7	1781.0 1792.3
				1721.5/1720.9
	ν CD		ν 2	2183.5/2179.7	2166.9 2138.8
	ν OH		ν 1	3615.4/3617.4	3624.8 3599.7

a

TABLE IX :

 IX For half-deuterated cis-formic acid (cis-DCOOH), fundamental vibrational energies (in cm -1

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE ASDOI:10.1063/1.5135571

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

Acknowledgments

The research was partly funded by the Natural Sciences and Engineering Research Council of Canada (NSERC Discovery Grant to A.B.). A.A. acknowledges the financial support of the Commission de Classement des Crédits Internationaux (CCCI) from ULB and the Lucia DeBrouckere fund from the chemistry department of ULB. Most of the computational resources have been provided by the Shared ICT Services Centre, Université Libre de Bruxelles. Test runs were performed on the computers of the Consortium des Équipements de

and help the experimentalists with assignments.

V. SUPPLEMENTARY MATERIAL

See Supplementary Material for the B3LYP/aug-cc-pVTZ harmonic fundamental transition frequencies of the deuterated cis-and trans-isomers of formic acid (DCOOD, HCOOD, and DCOOH); optimized B3LYP/aug-cc-pVTZ equilibrium geometries of the trans-and cis-isomers; MCTDH overtone and combination band energies for fully-deuterated transformic acid (trans-DCOOD); experimental (in Ne and Ar), MCTDH, and CC-VSCF overtone and combination band energies for half-deuterated trans-formic acid (trans-HCOOD); fundamental vibrational transition frequencies of half deuterated trans-formic acid (trans-HCOOD) as measured in the gas-phase and using rare-gas matrix isolation (Ne, Ar, Kr, and Xe); MCTDH overtone and combination band energies for half-deuterated trans-formic acid (trans-HCOOD); fundamental vibrational transition frequencies of half deuterated transformic acid (trans-DCOOH) as measured in the gas-phase and using rare-gas matrix isolation (Ar); MCTDH overtone and combination band energies for half-deuterated trans-formic acid (trans-DCOOH).

Calcul Intensif (CECI) on Vega (ULB) and Dragon2 (UMons).

1 P. Khare, N. Kumar, K. Kumari, and S. Srivastava, Rev. Geophys. 37, 227 (1999).

2 N. Kitadai and S. Maruyama, Geosci. Front. 9, 1117 (2018).


est destinée au dépôt

  Specific blockade of morphine-and cocaine-induced reinforcing effects in conditioned place preference by nitrous oxide in mice. Nadia Benturquia, Stéphanie Le Guen, Corinne Canestrelli, Vincent Lagente, Gabriella Apiou, Bernard P. Roques, Florence Noble To cite this version: Nadia Benturquia, Stéphanie Le Guen, Corinne Canestrelli, Vincent Lagente, Gabriella Apiou, et al.. Specific blockade of morphine-and cocaine-induced reinforcing effects in conditioned place preference by nitrous oxide in mice.. Neuroscience, 2007, 149 (3), pp.477-86. 10.1016/j.neuroscience.2007.08.003. inserm-


Juan Carlos Gonzalez-Torres
email: juan-carlos.gonzalez@supergrid-institute.com

Gilney Damm

Valentin Costan

Abdelkrim Benchaib

Françoise Lamnabhi-Lagarrigue

Transient stability of power systems with embedded VSC-HVDC links: Stability margins analysis and Control

This paper investigates the impact of embedded Voltage Source Converter-based High Voltage Direct Current (VSC-HVDC) links on AC grids transient stability. Firstly, using Transient Energy Functions (TEF), it is demonstrated that VSC-HVDC links controlled to track constant power references, do not inherently improve transient stability of the surrounding AC grid as an AC line naturally does. Then, a control law using the feedback linearization technique on a simple but representative power system is derived. The control law highlights and combines the three main actions the VSC-HVDC link can offer to enhance rotor angle stability: fast power reallocation, injection of synchronising power and injection of damping power. The control law is implemented and validated in EMT simulation. It is then shown that an HVDC link can assure the synchronisation of different AC areas even if no AC transmission lines interconnect them. Through another case study, it is shown how the HVDC link can help to share dynamic frequency reserves in order to not jeopardise the stability of the system. A last example investigates the effect of a DC fault on AC transient stability and how the control can help improving the system response.

Introduction

The structure and operation of electric transmission systems is evolving faster than ever. The introduction of renewable energies in centralised facilities often located far from end-users, as well as the disconnection of conventional generation units, pose new challenges for system operators. Because of these trends, the transport of bulk-power over long distances has become a challenging necessity. Indeed, if large amounts of power are transported over large AC corridors, the security of the system can be at risk [START_REF] Kundur | Power system stability and control[END_REF]. A possible solution to assure the system security under these conditions is the reinforcement of the transmission system. To this end, the High Voltage Direct Current (HVDC) technology has raised not only as a cost-effective solution for long distance transmission but also presents other technical advantages over AC lines. When a system is reinforced with an HVDC line based on Voltage Source Converters (VSC-HVDC), the fast and independent control of active and reactive power gives the possibility to use the VSC-HVDC link to enhance the stability of the surrounding AC grid. In fact, dynamic stability studies have been performed over many decades showing that power modulation of HVDC links can enhance rotor angle stability. The injection of active and reactive power by HVDC systems allows to damp power oscillations as well as to counteract power imbalances, which can be the source of transient stability problems in a multi-machine system when a disturbance occurs. Actually, as defined in [START_REF] Entso | ENTSO-E Draft Network Code on High Voltage Direct Current Connections and DC connected Power Park Modules[END_REF], it is already expected for the HVDC systems to participate on the stability enhancement of the surrounding system. This can be done by means of supplementary controllers such as the one proposed in this paper.

Power system transient stability stands for the ability of the system to keep synchronism after being subjected to a large disturbance such as a three-phase fault [START_REF] Kundur | Application of power system stabilizers for enhancement of overall system stability[END_REF]. Its assessment requires the analysis of the systems described in their non-linear nature, which can be complex, especially when there is a large number of dynamic elements. The main tasks of the Transient Stability Assessment (TSA) consist in:

• Providing transient stability limits: approximations of the Region Of Attraction (ROA) of the power system, i.e. how far the system can move away from its equilibrium point without losing synchronism.

• Providing transient stability margins: indicators of how close a system is to the transient stability limit following a given disturbance.

In power system engineering, the most common practice is to quantify these limits and margins in terms of the Critical Clearing Time (CCT) for a given disturbance. The CCT can be computed by iterative numerical solution of the systems' equations. Alternatively, the so-called Direct Methods, discussed in the classical literature of power systems [START_REF] Pai | Energy function analysis for power system stability[END_REF][START_REF] Machowski | Power system stability enhancement by WAMS-based supplementary control of multi-terminal HVDC networks[END_REF][START_REF] Padiyar | Structure Preserving Energy Functions in Power Systems: Theory and Applications[END_REF][START_REF] Anderson | Power system control and stability[END_REF], give the possibility to analytically approximate the region of attraction of a system. These methods are based on the use of the Transient Energy Functions (TEF), and they aim to provide information about the stability limits and margins in terms of energy. Although conservative, these methods can assess stability limits independently from the experienced disturbance and without solution of the equations (which is not the case for the CCT). Therefore, to have more general conclusions, Direct Methods are used in this paper to compare the stability margins of an AC transmission system against its equivalent AC/DC system. The comparison of transient stability margins is used to highlight the need of HVDC power modulation through supplementary controls in AC/DC grids.

The study of supplementary controllers allowing to modulate the control variables of HVDC links for AC stability enhancement started almost at the same time as the HVDC technology itself. A survey on controllers for the particular use of power oscillation damping can be found in [START_REF] Elizondo | Interarea oscillation damping control using high-voltage dc transmission: A survey[END_REF]. Supplementary controllers of HVDC links found in literature mainly differ depending on:

• HVDC technology (LCC or VSC),

• measured feedback signals (phase, voltages, etc...),

• control objectives (damp oscillations, security improvement, voltage stabilisation, etc.. ),

• control input (active or reactive power, DC voltage, AC or DC current, etc..),

• applied control technique, • the use (or not) of Wide Area Measurement Systems (WAMS),

• coordination of the HVDC link with other FACTS devices, • topology of the AC/DC grid (embedded HVDC links or links connecting two asynchronous areas) For decades, the transient stability problem has been treated via the control of the excitation of generators using the Automatic Voltage Regulators (AVR) and Power System Stabilisers (PSS). For this reason, in early development of point to point HVDC links, a natural choice was to equip them with PSSs aiming the stability of the surrounding AC grid [START_REF] Cresap | Small-signal modulation of the Pacific HVDC intertie[END_REF]. Another intuitive control scheme consisted in modulating the HVDC power proportionally to the speed deviation of a nearby machine [START_REF] Pai | Transient stability analysis of multimachine AC/DC power systems via energy-function method[END_REF]. This early development led to the implementation of controllers on existing projects based on measurements done near one HVDC station [START_REF] Cresap | Small-signal modulation of the Pacific HVDC intertie[END_REF][START_REF] Hao | Design of modulation controller to damp power oscillations of parallel AC line in the Tianshenqiao to Guangdong HVDC transmission[END_REF].

The introduction of WAMS allowed the implementation of control schemes using feedback signals measured at different locations of the power grid. For example, in [START_REF] Latorre | Application of Control Lyapunov Functions to Voltage Source Converters-based High Voltage Direct Current for Improving Transient Stability[END_REF][START_REF] Latorre | Active and reactive power control of a VSC-HVDC[END_REF] using Lyapunov's theory, it was proven that transient stability is enhanced by modulating active power proportionally to the difference of frequency measured at the PCC of both converters (also analysed in [START_REF] Sigrist | A fundamental study on the impact of HVDC lines on transient stability of power systems[END_REF][START_REF] Gonzalez-Torres | Power system stability enhancement via vsc-hvdc control using remote signals: Application on the nordic 44-bus test system[END_REF]). This control scheme was recently implemented in the western North American Power System on a real time simulation platform with actual measurements from PMUs [START_REF] Schoenwald | Design and implementation of a wide-area damping controller using high voltage DC modulation and synchrophasor feedback[END_REF]. This results have been upgraded and extrapolated for the proposition of supplementary controllers of Multi-Terminal DC grids in [START_REF] Renedo | Active power control strategies for transient stability enhancement of AC/DC grids with VSC-HVDC multi-terminal systems[END_REF]. Regarding the use of voltage phase angles, two industrial applications have been reported to modulate the power as a proportional function of the angle difference between inverter and rectifier sides [START_REF] Henry | Influence of embedded HVDC transmission on system security and AC network performance[END_REF][START_REF] Francos | INELFE-Europe's first integrated onshore HVDC interconnection[END_REF]. Although this kind of control has been tested on real applications, no dynamic analysis has been done. Also, regarding the use of angles as feedback signal, recently a PSSlike control using angle difference as feedback signal was proposed in [20] for designing an oscillation damping control.

The implementation of real WAMS in different power systems [START_REF] Aminifar | Synchrophasor measurement technology in power systems: Panorama and stateof-the-art[END_REF] also gave value to the research done regarding advanced controllers, with more feedback signals and more control inputs. For instance, as for the coordination of different HVDC links, a centralised scheme using a model predictive control was presented in [START_REF] Fuchs | Stabilization of large power systems using VSC-HVDC and model predictive control[END_REF]. There, it is proposed to use speed measurements of different machines. A compilation of other results on supplementary controllers of HVDC links using WAMS can be found in [START_REF] Li | Interconnected Power Systems: Wide-Area Dynamic Monitoring and Control Applications[END_REF].

Recently, promising tracks using the inner loops of the converter for improving transient stability have been proposed. In [START_REF] Aouini | Synchronverter-based emulation and control of HVDC transmission[END_REF], the authors show that controlling the inverter (rectifier) station in order to emulate the behaviour of a synchronous motor (generator) has a favourable impact on transient stability. Furthermore, in [START_REF] Arioua | Robust grid-oriented control of high voltage DC links embedded in an AC transmission system[END_REF][START_REF] Arioua | Multivariable control with grid objectives of an HVDC link embedded in a large-scale AC grid[END_REF] the authors design robust inner controllers of the converters taking into account the dynamics of the surrounding AC grid. Such control techniques can be achieved only based on local measurements.

In the cited works above, the advantage of synchronising power injections (proportional to phase angle deviations) is not highlighted. To the authors' best knowledge, it has not been shown how the synchronisation of different regions is possible through the control of HVDC links even when no AC links are connecting those grids. The present work brings answers to these points, and in addition, analyses two different events poorly reported in the literature. The first one is the uncoordinated poor fast frequency response, and the second is the study of the impact of DC faults on the overall system stability.

The present paper is organised as follows. In section 2, a modelling technique allowing to represent mixed AC/DC systems as a set of ordinary differential equations is presented. In section 3, a comparison between an AC system and an equivalent AC/DC system is performed using Direct Methods. This comparison shows that HVDC systems with constant power references do not intrinsically increase transient stability margins as an AC line does. Therefore, the need of a supplementary control loop for the HVDC link to enhance the AC stability is highlighted. In section 4, feedback linearization technique is used to derive a control law, and its physical meaning is studied. The proposed control law is based on three actions: damping power injection, synchronising power injection and feed-forward actions. In section 5, EMT simulations using detailed models are performed. The advantages of the proposed control are emphasised, through the analysis of different contingencies, such as three phase-faults, uncoordinated fast frequency response and DC faults.

AC/DC Power systems modelling for transient stability studies

In this section, the Internal Node Model (INM) of a multimachine/multi-converter power system is described. The INM allows to represent the system by a set of ordinary differential equations through reduction of the admittance matrix of the system [START_REF] Ghandhari | Control Lyapunov functions for controllable series devices[END_REF][START_REF] Fan | Coordinated control of dc grid and offshore wind farms to improve rotor-angle stability[END_REF]. For the derivation of the model, consider power system composed of n generators and m HVDC stations interconnected via a transmission network described by the nodal admittance matrix Y.

δi = ω i M i ωi = -K D,i ω i + P m,i -P e,i (1)
where for the i -th generator (∀i = 1, ... , n), δ i represents the rotor angle deviation with respect to a synchronous reference, ω i is the rotor speed deviation with respect to the synchronous speed ωs, M i is the inertia coefficient, P m,i is the mechanical input power, P e,i is the electrical output power. In the classical model, the angle of the internal voltage E d,i ∠δ i represents also the rotor angle deviation δ i of a synchronous generator with one pair of poles.

Regarding the converters, their dynamics are much faster when compared with the electromechanical dynamics of the generators. For the dynamics of interest, converters can be considered as complex current sources ( Īhvdc,k ) with independent active power (P hvdc,k ) and reactive power (Q hvdc,k ) injections. This assumption is particularly valid when inner controllers designed for decoupling the AC side and the DC side are used [START_REF] Zama | High dynamics control for MMC based on exact discrete-time model with experimental validation[END_REF][START_REF] Chen | Control induced explicit time-scale separation to attain DC voltage stability for a VSC-HVDC terminal[END_REF][START_REF] Chen | Control-induced time-scale separation for multiterminal high-voltage direct current systems using droop control[END_REF][START_REF] González-Torres | State feedback control assisted by a gain scheduling scheme for three-level NPC VSC-HVDC transmission systems[END_REF][START_REF] González-Torres | Long-length horizons dynamic matrix predictive control for a MMC inverter[END_REF]. The controlled active and reactive power injections can be divided in two components as follows:

P hvdc,k = P hvdc0,k + ∆P hvdc,k Q hvdc,k = Q hvdc0,k + ∆Q hvdc,k (2)
where P hvdc0,k and Q hvdc0,k are the set-points determined by the operator of the system and ∆P hvdc,k and ∆Q hvdc,k are the supplementary control references of the converters.

The impact of the active and reactive power injections on the generators' dynamics is given by the development of the electrical output power term (P ei ) in Equation [START_REF] Kundur | Power system stability and control[END_REF]. Using Kirchoff's law, we have that the relation between currents and voltages through the admittance matrix including the transient reactances of the synchronous generators Y bus is given by:

Ig I hvdc = Y bus Y A Y B Y C Y D E d V hvdc (3)
where

Y A ,Y B , Y C and Y D are of dimensions (n × n), (n × m), (m × n) and (m × m) respectively.
Ig and E d respectively contain the injected currents and internal voltages of generators. I hvdc and V hvdc are the injected currents and voltages at AC bus-bars of the converters [START_REF] Pai | Transient stability analysis of multimachine AC/DC power systems via energy-function method[END_REF]. We can solve for the generator's currents as follows:

Ig = Y SN E d + K hvdc I hvdc (4)
with

Y SN = (Y A -Y B Y -1 D Y C ) and K hvdc = Y B Y -1 D .
Matrix Y SN is the matrix relating all generators voltages and currents, K hvdc describes the effect of HVDC injected currents on the generator. Elements of Y SN and K hvdc are:

ȲSNij = G ij + jB ij and Khvdc,ik = D hvdc ik + jF hvdc ik (5)
G ij and B ij are the equivalent conductance and susceptance between generators i and j. D hvdc ik and F hvdc ik are the real and imaginary parts of the element K hvdc,ik which denote the effect of the injected current of the k -th converter on the i -th generator. Then, the expression of the active power supplied by the i -th machine P e,i in Equation (1), is developed as:

P ei =   Ē di n j=1 Ȳ * SNij Ē * dj -Ē di m k=1 K * hvdc ik Ī * hvdc k   (6)
Substituting ( 5) into (6), the electric power supplied by the i -th generator is derived as:

P e,i = n j=1 E d,i E d,j (G ij cos δ ij + B ij sin δ ij ) + m k=1 Ē d,i Vhvdc,k (D hvdc ik + jF hvdc ik )(P hvdc,k + jQ hvdc,k ) (7) with δ ij = δ i -δ j .
The complex ratio of voltages ( Ē di / Vhvdc k ) can be approximated to constant values [START_REF] Pai | Transient stability analysis of multimachine AC/DC power systems via energy-function method[END_REF]. Following this assumption, and applying Equations (2) into Equation ( 7), the electrical power can be expressed as follows:

P e,i = P gg,i + ∆P gv,i (8)
with

P gg,i = n j=1 E d,i E d,j (G ij cos δ ij + B ij sin δ ij ) + m k=1 a hvdc ik P hvdc0,k + m k=1 b hvdc ik Q hvdc0,k ∆P gv,i = m k=1 a hvdc ik ∆P hvdc,k + m k=1 b hvdc ik ∆Q hvdc,k (9)
where a hvdc ik and b hvdc ik are called "distribution factors" and they are computed as follows:

a hvdc ik = E d,i V hvdc,k F hvdc ik cos(δ i -θ hvdc,k ) + D hvdc ik sin(δ i -θ hvdc,k ) b hvdc ik = E d,i V hvdc,k D hvdc ik cos(δ i -θ hvdc,k ) -F hvdc ik sin(δ i -θ hvdc,k ) (10)
The same procedure can be performed in order to include loads with constant power characteristic [START_REF] Pai | Energy function analysis for power system stability[END_REF] as further done in this paper. Finally, Equations (1) and 9 can be rewritten to obtain the state space model of the AC power system with HVDC power injections used in this paper for the derivation of the control law, which is given by:

ẋ = δ ω = f (δ, ω) + [B P B Q ] u (11)
where

δ = [δ 1 . . . δn] T , ω = [ω 1 . . . ωn] T u = [∆P hvdc,1 . . . ∆P hvdc,k ∆Q hvdc,1 . . . ∆Q hvdc,k ] T f (δ, ω) =           ω 1 . . . ωn 1 M1 P m,1 -K D,1 ω i -P gg,1 (δ) . . . 1 Mn Pm,n -K D,n ωn -Pgg,n(δ)           Bp =        0 n×k a hvdc 11 M 1 • • • a hvdc 1k M 1 . . . . . . . . . a hvdc n1 Mn • • • a hvdc nk Mn        , B Q =        0 n×k b hvdc 11 M 1 • • • b hvdc 1k M 1 . . . . . . . . . b hvdc n1 Mn • • • b hvdc nk Mn       

Impact of embedded HVDC on AC grid stability

In this section, Direct Methods are used to assess the transient stability of the Single Machine Infinite Bus (SMIB) system in Figure 2.

The case corresponds to an AC system that has been reinforced by two different means: the system in figure 2a has been reinforced with an AC line, while in the case represented in Figure 2b, an HVDC link has been chosen for the reinforcement. All the AC lines are assumed to be loss-less lines of reactance X 12 . The internal reactance of the generator is assumed negligible against the reactance of the AC transmission line (i.e. X d1 << X 12 ). Using the methodology proposed in Section 2, the dynamics of both systems are described by:

• a) SMIB AC system:

M 1 ω1 + k D,1 ω 1 = P m,1 - Pe,1(δ1) (2Pmax sin δ 1 ) (12)
• b) SMIB AC/DC system:

M 1 ω1 + k D,1 ω 1 = P m,1 - Pe,1(δ1) Pmax sin δ 1 + P hvdc,1 (13)
with Pmax = E d,1 V 2 /X 12 .

The Equal Area Criterion

The power-angle (P-δ) curve of both systems is depicted in Figure 3. This curve represents the power delivered by the generator as function of the rotor angle. By means of the Equal Area Criterion (EAC), the stability margins of both cases can be analysed [START_REF] Gonzalez-Torres | HVDC protection criteria for transient stability of AC systems with embedded HVDC links[END_REF].

Fig. 3: Equal area criterion comparison

In this comparison, every transmission line (HVAC or HVDC) in both systems transmits the same amount of power in the pre-fault steady-state situation. Therefore, the constant active power reference for the HVDC link is set to be equal to the active power for one AC line in steady state (i.e. P hvdc1 = Pmax sin δ s ). This assumption allows to have the same equilibrium points in both systems (the same angle in steady state). In the plotted region systems have two equilibrium points, one being stable and the other unstable [START_REF] Kundur | Power system stability and control[END_REF]. It is possible to note that, the stable equilibrium point (SEP) of both systems is the same (δ 1 = δ s ), as well as their unstable equilibrium point (UEP) (δ 1 = δ u ).

Let us assume that a temporary three-phase fault occurs in bus 1. When the fault occurs, the power delivered by the generator suddenly changes from point 1 to 2 in Figure 3 (it goes to zero). During the fault, both systems accelerate following the same trajectory as the machine's delivered power is zero until the fault disappears. Due to this acceleration, the angle increases from δs to δ cl . During the timespan between fault occurrence and fault disappearance, the system acquires a kinetic energy proportional to surface 1-2-3-4. For the described fault, the acquired energy of both systems is the same (i.e.

A AC 1 = A DC 1 )
. When the fault is cleared, the systems get back to their pre-fault power-angle curve (trajectories go from points 3 to 5 in Figure 3). The angle continues to deviate until the kinetic energy absorbed during the fault is expended by transferring it to the grid [START_REF] Kundur | Power system stability and control[END_REF]. The systems move from 5 to 6, such that area

A AC 2 (A DC 2 ) is equal to area A AC 1 (A DC 1 )
. For the studied fault, it can be observed that the maximum angle reached during the first swing in the AC/DC SMIB is larger than the one in the AC SMIB (i.e. δ DC max < δ AC max ). This shows that the AC/DC SMIB is closer to instability for the same fault. Furthermore, after fault clearing, the available energy for decelerating the system is proportional to the surface 4-5-8, it can be easily observed that in the AC/DC SMIB, this surface is smaller than the surface in the AC SMIB case. This means that transient stability margins are larger in the AC SMIB than in the AC/DC SMIB. This phenomenon is due to the nature of the AC lines, their synchronizing power injections (AC lines transmit power as function of the transport angle) contribute to decrease the separation of areas. In the AC/DC SMIB case, as the HVDC link is set in constant power mode, its contribution against angle deviations is nonexistent.

Comparison via Energy functions

For more general conclusions, the transient stability of both systems can be compared by means of the Transient Energy Functions (TEF). These functions represent the sum of the kinetic and potential energy of the system [35]. For a SMIB system represented by Equation [START_REF] Hao | Design of modulation controller to damp power oscillations of parallel AC line in the Tianshenqiao to Guangdong HVDC transmission[END_REF], the TEF can be found by computing the first integral of motion, as follows:

W (ω 1 , δ 1 ) = W KE 1 2 M 1 ω 2 1 + W P E δ1 ∫ δ s
-[P m,1 -P e,1 (δ 1 )] dδ 1 [START_REF] Sigrist | A fundamental study on the impact of HVDC lines on transient stability of power systems[END_REF] where the total energy W is the sum of the kinetic and potential energy noted W KE and W P E respectively. From Equation ( 14), it must be noted that there is a relation between the Equal Area Criterion and the potential energy function W P E . Indeed, W P E represents the surface 1-6-7 on the equal area criterion in Figure 3.

The transient energy functions of the described systems under comparison are the following:

• a) AC SMIB system:

W ac = 1 2 M 1 ω 2 1 - W ac P E (δ1) P m,1 (δ 1 -δ s ) -2Pmax(cos δ 1 -cos δ s ) (15)
• b) AC/DC SMIB system:

W dc = 1 2 M 1 ω 2 1 - W dc P E (δ1) (P m,1 -P hvdc,1 ) (δ 1 -δ s ) -Pmax(cos δ 1 -cos δ s ) (16)
where W ac P E and W dc P E stand for the potential energy of the AC SMIB and the AC/DC SMIB respectively. The evaluation of the time derivative of the total energy (W ac and W dc ) gives an estimation on how the total energy decays when the system is subjected to a disturbance. The computation of the time derivative of Equations ( 15) and ( 16) yields (see [START_REF] Latorre | Application of Control Lyapunov Functions to Voltage Source Converters-based High Voltage Direct Current for Improving Transient Stability[END_REF] for more details):

Ẇ ac = Ẇ dc = -k D,1 ω 2 1 (17)
This equation shows that the energy decays at a rate proportional to the damping coefficient (k D,1 > 0 ) and the square of the speed deviation. It is also found that the rate of change of energy in both systems is the same. Therefore, for a more convenient comparison it is proposed to analyse the potential energy on both systems. For this purpose, the potential energy functions W ac P E and W dc P E are plotted alongside the power-angle curve of each system in Figure 4. The evolution of the potential energy as function of the angle deviation shows that there is a local minimum in the SEP. Besides, as described in [START_REF] Padiyar | Structure Preserving Energy Functions in Power Systems: Theory and Applications[END_REF], the Potential Energy Boundary Method states that there is a local maximum in the UEP (δ 1 = δ u , ω 1 = 0) as it is confirmed in Figure 4. According to the the Closest Unstable Equilibrium Point method in [35], an approximation of the stability boundary can be made by finding the equi-potential contour crossing the closest UEP. The potential energy value at this point is defined as the critical energy V cr . In both systems under study, V cr is equal to the potential energy at the UEP. The comparison of the local maximum values of potential energy of both systems shows that the critical energy is bigger in the AC SMIB than in the AC/DC SMIB (V cr ac > V cr dc ). This means that for the systems under study the transient stability margins are bigger in the AC SMIB case.

The same comparison can be done for all stable equilibrium points. In order to do so, the critical energy V cr as function of a given equilibrium point (δ 1 = δ s , ω 1 = 0) is expressed as follows:

• a) Critical energy for the AC SMIB system:

V cr ac (δ s ) = W ac (π -δ s , 0) = P m,1 (2δ s -π) + 4Pmax cos δ s (18)
• b) Critical energy for the AC/DC SMIB (the power reference of the HVDC is constant and it is equal to P hvdc,1 = Pmax sin δ s ):

V cr dc (δ s ) =W dc (π -δ s , 0) = (P m,1 -Pmax sin δ s ) (2δ s -π) + 2Pmax cos δ s (19)
Dividing Equations ( 18) and 19 by Pmax and eliminating the common terms, it can be shown that V cr ac > V cr dc when

2 cos δ s > sin δ s • (π -2 δ s ) (20)
which is true for all stable equilibrium points of interest (i.e. 0 < δ s < π/2). Consequently it can be stated that, in all the stable equilibrium points (0 < δ s < π/2), the critical energy of the AC SMIB system is higher than the critical energy of the AC/DC SMIB. Further analysis can show that, for a given equilibrium point (0, δ s ), all positive values of P hvdc1 will result in lower values of critical energy than its equivalent AC SMIB system. The vector field and the region of attraction of both systems are obtained by numerical integration and are depicted in the phase plane in Figure 5. The SEP is (δ 1 = δ s =30 o , ω 1 =0), and the region of attraction is delimited by the thick line. Some trajectories starting at different initial points (or points just after fault clearing) are also depicted. It can be observed that, when the post-fault trajectory starts inside the region of attraction, its energy is lower than the critical energy (W < V cr ) and the system remains stable, otherwise the system will diverge. It can be also seen that, as damping has been neglected (K D,1 ), those trajectories continue to oscillate. As suggested by Equation ( 17) if the damping is neglected, the total energy (W ac and W dc ) is constant after fault clearing (equiscalar contour). On the contrary, if the damping coefficient is positive, the energy dissipation will be positive leading the system to reach back the SEP.

Eigenvalue analysis

If after the first swing, the described SMIB system remain stable, and the damping coefficient is positive, the systems undergo oscillations until they reach their equilibrium point. These oscillations in power systems are generally analysed using the Lyapunov's first method. Linearization of the system around its equilibrium point can supply information about the frequency and damping factor of such oscillations. Linearizing Equations ( 12) and ( 13) around an stable equilibrium point, the eigenvalues of the SMIB system can be found. The damping factor ξ and the oscillation frequency f are found to be:

ξ = K D,1 2 √ KsM 1 , f = Ks M 1 - K D,1 2M 1 2 ( 21
)
where Ks is known in the literature as the synchronising power constant and, in this case, represents the slope of the power-angle curve.

For both systems it can be calculated as follows:

• a) SMIB AC system

Ks = K AC s = 2Pmax cos δ s (22)
• b) SMIB AC/DC system (the power reference of the HVDC is constant and equal to P hvdc,1 = Pmax sin δ s )

Ks = K DC s = Pmax cos δ s (23)
The synchronising power constant from case (a) K AC s is bigger than the one in case (b) K DC s for all the stable values of δ s . Consequently, the damping coefficient is higher for the hybrid AC/DC case than in the pure AC system.

It is possible to conclude, that HVDC links with constant references can damp power oscillations in a mixed AC/DC system, with respect to their equivalent AC system; as it has been actually observed by power system operators (e.g. [START_REF] Breuer | Application of HVDC for large power system interconnections[END_REF]). However it doesn't mean that transient stability is enhanced, as it was previously shown in this Section.

Linearizing feedback control

In order to enhance the stability margins of the AC network, the active power references of the VSC-HVDC are modulated based on feedback linearization techniques, such as to actively support the AC grid. The idea behind the feedback linearisation technique is to compensate the nonlinearities of the system by a control signal, produced by feedback of measured variables and knowledge of the system [START_REF] Slotine | Applied nonlinear control[END_REF]. The obtained system, with reduced nonlinearities is further controlled applying linear control techniques for deciding the dynamic behaviour of the system. The proposed control is used to derive a supplementary control scheme which computes the necessary active power modulation for the power-mode converter allowing to improve the rotor angle stability of the surrounding AC grid, damp the electro-mechanical oscillations and to ensure the synchronisation of two AC interconnected networks.

E d,1 1 T 1 3 4 T 2 2 E d,2
P load,1 + jQ load,1 P load,2 + jQ load,2

P hvdc,1 + jQ hvdc,1 -P hvdc,1 +jQ hvdc,2 MMC1 MMC2 X d,1 X d,2 2X ac 2X ac G 1 G 2 1
Fig. 6: 2-Area 2-Machines AC/DC power system Following the procedure described in Section 2, and neglecting damping effects, the dynamic equations of the power system depicted in Figure 6 are given by:

    ω1 ω2 δ1 δ2     =               1 M 1 P m1 -E d1 E d2 Xeq sin(δ 1 -δ 2 ) -a L 11 P L1 -a L 12 P L2 -(a hvdc 11 + a hvdc 12 )P hvdc0,1 1
M 2 P m2 -E d1 E d2 Xeq sin(δ 2 -δ 1 ) -a L 21 P L1 -a L 22 P L2 + (a hvdc 21 + a hvdc 22 )P hvdc0,1 ω 1 ω 2               +     1 M1 (a hvdc 11 + a hvdc 12 ) -1 M2 (a hvdc 21 + a hvdc 22 ) 0 0     ∆P hvdc1 (24)
where ∆P hvdc1 stands for the control input and it represents the supplementary power injection of the VSC-HVDC link and coefficients a L il are the distribution factors expressing the impact of the l -th load on the i -th generator. The equivalent reactance is the sum of the reactances between generator G1 and G2 (i.e. Xeq = X d1 + X T 1 + Xac + X T 2 + X d2 ). Generators are said to be in synchronous equilibrium if all angular distances δ i -δ j are bounded and all frequencies are identical ω i = ω j , therefore the objective of the control is to make ω 1 -ω 2 = 0. The error between the two regions can be modelled following the change of variables ω 1 -ω 2 = ω 12 and δ 1 -δ 2 = δ 12 as:

δ12 = - M 1 + M 2 M 1 M 2 E d1 E d2 Xeq sin(δ 12 ) + 1 M 1 P m1 - 1 M 2 P m2 -( a L 11 M 1 - a L 21 M 2 )P L1 + ( a L 12 M 1 - a L 22 M 2 )P L2 + ( a hvdc 1 M 1 + a hvdc 2 M 2 )P hvdc0,1 + ( a hvdc 1 M 1 + a hvdc 2 M 2 )∆P hvdc0,1 (25)
with δ12 = ω 12 . It can be observed that using Equation (25), the system is now written in strict-feedback form (see [START_REF] Khalil | Nonlinear Systems[END_REF]), where:

δ12 = f (δ 12 ) + c + b u (26)
with u the control input, f (δ 12 ) is the nonlinear function of the angle difference, and c is a constant making its derivative equal to zero in steady state (i.e. δ12 = δ12 = 0). The linearizing feedback control law can be computed as follows:

u = u lf + v = -b -1 (f (δ 12 ) + c + v) (27)
where u lf is the linearizing feedback part of the control law and v introduces the linear dynamics that are to be imposed to the system. For the system represented in Equation ( 25) the linearizing control law is given by:

u = -P hvdc0,1 - M1M2 a hvdc 1 M1+a hvdc 2 M2 M1+M2 M1M2 E d1 E d2 Xeq sin δ 12 ) + Pm1 M1 -Pm2 M2 -( a L 11 M1 - a L 21 M2 )P L1 + ( a L 12 M1 - a L 22 M2 )P L2 -v (28
) Note that if the internal reactance of the machines and the transformer are small compared to the AC line (X d1 << Xac), then a L 21 = a L 12 = 0 and a hvdc 2 = a hvdc 2 = 1, and the system is found to be described by:

δ12 = Pm1-P L1 M1 -Pm2-P L1 M2 -M1+M2 M1M2 (P hvdc0,1 + E d1 E d2
Xac sin(δ 12 ) + ∆P hvdc1 ) (29) Then the linearizing control law is:

u = -P hvdc0,1 -E d1 E d2 Xeq sin(δ 12 ) + M1M2 M1+M2 ( Pm1-P L1 M1 -Pm2-P L1 M2 ) -M1M2 M1+M2 v (30)
The first part of the control law compensates the nonlinearities of the original system, while the term v is used to define the desired closed loop dynamics. Then using linear control techniques, v can be designed to impose the electromechanical mode between two regions [START_REF] Gonzalez-Torres | Procédé de commande d'un lien de transmission électrique incluant une ligne haute tension continu[END_REF]:

v = -K δ -Kω δ12 -δ 12 ω 12 (31)
where K δ and Kω are control gains. Injecting the linear control law [START_REF] Chen | Control-induced time-scale separation for multiterminal high-voltage direct current systems using droop control[END_REF] and the feedback linearizing control law (28) into the system [START_REF] Arioua | Robust grid-oriented control of high voltage DC links embedded in an AC transmission system[END_REF], the dynamics of the closed-loop are found to be:

δ12 = -Kω δ12 -K δ ( δ12 -δ 12 ) (32)
with δ12 = ω 12 . The control gains can be designed using the pole placement technique. For this case, the gains can be chosen as function of the desired natural frequency (ωn) and damping coefficient (ξ) of the controlled system, as follows:

K δ = ω 2 n and Kω = 2ξωn (33)
Since the system has been linearized by the control action, it can be concluded that the region of attraction of the controlled system is theoretically determined only by the power limits (rated power) of the VSC-HVDC link.

Physical interpretation of the control actions

The control law in ( 30) is composed of three main actions:

• The first action (noted u lf in Equation ( 27)) allows to compensate for nonlinearities of the system. This action can also work as a feedforward action that compensates disturbances, especially if some of those values are measured (e.g. the power flowing through the parallel AC transmission line). Indeed, from the control law [START_REF] Ghandhari | Control Lyapunov functions for controllable series devices[END_REF] it can be remarked that the term f (δ 12 ) corresponds to the actual power flowing through the AC corridor, therefore the term c corresponds to the scheduled power to flow through such AC corridor. This property can highly simplify the practical implementation of this control action.

• An action injecting power proportionally to the speed deviation of both regions (see Equation ( 31)), contributing to the damping torque of the system. • An action injecting power proportionally to the transport angle between both regions (see Equation ( 31)), contributing to the synchronising torque. As it will be shown in the examples, this action allows to synchronise both regions even if the parallel AC lines are tripped.

Test of the proposed control under different disturbances

In this section, the proposed control scheme is applied to two different systems using the electromagnetic transient software EMTP-RV. The interest of using an EMT software is to take into account the dynamics of the Phasor Measurement Units (PMU), detailed controls of the converters, as well as studying the behaviour of the control in case of a DC fault. The HVDC stations in the link are two 500MVA Half-bridge Modular Multilevel Converters linked by a 200km wide-band DC cable model, the implemented EMT models are described in [START_REF] Shinoda | Virtual capacitor control: mitigation of DC voltage fluctuations in MMC-Based HVDC Systems[END_REF][START_REF] Loume | DC cable modelling and High Voltage Direct Current grid grounding system[END_REF].

In the first two examples of this section, the simulated power system corresponds to a two-area power system connected by two 200km AC lines in parallel with a VSC-HVDC link as shown in Figure 6. In the simulated system, two 900MVA generators equipped with AVR and PSS are implemented. Their parameters are the ones described in [START_REF] Kundur | Power system stability and control[END_REF]. The load in the first area is 300MW, while in the second area the load is 800MW, so the power flowing through the AC/DC corridor is 500MW. The implemented supplementary control is depicted in Figure 7. By means of the PMUs, phase angle (δ i ) and electrical pulsation (ω i ) at the PCC of both converters are measured. The communication system between the PMUs and the controller is represented by a delay. In this article it was considered a delay of 60ms [START_REF] Cui | Modeling and simulation of dynamic communication latency and data aggregation for wide-area applications[END_REF]. The reference angle between both stations δ12 and the scheduled power exchanged are signals received from a higher level controller. The first case study corresponds to a 250ms 3-phase fault occurring in the middle of the AC line at t = 1s. Three different control strategies are compared:

• CR (Constant references): The power reference of the converter is maintained constant.

• POD (Power oscillation damping): The control proposed in [START_REF] Latorre | Active and reactive power control of a VSC-HVDC[END_REF] is implemented. This controller modulates the power reference as a function of the frequency deviation in both PCCs of the HVDC link (∆P hvdc = K(ω 1 -ω 2 )).

• NLF (Nonlinear Feedback): The proposed control is implemented and tuned to obtain an inter-area oscillation with a natural frequency of ωn = 0.5Hz and a damping factor of ξ = 0.5. The gains are K δ = 1.7MW/deg and Kω = 197MW/Hz

In Figure 8, it is seen that in the CR case, angle deviations during the first swing are higher compared with the other strategies. This means that the AC system is closer to lose transient stability (or power system split). With the CR strategy, it is also clearly observed that inter-area oscillations are poorly damped. When the POD strategy is implemented, it can be seen that the VSC-HVDC link is able to damp power oscillations via the modulation of the active power. It is also observed that the angle deviation during the first swing has been decreased. In the post-fault equilibrium point, the power reference of the HVDC goes back to its initial reference. In other words, (a) Frequency difference between both areas (b) Transport angle between both areas (c) Power transferred through the HVDC link Fig. 8: Comparison of different control schemes after a 250ms threephase fault on one AC line followed by the line tripping the power flowing through the tripped AC line before the fault, is reallocated on the remaining AC line. This is reflected in the increment of the angle in the post-fault steady-state compared with the pre-fault value (35 o before the fault against 58 o after the fault). The NFL control, by taking into account the disturbance in the system, will rapidly compensate the power imbalance due to the fault. This rapid action will limit the angle deviation between regions during the first swing. After the fault, if power is available, the power of the tripped AC line is reallocated to the HVDC link. This action will take the transport angle near its initial condition. As a consequence, stability margins of the post-fault configuration are similar as in the initial configuration.

In order to quantify the transient stability margins for each strategy, the Critical Clearing Time for two different fault locations are computed by iterative simulations and are summarised in Table 1. In every case, the fault is cleared by the tripping of the faulty line. It can be observed that both supplementary controllers (POD and NLF) increase the stability limits with respect to the CR strategy. Besides, it is also shown that the stability margins are higher when the NLF strategy is implemented. In this scenario, the disconnection of the whole group of AC transmission lines connecting both areas after a fault is studied. Let us consider that after a fault on the AC transmission lines, both of them are tripped, then no AC lines are interconnecting both areas. In order to preserve the power exchanged between areas, a standard solution is to reallocate the power transmitted by the AC lines into the VSC-HVDC link (similar to the fast ramp strategy studied in [43]). Power reallocation would assure the right power exchange between both regions, however the split of the system won't be avoided. If both regions are not completely synchronised, the reconnection of the AC line cannot be done. In this case study, two strategies are compared: POD and NFL. Figure 9 shows that, after a few seconds, the POD control is able to ensure the right power exchange between regions and takes the value of the speed difference ω 12 to a stable value different from zero. This means that synchronism of both areas is not achieved. If the system are not synchronised it is not possible to reconnect the AC lines (a reconnection attempt is shown at t=10s).

The implementation of the NFL control shows that synchronisation of both areas can be achieved even if no physical AC line is connected. The proposed control scheme allows to automatically replace the disconnected group of lines while damping power oscillations, assuring power exchange and maintaining synchronism by emulating the behaviour of a modifiable transmission line connecting both areas. Besides, it can be seen that if the AC group of lines is reconnected later, at t=10s, the system remains stable and goes back to its initial equilibrium point. Synchronisation is achieved when ω 12 = 0. If both areas are synchronised, the AC transmission line can be reconnected, even if the frequency of each region is not stabilised at 1pu.

Split of the system due to uncoordinated Frequency Response

The system used in the last two cases correspond to a modified version of the Klein-Rogers-Kundur's two-area power system [START_REF] Kundur | Power system stability and control[END_REF]. In this modified version the two areas are interconnected by a group of AC transmission lines and a VSC-HVDC link in bipolar configuration. In addition, power electronics interfaced power sources are connected in the first region in order to provide fast frequency response (FFR). The procedure described in Section 4 is applied to this power system, and it is chosen to derive a control law that linearizes the input with respect to the error of frequency and angle between two regions given by:

x = δ COI1 -δ COI2 ω COI1 -ω COI2 (34)
where δ COI1 and δ COI2 are the angles of the centre of inertia (COI) of Area 1 and Area 2 respectively, ω COI1 and ω COI2 are the frequencies of the COI of Areas 1 and 2.

δ COIq = 1 M T p i=1 M i δ i , ω COIq = 1 M T p i=1 M i ω i M T = p i=1 M i ( 35
)
where p is the number of generators of synchronous generators in the Area q of interest.

Poorly coordinated frequency response:

A result of the increase in the penetration of renewable energies in a grid is the loss of inertia. The lack of inertia in the grid will lead to a more volatile power system. Controllable devices such as HVDC transmission and Power electronics interface sources must take part in the stabilisation of the grid. As reported in [START_REF] Eriksson | Synthetic inertia versus fast frequency response: a definition[END_REF], one possibility to do so, is to make the renewable sources participate into frequency control by supplying fast frequency response. While the primary response of synchronous machines acts some seconds after the disturbance, the power electronics interfaced sources can supply a faster response. However, as mentioned in [START_REF] Eriksson | Synthetic inertia versus fast frequency response: a definition[END_REF], this fast power injections can jeopardise system transient stability. In the case presented in this section, it is shown how those fast power injections may lead to system loss of synchronism. The evaluated contingency in this case study is the loss of generator G5 in Fig. 10. In order to describe the phenomena, two cases are compared. In both cases, the point to point VSC-HVDC system is controlled to track constant power references.

Stable case: No fast frequency response is injected by the renewable sources after generator tripping, only the primary reserve of the generators participates on frequency control. In Figure 11a, the frequency of the COI of each region as well as the COI of the whole system are represented. It can be seen that after generator tripping the global frequency drops, and frequency of each areas oscillates around the global frequency. The system remains in synchronism and some seconds after the contingency, frequencies reach the same steady state value.

Unstable case: Fast frequency response is injected. While generator tripping occurs in zone two the fast frequency response is injected in zone one, so the system will naturally transport the injected power from zone 1 to zone 2 by increasing the transport angle between regions. FFR power injections will limit the drop of the frequency on Area 1, while Area 2 has to wait the transport angle to become larger in order to receive that increase of power. As shown in Figure 11a, during the first hundred milliseconds, the FFR limits the frequency drop of Area 1 but not in Area 2, this will lead to higher frequency difference between regions, thus a transport angle beyond the stable limit and consequently to system split.

The HVDC link is a suitable mechanism to share the fast reserves while maintaining the synchronism of the system and damping power oscillations. As done before, two controllers are studied, the POD and the NFL control. Figure 12 shows the system response with different controllers. It can be observed that with both strategies, the synchronisation is assured and frequency oscillations are well damped while the FFR response is correctly delivered by the power electronic interfaced source.

Comparing those results with the system response without supplementary controllers, it is observed that the stable frequency in post-fault operation is higher than in the case where no FFR is supplied (0.996 p.u. vs 0.994 p.u.). When no FFR is injected, the minimum value of the global frequency is around 0.992 p.u. during the first swing. This minimal value is improved when controllers are implemented. Since the proposed control takes the angle difference as a feedback signal, in the equilibrium point after the disturbance, the VSC-HVDC link transmits more power than in the initial equilibrium point (250 MW vs. 200 MW). This means that the angle difference between regions in the post fault equilibrium point is smaller when the NFL control is implemented, so stability margins of the post-fault situation are bigger when the NFL control is implemented. The impact of DC faults on transient stability as well as strategies to minimise this impact, have been scarcely addressed in the literature [START_REF] Tang | Impacts of three MMC-HVDC configurations on AC system stability under DC line faults[END_REF]. Hence, this study case aims to show how transient stability can be jeopardise in case a fault occurs in the DC side of the transmission line. In addition, we aim to show how the modulation of active and reactive powers can be used to enhance transient stability in such cases. To do so, we can consider that a temporary pole-to-ground fault occurs on DC side of the negative pole of the first station. All stations in this case are half-bridge Modular Multilevel Converters. The protection strategy used for treating the DC fault is the Converter-Breaker strategy presented in [START_REF] Loume | A multi-vendor protection strategy for HVDC grids based on low-speed DC circuit breakers[END_REF]. The different steps of the protection strategy shown in Figures 14 and13 When the fault occurs (t 0 ), the current flowing through the station increases until the fault is detected and the DC circuit breakers (DCCB) disconnect the faulted cable, and hence suppress the current contribution to the fault. It should be noted that after fault and before DCCB breaker opening, MMC converters are rapidly blocked and the AC contribution to the fault passes through the freewheeling diodes [START_REF] Bertinato | Development of a protection strategy for future DC networks based on low-speed DC circuit breakers[END_REF]. When the converter is blocked it exhibits the behaviour of an uncontrolled rectifier. When the faulted cable is isolated, the MMC can be deblocked (t 1 ), therefore the faulted pole can operate in STATCOM mode. As it can be seen in Figure 13d, while the MMC is blocked (between t 0 and t 1 ), reactive power is zero, while after deblocking the reactive power control is recovered. Before reconnection, the overhead line requires deionization, which in this simulation lasts 350ms (t 2 ) after DC circuit breaker opening [START_REF] Tang | Impacts of three MMC-HVDC configurations on AC system stability under DC line faults[END_REF]. After this time, the cable can be reconnected and the control of the DC voltage is recovered. Then, the system is ready to return to operation (t 3 ) and power is smoothly restored to its initial value.

For this example power references are set to 300 MW for the VSC-HVDC link and the total exchanged power between 2 regions is 500 MW. Figures 13 and14 show the response of the DC and AC systems respectively, in case of a DC fault. Three control strategies are compare:

• Constant references: Active and reactive power references are maintained constant • NFL-Q ref constant: Active power references are modulated using the proposed non-linear feedback control • NFL-∆Q hvdc : Reactive power injections are taken into account in the linearizing feedback law When no supplementary control is implemented, the temporary power interruption in one pole leads to the loss of transient stability. In the case where NFL is implemented, the linearizing feedback (a) Speed deviation between zones (b) Angle difference between zones Fig. 14: AC System response after a DC fault and implementation of the converter breaker strategy [START_REF] Loume | A multi-vendor protection strategy for HVDC grids based on low-speed DC circuit breakers[END_REF] action of the controller includes the information about the power references (P hvdc0 ), thus it is possible to update the value of the static references including the loss of power of the faulted pole with a negligible communication delay. In the second and third case, this action compensates the lost pole with the power available in the healthy pole as it can be seen in Figure 13a. The healthy pole goes from its initial value (150MW) to its maximum power (around 220MW). From the AC grid point-of-view, the power lost is only of around 80 MW. This power compensation will reduce the effect of the contingency on transient stability.

When the MMC station is deblocked after fault clearing, it is possible to work in STATCOM mode before recovering the active power injections. As shown in Equations ( 11) reactive power injections of the MMC stations have an impact on the active power of the synchronous machines. Therefore when the STATCOM operation begins, then reactive power injections can be used as a control input for enhancing transient stability in case of DC fault. In Figures 14 and13c reactive power injections of the negative and positive poles are shown. While the healthy pole compensates the power loss (between t 1 and t 2 in Figure 13a), the inner current control gives priority to the active power injections, this is why reactive power goes to zero (Figure 13c). However, after deblocking, reactive power injections of the faulted pole can be modulated for reducing the impact of the active power loss. Figure 14b shows that angle difference between regions are slightly reduced when reactive power is modulated. The impact of reactive power injections depends on the distribution factors determined by the topological situation of the VSC-HVDC link.

Conclusions

In this article, the impact of embedded VSC-HVDC links on the surrounding AC grid is discussed and a supplementary control for stability purposes has been developed. The proposed control has been tested for different contingencies. Important conclusions should be taken from this article:

• Transient stability margins are lower when a system is reinforced with VSC-HVDC lines instead of AC lines. Indeed, VSC-HVDC links controlled to track constant power references, do not inherently improve transient stability of the surrounding AC grid as an AC line naturally does. Hence, supplementary controllers on VSC-HVDC links are needed for supporting AC transient stability.

• Through the proposition of the control, three main actions for supporting transient stability were identified: injection of synchronising power, injection of damping power and fast compensation of power disturbances.

• Synchronisation of two power grids can be achieved through the proposed control even when no AC link is interconnecting both grids.

• Poor Fast Frequency Response can destabilise the grid. Through the proposed control, the VSC-HVDC link can help the system to properly share those reserves.

• The effect of a DC fault on transient stability has been analysed. The importance of fast power reallocation as well as reactive power injections in case of a fault on the DC side of the converter has been highlighted.

Fig. 1 :

 1 Fig. 1: Multi-machine/multi-converter power system

Fig. 2 :

 2 Fig. 2: Comparison between: (a) AC SMIB and (b) AC/DC SMIB

Fig. 4 :

 4 Fig. 4: Potential energy and power-angle curve comparison.

Fig. 5 :

 5 Fig. 5: Comparison of the regions of attraction. Characteristics of the system: δ s = 30, P m,1 = 1pu , Pmax = 1pu, P hvdc,1 = 0.5pu, k D,1 = 0, M = 0.12MW s 2 /MVA rad (i.e. H 1 = 4s ).

Fig. 7 :

 7 Fig. 7: Control implementation

5 . 2

 52 Synchronisation in case of the disconnection of the AC parallel lines

( a )Fig. 9 :

 a9 Fig. 9: Comparison of different control schemes when the AC transmission line is disconnected.

Fig. 10 :

 10 Fig. 10: Modified version of the Klein-Rogers-Kundur's two-area power system with power electronics interfaced energy sources

Fig. 11 :

 11 Fig. 11: System split caused by a poorly coordinated fast frequency response.

Fig. 12 :

 12 Fig. 12: Comparison of controllers when generator 5 is tripped and dynamic frequency response is applied

  are summarised in the following table.

Fig. 13 :

 13 Fig.13: Active and reactive power injections during execution of the converter breaker strategy triggered by a DC fault[START_REF] Loume | A multi-vendor protection strategy for HVDC grids based on low-speed DC circuit breakers[END_REF]

Table 1

 1 Critical clearing times for different fault locations

Table 2

 2 Description of the protection sequence. Temporary Pole-to-ground fault on negative pole t 1 1.1 MMC deblocking and start of STATCOM operation t 2 1.35 DC Line reconnection and DC voltage recovery t 3 1.6 Beginning of active power restoration

	tn t [s] Description
	t 0	1

c The Institution of Engineering and Technology 2020

IET Research Journals, pp. 1-12

IET Research Journals, pp. 1-12 c The Institution of Engineering and Technology 2020

Acknowledgement

This work is supported by the French Government under the program Investissements d'Avenir (ANE-ITE-002-01).