S Speziale
email: speziale@gfz-potsdam.de

J Immoor

A Ermakov

S Merkel

H Marquardt

H P Liermann

The equation of state of TaC0.99 by X-ray diffraction in radial scattering geometry to 32 GPa and 1073 K

We have performed in situ synchrotron X-ray diffraction experiments on TaC0.99 compressed in a diamond anvil cell along 3 isothermal paths to maximum pressure (P) -temperature (T) conditions of 38.8 GPa at 1073 K. By combining measurements performed in axial diffraction geometry at 296 K and in radial geometry at 673 K and 1073 K, we place constraints on the pressure-volume-temperature (P-V-T) equation of state of TaC in a wide range of conditions. A fit of the Birch-Murnaghan equation to the measurements performed in axial geometry at ambient temperature yields a value of the isothermal bulk modulus at ambient conditions 𝐾 𝑇0 = 305  5 (1) GPa and its pressure derivative (𝜕𝐾 𝑇 /𝜕𝑃) 𝑇0 = 6.1  0.5. The fit of the Birch-Murnaghan-

Debye model to our complete P-V-T dataset allows us to constrain the Grüneisen parameter at ambient pressure 𝛾 0 = 𝑉(𝜕𝑃/𝜕𝐸) 𝑉0 to the value of 1.2  0.1.

I. INTRODUCTION

Tantalum carbide TaCx (0.6 < x < 1) is a B1 structured ultra-high temperature ceramic (UHTC) [START_REF] Gasch | Handbook of Ceramic Composites[END_REF]. UHTC materials have the potential to be used in different technological applications because of their combined high mechanical strength, extremely good thermal stability as well as its resistance to harsh chemical environments.

Transition metal carbides are in general among the ceramic materials with the best characteristics for high temperature applications where enhanced refractory and mechanical properties are required [START_REF] Upadhya | [END_REF]. Tantalum carbide has an extremely high melting temperature, low electrical conductivity and unusual mechanical behavior with respect to the other B1 structured (space group 𝐹𝑚3 ̅ 𝑚) monocarbides of the group IV and V transition metals (TMC) [START_REF] Toth | Transition Metal Carbides and Nitrides[END_REF][START_REF] Williams | [END_REF].

The structural, thermodynamic and transport properties of TaCx (in the following we will use TaC except when discussing the effect of non-stoichiometry on specific physical properties) have been the subject of many experimental and computational investigations and part of the existing results are tabulated in several review publications together with the properties of the other TMC . In particular, the elastic properties of TaC have been the subject of more than ten different experimental and computational studies [6,15,16,[28][29][30][31][32][33][34][35][36][37][38][39] exploring the effect of pressure or temperature and non-stoichiometry on the tensor and the aggregate elastic moduli of single crystals or polycrystals.

While several computational studies provide models of the isothermal compression behavior of TaC, only two experimental diamond-anvil cell studies investigate its compression behavior to the multi-GPa stress range at ambient temperature [32,37]. The first study focuses on the response of micron-sized polycrystalline TaC0.98 to compression under different stress conditions, and the second explores the behavior of nanocrystalline TaC under strongly nonhydrostatic stress.

To date there are no data on the combined effect of high-pressures, stress andtemperatures on the structure and the thermoelastic properties of tantalum carbide. With this study we give a first account of the compression behavior of TaC0.99 compressed to more than 40 GPa at ambient temperature and to a maximum of 38.8 GPa at 1073 K, and we place quantitative constraints on the pressure-volume-temperature (PVT) equation of state of this material.

Here we use a new approach to place quantitative constraints on the thermoelastic properties of TaC0.99. We combine the results of static compression experiments in a graphite resistive heated diamond-anvil cell using in situ synchrotron X-ray diffraction in radial scattering geometry (an established technique [40][41][42][43][44][45][46] yet never used before in studies dedicated to constrain the pressure-volume-temperature equation of state) with those from axial geometry at room temperatures. The use of high-temperature highpressure radial diffraction geometry in combination with a theory of non-hydrostatic lattice strains allows us to quantitatively constrain the departure from hydrostatic stress conditions, which is not straightforward in experiments performed in axial geometry. This is particularly relevant for stiff samples such as TaC when they bridge between the culets of the two diamond anvils under combined high pressure and high temperature. In 4 addition the experimental approach used here allows us to measure high temperature volume compression along strictly isothermal compression paths at moderately high temperatures, which is difficult with laser-heated diamond anvil-cell experiments, thus producing the best suited dataset for P-V-T equation of state fitting.

II. MATERIALS AND METHODS

A. Sample material

The starting material for all experiments presented here is a powder of TaC0.99 (with nominal stoichiometric TaC composition) produced by American Elements©. We have determined the chemical composition of our starting material by measuring its unit-cell volume by X-ray diffraction at the German Research Centre for Geosciences (GFZ) using a STOE Stadi P diffractometer with Cu K radiation monochromatized by a Ge(111) crystal. Two separate amounts of sample were mixed with 10 wt% Si standard (NIST 640d) as an internal calibrant. The average of the two measurements yields a unit-cell parameter 𝑎 0 = 4.4555  0.0007 Å. This value corresponds to a composition of TaC0.99 based on the calibration by Bowman [8]. The unit cell parameter of our sample is equivalent within uncertainty to the results of neutron diffraction measurements of a disordered stoichiometric TaC [9] and 3 larger (rather than smaller, as expected) than a more recent neutron diffraction study of nominally stoichiometric TaC [10].

B. High-pressure synchrotron X-ray diffraction in axial geometry

In these experiments the axis of the diamond anvils is parallel to the incident X-ray beam. 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.5115350

TaC0.99 was compressed quasi-hydrostatically in a symmetric piston-cylinder type diamond anvil cell [START_REF] Mao | Advances in high pressure research in condensed matter: proceedings of the international conference on condensed matter under high pressures[END_REF] with 300 m culet diamonds. Fine sample powder (with less than 1m average grain-size) was loaded in a 150 m wide cylindrical chamber drilled in a 250 m thick rhenium disk pre-indented to a thickness of 30 m. Neon gas was loaded as a pressure transmitting medium using the gas-loading facility from Core Labs located at PETRA III and operated by the Extreme Condition Science Infrastructure. A ruby sphere and a  5 m speckle of gently packed gold powder were used as internal pressure standards [START_REF] Mao | [END_REF][START_REF] Fei | Proc. Natl. Acad. Sci[END_REF]. The diamond anvil cell was compressed remotely with a gasmembrane system. X-ray diffraction images were collected at the general purpose experiment table of Extreme Conditions Beamline (ECB, P02. The images were analyzed with the Fit2d software package [START_REF] Hammersley | [END_REF]. Le Bail fit [51] of the integrated spectra were performed with the software MAUD [52] to determine the unit cell parameter a (TaC crystallizes in the cubic system in the space group 𝐹𝑚3 ̅ 𝑚). A typical X-ray diffraction image and the result of the Le Bail fit are shown in Figure 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 AS DOI: 10.1063/1.5115350

The conditions and the results of the room temperature compression experiments are reported in Tables I andII.

C. Simultaneous high-pressure high-temperature synchrotron X-ray diffraction in radial scattering geometry Two series of measurements were performed in radial scattering geometry at the general purpose experiment table of the ECB. In these experiments the axis of the diamond anvils is perpendicular to the incident X-ray beam, passing through an X-ray transparent sample chamber, made of a mixture of amorphous boron and epoxy.

The setup for these experiments is based on a modified Mao-Bell piston-cylinder type diamond anvil cell equipped with a graphite resistive sandwich heater for hightemperature experiments [53]. Pressure in the diamond anvil cell is controlled remotely by a gas membrane system. For temperature measurements two type-R thermocouples are placed in contact with the diamond at short distance from the culets. The diamond anvil cell is placed in a vacuum vessel (capable of maintaining 10 -4 mbar pressure during the experiments) to protect diamonds, the metallic parts of the diamond anvil cell and heater setup from oxidation [54].

TaC0.99 powder with a grain size < II.D). The ceramic gaskets are supported by a kapton ring. In all the experiments the diameter of the diamond culets was 300 m and no pressure transmitting medium was employed.

In the two high-temperature experimental runs a 5-10 m wide fragment of less than 5 m thick Au foil was placed above the sample as a pressure calibrant. The two hightemperature experimental runs were performed using a monochromatic X-ray beam of 0.4847 Å wavelength with an identical focus spot-size as for the axial diffraction experiment but slightly different sample to detector distance of 405 mm based on the same calibration approach as for the axial diffraction experiment.

The exposure time for all images in the radial X-ray diffraction experiments was 4s. The images were processed using Fit2d [START_REF] Hammersley | [END_REF]. Each image was integrated in 5 o azimuthal (𝜓) sectors, for a total of 72 patterns per image. The sets of patterns were further analyzed with MAUD software package [55] to determine structural parameters.

In ideal conditions, the geometry of the diamond anvil cell imposes a cylindrical symmetry of stress and strain at the center of the sample chamber (i.e. along the culets' axial direction). The stress state far from the center of the culets is more complex [56].

For this reason at each pressure step the optimized position was determined by X-ray absorption scans along the vertical and horizontal direction. An additional diffraction image of gold was collected at each compression step for pressure calibration at high temperatures [START_REF] Fei | Proc. Natl. Acad. Sci[END_REF].

In both the 673 K and 1073 K runs the sample was compressed at ambient temperature to 2.5  1 GPa and then progressively heated to the target temperature at an average rate of 3.2 °C/min and 4.6 °C/min, respectively. A total of 47 X-ray diffraction images were 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.5115350 collected along two isotherms at 673 K and 1073 K. The parameters of all the experimental runs are reported in Table I.

D. Data analysis of radial diffraction data with MAUD

All radial X-ray diffraction images were analyzed using MAUD software package.

MAUD determines both crystal structure parameters, and texture by combining Rietveld full-spectrum fitting and a selection of models to describe crystallite sizes, microstrains, macroscopic stress and the crystals' preferred orientation [53,57].

In radial X-ray diffraction experiments the sample is subject to strongly non-hydrostatic macroscopic stress because it is loaded without any soft pressure transmitting medium.

The X-ray diffraction images display anisotropic strains caused by the uniaxial stress applied to the sample. This results in non-circular Debye diffraction rings, with their ellipticity being related to the amount of deviatoric stress and the elastic moduli of the sample material. Specialized theories to describe non-hydrostatic crystal strains have been developed to analyze the X-ray diffraction images collected from samples compressed in the diamond anvil cell and other high-pressure devices [57][58][59]. In our analysis here we are interested in the study of the (volume) compression of TaC, and we correct for the effect of nonhydrostaticity by using the theory developed by Singh [58].

Here we assume that the TaC0.99 polycrystalline aggregate sample behaves as an elastically isotropic material. Following the formulation presented by Singh [58] for the isotropic aggregate case:

𝑑 𝑑 (ℎ𝑘𝑙) = 𝑑 𝑝 (ℎ𝑘𝑙)[1 + (1 -3𝑐𝑜𝑠 2 𝜓)𝑄], (1)
where 𝑑 𝑑 (ℎ𝑘𝑙) is the d-spacing of the (hkl) lattice planes family as determined in the spectrum integrated at an angle 𝜓 measured anticlockwise from the compression direction. 𝑑 𝑝 (ℎ𝑘𝑙) is the d-spacing under hydrostatic stress (corresponding to pressure 𝑝),

and 𝑄 is a measure of the strain due to the deviatoric component of the stress. For isotropic aggregates of elastically anisotropic materials, and assuming that stress is homogeneous in the aggregate Singh [58] derives:

〈𝑄〉 = 𝑡 (1 + 𝜈) 3𝐸 = 𝑡/6𝐺 ⁄ , ( 2
)
Where 〈𝑄〉 is the isotropic aggregate 𝑄, averaged on all the diffracting lattice families (see equation ( 1)), 𝑡 = 𝜎 3 -𝜎 1 is the difference between the normal stress along the compression axis and perpendicular to it, 𝜈 is the average Poisson's ratio, 𝐸 is the Young's modulus of the isotropic aggregate sample, and 𝐺is the isotropic aggregate shear modulus. The model of Singh [58] is incomplete for aggregates subject to plastic deformation and with large density of defects [60]. However, equation (1) has been shown to be useful to describe such experimental data, extract average stress levels [61],

and model plastic behavior of materials using more advanced methods [62].

The unit cell parameter a corresponding to the hydrostatic stress component is determined from the spectrum at 𝜓 = 54.7 degrees. The implementation of the model available in MAUD allowed us to fix the values of 𝜈 and 𝐸 and refine the values of 𝜎 3

and 𝜎 1 to fit the observed d-spacings. We calculated the value of 〈𝑄〉 by using equation [START_REF] Upadhya | [END_REF]. This can finally be easily converted (see equation ( 1)) into the non-hydrostatic stress contribution to the uncertainty on the unit cell parameter. The uncertainties due nonhydrostatic conditions are in average of the order of 10 -2 Å, substantially larger than the 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. In both the experiments performed at 673 K and 1073 K pressure was determined by fitting the unit cell parameter of gold, loaded together with the sample in the sample chamber. The analysis of X-ray diffraction of gold was performed similarly to that of the sample material itself. The stress/lattice strain analysis generally gave results of deviatoric stress comparable with those of the TaC0.99 sample. In some few cases the fit of stress and lattice strain of gold did not converge to consistent results due to lower than average gold diffraction signal.

PLEASE CITE THIS ARTICLE AS

The 3 X-ray diffraction images collected at 6.0, 6.8 and 8.2 GPa at 673 K and the 5 ones collected at 1073 K from 32.9 to 38.8 GPa show a deviation from the ideal geometry, with a maximum of lattice strains in a direction tilted with respect to the diamond anvil axis (Figure 3). Non-ideal stress conditions are typically due to gasket partial failure at the highest pressures (as in the case of the 1073 K isotherm) or, less frequently, at the lowest pressures, at the first contact between the diamond culets and the ceramic gasket (as in the case of the 673 K isotherm). Nevertheless, we performed a full analysis also of these images and we will discuss the results below.

All parameters determined by Rietveld analysis of the diffraction images of TaC0.99 in the two radial X-ray diffraction experimental runs are reported in Table II.

E. P-V-T equation of state analysis

A third order Birch-Murnaghan P-V equation of state [63] was fitted to the unit cell volumes of TaC0.99 measured in the experiments performed in axial geometry and quasihydrostatic stress. Pressures were determined from the unit cell volume of gold loaded in the diamond anvil cell together with the sample [START_REF] Fei | Proc. Natl. Acad. Sci[END_REF]. The pressure -volume data were weighted by the pressure uncertainties on the measurements of both gold and TaC0.99 propagated from the uncertainties of their respective experimentally determined unit cell volumes.

The three coefficients of the third order Birch-Murnaghan equation are the unit cell volume at ambient conditions 𝑉 0 , the isothermal bulk modulus at ambient conditions 𝐾 𝑇0

and its isothermal pressure derivative at ambient conditions(𝜕𝐾 𝑇 /𝜕𝑃) 𝑇0 . We fixed the value of 𝑉 0 to the one we measured at the GFZ (see above) and we refined the values of 𝐾 𝑇0 and (𝜕𝐾 𝑇 /𝜕𝑃) 𝑇0 . The individual data points used in the fit are reported in Table III ("TaC-Ne"). In the fit we did not consider all the data points collected in decompression, however we will mention, in the discussion, their effect on the fit results. After the first decompression step the piston and cylinder of the diamond anvil cell were jammed. We observed a sudden decrease of sample (and Au standard) unit cell volume between each of the decompression points. The decompression volumepressure data do not plot on the compression dataset, even though they overlap within reciprocal pressure uncertainties. Our interpretation is that, due to the sudden stress release, the stress field applied to the sample was different with respect to that generated by the smooth stress increase in compression, and the anomalous stress condition affected the pressure determination. The last data point, collected at 1 bar after complete pressure release, 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.5115350 overlaps with the ambient conditions X-ray diffraction measurement and refinements performed at the GFZ (see Figure 4).

The pressure volume data of TaC0.99 collected in radial diffraction geometry along the two high temperature isotherms at 673 K and 1073 K were used to fit the parameters of the Birch-Murnaghan-Debye P-V-T equation of state (e.g. [START_REF] Duffy | Ultrahigh-Pressure Mineralogy: Physics and Chemistry of the Earth's Deep Interior[END_REF]). Pressures were determined from the measured unit cell volume of gold loaded in the diamond anvil cell together with the sample. Error propagation of the uncertainties on sample and gold experimental unit cell volumes was performed using the same strategy as for the room temperature quasi-hydrostatic experimental dataset. We additionally weighted the data of the high temperature runs (performed without soft pressure transmitting medium) by propagating the additional volume strain caused to the uniaxial stress present in the sample chamber (see section about Qlattice strains analysis). The uncertainty on temperature in the experimental setup used in this study is estimated to be  20 K at 673 K and  60 K at 1073 K [53].

The coefficients of the Birch-Murnaghan-Debye equation are 𝑉 0 , 𝐾 𝑇0 and (𝜕𝐾 𝑇 /𝜕𝑃) 𝑇0

for the reference isotherm, which we fixed to be the ambient temperature isotherm. Three thermal parameters were refined: the Grüneisen parameter at ambient conditions 𝛾 0 , its logarithmic volume derivative 𝑞 = (∂ln 𝛾 𝜕 ln 𝑉 ⁄ ) 𝑇 , and the Debye temperature at ambient conditions 𝜃 0 . Due to our limited dataset we fixed the initial unit cell volume of the reference isotherm (𝑉 0 ) to our own measurement result and the value of 𝜃 0 to 567  10 K which is the average of the available literature data [19,24,30,31,36,39] excluding the two extreme values [34,37]. We thus refined the value of 𝑞 and 𝛾 0 , for which the extant published results are in strong disagreement [17,31,36].

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.5115350

The whole dataset used to fit the parameters of the equation of state of TaC0.99 are reported in Table II.

III. RESULTS

A. Isothermal equation of state of TaC0.99 at 296 K The unit-cell volumes determined in the quasi-hydrostatic experiment are reported in Table II. The compression curve is plotted in Figure 4. The third order Birch-Murnaghan fit of our data yields 𝐾 0𝑇 = 305  5 GPa and (𝜕𝐾 𝑇 /𝜕𝑃) 𝑇0 = 6.1  0.5, where 𝐾 0𝑇 is the initial bulk modulus and (𝜕𝐾 𝑇 /𝜕𝑃) 𝑇0 is its initial pressure derivative. We performed our fit by fixing 𝑉 0 (the initial unit-cell volume) to the value of 88.478 Å 3 based on our measurement at ambient pressure. The root mean square (RMS) misfit of the pressure dataset is 0.29 GPa. We tested the effect of including the three decompression data (these collected at 38.9, 32.0 and 18.8 GPa) in the fitted dataset. The fitted coefficients are insensitive, within estimated uncertainties, to the inclusion of the decompression data.

However, the RMS misfit increases by more than 20% by adding these three data points, and we decided to exclude them in our final fit (see Table II). The model isothermal compression curve is plotted as a solid line in Figure 4.

B. P-V-T equation of state of TaC0.99

The high-temperature unit-cell volumes dataset (Table II) combined with the parameters of the ambient temperature isotherm were used to fit the thermal coefficients of a Birch-Murnaghan-Debye equation of state. We did not include in the dataset the data collected at 6.0, 6.8 and 8.2 GPa and 673 K and the 5 ones collected at 1073 K from 32.9 to 38.8

GPa (see also discussion). In the fit we fixed the value of 𝜃 0 to the value of 570  10 K averaged between the available literature data (see methods above). The fit yields an initial value of the Grüneisen parameter 𝛾 0 = 1.2  0.1. The value of the logarithmic pressure derivative of the Grüneisen parameter (𝑞) is poorly constrained by the current dataset and the best fit value is 𝑞 = 1  2. The RMS pressure misfit of the full high temperature dataset is equal to 0.46 GPa. The complete set of coefficient of the pressurevolume-temperature equation of state of TaC0.99 based on our experimental results is reported in Table III. The experimental unit-cell volumes and the model isotherms at 296 K, 673 K and 1073 K are plotted in Figure 5.

IV. DISCUSSION

A. Isothermal equation of state of TaC0.99 at 296 K The values of 𝐾 0𝑇 = 305  5 and (𝜕𝐾 𝑇 /𝜕𝑃) 𝑇0 = 6.1  0.5 constrained by our new quasihydrostatic dataset compare well with the available computational and experimental results for TaCx (Table IV).

Several ab initio DFT computational studies have investigated the elastic properties of

TaC [6,13,15,18,[33][34][35][36][37][38][39]. The values of the ambient pressure bulk modulus are generally similar in all the studies when the same approximation is used for the exchange-correlation electronic interaction. The average for the studies using the generalized gradient approximation (GGA) is 315  16 GPa and overlaps with our result within reciprocal 1 uncertainties. One of the computational studies [33] combines DFT calculations and experimental ultrasonic interferometry data, for only two of the three elastic stiffness coefficients of TaC. Using the elastic anisotropy expressed as the Zener 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.5115350

anisotropy ratio 𝐴 = 2𝐶 44 /(𝐶 11 -𝐶 12 ) [START_REF] Zener | Elasticity and Anelasticity in Metals[END_REF], from their calculations we can estimate the value of the experimental bulk modulus 𝐾 0𝑆 = 310 GPa, which agrees with our result within 1. Here we point out that the difference between the isothermal bulk modulus 𝐾 0𝑇 and the isentropic bulk modulus 𝐾 0𝑆 expressed by the thermodynamic relationship 𝐾 0𝑆 = 𝐾 0𝑇 (1 + 𝛼𝛾𝑇)(where 𝛼 is the thermal expansion coefficient, and 𝛾 is the Grüneisen parameter) in the case of TaC is of the order of 3 ± 1 GPa, which is comparable with the uncertainties of all the existing experimental studies (Table IV).

Our best fit 𝐾 0𝑇 is at the lower end of the values of isothermal and isentropic bulk modulus reported in the literature except for a single-crystal ultrasonics study of TaC0.90 [29] and a neutron inelastic scattering study of single-crystal TaC [17] which yield 𝐾 0𝑆 = 217  12 GPa and 𝐾 0𝑆 = 283 GPa, respectively. The three studies performed on TaC0.99, present values of the isentropic bulk modulus 𝐾 0𝑆 between 345 GPa, and 317 GPa [28,30,[START_REF] Bukatov | [END_REF] all calculated by extrapolation to full density of ultrasonic experimental studies of porous ceramic materials. If we consider all the available data from elasticity studies including a study of TaC0.98 [31] the average value of 335  13 GPa is compatible with our fitted value only at 3 uncertainty level.

Only three experimental and three computational studies investigated the pressure dependence of the elastic moduli of TaCx [16, 31-33, 36, 37]. The pressure derivatives (𝜕𝐾 𝑇 /𝜕𝑃) 𝑇0 based on experimental studies [31,32,37] are equal to 4.97  0.27, 4.0  0.4 and 4 respectively, and are substantially different from our result. However, the pressure derivative presented in [31] is calculated from ultrasonics pulse-echo measurements of a compressed porous sample without direct measurements of the sample volume. The value reported by [37] is based on a high-pressure X-ray diffraction study of nano-particle TaC, which could have a different elastic behavior with respect to the bulk material (see for instance the case of MgO [67]), and it has been fixed in the data analysis. The results of [32] are based on a similar approach as the one used in the present study and the large disagreement deserves a more detailed explanation. Liermann and coauthors [32] performed high-pressure X-ray powder diffraction of TaC0.98 compressed in the diamond anvil cell using Al as a pressure transmitting medium. Even though Al is a metal with low shear modulus at ambient conditions, it represents a strong pressure medium if compared to Ne at pressures of tens of GPa. The use of strong pressure media unavoidably produces deviatoric stress in the sample chamber. In axial X-ray diffraction geometry the crystallites which contribute to the measured signal are those with diffracting vectors (the directions normal to the diffracting lattice planes) close to perpendicular to the diamond axes, and thus subject to deviatoric extension (relative to the ideal hydrostatic strain) [58]. This causes underestimation of volume compression and overestimation of the material bulk modulus [67,[START_REF] Singh | High Temp[END_REF]. In addition to this, the value of bulk modulus and pressure derivative at ambient conditions obtained by fitting high-pressure X-ray diffraction results suffer from trade-offs between the two coefficients of the fitted equation of state [START_REF] Bass | [END_REF]. In order to better compare the results of the two studies we have plotted the difference between the P-V isotherms calculated using the two sets of parameters and a common starting volume (Figure 6). This procedure corresponds to comparing pressures along the two isotherms at the same values of compression. The difference between the two compression curves is within three times the average uncertainty of our data if we include the error propagated from the unit-cell volume uncertainty. If we consider that the experimental uncertainties on pressure and unit-cell 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.5115350

volume in [32] are of the same order of magnitude as the one in our dataset, the disagreement between the two isothermal equations of state is within 2 reciprocal uncertainties. The similarity between the compression curves of TaC0.99 and TaC0.98 suggests that, within experimental uncertainties, the isothermal compressibility of TaCx at 300 K is only marginally sensitive to low levels of nonstoichiometry (0.98  x < 1).

The pressure derivative of the bulk modulus from the four ab initio computational studies with local density approximation (LDA) [33], generalized gradient approximation (GGA) [36] and with both approximations [16,37] is in all cases substantially smaller (between 29% and 42%) than our result. However, this is partially compensated by their larger bulk modulus at room pressure. The discrepancy between our experimental results and ab initio computations could be due to the presence, in our sample material, of a complex defect structure in addition to slight non-stoichiometry, which includes dislocations, grain boundaries and microstrains associated to them. Indeed it has been observed that carbondefective TaCx ceramics show higher mechanical strength than stoichiometric TaC [START_REF] Toth | Transition Metal Carbides and Nitrides[END_REF].

The characterization of defects states and their energetics in TaC has been subject of several computational and theoretical studies [70,71]. A recent study has shown that increasing C vacancies concentration reduces the elastic moduli of TaCx [72]. However, it is not yet clear what effect this deficiency has on the pressure dependence of bulk modulus. The isothermal compression curves of all the existing high-pressure studies of TaCx [16, 31-33, 36, 37] are calculated with the same procedure as for [32] (see above)

and compared with ours in Figure 6.

B. P-V-T equation of state of TaC0.99

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.5115350

Our full high-pressure / high-temperature dataset places a strong constrain on the value of 𝛾 0 = 1.2  0.1. In the framework of the Birch-Murnaghan-Debye model this coefficient is the thermodynamic Grüneisen parameter at ambient conditions. This value is obtained by fitting both our full dataset, or after removing the 3 data points at 6.0, 6.8 and 8.2 GPa and 673 K and the 5 data for 32.9, 34.9 36.1, 37.9 and 38.8 GPa and 1073 K due to their non-ideal geometry (see section 2. Materials and Methods). Our best-fit value of 𝛾 0 compares very well with the value of the elastic Grüneisen parameter 𝛾 𝑒𝑙 = 1.22 determined for polycrystalline TaC0.98 by ultrasonic techniques [31]. However, the elastic Grüneisen parameter is the average value of the volume dependence of the acoustic phonon frequencies at the Brillouin zone center, a parameter which is different from the thermodynamic Grüneisen parameter 𝛾 𝑡ℎ = 𝛼𝐾 𝑆 (𝜌𝐶 𝑝 ) ⁄ (where 𝛼 is the volume thermal expansion coefficient, 𝐾 𝑆 is the isentropic bulk modulus, 𝜌 is density and 𝐶 𝑝 is the specific heat at constant pressure) which controls the thermal pressure contribution in the P-V-T equation of state. The value of the thermodynamic Grüneisen parameter at ambient conditions calculated from the published thermoelastic properties of TaCx (Table V) ranges between 1.47 and 2.34, that is 30% to 90% larger than the best-fit value to our experimental dataset. The reason of this difference is clearly due to the limitations of the Birch-Murnaghan-Debye model to describe the thermal energy contribution in TaC0.99.

The high temperature study of the elastic properties of TaC0.99 by Jun and Shaffer [30] shows that the thermodynamic Grüneisen parameter has very strong temperature dependence below 600 K at ambient pressure with a decrease from 1.48 at 300 K to 1.2 at 600 K, and then it remains almost constant at the value of 1.2 at temperatures as high as 1500 K. On the opposite, Peng et al. [36] who use the Debye model and a quasi-harmonic approximation obtain 𝛾 slowly increasing with temperature at ambient pressure (Figure 7).

The quasi-harmonic approximation underlying our equation of state assumes that 𝛾 is not explicitly temperature dependent (it is indeed temperature dependent through its volume dependence). This approximation seems to hold for TaC0.99 at temperatures above 600 K.

The result of the fit to our 673 K and 1073 K isotherms constrains an effective hightemperature value of the ambient pressure thermal Grüneisen parameter, which is more consistent with the results of high temperature ultrasonics measurements [30] rather than room temperature ultrasonics [30,31] or computations [36].

The value of 𝑞, the logarithmic volume derivative of the Grüneisen parameter, is poorly constrained by our high-temperature dataset. The best-fit value of 𝑞 is strongly affected by the highest pressures data along the 1073 K isotherm. The fit to our full dataset yields 𝑞 = -0.3  1.4, while excluding the non-ideal data points at 32.9, 34.9 36.1, 37.9 and 38.8

GPa at 1073 K from the fit yields 𝑞 = 1  2 (see also results). The average RMS misfit of the full dataset (50 data points) is 0.46 GPa, that of the non-ideal 5 data points is 0.59

GPa while that of all the others is 0.44 GPa. The non-ideal data points have a 34% higher disagreement with the fit than all the rest of the dataset. We interpret this as a systematic bias probably connected to the non-ideal geometry in these measurements, and we consider the result of the analysis of the reduced dataset as our best fit result. The fit results are instead unaffected by including or excluding the 3 non-ideal data points at 6.0, logarithmic volume sensitivity of 𝛾 in the high temperature regime (T > 600 K). This corresponds to:

𝛾(𝑉) ≅ 𝛾 0 𝐻𝑇 (𝑉 𝑉 0 ⁄ ), (4)
where 𝛾 0 𝐻𝑇 is the high-temperature, ambient pressure value of 𝛾 = 1.2  0.1 based on our fit (see discussion above), 𝑉 represents volume and 𝑉 0 the volume at ambient conditions.

In conclusion, our P-V-T model is able to successfully describe the pressurevolumetemperature relationship for TaC0.99 at 296 K and in the high-temperature range between 673 and 1073 K at pressures as high as 35 GPa, while we expect that it fails to describe the high-pressure behavior of TaC0.99 in the 300-600 K temperature regime, especially at low pressures. The parameters for the high-temperature equation of state of TaC0.99 are reported in Table III. The extant published results for TaCx are reported for comparison in Tables IV and V.

V. CONCLUSIONS

Our extensive experimental investigation gives new information about the compression behavior of TaC0.99 under very high pressures and temperatures in the diamond anvil cell.

a) We have constrained the parameters of the P-V-T equation of state of TaC0.99 to 1073 K and 32 GPa. The ambient pressure derivative of the bulk modulus is substantially higher (50%) with respect to previous studies at the same composition. The thermodynamic Grüneisen parameter refined from our results is lower than the available data at ambient conditions, and it seems to be consistent with the experimental high temperature limit. The approach used here of determining isothermal EOS data from xray diffraction data collected in the radial diffraction geometry under nonhydrostatic conditions is novel and shows the potential of this technique, in particular the use of ceramic gaskets that are more stable at higher temperatures and the lack of a hydrostatic compression medium. Thus, the technique is ideally suited to characterize the thermalelastic properties of hard materials such as transition metal carbides. indicates experimental studies, "DFT" indicates computational studies using density functional theory.

Figure 7. Temperature dependence TaC0.99 thermodynamic Grüneisen parameter at ambient temperature. The results from our P-V-T are compared with the other available results from experimental and computational studies.

Tables titles

Table I. Experimental condition of the four X-ray diffraction datasets. "Cold compr." is the maximum sample pressure before starting heating.

Table II. Summary of the all the experimental results from all the three runs. "stdev P + V-TaC" is the complete uncertainty on pressure including the propagation of the uncertainty on the unit-cell volume of TaC0.99. "<Q>" is the isotropic aggregate value of Q, a measure of the strain due to the deviatoric component of the stress applied to the sample (see main text). "PV fit" indicates whether the data points were used or not in the isothermal (296 K) P-V equation of state fit. "PVT fit" indicates whether the data points were used or not in the P-V-T equation of state fit.

2 )

 2 of the 3 rd generation light source PETRA III at DESY, Hamburg, Germany. The monochromatic incident Xray radiation had a wavelength of 0.2966 Å. The beam was focused to a spot size of 3m (vertical) by 8 m (horizontal) full width at half maximum by means of Compound Refractive Lenses (CRL). X-ray diffraction images were collected with a Perkin Elmer XRD1621 flat panel detector at 420 mm distance from the sample. The sample to detector distance and tilting angles of the detector were calibrated using a CeO2 standard from the National Institute of Standards and Technology (NIST, 674b). Twenty X-ray diffraction images were collected at pressures between 5.1 GPa and 43.6 GPa, including 3 images collected during decompression. The exposure time for each image was 20 s.

  Figure 2(a). The unrolled sequence of spectra integrated at 5 degree azimuthal-angle

6.8 and 8 . 2

 82 GPa along the 673 K isotherm. Based on the overall fit of the model isotherms at 673 K and 1073 K, a value of 𝑞 close to unity is a good approximation of the 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.5115350

Figure 1 .

 1 Figure 1. X-ray diffraction of TaC0.99 compressed in Ne pressure transmitting medium to

Figure 2 .

 2 Figure 2. Radial X-ray diffraction of TaC0.99 compressed without pressure transmitting

Figure 3 .

 3 Figure 3. Example of radial X-ray diffraction measurement in non-ideal sample

Figure 4 .

 4 Figure 4. Isothermal compression of TaC0.99 at ambient temperature (296 K). The dataset

Figure 5 .

 5 Figure 5. P-V-T compression behavior of TaC0.99. The data used for fitting are colored as

Figure 6 .

 6 Figure 6. Comparison of the 296 K P-V equation of state of TaC0.99 with the other

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/1.5115350

  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.

	1 m was loaded in 50 m diameter cylindrical
	chambers laser-drilled in amorphous boron -epoxy gaskets of 380 m diameter and 50
	m thickness. With such a narrow sample-chamber the entire sample remains as close as
	possible to the center of the diamond culets, where stress gradients are low and where

stress and strain fields are close to ideal axial symmetry that is required in the interpretation of X-ray diffraction experiments in radial scattering geometry (see section

Table III .

 III Best fit parameters of the P-V-T equation of state of TaC0.99 based on our experimental data.TableIV. Comparison of our best-fit thermoelastic parameters for TaC0.99 and all the available literature data. "E" is the isotropic aggregate Young's modulus. "Poisson's" is the isotropic aggregate Poisson's ratio. All the parameters refer to ambient conditions.  The elastic coefficients are not presented in the original paper. Parameter fixed in the analysis of the data.

*

Table V .

 V Summary of the available results about the thermal properties of TaCx from the literature. 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

	642	Table I
		Experiment	TaC-Ne	TaC-HT1	TaC-HT2
		Diffraction geometry	Axial	Radial	Radial
		X-ray Wavelength (Å)	0.2966	0.4847	0.4847
		DAC type	Symmetric P-C	modif. Mao-Bell	modif. Mao-Bell
		diamond anvil culet (m)	300	300	300
		P transmitting Medium	Ne	no	no
		P calibrant	Au, ruby	Au	Au
		Temperature (K)	296	673	1073
		N data points	20	22	25
		Max P (GPa)	43.6	33.6	38.8
		Exp time (s)	20	4	4
		Max membrane P (bar)	28.4	58.9	65.0
		Exp duration (h)	3	5	6
		Cold compr. (GPa)		2.4	2.6
	643

CITE THIS ARTICLE AS DOI: 10.1063/1.5115350

  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

	644	Table II
						stdev P + V
		Experiment	T (K)	P (GPa)	V (Å 3 )	TaC (GPa)	<Q>	PV fit	PVT fit
		TaC-Ne	296	5.10.1 87.060.01	0.19	--	yes	yes
		TaC-Ne	296	6.40.2 86.710.02	0.31	--	yes	yes
		TaC-Ne	296	10.30.1 85.790.02	0.22	--	yes	yes
		TaC-Ne	296	11.60.1 85.470.02	0.20	--	yes	yes
		TaC-Ne	296	14.50.3 84.810.02	0.35	--	yes	yes
		TaC-Ne	296	17.30.5 84.210.04	0.67	--	yes	yes
		TaC-Ne	296	19.70.6 83.820.04	0.82	--	yes	yes
		TaC-Ne	296	22.20.6 83.350.04	0.84	--	yes	yes
		TaC-Ne	296	24.60.5 82.910.04	0.74	--	yes	yes
		TaC-Ne	296	27.10.4 82.430.04	0.66	--	yes	yes
		TaC-Ne	296	29.50.5 82.020.04	0.76	--	yes	yes
		TaC-Ne	296	31.60.5 81.640.01	0.60	--	yes	yes
		TaC-Ne	296	34.10.6 81.210.01	0.64	--	yes	yes
		TaC-Ne	296	36.40.8 80.830.01	0.82	--	yes	yes
		TaC-Ne	296	38.40.6 80.480.01	0.66	--	yes	yes
		TaC-Ne	296	40.80.8 80.070.01	0.81	--	yes	yes
		TaC-Ne	296	43.60.7 79.630.01	0.77	--	yes	yes
		TaC-Ne	296	38.90.8 80.470.01	0.83	--	no	yes
		TaC-Ne	296	32.00.7 81.720.02	0.68	--	no	yes
		TaC-Ne	296	18.81.0 84.110.04	1.06	--	no	yes
		TaC-Ne	296	0.0001 88.4780.003	0.001	--	yes	yes
		TaC-HT1	67310	7.40.3 86.950.10	0.68	0.0009	no	yes
		TaC-HT1	67410	6.00.3 87.100.11	0.70	0.0010	no	no
		TaC-HT1	67310	6.80.3 86.890.09	0.64	0.0011	no	no
		TaC-HT1	67310	8.20.3 86.680.09	0.66	0.0012	no	no
		TaC-HT1	67310	9.40.3 86.460.11	0.75	0.0014	no	yes
		TaC-HT1	68210	11.00.3 86.170.10	0.72	0.0015	no	yes
		TaC-HT1	67310	12.60.3 85.710.11	0.78	0.0018	no	yes
		TaC-HT1	67310	15.30.4 85.130.10	0.81	0.0022	no	yes
		TaC-HT1	67310	17.40.4 84.610.11	0.89	0.0024	no	yes

CITE THIS ARTICLE AS DOI: 10.1063/1.5115350

  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.

	Table III
	This study	Parameter 1 stdev
	V0 (Å 3 )	88.478	0.020
	KT0 (GPa)	305	5
	(dKT /dP) T0	6.1	0.5
	 th0	1.2	0.1
	Q	1	2
	0 (K)	567	10

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/1.5115350

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.5115350

ACKNOWLEDGMENTS

Part of this research was supported through the projects "GeoMaX" funded under the Emmy-Noether Program of the German Science Foundation (MA4534/3-1) as well as grant MA4534/4-1. HM acknowledges support from the Bavarian Academy of Sciences.

We acknowledge technical assistance by A. Ehnes and I. Schwark from the sample environment group at PETRA III within the framework of the Extreme Conditions Science Infrastructure ECSI). We thank J. Mainberger for preparing the carbides samples for the carbon analysis and Martin Harms for proofreading and reformatting the reference list.


Zohra Kader

Antoine Girard

Reference Tracking for Linear Time Invariant Systems with a Relay Control *

Zohra Kader 1 and Antoine Girard 1

Abstract-This paper investigates the problem of tracking a reference trajectory for Linear Time Invariant (LTI) systems using a relay feedback. First of all, the problem of a reference tracking is reformulated as a stabilization problem of the tracking error. Then, under some conditions on the reference inputs a design approach of a relay feedback controller that guarantees a perfect tracking of the reference trajectory is provided. Moreover, conditions in LMI form are proposed in order to synthesize a relay feedback and provide an estimation of the domain of attraction. Secondly, the case where the reference inputs can take the same values as the control input is considered. In this context, achieving a perfect tracking is no longer possible. Thus, using the comparison principle, a relay feedback controller that guarantees a practical stabilization of the tracking error is proposed. An estimation of the domain of attraction along with the bound of the tracking error are equally provided. Finally, computer simulations show the efficiency of the developed method over a numerical example.

I. INTRODUCTION

Relay feedback controllers are omnipresent in different application fields [START_REF] Gonc ¸alves | Global stability of relay feedback systems[END_REF], [START_REF] Johansson | Relay feedback and multivariable control[END_REF], [START_REF] Í¡sypkin | Relay control systems[END_REF], [START_REF] Utkin | Sliding mode control in electromechanical systems[END_REF]. They are known to be simple, robust, and able to emulate locally the behavior of a linear static feedback [START_REF] Edwards | Sliding mode control: theory and applications[END_REF]. These properties make them an interesting substitute to continuous controllers -see for instance [START_REF] Gonc ¸alves | Global stability of relay feedback systems[END_REF], [START_REF] Hetel | Variable structure control with generalized relays: A simple convex optimization approach[END_REF]. Despite these advantages, relay feedback controllers have been shown to exhibit some difficulties and undesirable phenomena [START_REF] Johansson | Limit cycles with chattering in relay feedback systems[END_REF], [START_REF] Johansson | Fast switches in relay feedback systems[END_REF]. Indeed, it has been shown in the literature that relay feedback controllers present sliding modes, chattering and limit cycles. The presence of these complex behaviors make their study very challenging. However, these phenomena must not be neglected and should be considered in the design procedures. In particular, for systems with sliding modes the notion of system's solution must be reviewed in order to take into account the dynamics obtained by fast switching [START_REF] Cortes | Discontinuous dynamical systems[END_REF], [START_REF] Filippov | Differential equations with discontinuous righthand sides: control systems[END_REF]. Up to now, the investigations on relay feedback controllers design are limited to the stabilization problem to an equilibrium point. Frequency domain methods [START_REF] Boiko | Discontinuous control systems: frequency-domain analysis and design[END_REF] and LMI approaches [START_REF] Polyakov | Practical stabilization via relay delayed control[END_REF], [START_REF]On practical stabilization of systems with delayed relay control[END_REF] have been used for this purpose.

Here, we are interested on designing relay feedback controller that enforces the trajectories of an LTI system to track a reference trajectory. The proposed approach is inspired from the results in [START_REF] Hetel | Variable structure control with generalized relays: A simple convex optimization approach[END_REF] and [START_REF] Kader | Stabilization by a relay control using non-quadratic lyapunov functions[END_REF]. First, we reformulate the {zohra.kader ,antoine.girard}@l2s.centralesupelec.fr.

problem of a reference tracking as a stabilization problem of the tracking error. Then, under some conditions on the reference inputs a design approach of a relay feedback controller that guarantees a perfect tracking of the reference trajectory together with an estimation of the domain of attraction is provided. Moreover, an LMI approach that allows the design of a stabilizing relay feedback while providing an estimation of the domain of attraction is proposed. Secondly, the case where the reference inputs can take the same values as the relay feedback is considered. In this context, using the comparison principle, the designed relay feedback guarantees a practical stabilization of the tracking error. An estimation of the domain of attraction and the bound of the tracking error are equally provided. Finally, computer simulations are performed in order to show the efficiency of the developed method.

The paper is organized as follows: Section II presents the system description, the stability notions, and the problem under study. Under some restrictions on the reference inputs, a design approach of a relay feedback controller that ensures the perfect tracking of the reference trajectory is provided in Section III. In addition, a constructive method based on LMIs is proposed in the same section. In Section IV, using the comparison principle a relay feedback controller that ensures the boundedness of the closed-loop system trajectories is provided. A numerical example is given to illustrate the efficiency and the limits of the presented method. Finally, concluding remarks are provided in the last section.

A. Notations

In this paper we use the notation R + to refer to the interval [0, ∞). The transpose of a matrix M is denoted by M T and if the matrix is symmetric the symmetric elements are denoted by * . The notation M 0 (resp. M 0) means that the matrix M is positive (resp. negative) semi-definite, and the notation M 0 (resp. M ≺ 0) means that it is positive definite (resp. negative definite). M (i) refers to the i-th row of a matrix (or a vector).

The identity matrix is denoted by I and both notations eig min (M ) and eig max (M ) are used to refer to the minimum and maximum eigenvalue respectively of a matrix M . For a positive definite matrix M ∈ R n×n and a positive scalar γ, we denote by E(M, γ) the ellipsoid

E(M, γ) = {x ∈ R n : x T M x ≤ γ}, (1)
and for all positive scalar r, we denote by B(0, r) the ball

B(0, r) = E(I, r) = {x ∈ R n : x ≤ r}, (2)
where . is used to denote the Euclidean norm for a vector and the associated norm for a matrix. For a given set S, the notation Conv{S} indicates the convex hull of the set, int{S} its interior and S its closure and finally the closed convex hull of the set S is denoted by Conv{S}. The minimum argument of a given function f : S -→ R such that the set S ⊂ R is a finite set of vectors is denoted by

arg min f = {y ∈ S : f (y) ≤ f (z), ∀z ∈ S}. (3)
For a positive integer N , we denote by I N the set {1, . . . , N }. By ψ N we denote the unit simplex

ψN = α = (α1, . . . , αN ) T ∈ R n : N i=1 αi = 1, αi ≥ 0, i ∈ IN . (4)

II. SYSTEM DESCRIPTION AND PROBLEM STATEMENT

A. System description Consider the linear system

ẋ = Ax + Bu, (5)
with x ∈ R n and an input u which takes values in the set

V = {v 1 , . . . , v N } ⊂ R m . A ∈ R n×n
, and B ∈ R n×m are the matrices describing the system. This paper deals with the design of a relay feedback controller that enforces the system state to asymptotically track the reference signal x r which is given by the reference model

ẋr = Ax r + Bu r , (6)
where u r is a piece-wise constant and right continuous function from R + to the finite set of vectors V r = {v r 1 , v r 2 , . . . , v r Ñ }, with a finite number of discontinuities on every bounded interval of R + , A and B are defined in [START_REF] Cortes | Discontinuous dynamical systems[END_REF]. Note that the number of the reference input vectors Ñ > 0 can be equal to N .

Given the desired reference x r from (6) and the system model ( 5), the open loop tracking error system is

ė = Ae + Bu -Bu r , u ∈ V, u r ∈ V r , e := x -x r . (7)
Due to the formulation [START_REF] Filippov | Differential equations with discontinuous righthand sides: control systems[END_REF], the tracking problem can be restated as the design of a relay feedback controller that asymptotically stabilizes the error system to the origin. To this end, in the sequel we assume that: A-1 The pair (A, B) is stabilizable, which means that there exists a matrix K such that the closed-loop matrix

A cl = A + BK is Hurwitz. A-2
The set Conv{V} is nonempty, and the null vector is contained in its interior (0 ∈ int{Conv{V}}). Here, we consider a controller given by u(e) ∈ arg min

v∈V e T Sv, (8)
where the matrix S ∈ R n×m characterizes the switching hyperplanes. The formulation (8) can be interpreted as networked control systems with quantization [START_REF] Nesic | A unified framework for design and analysis of networked and quantized control systems[END_REF]. Moreover, the formulation of the controller (8) encompasses the classical sign function in the classical relay feedback. Therefore, if

V = {v 1 , v 2 } = {-v, v} with v > 0 then we get u(e) = -vsign(S T e) ∈    v if S T e < 0, {-v, v} if S T e = 0, v if S T e > 0. (9)
The closed-loop system is modeled by a differential equation with a discontinuous right hand-side. Consequently, to study the stability of the system ( 7), (8) we will consider the Filippov solutions of differential inclusions. Next, we provide the definitions of differential inclusions and their solutions [START_REF] Aubin | Differential inclusions: set-valued maps and viability theory[END_REF], [START_REF] Bacciotti | Liapunov functions and stability in control theory[END_REF], [START_REF] Cortes | Discontinuous dynamical systems[END_REF], and [START_REF] Filippov | Differential equations with discontinuous righthand sides: control systems[END_REF].

B. Solution concept and stability notions

The interconnection ( 7), ( 8) is the closed-loop system modeled by a discontinuous differential equation of the form

ė = Ae + Bu(e) -Bu r = f (t, e). ( 10
)
Therefore, to the discontinuous closed-loop system [START_REF] Johansson | Relay feedback and multivariable control[END_REF] we associate the differential inclusion

ė ∈ F(t, e), (11)
with F(t, e) the set-valued map which can be computed from the differential equation with a discontinuous right hand-side using the construction given by :

F (t, e) = δ>0 µ(S)=0 Conv{f (t, B(e, δ))\ S}, e ∈ R n , t ∈ R + , ( 12
)
where Conv is the closed convex hull, B(e, δ) is the open ball centered on e with radius δ, and S is a set of measure zero with µ(S) its measure in the sens of Lebesgue -see for instance [START_REF] Bacciotti | Liapunov functions and stability in control theory[END_REF], [START_REF] Cortes | Discontinuous dynamical systems[END_REF], [START_REF] Filippov | Differential equations with discontinuous righthand sides: control systems[END_REF], [START_REF] Poznyak | Attractive Ellipsoids in Robust Control[END_REF], and the references therein. The closed-loop system is then modeled by a differential inclusion for which the notion of a solution was defined in [START_REF] Filippov | Differential equations with discontinuous righthand sides: control systems[END_REF], and recalled in the following. Definition 1: (Filippov solution) Consider the closed-loop system [START_REF] Johansson | Relay feedback and multivariable control[END_REF] and its associated differential inclusion [START_REF] Johansson | Limit cycles with chattering in relay feedback systems[END_REF]. A Filippov solution of the discontinuous systems [START_REF] Filippov | Differential equations with discontinuous righthand sides: control systems[END_REF], [START_REF] Gonc ¸alves | Global stability of relay feedback systems[END_REF]

over the interval [t a , t b ] ⊂ [0, ∞) is an absolutely continuous mapping y : [t a , t b ] -→ R n satisfying: ẏ(t) ∈ F(t, y(t)), for almost all t ∈ [t a , t b ], (13)
with F(t, e) given by [START_REF] Johansson | Fast switches in relay feedback systems[END_REF].

A differential inclusion has at least one solution if the set valued map F(t, e) is locally bounded and takes nonempty, compact and convex values [START_REF] Aubin | Differential inclusions: set-valued maps and viability theory[END_REF], [START_REF] Bacciotti | Liapunov functions and stability in control theory[END_REF], [START_REF] Cortes | Discontinuous dynamical systems[END_REF], [START_REF] Filippov | Differential equations with discontinuous righthand sides: control systems[END_REF]. Definition 2 (Equilibrium point): e eq is said to be an equilibrium point of the differential inclusion [START_REF] Johansson | Limit cycles with chattering in relay feedback systems[END_REF] if 0 ∈ F(t, e eq ) for all t ≥ t 0 ≥ 0. Hereafter the notions of stability that we will use are introduced.

Definition 3 (Local uniform exponential stability): The equilibrium point e eq = 0 of the differential inclusion ( 11) is said to be locally uniformly exponentially stable with a decay rate α, if there exist a set D ⊂ R n , 0 ∈ int{D}, and strictly positive scalars κ independent of t 0 , and α such that for all Filippov solutions e(t) of ( 11) with e(t 0 ) ∈ D, e(t) ≤ κexp(-α(t -t 0 )) e(t 0 ) , ∀t ≥ t 0 ≥ 0. (14) A set D satisfying these properties is usually called an estimation of the domain of attraction.

We recall that sufficient conditions for the local uniform exponential stability of Filippov solutions in the case of systems modeled by differential equations with a discontinuous right hand-side ė = f (t, e), are given by the existence of a strict Lyapunov function V , V (t, 0) = 0, V (t, e) ≥ 0, ∀e = 0, t ≥ t 0 ≥ 0 and continuous positive definite functions W 1 and W 2 (i.e. W i (0) = 0 and W i (e) > 0 for all e = 0 and i ∈ {1, 2}) (see [START_REF] Filippov | Differential equations with discontinuous righthand sides: control systems[END_REF], Chapter 3, page 153, Theorem 1), such that

W 1 (e) ≤ V (t, e) ≤ W 2 (e), (15)
and

sup y∈F (t,e) ∂V ∂t + ∂V ∂e y ≤ -2αV (t, e), ∀e ∈ D\{0}, ∀t ≥ 0, (16)
for some domain D such that 0 ∈ D and some positive scalar α.

We adapt as follows the concept of practical stability from [START_REF] Retchkiman | Stabilization of nonlinear systems using vector lyapunov functions[END_REF] to the context under study.

Definition 4 (Practical stability): System ( 11) is said to be

• Practically stable, if given (R 1 , R 2 ) with 0 < R 1 < R 2
then e(t 0 ) < R 1 ⇒ e(t) < R 2 , ∀t ≥ t 0 ≥ 0; • Uniformly practically stable, if it is practically stable for all t 0 ≥ 0. For more details about practical and exponential stability, one can refer to [START_REF] Filippov | Differential equations with discontinuous righthand sides: control systems[END_REF], [START_REF] Khalil | Nonlinear systems[END_REF].

III. CASE 1. PERFECT REFERENCE TRACKING

A. Control design

In this section, we consider the case where a perfect tracking of the reference trajectory can be achieved by a relay feedback. To this end, let us make the following additional assumption: A-3 The set of reference inputs Conv{V r } is contained in the interior of the convex set Conv{V} i.e., Conv{V r } ⊂ int{Conv{V}}.

Theorem 1: Assume that A.1 and A.3 hold. Then there exist a matrix S such that the closed-loop system (7), ( 8) is locally uniformly exponentially stable with a decay rate α in a domain that contains the origin.

The control design approach proposed in Theorem 1 can be used to provide a constructive method of relay feedback controllers that enforce system (5) to track the reference trajectory [START_REF] Edwards | Sliding mode control: theory and applications[END_REF]. In the next section, using a quadratic Lyapunov function we propose LMIs (Linear Matrix Inequalities) criteria for the design of relay controllers while providing an estimation of the domain of attraction.

B. LMIs solution

The result of Theorem 1 has a qualitative nature. In practice, it is useful to find a constructive procedure which provides a switching law which ensures the stability of the system [START_REF] Filippov | Differential equations with discontinuous righthand sides: control systems[END_REF] together with an estimation of the domain of attraction. Here, we provide a numerical approach to deal with this problem. More precisely, an LMI solution is proposed hereafter. In order to express the result note that for any set V there exists a finite number n h of vectors

h i ∈ R 1×m , i ∈ I n h such that Conv{V} = {u ∈ R m : h i u ≤ 1, i ∈ I n h }. ( 17
)
Theorem 2: Assume that A.1 and A.3 hold. Consider the linear closed-loop system ( 7), ( 8) and a tuning parameter δ > 0. If there exist a positive definite matrix Q ∈ R n×n and a matrix Y ∈ R m×n such that the LMIs

QA T + AQ + Y T B + BY ≤ -2δQ (18)
1 hi,j Y Y T hT i,j Q ≥ 0, ∀i ∈ I n h , ∀j ∈ I Ñ , (19)
with hi,j = hi 1-hiv r j , are feasible then the origin of system (7) with the switching law u(e) ∈ arg min

v∈V e T Q -1 Bv, (20)
is locally uniformly exponentially stable with a decay rate δ. Remark 1: The feasibility of the set of LMIs ( 18)-( 19) allows the design of the matrix S = M B that characterizes the switching hyperplane and provides an estimation E(M, 1) of the domain of attraction. More precisely, any solution of the closed-loop system [START_REF] Filippov | Differential equations with discontinuous righthand sides: control systems[END_REF], [START_REF] Gonc ¸alves | Global stability of relay feedback systems[END_REF] starting in the domain of attraction E(M, 1) converges to the origin exponentially with a decay rate δ. In this case, the trajectories of system (5) will converge to (or track) the trajectories of the reference system [START_REF] Edwards | Sliding mode control: theory and applications[END_REF]. In this section, we have shown that if the reference input u r satisfies Assumption A.3 then we are able to design a relay controller ensuring that the tracking error will tend to zero as t tends to infinity. When Assumption A.3 does not hold, only practical stabilization of the tracking error is possible. This case is considered in the next section.

IV. CASE 2. PRACTICAL STABILIZATION OF THE TRACKING ERROR

Here, we are interested in the case where Assumption A.3 does not hold i.e, the reference imputs u r can take values in Conv{V}. In this case, one can only achieve practical stability of the error system (7) by the relay controller [START_REF] Gonc ¸alves | Global stability of relay feedback systems[END_REF] i.e., one is interested on specific trajectory bounds of ( 7), [START_REF] Gonc ¸alves | Global stability of relay feedback systems[END_REF]. In order to provide stabilization criteria of system [START_REF] Filippov | Differential equations with discontinuous righthand sides: control systems[END_REF], we use the so-called comparison Lemma [START_REF] Khalil | Nonlinear systems[END_REF], [START_REF] Retchkiman | Stabilization of nonlinear systems using vector lyapunov functions[END_REF]. In order to build the comparison system, we consider that Assumptions A.1 and A.2 hold. Therefore, thanks to the fact that the pair (A, B) is stabilizable, there exists a gain K such that A cl = A + BK is Hurwitz. Furthermore, for all δ > 0 there exists a matrix M = M T 0 verifying

A T cl M + M A cl -2δM. ( 21
)
Moreover, since the set Conv{V} is nonempty and the null vector is contained in its interior (0 ∈ int{Conv{V}}), then there exists a neighbourhood of the origin E(M, γ) with γ > 0 such that for all e ∈ E(M, γ) we have

Ke ∈ Conv{V}. ( 22
)
Thus, for all e ∈ E(M, γ) there exist positive scalars α j (e), j ∈ I N , such that

η = -λη + πu T r u r , η(t 0 ) ≥ 0, ( 24
)
where π is a positive real satisfying π > ξmax 2δ with δ > 0 is defined in [START_REF] Utkin | Sliding mode control in electromechanical systems[END_REF] and ξ max > 0 satisfies the inequality e T M BB T M e ≤ ξ max e T M e for all e ∈ R n , and λ = 2δ -π -1 ξ max > 0. Moreover, for all t ≥ t 0 , η(t) ≤ ε ≤ γ where ε > 0 is such that KE(M, ε) ⊂ Conv{V}. The comparison system (24) is utilized in order to prove the uniform practical stability of the closed-loop system [START_REF] Filippov | Differential equations with discontinuous righthand sides: control systems[END_REF], [START_REF] Gonc ¸alves | Global stability of relay feedback systems[END_REF]. This result is given in the following theorem.

Theorem 3: Consider system ( 7), (8) under Assumptions A.1, A.2 and A.4. Then, there exist a positive scalar γ and matrices M = M T 0 and S = M B such that system (7), ( 8) is uniformly practically stable with R 1 = η(t0) eigmin(M ) and R 2 = ε eigmin(M ) . Remark 2: One may remark that: 1) the result in Theorem 3 means that the existence of a static feedback controller ũ = Ke and a comparison system given by Assumption A.4 implies the existence of a relay feedback that stabilizes system [START_REF] Filippov | Differential equations with discontinuous righthand sides: control systems[END_REF]. The designed relay feedback guarantees that the closedloop system has the same stability properties as the comparison system (24). 2) the result of Theorem 3 is constructive. Indeed, if the inequality ( 21) is satisfied and if there exists a comparison system defined as in Assumption A.4 then the closed-loop system (7), ( 8) is uniformly practically stable with a switching hyperplane given by S = M B. Inequality ( 21) can be easily reformulated as an LMI design condition [START_REF] Boyd | Linear matrix inequalities in system and control theory[END_REF]: there exist Q = Q T 0 and a matrix Y ∈ R m×n such that

QA T + AQ + Y T B + BY ≤ -2δQ (25)
with Q = M -1 and K = Y M . Then, the switching hyperplane is characterized by S = M B. - V. ILLUSTRATIVE EXAMPLE Consider the linear system [START_REF] Filippov | Differential equations with discontinuous righthand sides: control systems[END_REF] with u ∈ V = {-v, v} = {-1, 1}, and matrices A = -1 0 0 0.1 and B = 1 1 .

A. Perfect tracking

First we consider that the reference input u r takes values inside the set V r = { -v 2 , v 2 } and is shown in Figure 3. Considering a decay rate δ = 1.5, the set of LMIs ( 18)- T and x r (0) = 0 0 T and are reported in Figures 12345.

As we can see from Figures 1 and2, the tracking error trajectories starting in the domain of attraction E(Q -1 , 1) contained in the convex Conv{V} (in red line) converge to the origin. From Figures 4 and5 we can observe that the trajectory of the closed-loop system converge to the reference trajectory. These observations are consistent with the theoretical development in Section III

B. Non-perfect tracking

Here, we consider that the reference input u r takes values in the set V r = {-1, 0, 1} and is shown in Figure 7. Considering a decay rate δ = 0.5, the LMI (25) is found feasible for The comparison system (24) is such that ξ max = 2.0777, η(0) = 1.5, π = 3, λ = 0.3074 and u r is shown in Figure 6. Thus, under this conditions η(t) ≤ 5 for all t ≥ 0. Computer simulations of the tracking error evolution are performed for initial states x(0) = -1.53 .02 T and x r (0) = 0 0 T and are reported in Figure 7-10. From Figures 7 and 10 we can observe that the trajectories of the tracking errors initialized such that V (e(0)) ≤ 1.5 remains bounded as t tends to infinity. The state trajectories get and remain closer to the reference states as it can be seen from Figures 8 and9. These observations are consistent with the theoretical result shown in Section IV.

VI. CONCLUSION

In this paper the problem of tracking a reference trajectory for Linear Time Invariant (LTI) systems using a relay feedback has been investigated. First, the problem of a reference tracking has been reformulated as a stabilization problem of the tracking error. Then, under some conditions on the reference inputs a design approach of a relay feedback controller that guarantees a perfect tracking of the reference trajectory together with an estimation of the domain of attraction has been provided. Moreover, conditions in LMI form have been proposed to compute a relay feedback and an estimation of the domain of attraction. Secondly, the case where the reference inputs can take the same values as the relay feedback has been considered. In this context, using the comparison principle, it has been shown that the designed relay feedback guarantees a practical stabilization of the tracking error. An estimation of the domain of attraction and the bound of the tracking error have been equally provided. Finally, computer simulations have been performed in order to show the efficiency of the developed method over a numerical example. In the future, the approach may be improved in various manners. For example, it is of interest to reduce the conservatism induced by Assumptions A-3 and A-4. In order to reduce the chattering, we can try to extend the approach to the case of min-switching strategy with mindwell time condition.

APPENDIX I PROOF OF THEOREM 1

Proof: Thanks to the fact that the pair (A, B) is stabilizable then there exists a gain K such that A cl = A + BK is Hurwitz. Furthermore, for all δ > 0 there exists a matrix M = M T 0 satisfying

A T cl M + M A cl -2δM. ( 26
)
Let us consider the quadratic Lyapunov function V (e) = e T M e. Our objective is to show that the differential inclusion (11) associated to the closed-loop system [START_REF] Filippov | Differential equations with discontinuous righthand sides: control systems[END_REF], [START_REF] Gonc ¸alves | Global stability of relay feedback systems[END_REF]

for some positive scalar α > 0 in some domain D ⊂ R n that will be determined later. From (31), and thanks to the fact that ψ * (z) is compact for all z ∈ R n , we have Therefore, there exists a neighbourhood of the origin E(M, γ) with γ > 0 such that for all e ∈ E(M, γ) we have Ke ∈ B(0, ε).

(37)

  us make the following assumption: A-4 There exists a function η : R + → R satisfying

Fig. 1 .

 1 Fig. 1. Evolution of the tracking error dynamics starting at e = [0.60.25] T .

Fig. 2 .

 2 Fig. 2. Estimation of the domain of attraction (red line); trajectory of the tracking error (black line); Switching hyperplane (blue line).

( 19 )

 19 are feasible for Q = 0.4070 0.1519 0.1519 0.0661 and Y = -0.2221 -0.1125 . Computer simulations of the tracking error evolution are performed for initial states x(0) = 0.6 0.25

Fig. 3 .Fig. 4 .Fig. 5 .

 345 Fig. 3. Reference input ur.

Fig. 6 . 2 Fig. 7 .Fig. 8 .

 6278 Fig. 6. Reference input ur.

Fig. 9 .Fig. 10 .

 910 Fig. 9. Evolution of the state variables x (2) and x r(2)

  α > 0 in a domain D ⊂ R n that we will determine later on. Now, let I * (z) be the set of index for any z ∈ R n such thatI * (z) = {i ∈ I N : z T M B(v j -v i ) ≥ 0, ∀j ∈ I N }. (28)To I * (z) we associate for all z ∈ R n the set ψ * (z) of vectors defined byψ * (z) = {β ∈ ψ N : β i = 0, ∀i ∈ I N \ I * (z)}.(29)Using (28) and (29), the set valued map F(t, e) in (12) satisfiesF(t, e) ⊆ F * (t, e),(30)withF * (t, e) = Conv i∈I * (e) {Ae + Bv i -Bu r } ={Ae + Bv(β) -Bu r : β ∈ ψ * (e)},(31)and v(β) = N i=1 β i v i . Thus, in order to show (27), it is sufficient to prove that sup y∈F * (t,e) ∂V ∂e y ≤ -2αV (e),

  sup y∈F * (t,e) ∂V ∂e y = sup β∈ψ * (e) ∂V ∂e (Ae + Bv(β) -Bu r ) = max β∈ψ * (e) ∂V ∂e (Ae + Bv(β) -Bu r ) . (33) Then, showing (32) is equivalent to demonstrating that for some α > 0 max β∈ψ * (e) ∂V ∂e (Ae + Bv(β) -Bu r ) ≤ -2αV (e), (34) in a domain D ⊂ R n to be determined below. From inequality (26), we obtain ∂V ∂e (A cl e) ≤ -2δV (e), ∀e ∈ R n . (35) Since Conv{V r } ⊂ int{Conv{V}} then there exists ε > 0 such that for all c(e) ∈ B(0, ε) and v r ∈ V r we have v r + c(e) ∈ Conv{V}. (36)

* This work has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 725144).

Thus, From (36) and (37), for all e ∈ E(M, γ) we have v r + Ke ∈ Conv{V}.

(38)

Therefore, for all e ∈ E(M, γ) there exist positive scalars α j (e), j ∈ I N , such that N j=1 α j (e) = 1 and

From now on, we consider the case where (35) and (39) are verified (i.e. for all e ∈ E(M, γ)). From (28), for all i ∈ I * (e) we have

Then, for any β ∈ ψ * (e) we have

Multiplying this last inequality by α j (e) for j ∈ I N (α j (e) defined in (39)) and summing the N elements, we obtain

for all v r ∈ V r . Adding this to the left part of (34), for all u r = v r ∈ V r and for all e ∈ E(M, γ) we get

From the last inequality and taking α = δ, (34) is verified for all u r = v r ∈ V r and for all e ∈ E(M, γ). Therefore, (26) is verified for all u r = v r ∈ V r and for all e ∈ E(M, γ).

APPENDIX II PROOF OF THEOREM 2

Proof: We want to prove that if LMIs ( 18)-( 19) are feasible then the closed-loop system (7), ( 8) is locally uniformly exponentially stable in some domain D. Note that LMI [START_REF] Poznyak | Attractive Ellipsoids in Robust Control[END_REF] guarantees that the function

with A cl = A + BK and K = Y M . Applying the Shur complement and taking into account the fact that

with hi,j = hi 1-hiv r j . Let us remark that the minimum of V along the hyperplane hi,j Ke = 1 is given by min hi,j Ke=1

e T M e = ( hi,j KQK T hT i,j ) -1 . From (45) we have that for all e ∈ E(M, 1), hi,j Ke ≤ 1, ∀i ∈ I n h , ∀j ∈ I Ñ . By the definition of hi,j , we obtain

Thus, the relation u r + Ke ∈ Conv{V} holds. Therefore, from (46), ( 44) is satisfied for all e ∈ E(M, 1). From now on, using the same argument as in Theorem (1) one is able to show that sup y∈F (t,e)

∂V ∂e y ≤ -2δV (e), ∀e ∈ E(M, 1).

APPENDIX III PROOF OF THEOREM 3 Proof: First, let us define for any z ∈ R n the set of index I * (z) such that

To I * (z) we associate for all z ∈ R n the set ψ * (z) of vectors defined by

Using ( 48) and ( 49), the set valued map F(t, e) in ( 12) satisfies

with

From (51), and using the fact that the set ψ * (z) is compact for all z ∈ R n , we have

From ( 48), for all i ∈ I * (e) we have

Then, for any β ∈ ψ * (e) we have

Multiplying this last inequality by α j (e) for j ∈ I N (α j (e) defined in (23)) and summing the N elements, we obtain

Adding this to the left part of (54), we get

Now, using (52), we obtain

Recall that for any positive number θ (63)

From Assumption A.4 and by the comparison principle [START_REF] Khalil | Nonlinear systems[END_REF] taking V (e(t 0 )) ≤ η(t 0 ) i.e., e(t 0 ) ≤ R 1 leads to

Thus, e(t) ≤ R 2 , ∀t ≥ t 0 ≥ 0,

which ends the proof.