Roman Matkovskyy
email: roman.matkovskyy@rennes-sb.com

Centralized and Decentralized Bitcoin Markets: Euro vs USD vs GBP

Keywords: JEL classification: G1, C14, C22 Bitcoin, centralized exchange, decentralized exchange, volatility, dependency, ARMA-GARCH copula JEL classification: G1, C14, C22 Bitcoin, centralized exchange, decentralized exchange, volatility, dependency, ARMA-GARCH copula

In this study, I compared the euro, U.S. dollar, and British pound sterling (GBP) centralized and decentralized bitcoin cryptocurrency markets in terms of return volatility and interdependency. This comparison showed the decentralized bitcoin market has higher volatility and the centralized markets have higher tail dependence regarding returns. The volatility analysis results are contrary to the established leverage reasons that market drops cause volatility. The results demonstrate a higher left tail dependence is in line with the general pattern in "traditional" financial markets which more extreme dependent in downturns. It was also shown trade volume increases as prices decrease, demonstrating participants' lack of confidence and consensus in a price-jump period.

Introduction

Bitcoin, a communication protocol used in virtual currency transactions [START_REF] Böhme | Bitcoin: Economics, Technology, and Governance[END_REF], has recently entered into and assumed a place in financial markets and portfolio management (Dyhrber, 2016). In December 2017, bitcoin was approved by the U.S. Commodities and Futures Trading Commission (CFTC), which made the Chicago Board Options Exchange (CBOE) and CME bitcoin futures contract launches possible. This approval opened up bitcoin trading to a larger group of investors and traders and thus allowed more international participants to become involved in this market. 1 Despite bitcoin being a homogeneous and identical virtual good across all online markets in which it is traded, its prices behave differently across these markets due to their characteristics' variations [START_REF] Pieters | Financial Regulations and Price Inconsistencies across Bitcoin Markets[END_REF], and high fluctuations in bitcoin prices have resulted in periods of high volatility.

There are two main bitcoin markets: centralized and decentralized markets (exchanges).

Centralized bitcoin exchanges include intermediaries such as companies that act as a proxy in order to facilitate trading. Unlike centralized exchanges, the decentralized market is a personto-person (P2P) bitcoin trading site that allows users to post advertisements indicating exchange rates and payment methods for buying or selling bitcoins without disclosing their identities. Trading in the decentralized bitcoin markets has disadvantages, however, because the centralized exchanges offer significantly better trading functionalities. But the major centralized exchanges are not immune to difficulties and have experienced delays and technical difficulties in transactions. For instance, Kraken and other centralized exchanges have halted deposits and withdrawals without warning at various times.

Recently, various aspects of cryptocurrencies, particularly of bitcoins, have drawn increasing research attention. For instance, [START_REF] Luther | Bitcoin and the bailout[END_REF] focus on bitcoin-related apps found downloads of these apps depended on crises occurring in banking systems. Bouri et al. (2017a) and Bouri et al. (2017b) studied the issues of hedging and safe-haven properties of cryptocurrencies. [START_REF] Urquhart | Price clustering in Bitcoin[END_REF] documented that price and volume have a significant positive relationship with price clustering in bitcoin, and [START_REF] Balcilar | Can volume predict Bitcoin returns and volatility? A quantiles-based approach[END_REF] showed that volume cannot be used to predict the volatility of bitcoin returns at any point of the conditional distribution. [START_REF] Katsiampa | Volatility estimation for Bitcoin: A comparison of GARCH models[END_REF], based on goodness-of-fit tests, showed that the AR-CGARCH model can better explain bitcoin price volatility. Nadarajah and Chu (2017) and [START_REF] Urquhart | The inefficiency of Bitcoin[END_REF] demonstrated the high inefficiency of bitcoin, and [START_REF] Charfeddine | Are shocks on the returns and volatility of cryptocurrencies really persistent?[END_REF] also confirmed the inefficiency of the bitcoin market. Phillip et al. (2018) showed in their study of 224 cryptocurrencies that they have some common features (i.e., leverage effects and Student-t error distributions). [START_REF] Blau | Price dynamics and speculative trading in Bitcoin[END_REF] demonstrated that speculative trading is not associated with higher levels of volatility in bitcoin, and in analyzing bitcoin volatility, [START_REF] Conrad | Long-and short-term cryptocurrency volatility components: A GARCH-MIDAS analysis[END_REF] found it is closely linked to global economic activity. These papers, however, mainly focused on centralized U.S. dollar bitcoin markets. The empirical comparison of return volatility and the correlation of the centralized and decentralized bitcoin markets using different currencies has been largely ignored. Such research could be potentially useful in reducing risks in investment decision making. Thus, in this study, I fill this gap by comparatively analyzing volatility and correlation among the main euro, U.S. dollar, and British pound bitcoin exchanges2 as well as between the centralized (e.g., Kraken, GDAX-CoinsBank,3 bitcoin.de, Bitstamp, Rock Trading, and the exchangebacked peer-to-peer marketplace Coinfloor) and decentralized (P2P) exchanges (e.g., LocalBitcoins) to contribute to a better understanding of the cryptocurrency market's behavior.

For this purpose, I analyzed dependences across the selected markets by applying tail copulas and modeling the dynamics of the average values and variances of the time series using multivariate models. These included the generalized autoregressive conditional heteroskedasticity (GARCH) model (on the residuals from the autoregressive-movingaverage [ARMA] models), a multivariate exponentially weighted moving average (EWMA) for the covariance matrix, and a multivariate model based on copulas that uses the residuals of the ARMA-GARCH process.

The findings provide the first empirical evidence on the dynamics of return volatility and interdependence of the euro, U.S. dollar, and British pound centralized and decentralized bitcoin markets. These were analyzed jointly and, thus, contribute to a better understanding of behavior of the bitcoin centralized and decentralized markets that use different currencies.

The results can help shed light on timely issues researchers, practitioners, and policymakers face in terms of risk reduction in bitcoin investment decision making.

The results extend the findings of [START_REF] Blau | Price dynamics and speculative trading in Bitcoin[END_REF], [START_REF] Katsiampa | Volatility estimation for Bitcoin: A comparison of GARCH models[END_REF], and [START_REF] Klein | Bitcoin is not the New Gold -A comparison of volatility, correlation, and portfolio performance[END_REF] by showing decentralized markets have higher volatility and centralized markets have higher tail dependence in terms of returns. The volatility analysis results are contrary to the accepted leverage reason that market drops cause volatility. The results demonstrate the decentralized bitcoin markets are more interdependent, but a price setting in the P2P bitcoin markets is not correlated over time or across markets. Correlation between markets is not stable over time and is higher for the centralized U.S. dollar and British pound exchanges.

The estimates show the centralized bitcoin markets demonstrate a higher left tail dependence that is in line with the general pattern of the "traditional" financial markets, which are more extreme dependent in downturns.

Analysis also revealed trade volume increases as prices decrease, demonstrating participants' lack of confidence and consensus in the price-jump period. The existing pattern can also mean that based on their sentiments, investors may overreact because they do not know how the precise technique to estimate the bitcoin price.

Based on this thinking, the obtained results are potentially useful to researchers, practitioners, and participants in the highly volatile bitcoin market to make better investment and risk management decisions.

The structure of this study is as follows. The models used are presented in Section 2, the model results are discussed in Section 3, and the study conclusions are offered in Section 4. The appendix furnishes supplementary material.

Methodology

The GARCH framework is frequently used as an efficient instrument to measure the interdependence of returns and volatilities across stock markets, and in a growing stratum of literature, these empirical models and measurements are also applied to bitcoin volatility analyses. [START_REF] Katsiampa | Volatility estimation for Bitcoin: A comparison of GARCH models[END_REF] Some recent studies on bitcoin have focused on distributions, which allow for asymmetry as well as fat tails. For example, [START_REF] Bouri | Bitcoin and global financial stress: A copula-based approach to dependence and causality in the quantiles[END_REF] applied a copula-based approach to document a right tail dependence between the global financial stress index and bitcoin returns. They also showed limited directional predictability from the global financial stress index regarding bitcoin returns in the medium term.

In this study, I used well-established models (i.e., tail copulas, GARCH models [on the residuals from the ARMA models], a multivariate EWMA for the covariance matrix, and a multivariate model based on copulas that uses the residuals of the ARMA-GARCH process). I start with a kernel density estimation to depict the distributions. The determination of the optimal bandwidth was based on choosing the one that minimizes an optimality criterion (e.g., the mean integrated squared error).

The multivariate EWMA models for the covariance matrix were then built to show the volatility of returns across different markets. Among the attractive properties EWMA has over GARCH modelling is in the greater weight put upon more recent observations. The EWMA model was proposed in [START_REF] Riskmetrics | Riskmetrics Technical Document[END_REF], where the variances and covariances are as in the IGARCH-type models. In general, this model is

Σ = Σ + 1 - ′
, where are residual from the mean equation, and λ (0<λ<1) is the persistence parameter, which equals 0.96 for daily observations [START_REF] Riskmetrics | Riskmetrics Technical Document[END_REF].

Correlation was tested by different measurements, namely tail dependence, tail copulas, and Spearman's rho. The well-established tail dependence coefficients (TDC) based on the concept of copulae were used as the scalar measures for correlation in tails (see, for instance, [START_REF] Frahm | Estimating the tail-dependence coefficient: Properties and pitfalls[END_REF][START_REF] Schmidt | Non-parametric Estimation of Tail Dependence[END_REF]. They showed how extreme values in one series are likely to be accompanied by equally large or low extreme values in other series. A tail copula is defined as a function that explains the dependence structure of joint distributions in lower (left) or upper (right) tails. Let F be a distribution function with corresponding copula C. Then, the function Λ : ℝ → ℝ is called a lower (left) tail copula associated with F, if the following limit exists everywhere on ℝ :

Λ , : = lim → tC x # ⁄ , y # ⁄ (1)
Analogously, the function Λ & : ℝ → ℝ is defined as an upper tail copula if the following limit exists on ℝ :

' & , : = ()* → #+ ̅ # ⁄ , # ⁄ (2)
According to [START_REF] Schmidt | Non-parametric Estimation of Tail Dependence[END_REF],

' , = ()* → -. ≤ /#|2 ≤ |# = ()* → #-3 ≤ 4 /# , 5 ≤ 6 /#
(3) The tail dependence coefficients were estimated in a non-parametrical way as in [START_REF] Schmidt | Non-parametric Estimation of Tail Dependence[END_REF] (see Appendix 1). This approach establishes the asymptotic normality and strong consistency, and received estimates provide the general pattern of interdependence. Its application derives a smaller finite sample bias than other approaches such as those of [START_REF] Einmahl | Multivariate extremes estimation[END_REF] and [START_REF] Fermanian | Weak convergence of empirical copula processes[END_REF], to name a few.

' & , = ()* → #-3 > 4 1 -/# , 5 > 6 1 -/# (
Tail dependence was calculated by applying ARMA-GARCH-copula models (i.e., by fitting and testing the Gaussian, t, and Gumbel copulas to residuals from the ARMA-GARCH models in the way as discussed in [START_REF] Cherubini | Copula Methods in Finance[END_REF]). The copula approach has several advantages (see [START_REF] Cherubini | Copula Methods in Finance[END_REF][START_REF] Frees | Understanding Relationships Using Copulae[END_REF]Nelson, 2006). For example, different marginal distributions can be fitted and different types of dependence structures can be analyzed without imposing restrictions that increases a model flexibility.

The Gaussian copula [START_REF] Lee | Generalized econometric models with selectivity[END_REF]) is as follows:

+ ? , ? ; A = Φ C Φ ? , Φ ? ; A = D D E F G H G I J KH L G × N O G FO G F G P QRQ# J KH L H (5)
where S is the CDF of the standard normal distribution, Φ C ? , ? is the standard normal distribution, and A is the correlation parameter in the interval (-1, 1).

Student t copula is defined as [START_REF] Genz | Computation of Multivariate Normal and t Probabilities[END_REF]:

+ ? , ? ; A , A = D D E F G G H G I TG KH L G × N1 + O G F G O G U F G G P F H / QRQ# T H KH L H (6)
where # F H ? is the inverse of the CDF of the standard univariate t-distribution with A degrees of freedom. A controls the heaviness of the tails (i.e., increasing the value of A and decreasing the tendency to exhibit extreme co-movements: if A < 3, then the variance does not exist, and if A < 5, the fourth moment does not exist).

The [START_REF] Gumbel | Bivariate Exponential Distributions[END_REF] copula is also an Archimedean copula and is given by + ? , ? ; A = exp -? [ F + ? [ F /F , (

where ? [ \ = -log ? [ \ .

Maximum likelihood was used to derive the estimates. Finally, I tested whether the correlation is stable over time by means of Spearman's rho.

Data and empirical results

Due to data availability for centralized and decentralized bitcoin markets that use the euro, U.S. dollar, and British pound, I formed the following datasets. For the centralized bitcoin exchanges, I used data provided by GDAX/CoinsBank (for the euro, U.S. dollar, and British pound), Kraken (for the euro and U.S. dollar), Rock Trading (for the euro and U.S. dollar), Bitcoin.de (for the euro), Bitstamp (for the U.S. dollar), and the first exchange-backed peer-to-peer marketplace, Coinfloor (for the British pound). LocalBitcoins was used to form a dataset on a decentralized exchange (a P2P market). To calculate returns, I used the daily closing prices for bitcoin in euros, U.S. dollars, and British pounds from April 27, 2015 to August 9, 2018.4

The estimated kernel densities of the returns for the euro, U.S. dollar, and British pound bitcoin markets are presented in Figure 1 Figures 2 and3 show that volatility in the decentralized exchanges is significantly higher than in the centralized ones. One can observe that in the centralized markets, volatility increases as prices jump (including the exchange-backed peer-to-peer marketplace Coinfloor).

It is contrary to simple leverage reasons that market drops cause volatility. Compared to "traditional" institutional money, the significant increases in a price are not typically accompanied by jumps in volatility.

The results also show, that in the centralized bitcoin markets volatility increased after the announcement of trading in Bitcoin futures. This is in line with [START_REF] Corbet | Bitcoin Futures-What use are they?[END_REF] findings.

Regarding trade volume, it increased as prices decreased (see Figure A1 in the appendix), demonstrating a lack of confidence and consensus of a large number of participants in the price-jump period. Thus, volumes did not go with the trend.

Regarding correlation in tails, the tail dependence coefficients showed the following pattern (see Tables 1 and2): Trading also demonstrated the analogous pattern. Similar results were derived from the ARMA-GARCH copula models (see Figure A2 in the Appendix). Also, tail correlation within the same exchanges was higher than between currency pairs across different exchanges.

Several tail copulas studies, such as [START_REF] Ang | Asymmetric correlations of equity portfolios[END_REF], [START_REF] Giacomini | Inhomogeneous dependence modeling with time-varying copulae[END_REF][START_REF] Hu | Dependence patterns across financial markets: a mixed copula approach[END_REF], and [START_REF] Hong | Asymmetries in stock returns: statistical tests and economic evaluation[END_REF], demonstrated that "traditional" financial markets are more extreme dependent in downturns; thus, a lower tail coefficient is superior to an upper tail one.

Therefore, the results support this hypothesis also for the centralized cryptocurrency markets.

The tail dependence coefficients for the P2P bitcoin market were different from those obtained for the centralized markets; there was no significant tail dependence in bitcoin returns (see Table 2 here and Figure A3 in the Appendix). Because correlation is not linear and variance is not constant, Spearman's rho was applied to test stability of correlation over time between the pairs of the selected markets (see Figure 4). Figures 4 and5 show, in general, the correlation was not stable over time and was higher for the centralized exchanges. The currency pairs of the analyzed decentralized bitcoin market were less correlated, and this correlation was more volatile. In general, the estimates showed the correlation between the euro and U.S. dollar centralized markets was slightly higher, and the correlation between the currency pairs was higher within the same exchange.

The only exception was the U.S. dollar-British pound, where the correlation between Kraken

Conclusion

The performed analysis depicts the correlation between the currency pairs across the decentralized (LocalBitcoins) and six centralized bitcoin markets (GDAX, Kraken, Bitcoin.de, Bitstamp, Rock Trading, and the exchange-backed peer-to-peer marketplace, Coinfloor). The estimates derived reasonable results and were consistent with the data provided.

The results showed that the decentralized bitcoin exchange has higher volatility than the centralized markets. In the centralized markets, volatility increased as prices jumped, but that was not the case in the decentralized exchanges. The volatility analysis results are contrary to the conventional leverage reason that market drops cause volatility. Thus, bitcoin price increases cause an increase in risk because a greater dispersion of returns around the mean causes a larger drop in the compound return.

In addition, trade volume increased as price decreased. This discovered pattern can also mean that investors, based on their sentiments, may overreact because they do not know how to precisely estimate the bitcoin price.

The centralized bitcoin markets demonstrated a higher left tail dependence that is in line with the general pattern in the "traditional" financial markets that are more extreme dependent in downturns. Having long tails, the decentralized bitcoin market does not show a strong pattern of extreme (tail) dependency. Nevertheless, the centralized euro bitcoin markets have a lower average value of tail dependency, and the U.S. dollar exchanges have a superior dependency in terms of extreme low and high returns. The decentralized bitcoin market is more interdependent. Thus, a price setting in the P2P bitcoin markets is not correlated over time or across markets. Correlation between markets is not stable over time and is higher for the centralized exchanges.

Thus, the results support the existing findings regarding bitcoin volatility [START_REF] Blau | Price dynamics and speculative trading in Bitcoin[END_REF][START_REF] Katsiampa | Volatility estimation for Bitcoin: A comparison of GARCH models[END_REF][START_REF] Klein | Bitcoin is not the New Gold -A comparison of volatility, correlation, and portfolio performance[END_REF] and extend them with the comparative analysis of the centralized and decentralized bitcoin markets that use the euro, U.S. dollar, and British pound.

They also support [START_REF] Pieters | Financial Regulations and Price Inconsistencies across Bitcoin Markets[END_REF] findings that the bitcoin markets fail to achieve the law of one price (LOOP). Due to the identical nature of all bitcoin, this means that the analyzed bitcoin markets differ in their characteristics. In general, the results are in line with the hypothesis that bitcoin performs poorly as a unit of account and as a store of value and can be considered rather as a speculative asset [START_REF] Corbet | Bitcoin Futures-What use are they?[END_REF][START_REF] Yermack | Is Bitcoin a real currency? An economic appraisal Handbook of Digital Currency[END_REF].

Overall, these findings are significant and potentially useful to researchers, practitioners, and participants in the highly volatile bitcoin market for making better investment and risk management decisions. Future research can extend the datasets to cover other currencies; focus on higher-frequency bitcoin data, such as hourly (because the results can vary from the current ones); and use other dynamic models based on the concept of copulae, such as the generalized autoregressive score (GAS) framework. Note: the ARMA-GARCH models used for residuals calculation are the same as in the EMWA calculations (see Figure 3).

  proved the AR-CGARCH model can better explain bitcoin price volatility, and Chan et al. (2018) used pairwise GARCH models and constant conditional correlation models to test the dynamic hedging abilities of bitcoin. Klein et al. (2018) applied a BEKK-GARCH model to estimate time-varying conditional correlations between bitcoin and gold. Symitsi and Chalvatzis (2018) used an asymmetric multivariate VAR-GARCH model to study spillover effects between bitcoin and energy and technology companies. Moreover, Baur et al. (2018) used GARCH models to demonstrate that bitcoin shows distinctively different return volatility and correlation characteristics compared to other financial assets.

  4)where Λ : ℝ → ℝ is a lower tail copula, the function Λ & : ℝ → ℝ is an upper tail function, and C is a copula associated with a distribution function; ℝ ∶= 90, ∞< { ∞, ∞ }. It means that a lower (left) tail copula and an upper (right) tail copula describe the dependence structure of the extreme-value distribution, C.

Fig. 1 .Fig. 2 .Fig. 3 .

 123 Fig. 1. Kernel density of bitcoin returns in the euro (black), U.S. dollar, and British pound

Fig. 5 .

 5 Fig. 5. Correlation over time between the currencies, decentralized bitcoin exchange,

Fig. A2 .

 A2 Fig. A2. Tail dependence between the currencies, the centralized bitcoin exchanges, ARMA-

Table 1 .

 1 Lower (left) and upper (right) tails of the Bitcoin return dependency, centralized exchanges.Based on Table1, one can assume there is a significant tail correlation in bitcoin returns across the centralized markets, especially in the lower (left) tail. The results obtained for Rock

					Kraken (euro, U.S.
			GDAX		dollar)/GDAX (British	Bitcoin.de/ BitStamp/Coinfloor
	Bitcoin
						pound)
	currency
									U.S.	British
	market		U.S.	British		U.S.	British	Euro,
		Euro			Euro				dollar,	pound,
			dollar	pound		dollar	pound	Bitcoin.de
									BitStamp	Coinfloor
				Lower tail of the Bitcoin returns dependency
	Euro	1	0.882	0.88	1	0.828	0.759	1	0.62	0.59
	U.S.
			1	0.91		1	0.862		1	0.82
	dollar
	British
				1			1			1
	pound
				Upper tail of the Bitcoins returns dependency
	Euro	1	0.85	0.824	1	0.862	0.896	1	0.62	0.56
	U.S.
			1	0.85		1	0.897		1	0.82
	dollar
	British
				1			1			1
	pound

Table 2 .

 2 Lower (left) and upper (right) tails of the bitcoin return dependency, the

	decentralized exchange
	Bitcoin currency	Decentralized exchange, LocalBitcoins
	market	Euro	U.S. dollar	British pound
	Lower tail of the Bitcoin returns dependency
	Euro	1	0.03	0
	U.S. dollar		1	0.029
	British pound			1
	Upper tail of the Bitcoins returns dependency
	Euro	1	0.0	0.03
	U.S. dollar		1	0.029
	British pound			1

As of December

2017, the bitcoin peaked with an estimated transaction volume of USD $4,652,736,998 and market capitalization of USD $325,000,000,000(Blockchain.info, 2018).

China is excluded from the analysis because Beijing shut down all virtual exchanges.

GDAX is CoinsBank now.

Due to data availability, for the dataset Rock Trading-Coinfloor, the dataset covers the period of April 27, 2015 to January 3, 2018.

Appendix 1. Non-parametric estimators of tail dependence (for the detail discussion see [START_REF] Schmidt | Non-parametric Estimation of Tail Dependence[END_REF].

Empirical copula is defined by

where + _ is the empirical copula, 4 _ , and 6 _ are the empirical distribution functions corresponding to marginal distribution functions G, H.

Following the definition of tail copulae, (see Equation 1and Equation 2):

the first step of estimators (Genest et al., 1995;[START_REF] Schmidt | Non-parametric Estimation of Tail Dependence[END_REF] known as empirical tail copulae, are

∑ 1 {o pH q r_ hO stu o pG q r_ hk} _ \v (A.5) where w _ \ and w _ \ are the rank of independent and identically distributed (i.i.d.) random vectors X (j) and Y (j) , j=1,…, m; x ∈ {1, … , *}, x = x * → ∞ and x/* → ∞ as * → ∞.


Keywords: eq Equivalent component, ref Reference component, n Nominal value, Superscripts h h th harmonic component, fund Fundamental component

h-harmonic positive gain. Selective h-harmonic voltage distortion.

!

Proportional h-harmonic distortion gain. "

Proportional gain of the harmonic droop.

#,

Damping factor and cutoff frequency of lowpass filter.

Introduction

In an islanded microgrid, Distributed Generation (DG) units use the droop method as a decentralized control to share the active and reactive power in proportion to their rated power [START_REF] Kim | Mode adaptive droop control with virtual output impedances for an inverter-based flexible AC microgrid[END_REF]- [START_REF] Xin | Control of island AC microgrids using a fully distributed approach[END_REF]. The droop method has been proved to be an interesting strategy that enables the DG units to share fundamental real and reactive power without inter-unit communication, i.e. in autonomous way. However, system stability and proper load sharing are not always maintained when using the conventional droop. This is due to the mismatch in line impedances and the slow dynamic response of the droop. Therefore, several improvements were considered in literature to enhance the performance of the droop control [START_REF] Han | An improved droop control strategy for reactive power sharing in islanded microgrid[END_REF]- [START_REF] Bojoi | Enhanced power quality control strategy for single-phase inverters in distributed generation systems[END_REF].

Another issue in the islanded microgrid deals with harmonic power sharing when nonlinear loads are connected. In fact, the conventional droop control is designed to share the fundamental positive-sequence power components. The traditional way to mitigate harmonics is by using passive or active filters which increase system complexity and cost. As a consequence, new research studies apply different control proposals to reduce the harmonic distortion and manage the harmonic currents between DG units. For example, a virtual output impedance loop is introduced to autonomously share the harmonic currents [START_REF] Han | An improved droop control strategy for reactive power sharing in islanded microgrid[END_REF], [START_REF] He | An islanding Microgrid power sharing approach using enhanced virtual impedance control scheme[END_REF]- [START_REF] Han | An Enhanced Power Sharing Scheme for Voltage Unbalance and Harmonics Compensation in an Islanded AC Microgrid[END_REF]. It was seen that, in order to attenuate the effect of different feeder line impedances, a large virtual impedance value is designed. Therefore, the more voltage distortions exist. Another approach, uses a virtual capacitive loop to improve the voltage distortion and to share accurately the harmonic current [START_REF] Micallef | Reactive power sharing and voltage harmonic distortion compensation of droop controlled single phase islanded microgrids[END_REF], [START_REF] Micallef | Mitigation of Harmonics in Grid-Connected and Islanded Microgrids via Virtual Admittances and Impedances[END_REF]. The idea behind using a capacitive impedance value is to compensate the inductive line impedance by introducing a virtual capacitive element which is equal in magnitude but with opposite phase shift. However, these approaches are based on knowing the line impedance to tune the appropriate virtual impedance value. In [START_REF] Lee | Design of a new cooperative harmonic filtering strategy for distributed generation interface converters in an islanding network[END_REF], [START_REF] Zhong | Harmonic droop controller to reduce the voltage harmonics of inverters[END_REF], the authors propose a harmonic droop controller to damp the harmonic distortion at PCC and to share the harmonic current between converters. However, additional computation and system complexity restrict the use of such controllers. It is to be noted that, in a weak low-voltage (LV) microgrid, the damping of total harmonic distortion (THD) is highly affected by the microgrid line impedances. In [START_REF] He | A flexible harmonic control approach through voltage-controlled DG-grid interfacing converters[END_REF], a voltage controller method for each DG in a microgrid is proposed to compensate the voltage distortion by adjusting the DG equivalent harmonic impedance. It was seen that the harmonic sharing is sensitive to the grid impedance variation or disturbance in the main grid. To increase the reliability and efficiency of the system, new proposals based on consensus algorithm are presented recently in [START_REF] Zhang | Distributed Adaptive Virtual Impedance Control for Accurate Reactive Power Sharing Based on Consensus Control in Microgrids[END_REF], [START_REF] Kim | Consensus-based improved droop control for suppressing circulating current using adaptive virtual impedance in microgrids[END_REF]. The consensus method is a control technique that requires a minimal communication infrastructure to significantly enhance the performance of conventional droop and harmonic power sharing as well.

In this paper, the main objective is to damp the voltage distortion at the PCC when non-linear loads are connected. Furthermore, the distribution of harmonic currents between different DGs is managed according to the harmonic power rating settings. The overall controller is applied in one loop design. The operation principles of the proposed controller are listed as:

• Design a droop function for the harmonic power. This enables the DG units to proportionally share the harmonic power between each other according to their rating powers.

• Mitigate the voltage harmonic distortion generated by non-linear loads.

• Apply a control loop that is immune from grid disturbance and line impedance mismatch.

• Simplify the control design as possible. There is no need to know the feeder impedances in advance.

Power Sharing in Islanded Microgrid

As shown in Fig. 1, a microgrid formed by n-Distributed Generators (DGs) with linear and non-linear loads linked together by means of feeder line impedances. The harmonic (dq) extraction block is there inorder to measure the voltage distortion at the PCC and share the data to all DGs by using a communication link infrastructure. It is to be noted that, the communication linked is used only to improve the performance of the harmonic power sharing and to compensate the PCC voltage distortion. The only data needed by the proposed harmonic droop is the voltage distortion at the PCC. These data are dc-components (in synchronous dq-frame). This will simplify the type of communication by using low-bandwidth communication.

Each DG is represented by a Distributed Interfaced Converter (DIC) with an output filter and a local control block. The local control can be decomposed of a primary control loop which is defined here by the droop controller and an inner control loop inorder to regulate the output voltage waveform as regard to the reference voltage.

A) Droop Control for fundamental power sharing

The droop control for voltage source converters (VSC) is commonly used to achieve the power balance in islanded microgrid. It consists of a simple proportional control and the selection of its gains are commonly based on power sharing criteria. To demonstrate the droop control, two DG-sources with general output impedances are connected to common linear and nonlinear loads as in Fig. 2. The real and reactive power delivered through the output impedance can be written as:

i PCC PCC i PCC i i i i i i i i i PCC PCC i PCC i i i i i i i i VV V VV P Z Z Z VV V VV Q Z Z Z δ θ δ θ δ θ δ θ   = - +       = - -     (1.1)
where i V is the output voltage amplitude of DG i ( {1, 2}) i ∈ , PCC V
is the common bus voltage, i δ is the power angle between i V and PCC V . i Z and i θ are, respectively, the corresponding magnitude and phase of the equivalent impedance between DG i and the common AC bus.

Non-Linear Load

1 L 1 DG f i L i ch R PCC Linear Load ch L 2 L 2 R 1 R 1 o i 2 o i 1 1 Z θ ∠ 2 2 Z θ ∠ 2 2 V δ ∠ 1 1 V δ ∠ 2 DG 0 PCC V °∠ Fig. 2.
3 2 3 ( ) 2 i PCC i i i PCC i i PCC i VV P X V Q V V X δ = = - (1.2)
Thus, the real power is strongly dependent on the power angle i δ , while the reactive power is mainly subjected to the voltage amplitude i V. Then, the conventional droop controller can be expressed as:

0 0 i i i i i i m P V V n Q ω ω = - = - (1.3)
where 0 ω and 0

V are the output frequency and voltage at no load, i m and i n are the frequency and amplitude droop coefficients of DG i , respectively. However, the conventional droop controller is not always true since in Low-voltage (LV) microgrid, the feeder lines are mainly resistive. In this case, the droop control in its previous configuration may cause stability problem due to the coupling between the real and reactive power.

ii. Case 2: Resistive Impedance ( 0

) i θ = o
In the case of a LV microgrid, i R is much larger than i X ( )

i i R X 
. Thus, after performing some assumptions in equation (1.1), it will be rewritten as follow

3 ( ) 2 3 2 PCC i i PCC i i PCC i i i V P V V R VV Q R δ = - = - (1.4)
So here, the real power is controlled by controlling the voltage amplitude i V and the reactive power is controlled by the power angle i δ . Hence the droop control strategy takes the form

0 0 i i i i i i m Q V V n P ω ω = + = - (1.5)

B) Harmonic Power Sharing

When significant non-linear loads are connected to a microgrid, harmonic power sharing errors present due to different DG ratings and grid parameters. Inaccurate sharing of harmonic power can lead to a DG overcapacity, overheating, losses and distortion in current and voltage waveforms. Therefore, a control strategy is investigated to accurately share harmonic power and compensate the THD of the PCC voltage. It was shown that, DG units location can have significant impact on power quality. In fact, the location of DG plays a vital role in improving the voltage profile and reducing losses. The way to share harmonic power between different DG units is achieved by changing the effective line impedances. If the equivalent harmonic impedance is reduced properly, the PCC voltage quality can be improved [START_REF] He | A flexible harmonic control approach through voltage-controlled DG-grid interfacing converters[END_REF]. Fig. 3 shows a simplified circuit of two DG units connected to a non-linear load represented by a passive impedance h Load Z and a harmonic current source h I . In ideal case, when harmonic current passed, the DG unit is equivalent to a zero impedance. In Fig. 3 the DG unit harmonic impedance is considered based on the existence of internal output impedance.

To effectively change the physical line impedance, two ways were reported in literature. The first method is done by adding a virtual conductance [29]. In this way, harmonic sharing error may exist if the difference between the feeder line impedances, Line1 Line2 ,

h h Z Z is important.
The second method was proposed in [START_REF] Han | An improved droop control strategy for reactive power sharing in islanded microgrid[END_REF], [START_REF] He | An islanding Microgrid power sharing approach using enhanced virtual impedance control scheme[END_REF]- [START_REF] Han | An Enhanced Power Sharing Scheme for Voltage Unbalance and Harmonics Compensation in an Islanded AC Microgrid[END_REF], which is based on adding a virtual impedance in series to the harmonic line impedance. This virtual part is designed to be larger than the line impedance. This will insure equal impedances for the DG units to share equally the harmonic power. However, this method will increase the voltage distortion at the PCC.

Proposed Harmonic Droop Controller

It is clear that sharing the harmonic power and compensating the PCC voltage distortion cannot be done unless using several virtual impedance loops. The proposed harmonic droop is able to make the two tasks in one control loop. Furthermore, the proposed method is less sensitive to physical line impedance mismatch. Therefore, accurate harmonic sharing can be carried out with a simpler controller. A) The essence of the proposed harmonic droop

DG PCC i h h h i V G V = -⋅ (1.6)
Therefore, the equivalent harmonic impedance of DG i seen from the PCC point can be obtained as follow:

i i i i i i i h h h h h h h i h h h h h i I V V Z G V Z Z V Z I G = = + = = + (1.7)
where, DG i h I and i h DG Z are the output harmonic current and the harmonic impedance of DG i respectively.

Based on equation (1.7), it is obvious to interpret that a proper control of the positive gain h i G leads to a reduction in DG i line impedance by a factor of ( 1)

h i G +
. In this way, it is clear that the PCC harmonic voltage is minimized. To accurately share the harmonic current between different DGs, the proposed harmonic droop is designed to control the positive gain h i G as regard to the harmonic sharing criteria given for the different DG units in a microgrid.

B) Proposed Harmonic Droop characteristics

The overall block diagram representing the DG i connected to a microgrid across a feeder line impedance with the proposed harmonic droop controller is shown in Fig. 5. It is seen that a minimal communication link is used to send the harmonic voltages at the PCC to the harmonic droop controller. The harmonic droop uses that information to generate the appropriate reference voltage DG _ref ( )

i h V . The
proposed harmonic droop works in collaboration with the fundamental droop control to generate the total reference voltage of DG i unit. The inner controller is there to insure the measure voltage follows perfectly the reference. 5. The overall proposed control strategy with the harmonic droop controller for DG interfaced converter.

V 7 PCC V - + - + - - 1 s Common Control Part max HD h i Q _rated h i Q h K i b h i G PCC . -h h i G V DG _ref i h V PCC h V h HD Selective h HD Fig.

Based on the proposed harmonic technique (

) 5, the positive gain ( )

h h i i Q G - presented in Fig.
h i
G is defined as function of the harmonic power rating _rated ( )

h i
Q and the measured harmonic power ( )

h i
Q . The characteristics of the harmonic droop controller is identical to the conventional droop where the respective parameters are designed in similar way. The harmonic droop coefficient ( ) i b for all DG units is defined as:

1 1_rated 2 2_rated _rated h h h n n b Q b Q b Q = = = L (1.8)
where _rated [START_REF] Kim | Mode adaptive droop control with virtual output impedances for an inverter-based flexible AC microgrid[END_REF] h i Q i n = K , is the rated harmonic capacity of each DG unit.

To maintain a proportional share of harmonic power between DG units without being affected by the mismatch in line impedances, a common control part is added in the proposed harmonic droop ( )

h h i i Q G - .
As seen in Fig. 5, a selective harmonic distortion (HD h ) is compared with an input setting value, HD max , the maximum allowed harmonic distortion. At steady state condition, the integrator input value will reach zero. Since the common control part is designed to be the same for all DG units, then we obtain the following:

1 1 1_rated 2 2 2 _rated _rated ( ) ( ) ( ) h h h h h h n n n b Q Q b Q Q b Q Q - = - = = - L (1.9)
And according to equation (1.8), a proportional share of harmonic power is satisfied. The following equation will be considered as:

1 1 2 2 h h h n n b Q b Q b Q = = = L (1.
10) The harmonic power drawn by each DG unit to the non-linear load is defined as:

DG DG 3 i i h fund h i Q V I = ⋅ ⋅ (1.11)
where, DG i fund V and DG i h I represent the fundamental rms-voltage at the output of DG i and the h-harmonic current traversing the DG i line impedance.

To set the harmonic droop coefficients, the following equation is used:

0 _rated / h i i b b Q = (1.12)
The gain b 0 is tuned to be high enough with a compromised strategy in order to adopt a proper harmonic sharing and ensure a stable system at the same time. The gain K h is selected to be the same for all DG units and defined as

0 max / h K b HD = (1.13)
Consequently, by adjusting the gain b 0 and the maximum harmonic distortion HD max , the proposed harmonic droop controller become easy to tune and flexible to work with. The proposed controller will dynamically change the gain h i G for each specific harmonic frequency in order to set the reference harmonic voltage for the DG unit to generate at the output.

The h-harmonic distortion of the PCC voltage can be evaluated using the following equation:

/ h fund h PCC PCC HD V V = (1.14)
where, fund PCC V and h PCC V represent the fundamental and the h-harmonic component of the voltage at PCC.

The harmonic component is extracted using real time dq-synchronous individual harmonic frame method as shown in Fig. 5.

After adding a harmonic reference voltage to the input of the DG inner controller, a saturation process should be avoided. For this reason, the integrator of the proposed harmonic droop should be limited to the inverter saturation boundaries when using a specific Modulation. In the case of sinusoidal PWM, the saturation limit is given as:

_ref DG 3 / 2( / 2) i dc V V ≤ (1.15)
After adding the harmonic reference voltage to the DG i inner controller, the reference rms-voltage, _ref DG i

V

, can be expressed as:

_ _ref _ref 2 2 DG DG ( ) ( ) i ref i i fund h DG V V V = + (1.16)
The h-harmonic reference voltage is written as:

DG _ref PCC DG _ref (max) max PCC . . . i i h h h i h h fund i V G V V G HD V = - = - (1.17)
Therefore, substituting equations (1.16) and (1.17) in the inequality (1.15), the following limiting values of the gain h i G is obtained:

2 2 3 dc DG 2 max DG ( 2) ( ) 0 i i fund h i fund V V G HD V - ≤ ≤ (1.18)

Simulation Results

The proposed harmonic droop controller is tested in MATLAB simulator for a microgrid with three DGs connected to a linear and non-linear load. The linear load is represented by a classical R-L impedance and the nonlinear is modeled by current sources that inject harmonics to the grid. The performance of the proposed harmonic droop is checked for 5 th and 7 th harmonic orders. The dynamic behavior of the harmonic current is limited to the bandwidth of the voltage controller of the DG unit.

The system parameters of the LV-microgrid considered in this test are given in Table I. The three DG units are provided with the proposed harmonic droop in addition to the fundamental droop controller. In the following, the proposed harmonic droop is tested in two different cases.

A) Equal harmonic power ratings of DG units

The three sources are designed to have equal power ratings for both fundamental and harmonics with different line impedances. In this case, we will analyze the difference before and after activating the proposed harmonic droop to illustrate its performance in balancing the draw of current between the three DG units. In Fig. 6, the fundamental power sharing is presented when P-V and Q-ω droop is adopted for LV microgrid. The results show an equal fundamental power sharing after using an improved droop controller as proposed in [START_REF] Zhong | Robust droop controller for accurate proportional load sharing among inverters operated in parallel[END_REF].

The 5 th and 7 th harmonic power sharing are shown in Fig. 7. In the first stage, when the proposed harmonic droop is turned OFF, the harmonic power sharing between DG units is not satisfied. After turning ON the harmonic droop controller in three DG units, the harmonic power is balanced between the units regardless of the mismatched in line impedances. It is to be noted that in [START_REF] Moussa | Voltage Harmonic Distortion Compensation with Non-Linear Load sharing of Droop Controlled Islanded Microgrid[END_REF], the harmonic droop control is tested on a microgrid with inductive line impedances. In this case P-ω and Q-V are adopted to share the fundamental power and according to Fig. 8, the harmonic power is shared equally between two DG units by applying the proposed harmonic droop.

Q 7 3 Q 7 2 Q 5 3 Q 5 1 Q 7 1 Q 5 2
Harmonic Droop stopped in DG3 Harmonic Droop stopped in DG3 Fig. 9. Simulation Results of Case A, after stopping Harmonic Droop in DG 3 : (a) 5 th Harmonic power for three DG units: , and

Q Q Q .
In Fig. 9 the harmonic droop is stopped in DG 3 , the unit that is near the load, where the other DG units still participate in sharing the harmonic droop. The idea behind this test is to examine the performance of the proposed droop when communication is lost in any DG unit. The DG unit that lose communication with the PCC is allowed to share only the fundamental power whereas the harmonic droop is turned OFF. It is shown in Fig. 9 that the harmonic power is perfectly distributed between DG 1 and DG 2 with no harmonic power drawn from DG 3 . Fig. 10 shows the output current waveforms generated by each DG in two stages. Fig. 10(a) represents the current waveforms when the proposed harmonic droop is running in the three DGs. While Fig. 10(b) represents the current waveforms when DG 3 stopped in sharing the harmonic power with other DG units. This part shows a sinusoidal waveform of current

Q 5 1 Q 5 2 Q 5 3 Q 7 1 Q 7 2 Q 7 3
Harmonic Droop start Proportional 7 th Harmonic power sharing

B) Different harmonic power ratings of DGs

The same microgrid configuration is adopted, but considering different harmonic power ratings for the DG units as shown in Table I. In this case, the efficiency of the proposed harmonic droop will be illustrated. All the DG units are supported by the proposed harmonic droop that will guarantee the proportional harmonic sharing between the units. The same simulation procedure is tested in this case. Fig. 11 shows the harmonic power sharing between the three sources. As represented in the figure, the two sources DG 1 and DG 2 are able to share equally the harmonic power which is half the quantity shared by DG 3 . This interpret the accuracy of the proposed harmonic droop in sharing the harmonic power. It was suggested to stop the harmonic droop in DG 3 , Fig. 12 shows the repartition of harmonic power between the DG units after stopping the harmonic droop in DG 3 . The 5 th and 7 th harmonic powers drawn by DG 3 are minimized as possible. Fig. 13 shows the output phase currents for each DG when the harmonic droop controller is applied to all DGs and after terminating the harmonic droop in DG 3 . Fig. 13( Phase currents (A) The FFT analysis is done to interpret the efficiency of the proposed harmonic droop controller in compensating the voltage distortion at PCC. Fig. 14 shows the study and the distribution of harmonic contents of the PCC output voltage before and after starting the harmonic droop. The total harmonic distortion (THD) drops from 1.22% to 0.42% after compensation.

Q 7 3 Q 7 2 Q 5 2 Q 7 1 Q 5 1 Q 5 3

Experimental Results

A microgrid with two DGs connected to a common R-L load and a non-linear full bridge rectifier load is constructed in the laboratory as shown in Fig. 15. The system parameters used in this test are given in Table III. The experiment test bench allows to examine the effectiveness and performance of the proposed harmonic droop controller in sharing the harmonic power and compensating the voltage distortions even in microgrid with more inductive feeder lines. For case 1, Fig. 16 shows the 5 th and 7 th harmonic Var for the two DG units before and after starting the harmonic droop. The same figure shows the output gains time response of the harmonic droop that is responsible on adding the reference harmonic voltage of each DG to regulate the harmonic output currents.

Q 5 2 Q 5 2 G 5 1 G 7 1 Q 7 2 Q 7 2 G 7 1
G Fig. 16. Experimental results for case 1: 5 th and 7 th Harmonic power repartition with output gain time response for both DGs. Fig. 17 verifies the operation of the harmonic droop controller under different harmonic ratings which is related to case 2. The distribution of the average harmonic power between the two DG units shows that a proportionally share of power occur according to the ratings given in table III. In order to demonstrate more the effect of harmonic contents on the voltage and currents signals, Fig. 18 shows the voltage and current waveforms for the two DGs before and after starting the harmonic droop. The figure indicates that more harmonic contents are added to the voltage of DG 2 in order to be able to absorb more harmonic current.

The compensation of the PCC voltage distortions is demonstrated by Fig. 19. The harmonic droop allows reducing the THD from 2.88% to 1.88%. The experimental test concerns only the reduction of 5 th and 7 th harmonic order contents. The 5 th and 7 th harmonic voltage components are widely reduced since, HD5 drops from 1.85% to 0.6% and HD7 goes down from 1.37% to 0.34%.

Conclusion

In this paper, a new harmonic droop controller is proposed to share the harmonic power between different DG units and to compensate the voltage distortions at PCC. The main contribution of this paper is that the proposed harmonic droop reduces the voltage distortion at the end users autonomously between different DGs in a one loop controller without the knowledge of line characteristics. The harmonic droop current sharing is unaffected by the mismatch in the line impedances. The proposed controller is efficient, as well, when mismatch in line impedances present. The overall control strategy is discussed deeply in section II and simulation results verify the effectiveness of the proposed controller for different case studies. The proposed method works in collaboration with the fundamental droop controller without any interaction between them. Experimental tests are applied to demonstrate the performance of the harmonic droop in real application. The proposed harmonic droop works in case of resistive or inductive microgrid as realized in the simulation and experimental tests. The results show the competence of the proposed harmonic droop in compensation as compared with other controllers.

7 Fig. 1 .

 71 Fig.1. Three-phase microgrid with parallel DGs, Linear and nonlinear loads.

Fig. 3 .

 3 Fig. 3. Equivalent circuit of two DGs at harmonic frequency range.

Figure 4 Fig. 4 .

 44 Figure4shows a simplified microgrid with two DG units supplying a common non-linear load. The load is modeled as a harmonic source current h I and a passive impedance Load h Z . Each DG unit is

Fig. 6 .

 6 Fig. 6. Simulation Results of Case A: (a) Fundamental Real Power for three DG units: P 1 , P 2 and P 3 ; (b) Fundamental Reactive power for three DG units: Q 1 , Q 2 and Q 3 .

Fig. 7 . 1 QFig. 8 .

 718 Fig. 7. Simulation Results of Case A: (a) 5 th Harmonic power for three DG units: 5 5 5 1 2 3 , and Q Q Q ; (b) 7 th Harmonic power for three DG units: 7 7 7 1 2 3 , and Q Q Q .

  b) 7 th Harmonic power for three DG units:

Fig. 11 .

 11 Fig. 11. Simulation Results of Case B: (a) 5 th Harmonic power for three DG units: 5 5 5 1 2 3 , and Q Q Q ; (b) 7 th Harmonic power for three DG units: 7 7 7 1 2 3 , and Q Q Q .

Fig. 12 .

 12 Fig. 12. Simulation Results of Case B, after stopping Harmonic Droop in DG 3 : (a) 5 th Harmonic power for three DG units: 5 5 5 1 2 3 , and Q Q Q ; (b) 7 th Harmonic power for three DG units: 7 7 7 1 2 3 , and Q Q Q .

Fig. 13 .Fig. 14 .

 1314 Fig. 13. Simulation Results of Case B: (a) Phase currents after harmonic compensation in all DG units; (b) Phase currents after harmonic compensation in DG 1 and DG 2 .

Fig. 15 .

 15 Fig. 15. Laboratory test bench with all equipment.

QFig. 17 .Fig. 18 .

 1718 Fig.[START_REF] Lee | Design of a new cooperative harmonic filtering strategy for distributed generation interface converters in an islanding network[END_REF]. Experimental results for case 2: 5 th and 7 th Harmonic power repartition with output gain time response for both DGs.

Fig. 19 .

 19 Fig.[START_REF] He | A flexible harmonic control approach through voltage-controlled DG-grid interfacing converters[END_REF]. Experimental results: PCC voltages with harmonic analysis before and after harmonic droop compensation.

TABLE I .

 I SIMULATION PARAMETERS FOR THREE DGS CONFIGURATION

					System Parameters
	Fundamental Line voltage (V)/frequency (Hz)	110V/60Hz
		Feeder impedance	1 DG Z (R 1 , L 1 )			(0.58 , 0.05 mH) Ω
		Feeder impedance	2 Z (R 2 , L 2 ) DG			(0.2 , 0.02 mH) Ω
		Feeder impedance	3 Z (R 3 , L 3 ) DG			(0.1 , 0.01 mH) Ω
	5 th harmonic current source 5 dq I	5 d = + I	j	I	5 q	6 -j2 (A)
	7 th harmonic current source 7 dq I	7 d = + I	j	I	7 q	-5 (A)
	Max. Harmonic distortion	max HD (%)	0.1%
				DG Droop Characteristics
	Control Parameters			DG 1			DG 2	DG 3
	i m	(rad/Var)	1.33 10 -×	4		1.33 10 -×	4	1.33 10 -×	4
	i n	(V/W)		2.5 10 -×	3		2.5 10 -×	3	2.5 10 -×	3
				Case A/Case B	Case A/Case B	Case A/Case B
	_rated (Var) h i Q		500/500			500/500	500/1000
	b i (1/Var)		0.06/0.06		0.06/0.06	0.06/0.03

  Fig. 10. Simulation Results of Case A: (a) Phase currents after harmonic compensation in all DG units; (b) Phase currents after harmonic compensation in DG 1 and DG 2 .

					6				Phase currents after starting the Harmonic droop in all DGs
			Phase currents (A)	-4 -2 0 2 4
					-6		17.01	17.02	17.03	17.04	17.05	17.06	17.07	17.08	17.09	17.1
												Time (seconds)
														(a)
					10				Phase currents after stopping the Harmonic droop in DG3
			Phase currents (A)	-8 -6 -4 -2 0 2 4 6 8			i DG 1 i DG 2	i DG 3
				-10	22	22.01	22.02	22.03	22.04	22.05	22.06	22.07	22.08	22.09	22.1
												Time (seconds)
														(b)
	5 th Harmonic Power	(Var)	200 300 400 500 600 700	6		7		8	9		10	Time (seconds) 11	12	13	DG i Proportional 5 th Harmonic 14 15 power sharing	3 16	and the
	harmonic contents are distributed equally between	DG i	1	,	2 i . DG (a)
	Harmonic Power	(Var)	300 400 500				Harmonic Droop start
	7 th		200
				6		7		8	9		10				11	12	13	14	15	16
												Time (seconds)
															(b)

TABLE II .

 II Repartition of Harmonic power between three DG units for the cases (A) and (B)

		Equal harmonic power ratings of DG units (Case A)	Different harmonic power ratings of DG units (Case B)
		Before Harmonic droop	After Harmonic droop	Before Harmonic	After Harmonic droop
		compensation		compensation		droop compensation	compensation
		5 th Har.	7 th Har.	5 th Har.	7 th Har.	5 th Har.		7 th Har.	5 th Har.	7 th Har.
		( % )	( ' )	( % )	( ' )	( % )		( ' )	( % )	( ' )
	DG1	356 Var	297 Var				356 Var		297 Var	314.25 Var	244.12 Var
	DG2	443 Var	337.5 Var	419 Var	325.5 Var	443 Var		337.5 Var	314.25 Var	244.12 Var
	DG3	458 Var	342 Var				458 Var		342 Var	628.5 Var	488.25 Var
		800
		600
		400
		200
		0
		14	15	16	17	18	19	20	21	22	23	24	25
		600
		400
		200
		0
		14	15	16	17	18	19	20	21	22	23	24	25

TABLE III .

 III SIMULATION PARAMETERS FOR TWO DGS CONFIGURATION

	1
				System Parameters
	Fundamental Line voltage (V)/frequency (Hz)		110V/60Hz
		Feeder impedance	1 DG Z (R 1 , L 1 )			(0.1 , 1 mH) Ω
		Feeder impedance	2 Z (R 2 , L 2 ) DG			(0.2 , 5 m H ) Ω
	Max. Harmonic distortion	max HD (%)				1%
			DG Droop Characteristics
					Case 1:				Case 2:
	Control	Equal harmonic power	Different harmonic
	parameter			ratings			power ratings
			DG 1			DG 2			DG 1	DG 2
	i m	(rad/W)	5 10 -×	5		5 10 -×	5		5 10 -×	5	5 10 -×	5
	i n	(V/Var)	6.7 10 -×	3	6.7 10 -×	3	6.7 10 -×	3	6.7 10 -×	3
	_rated (Var) h i Q	500			500			500	250
	b i (1/Var)	0.1			0.1			0.1	0.2


Romain Vandeputte

Dorra Khiri

Cécile Lafont

Laurent Cantrel

Florent Louis

Theoretical investigation of thermochemical properties of cesium borates species

Keywords: Cesium borates species, Ab initio calculations, Standard enthalpy of formation, Standard molar entropy, Heat capacity at constant pressure Vibrational frequencies

The coupled cluster single-double and perturbative triple (CCSD(T)) approach within highly correlated wave functions has been performed to better understand the thermochemical properties of gaseous cesium borates. Several corrections were added for core-valence correlation and relativistic effects. Structural parameters have been optimized with the B3LYP/aug-cc-pVTZ level of theory. The standard molar entropies and heat capacities at constant pressure have been also determined at the same level of theory. ∆ f H°2 98K have been obtained for the studied cesium borates species using a set of three reactions. The standard enthalpies at 298 K for CsBO 2 and Cs 2 B 2 O 4 are estimated to be (-689.0±3.8) and (-1578.5 ±7.0) kJ mol -1 , respectively. They are in good agreement with available literature data.

Numerical simulations show that in a nuclear context, CsBO 2 will be formed preferentially to CsOH in steam atmosphere for high temperature.

Introduction

The evaluation of the behaviour of fission products that can be released from the damaged fuel and transported through the reactor coolant system (RCS) to the reactor containment is a key safety issue of nuclear reactors in order to predict radiological consequences resulting from outside radioactive species. Among these fission products (FP), iodine and cesium, which can be released and transported outside in case of containment building leakages, are of particular interest due to their high radio toxicities at short term for iodine and at middle term for cesium due to mainly 131 I, 132 I, and 137 Cs isotopes. Boron, present in the RCS and for some reactors in the boron carbide control rod, can interact with cesium to form some cesium borates, [START_REF] Hontañón | The effect of chemical interactions on the transport of caesium in the presence of boron[END_REF][START_REF] Haste | Study of boron behaviour in the primary circuit of water reactors under severe accident conditions: A comparison of Phebus FPT3 results with other recent integral and separate-effects data[END_REF] which affect cesium and iodine transport through the RCS because it can prevent formation of CsI and its dimer, [START_REF] Gouëllo | 21st International Conference Nuclear Energy for New Europe[END_REF][START_REF] Kalilainen | Chemical reactions of fission product deposits and iodine transport in primary circuit conditions[END_REF] and promote gaseous iodine formation. Boron is expected to be under trivalent form (HBO 2 , B 2 O 3 , H 3 BO 3 …) under severe accident conditions. For very high temperature (> 2000 K), some other boron species like HBO and BO could exist [START_REF] Yetter | Kinetics of hightemperature B/O/H/C chemistry[END_REF], thus we have decided to make a complete study including possible boron species with valence lower than 3, possibly useful for other applications even if, it was marginal for nuclear context. In thermodynamic simulations, CsOB is obtained at trace levels (< 0.001%) for 2000 K. This paper is focused on the determination of thermochemical properties of cesium borates by computational chemistry tools because cesium borates are not well known and the literature is quite poor. These data are useful to well modelling chemical reactivity of FPs through the RCS and be able to calculate cesium chemical speciation with severe accident simulation software. This work will extent similar works already performed for other cesium species. [START_REF] Badawi | Molecular structures and thermodynamic properties of 12 gaseous cesium-containing species of nuclear safety interest: Cs 2 , CsH, CsO, Cs 2 O, CsX, and Cs 2 X 2 (X=OH, Cl, Br, and I)[END_REF][START_REF] Sudolská | Molecular structures and thermodynamic properties of monohydrated gaseous iodine compounds: Modelling for severe accident simulation[END_REF] The literature data [START_REF] Asano | Study of Vaporization Borosilicate Glasses[END_REF][START_REF] Cordfunke | Thermochemical data for reactor materials and fission products: The ECN database[END_REF][START_REF] Kulikov | An Electron Diffraction of the CsBO 2 Molecule[END_REF][START_REF] Seshadri | Infrared spectra of matrix isolated alkali metal metaborates[END_REF][12][START_REF] Lopatin | Thermodynamic properties and structure of gaseous metaborates[END_REF][START_REF] Büchler | Gaseous Metaborates. I. Mass-Spectrometric Study of the Vaporization of Lithium and Sodium Metaborates[END_REF][START_REF] Ezhov | Electron diffraction study of the molecular structure of cesium metaborate[END_REF][START_REF] Teghil | Matrix-Isolation IR Studies on Alkali Metal M(XO 2 ) Species. Bonding and Structure of Li(CO 2 ), M(ClO 2 ) and M(B02) Molecules[END_REF][START_REF] Konings | The vibrational spectrum of gaseous cesium metaborate[END_REF][18][19][START_REF] Nakajima | Thermodynamic study of gaseous CsBO 2 by Knudsen effusion mass spectrometry[END_REF][START_REF] Glushko | Termodinamicheskie svoistva individual'nykh veshchestv: Spravochnic (Thermodynamic Properties of Pure Substances: Handbook)[END_REF][START_REF] Makarov | Thermodynamic Characteristics of Gaseous Monomers of Alkali Metal Metaborates[END_REF][START_REF] Asano | Vaporization of glasses in the R 2 0-Cs 2 0-B 2 0 2 -SiO 2 (R = Li, Na, K and Rb) system[END_REF] concerning the molecular properties (geometric parameters and vibrational frequencies) and the thermochemical properties (∆ f H°2 98K , S°2 98K , C p = f(T)) for cesium borates have been collected in Tables 1-2 and are detailed here after. Our results will be discussed and respective formation of CsBO 2 and CsOH will be compared based on thermochemical equilibrium simulations.

Cesium metaborate (CsBO 2 ).

p i = kI i T/(σ i g i n i ),
where k is the pressure calibration constant, σ i is the relative ionization cross-section, g i is the gain of the electron multiplier and n i is the isotopic abundance ratio. k is obtained by an internal reference whereas σ CsBO2 was assumed to be equal to σ Cs(g) and σ Cs2(BO2)2 equal to 1.5σ Cs(g) . Enthalpies of formation were derived from these vapour pressure measurements using a second-law. Uncertainties exit concerning the ionization cross section of CsBO 2(g) , which has been assimilated to Cs (g) . The gas-transpiration method has been used by Cordfunke et al. [START_REF] Cordfunke | Vapour pressures of some caesium compounds I. CsBO 2[END_REF] to measure the vapour pressure of CsBO 2 over the solid and the liquid.

Cordfunke and Konings [START_REF] Cordfunke | Thermochemical data for reactor materials and fission products: The ECN database[END_REF] made a third-law analysis of the vapour pressure data of Gorokhov et al., [START_REF] Gorokhov | A mass spectrometric study of the vaprization of alkali metal metaborates[END_REF] and Biswas and Mukerji [START_REF] Mukerji | Vapor pressure of cesium metaborate[END_REF]. Geometrical parameters are reported in and H 3 BO 3 after heating up to 873 K. The samples were handled in an inert atmosphere glove box, the vapour pressure of CsBO 2 was deduced from the ionization current of Cs + assuming a negligible amount of the dimer in the temperature range of 792-915 K (< 1%). The modified integral method, which does not require such an ionization cross section is applied in this study to make a link between the ion current and the vapour pressure to avoid assumption on the ionization cross-section of Cs + . From the vapour pressure measurement, enthalpy of formation is derived as reported in Table 2.

In this paper, we have studied the thermochemical properties of 15 cesium metaborate species. The computational approach used in this work is detailed in section 2. The results and discussions are presented in section 3. a results are given for 11 B. b estimated value.

Table 2.

Literature data -thermochemical properties for cesium borates

Species ∆ f H° (kJ mol -1 ) S°2 98K (J mol -1 K -1 ) C p,298K (J mol -1 K -1 )

Computational methods

The equilibrium geometries of all studied species were optimized using the density functional theory (DFT) employed with the B3LYP [START_REF] Becke | Density-functional thermochemistry. III. The role of exact exchange[END_REF] functional and the augmented correlation consistent polarized triple zeta [START_REF] Dunning | Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen[END_REF] (aug-cc-pVTZ) basis set for H, B, and O atoms and with the aug-cc-pVTZ-PP basis set [START_REF] Hill | Gaussian basis sets for use in correlated molecular calculations. XI. Pseudopotential-based and all-electron relativistic basis sets for alkali metal (K-Fr) and alkaline earth (Ca-Ra) elements[END_REF] for the 9 valence electrons of the Cs atom (this basis set will be abbreviated as aug-cc-pVTZ) while the core electrons were described by the pseudopotential ECP46MDF [START_REF] Lim | All-electron and relativistic pseudopotential studies for the group 1 element polarizabilities from K to element 119[END_REF]. Vibrational frequencies and zero-point energies were determined at the same level of theory and were multiplied by the appropriate scaling factor (0.968) [START_REF] Johnson | NIST Computational Chemistry Comparison and Benchmark Database NIST Standard Reference Database Number 101[END_REF]. The stabilities of DFT wave functions were systematically checked [START_REF] Seeger | Self-consistent molecular orbital methods. XVIII. Constraints and stability in Hartree-Fock theory[END_REF][START_REF] Bauernschmitt | Stability analysis for solutions of the closed shell Kohn-Sham equation[END_REF]. All previous calculations were performed using the GAUSSIAN16 software package [34].

Electronic energies calculations were performed using the MOLPRO2015 [START_REF] Werner | version 2015, a package of ab initio programs[END_REF] software. The single and double coupled cluster theory with inclusion of a perturbative estimation for triple excitation [CCSD(T)] [START_REF] Čížek | On the Use of the Cluster Expansion and the Technique of Diagrams in Calculations of Correlation Effects in Atoms and Molecules[END_REF][START_REF] Purvis | A full coupled-cluster singles and doubles model: The inclusion of disconnected triples[END_REF][START_REF] Scuseria | An efficient reformulation of the closed-shell coupled cluster single and double excitation (CCSD) equations[END_REF][START_REF] Scuseria | Is coupled cluster singles and doubles (CCSD) more computationally intensive than quadratic configuration interaction (QCISD)?[END_REF] were performed at the DFT structures with the weighted corevalence basis set [START_REF] Peterson | Accurate correlation consistent basis sets for molecular core-valence correlation effects: The second row atoms Al-Ar, and the first row atoms B-Ne revisited[END_REF][START_REF] Deyonker | Systematically Convergent Correlation Consistent Basis Sets for Molecular Core-Valence Correlation Effects: The Third-Row Atoms Gallium through Krypton[END_REF] aug-cc-pwCVnZ (n = T, Q, 5) for H, B, and O atoms and the aug-cc-pwCVnZ-PP (n = T, Q, 5) for the cesium atom [START_REF] Hill | Gaussian basis sets for use in correlated molecular calculations. XI. Pseudopotential-based and all-electron relativistic basis sets for alkali metal (K-Fr) and alkaline earth (Ca-Ra) elements[END_REF]. We included all 9 electrons (4p 2 4p 6 5s 1 ) outside the relativistic ECP core of cesium atom in the valence electron correlation energy calculations for the cesium containing species. The CCSD(T) calculated potential energies were extrapolated to the complete basis limit(CBS) using the relation of Peterson et al. [START_REF] Peterson | Benchmark calculations with correlated molecular wave functions. IV. The classical barrier height of the H+H 2 →H 2 +H reaction[END_REF]

= + exp --1 + exp -n -1 (1)
with n= 3(awCVTZ), 4 (awCVQZ), and 5 (awCV5Z).

The CCSD(T)/CBS energies were then corrected with additional calculations. The contributions to the total energies due to core-valence correlation (CV) were calculated at the CCSD(T) level of theory using the aug-cc-pwCVTZ basis set. The scalar relativistic corrections (SR) were evaluated using the CISD/aug-cc-pVTZ level of theory (the mass-velocity and one-electron Darwin [START_REF] Davidson | Validity of first-order perturbation theory for relativistic energy corrections[END_REF] (MVD) corrections). The spin-orbit corrections (SOC)

were performed using the multi configuration interaction [START_REF] Knowles | An efficient method for the evaluation of coupling coefficients in configuration interaction calculations[END_REF] (MRCI) approach in conjunction with the aug-cc-pVTZ basis set. This correction is equal to zero for the studied species but we will consider it for reactions including OH. For OH, our calculated SOC (-0.830 kJ mol -1 ) is in excellent agreement with the values from either the NIST-JANAF thermochemical database [START_REF] Chase | NIST-JANAF Thermochemical Tables Fourth Edition Part I[END_REF] (-0.836 kJ mol -1 ) or the CASSCF/CASPT2/RASSI-SO calculation of Šulka et al. [START_REF] Šulka | A Theoretical Study of the X-Abstraction Reactions (X = H, Br, or I) from CH 2 IBr by OH Radicals: Implications for Atmospheric Chemistry[END_REF] (-0.833 kJ mol -1 ).

The enthalpy at 298 K has been calculated including the different corrections explained previously:

H 298K = E CBS + E ZPE +ddH 298K + E CV + E SR + E SO (2)
The ddH 298K correspond to the thermal contribution to enthalpy obtained using the THERMO script [47].

Results and Discussion

3.1.Structural properties and vibrational frequencies

The geometric structures of stable cesium borates species have been investigated and presented in Figures 12 ). A first study [START_REF] Vandeputte | Le rôle du bore sur la spéciation de l'iode dans le circuit primaire[END_REF] was performed to locate on the potential energy surface (PES) all conformers for each previous type of cesium borates. Among them, 15 species (Figure 1) were chosen because their relative standards Gibbs free energy at 300 K to the most stable structure were below 50 kJ mol -1 . Their geometric parameters computed at the B3LYP/aug-cc-pVTZ level of theory are given in the supplementary material (Tables S1-S3).

For To the best of our knowledge, there are no theoretical nor experimental studies on the characterization of structural properties of studied cesium borates species except the cesium metaborate CsBO 2 for which geometric parameters have been determined using gas phase electron diffraction [START_REF] Kulikov | An Electron Diffraction of the CsBO 2 Molecule[END_REF][START_REF] Ezhov | Electron diffraction study of the molecular structure of cesium metaborate[END_REF][START_REF] Komarov | Investigation of the molecular structure of alkali metal and thallium metaborates by the gas electron diffraction method: II. Average geometrical configuration of the cesium metaborate molecule[END_REF]. The main difference occurs in the type of the CsBO 2 structure.

Literature data found a non linear structure with a bond angle close to 130° while two different structures (one linear and one anchor-shaped) are obtained. It is worth mentioning that in our PES exploration, there are no minima corresponding to a non linear CsOBO structure; only a first-order saddle point was found.

Our calculated bond lengths r(B-O 1 ) and r(O 2 -B) compared well with experimental values obtained using gas phase electron diffraction [START_REF] Kulikov | An Electron Diffraction of the CsBO 2 Molecule[END_REF][START_REF] Ezhov | Electron diffraction study of the molecular structure of cesium metaborate[END_REF][START_REF] Komarov | Investigation of the molecular structure of alkali metal and thallium metaborates by the gas electron diffraction method: II. Average geometrical configuration of the cesium metaborate molecule[END_REF]. Concerning the bond length r(O 1 -Cs), the obtained values are 2.552 and 2.981 Å for CsOBO and anchor-shaped conformer, respectively. The analysis of literature data shows that r(O 1 -Cs) is about 2.7 Å, which is a median value between the two distinct structures.

The scaled vibrational frequencies of cesium borates species are gathered in Table 3. For CsOBO (CsBO 2 conf. 2), our calculated vibrational frequencies are in very good agreement with experimental results given in Table 1 with the exception of the low frequency mode. The standard enthalpies of formation at 298 K o f bo r o n -c o n t a in in g c o mpo u n d s u s ed in t h is w o r k h a v e b ee n d e t er min e d a t t h e s a me l ev e l o f t h e o r y a n d g iv e n in t h e s u ppl e me n t a r y ma t e r ia l (Ta bl e S5) w it h t h e a v a il a bl e l it e r a t u r e data of Cs, H, OH, H 2 , and H 2 O. Th e c a r t e s ia n c o o r d in a t es a n d s c a l e d v ibr a t io n a l f r e q u e n c ies o f bo r o n -c o n t a in in g s pe c ies a r e g iv in g w it h t h e r e a c t io n s u s e d t o r e t r ie v e t h e f H°2 98K o f t h es e c o mpo u n d s in t h e s u ppl e me n t a r y ma t e r ia l (Ta bl e s S6-S8). Fo r BOH, as we can see in Table S5, o u r c a l c u l a t e d f H°2 98K is in very good agreement with the CCSD(T)/CBS value obtained by Gr a n t e t Dix o n [START_REF] Grant | Heats of Formation and Bond Energies of the H (3-n ) BX n Compounds for (X = F, Cl, Br, I, NH 2[END_REF] w h il e t h e r e a r e u p t o 5 k J mo l -1 d if f e r e n c e w it h t h e o n e s calculated by Po l it zer et al. [START_REF] Politzer | Computational Determination of the Energetics of Boron Ignition/Combustion Reactions[END_REF] at either the B3PW91/6-311+G(2df,2pd) or the CBS-Q l e v e l s o f t h e o r y . Our calculated f H°2 98K for OBOH and BOBH are slightly different from literature data [START_REF] Chase | NIST-JANAF Thermochemical Tables Fourth Edition Part I[END_REF][START_REF] Politzer | Computational Determination of the Energetics of Boron Ignition/Combustion Reactions[END_REF] (differences are here greater than 10 kJ mol -1 ). These differences arise from the different methodology used in the calculations. The basis set and the functional used in the determination of the optimized structures and also the method used in the electronic energies calculations are different and our methodology has been validated in our previous works [START_REF] Fortin | Reactivity of Hydrogen Peroxide with Br and I Atoms[END_REF][START_REF] Khiri | Theoretical Study of the Reactions of H Atoms with CH 3 I and CH 2 I 2[END_REF] and it shows a very good agreement with experimental and theoretical literature values. The standard molar entropies (S°2 98K ) and heat capacities at constant pressure (C p = f (T))

were calculated over a wide range of temperature (300-1500 K) with the THERMO software [47] using the rigid-rotor harmonic oscillator approximation from optimized geometries, scaled vibrational frequencies, and moments of inertia at the B3LYP/aug-cc-pVTZ level of theory. Corresponding values are gathered in the supplementary material (Table S9). For CsBO 2 , our values for the linear conformer (S°2 98K = 316.16 and C p,298 = 64.14 J mol -1 K -1 ) are in good agreement with the literature values [START_REF] Cordfunke | Thermochemical data for reactor materials and fission products: The ECN database[END_REF]18,19] as shown in Table 2. Our values for CsOBO and anchor-shaped structure are slightly different from those estimated by Asano et al. [START_REF] Asano | Study of Vaporization Borosilicate Glasses[END_REF] (S°2 98K = 284.675 and C p,298 = 38.250 J mol -1 K -1 ).

Main dominant species in nuclear severe accident conditions

As we can see in Table 5, there are a large number of gaseous cesium borates, which have been identified. For a first ranking, species for which enthalpy of formation is positive are not expected to be formed. For others, we wish to know what the main species which could be formed in the RCS during a nuclear severe accident. Thermochemical equilibrium computations are performed based on thermochemical data determined in this work. The boundary conditions of the simulations are as follows: temperature range between 2000 and 500 K is explored, a pressure of 5 bar, the gas is steam and the cesium and boron are fixed respectively to 5 ×10 -6 and 2 ×10 -5 moles corresponding to respective molar ratios of 0.0025 and 0.0099. Results are reported in Table 6. In the temperature range where gaseous cesium borates are formed, the main species are CsBO 2 with the two conformers, CsBO 2 dimers are formed but in a lesser extent, around 1.7%. Cesium hydroxide is less stable than CsBO 2 . For lower temperature, CsBO 2 is under condensing form. These chemical equilibrium results are consistent with previous simulations of Minato [START_REF] Minato | Thermodynamic analysis of cesium and iodine behavior in severe light water reactor accidents[END_REF].

Table 6.

Percentages of cesium present in gas phase and condensed forms as a function of temperature. If we consider potential role of boron, it can promote gaseous iodine formation according to i) restricted formation of CsI and ii) reactivity between cesium iodide and boron species to form gaseous iodine. These two points are detailed here after.

The restricted formation of CsI can be expected if the equilibrium constants of the reactions I and II are compared.

CsOH + HI CsI + H 2 O (I) CsBO 2 + HI CsI + HBO 2 (II)

The equilibrium constant of reaction II is about four order de magnitude lower than reaction I at 1500 K. If we examine the reactivity between CsI and boron species, it appears that metaboric acid is the most stable boron species at high temperature [START_REF] White | Infrared Emission Spectrum of Gaseous HBO 2[END_REF] (typically 1500 K) whereas for lower temperature the most stable forms are H 3 BO 3 and in a less extent H 3 B 3 O 6 .

The main chemical reaction to consider is: CsI + HBO 2 CsBO 2 + HI (III)

For reaction III, the Gibbs free energy of reaction is negative for T higher than 1700 K thus the formation of HI by this chemical pathway is only possible at high temperature. Other possible chemical pathways at moderate temperature, around 1000 K, are:

Fig. 1 Fig. 2

 12 Fig. 1 Structures of cesium borates species (CsBO, CsBO 2 , CsB 2 O, CsB 2 O 2 , and CsB 2 O 4 )

  Asano et al. [8] heated, in a platinum Knudsen cell, synthesised glasses containing Na 2 O, B 2 O 3 , SiO 2 and Cs 2 O in the temperature range 845-1150 K and generated vapour pressure of about 4 10 -5 Pa. Pressure of CsBO 2(g) and Cs 2 B 2 O 4(g) were identified and measured by mass spectrometry. The intensities of the Cs + and CsBO 2 + ions were used for the determination of the partial pressure of CsBO 2(g) and Cs 2 (BO 2 ) 2 + intensity for the partial pressure of Cs 2 (BO 2 ) 2(g) but for some glasses, Cs 2 (BO 2 ) 2(g) could not be measured due to too low intensity of the ionization pic. Usual relationship is used to make a link between the measured ion intensity I i , of the vapour species i into the corresponding partial pressure p i , at the sample temperature T by the relation, as follows:

Table

Table 1 .

 1 Literature data -geometric parameters and vibrational frequencies for cesium metaborate (CsBO 2 ) a

	Geometric parameters	Vibrational frequencies (cm -1 )	Technique	Reference
	r(Cs-O) = 2.650 Å				Review	Cordfunke and Konings [9]
	r(B-O) = 1.270 Å
	90< θ(°)<143
	Recommended θ = 130°
	r(Cs-O) = (2.71 ± 0.04) Å				Gas-phase electron diffraction	Kulikov et al. [10]
	r(B-O) = 1.266 Å
	θ = 127°
		60 b 206 576 581	1077	1945	IR with isolation matrix (Ar, Kr)	Seshadri et al. [11]
	r(Cs-O) = (2.62 ± 0.03) Å				Gas-phase electron diffraction	Ezhov and Komarov (1978) [12]
	r(B-O) = (1.25 ± 0.01) Å
	θ = (135 ± 8.8)°
		610		1935	IR with gas phase	Büchler and Marram [14]
	r(Cs-O) = 2.69 Å				Gas-phase electron diffraction	Ezhov and Komarov (2000) [15]
	r(B-O) = 1.30 Å
	r(B-O) = 1.22 Å
	θ = 125°
		207 577 581	1077	1945	IR with isolation matrix (Ar)	Teghil et al. [16]
		233 588 632	1045	1933	IR with gas phase	Konings et al. [17]

  . They can be classed into two groups, the first one corresponds to the monomers (CsBO, CsBO 2 , CsB 2 O, CsB 2 O 2 , and CsB 2 O 4 ) and the second one corresponds to the dimers (Cs 2 B 2 O 2 and Cs 2 B 2 O 4

  CsBO, only two stable structures (linear CsBO and non linear CsOB) have been found on the PES. Two stable conformers have been localised on the PES of CsBO 2 species. One presents a linear structure Cs-O 1 -B-O 2 while the other one is anchor-shaped. The exploration of the PES of CsB 2 O shows that there are three stable structures to be considered. Among them, one corresponds to a linear Cs-B 1 -B 2 -O structure. Only one stable structure has been found for the CsB 2 O 2 compound. Two stable structures have been identified for the CsB 2 O 4 compounds. For the Cs 2 B 2 O 2 , three dimers have been characterized with cyclic structures and for the Cs 2 B 2 O 4 , two stable dimers have been found.

Table 3 .

 3 SR values are too weak, they are between -0.83 a n d 0.16 kJ mol -1 . The E SO are equal to zero except for reactions including OH radical ( E SO = -0.83 kJ mol -1 ). The most important contribution corresponds to the ∆E ZPE + d d H 298K .

	B3LYP/aug-cc-pVTZ scaled vibrational frequencies of cesium borates species

CBS , ∆E ZPE + d d H 298K , ∆E CV , ∆E SR , and ∆E SO ), the enthalpies at 298 K (∆ r H°2 98K ) a n d t h e s t a n d a r d e n t h a l py o f f o r ma t io n a r e c o l l ec t e d in Ta b l e 4 f o r t h e d if f e r e n t s t u d ie d r e a c t io n s o f Cs BO 2 a n d Cs 2 B 2 O 4 . Va l u e s c o r r es po n d in g t o o t h er cesium borates species are giving in the s u ppl e me n t a r y ma t e r ia l (Ta bl e S4). As s h o w n in Ta bl e 4, t h e in c l u s io n o f t h e d if f er e n t en er g y c o n t r ib u t io n s is v er y impo r t a n t t o o b t a in t h e a c c u r a t e v a l u e o f t h e s t a n d a r d e n t h a l py o f f o r ma t io n a t 298 K. Th e E

Table 4 .

 4 Standard reaction enthalpies at 298 K, standard enthalpy of formation at 298 K with different energy contributions for the various studied reactions of CsBO 2 and Cs 2 B 2 O 4 (kJ mol -1 )

	Reactions	∆E CBS	∆E ZPE +ddH 298K	∆E CV	∆E SR	∆E SO	∆ r H°2 98K	∆ f H°2 98K
			CsBO 2 conf. 1
	CsBO 2 + H → Cs + OBOH	-15.23	21.90	-0.13	0.15	0.00	6.69	-692.85
	2 CsBO 2 + H 2 O → 2 Cs + B 2 O 4 H + OH	816.38	-0.40	1.46	-0.83	-0.83	815.95	-692.60
	2 CsBO 2 + H → 2 Cs + B 2 O 4 H	291.17	27.26	0.56	0.05	0.00	319.04	-692.81
			CsBO 2 conf. 2
	CsBO 2 + H → Cs + OBOH	-17.69	20.78	-0.22	0.07	0.00	2.93	-689.10
	2 CsBO 2 + H 2 O → 2 Cs + B 2 O 4 H + OH	811.45	-2.64	1.29	-0.83	-0.83	808.43	-688.85
	2 CsBO 2 + H → 2 Cs + B 2 O 4 H	286.23	25.02	0.39	-0.11	0.00	311.53	-689.05
			Cs 2 B 2 O 4 dim. 1
	Cs 2 B 2 O 4 + H → 2 Cs + B 2 O 4 H	625.88	16.91	-0.44	0.11	0.00	-642.46	-1578.80
	Cs 2 B 2 O 4 + H 2 O → 2 Cs + B 2 O 4 H + OH	1014.74	-5.44	1.28	-0.83	-0.83	1008.93	-1578.19
	Cs 2 B 2 O 4 + H 2 → 2 Cs + 2 OBOH	489.53	22.22	0.39	-0.11	0.00	512.03	-1578.61
			Cs 2 B 2 O 4 dim. 2
	Cs 2 B 2 O 4 + H → 2 Cs + B 2 O 4 H	583.11	17.10	-0.31	0.16	0.00	600.06	-1536.40
	Cs 2 B 2 O 4 + H 2 O → 2 Cs + B 2 O 4 H + OH	971.98	-5.24	1.40	-0.77	-0.83	966.54	-1535.80
	Cs 2 B 2 O 4 + H 2 → 2 Cs + 2 OBOH	446.77	22.42	0.51	-0.06	0.00	469.64	-1536.21

Table 5

 5 H°0 K = -1577 ± 41.1 kJ mol -1 ). We have to notice that for the values of the standard enthalpies of formation a v a il a bl e in t h e l it er a t u r e f o r bo t h CsBO 2 and Cs 2 B 2 O 4 , the experimental uncertainties were very high.

	gathers the calculated standard enthalpies of formation and their corresponding
	uncertainties at 298 K of studied cesium borates species. Our calculated standard enthalpies of
	formation at 298 K for the two stable conformers of CsBO 2 ((-689.0± 3.8) a n d (-692.7± 3.8)
	kJ mol -1 for CsOBO and anchor-shaped structure, respectively) are in very good agreement
	with the experimental value given by Asano et al. [8] (-690.1 ± 20.3) kJ mol -1 using mass
	spectrometric Knudsen effusion method. It compares well also with other available literature
	data as shown in Table 1. Our calculated f H°0 K ((-687.6± 3.8) a n d (-690.2± 3.8) kJ mol -1
	for CsOBO and anchor-shaped structure, respectively) are close to those obtained by Makarov

Table 5 .

 5

	Calculated standard enthalpies of formation at 298 K a n d a s s o c ia t e d u n c e r t a in t ie s
	o f c es iu m b o r a t e s species
	Species	∆ f H°2 98K (kJ mol -1 )	∆∆ f H°2 98K (kJ mol -1 )
	CsBO	-197.8	4.3
	CsOB	-218.0	4.3
	CsBO 2 conf.1	-692.7	3.8
	CsBO 2 conf.2	-689.0	3.8
	CsB 2 O conf. 1	142.6	6.5
	CsB 2 O conf. 3	161.5	6.5
	CsB 2 O conf. 2	175.5	6.5
	CsB 2 O 2	-447.4	6.5
	CsB 2 O 4 conf. 1	-1188.0	6.5
	CsB 2 O 4 conf. 2	-946.4	6.5
	Cs 2 B 2 O 2 dim.1	-589.0	7.0
	Cs 2 B 2 O 2 dim. 2	-588.0	7.0
	Cs 2 B 2 O 2 dim.2	-584.5	7.0
	Cs 2 B 2 O 4 dim. 1	-1578.5	7.0
	Cs 2 B 2 O 4 dim. 2	-1536.1	7.0
	3.3. Standard molar entropies at 298 K and heat capacities at constant pressure

  CsOH thermodynamic data values from Cordfunke and Konings[START_REF] Cordfunke | Thermochemical data for reactor materials and fission products: The ECN database[END_REF] and b ∆ f H° (CsBO 2 cond) = -976.8 kJ.mol -1 at 298 K[START_REF] Cordfunke | The thermodynamic properties of cesium metaborate CsBO 2 from 5 to 1000 K[END_REF]

	Temperature (K)	CsBO 2 conf. 1	CsBO 2 conf. 2	Cs 2 B 2 O 4 dim. 1	Cs 2 B 2 O 4 dim. 2	CsOH a	CsBO 2 cond. b
	2000	28.6	68.9	0.02	0.05	2.4	0.
	1500	32.9	65.1	0.84	0.77	0.43	0.
	1000	0.24	0.34	0.05	0.01	0.03	99.3
	500	0.	0.	0.	0.	0.	100.0

a

ACKNOWLEDGMENTS

We t h a n k t h e Ce n t r e d e Re s s o u r c es In f o r ma t iq u e s (CRI) o f t h e Un iv er s it y o f Lil l e a n d t h e Ce n t r e Rég io n a l In f o r ma t iq u e e t d 'Appl ic a t io n s Nu mér iq u e s d e No r ma n d ie (CRIANN) f o r pr o v id in g c o mpu t er t ime f o r pa r t o f t h e t h e o r e t ic a l c a l c u l a t io n s . Th is w o r k w a s pa r t o f t h e LABEX Ca PPA (Ch e mic a l a n d Ph y s ic a l Pr o pe r t ie s o f t h e At mo s ph er e), w h ic h is f u n d ed by t h e Fr e n c h Na t io n a l Res e a r c h Ag en c y (ANR) t h r o u g h t h e PIA (Pr o g r a mme d 'In v e s t is s e me n t d 'Av en ir ) u n d e r c o n t r a c t "ANR-11-LABX-0005-01 a n d . a l s o s u ppo r t ed b y t h e Re g io n a l Co u n c il "No r d -Pa s d e Ca l a is " a n d t h e "Eu r o pea n Fu n d s f o r Reg io n a l Ec o n o mic De v e l o pme n t ". We t h a n k a l s o Br u n o PIAR f r o m IRSN f o r f r u it f u l d is c u s s io n .

CsI + BO 2 CsBO 2 + I (IV)

where Gibbs free energies are negative in the temperature range 2000-500 K. The reaction IV should be controlled by BO 2 formation, which is limited with respect to HBO 2 .

The last reaction V corresponds to reactivity between boric acid and cesium iodide:

with negative Gibbs free energies for T> 1300 K.

All these considerations do not include any kinetics. It may impact the iodine speciation, which be shared between CsI and gaseous iodine.

Conclusions

The thermochemical properties at 298 K are predicted for cesium borates species at the CCSD(T)/CBS(T,Q,5)//B3LYP/aug-cc-pVTZ level of theory. The standard enthalpies of formation at 298 K have been calculated using a set of three reactions. The standard molar entropies and heat capacities at constant pressure have been calculated at the B3LYP/aug-cc-pVTZ level of theory.

Notes

The authors declare no competing financial interest.