Vladimir I Korobov

J.-Ph Karr

Rovibrational spin-averaged transitions in the hydrogen molecular ions

We reconsider the calculation of rovibrational transition frequencies in hydrogen molecular ions. Some previously neglected contributions, such as the deuteron polarizability, are included into consideration in comparison with our previous work. In particular, one-loop and two-loop QED corrections at mα 7 and mα 8 orders are recalculated in the framework of the adiabatic approximation, with systematic inclusion of corrections associated with vibrational motion. Improved theoretical transitions frequencies are obtained and found to be in very good agreement with recent highprecision spectroscopy experiments in HD + . New values for the mp/me and m d /mp mass ratios are determined.

I. INTRODUCTION

Considerable progress has been recently achieved in the spectroscopy of hydrogen molecular ions. Several rovibrational transitions in HD + have been measured with relative uncertainties in the 10 -11 -10 -12 range through Dopplerfree spectroscopy of ultracold trapped ions in the Lamb-Dicke regime [START_REF] Alighanbari | Precise test of quantum electrodynamics and determination of fundamental constants with HD + ions[END_REF][START_REF] Patra | Proton-electron mass ratio from laser spectroscopy of HD + at the part-per-trillion level[END_REF][START_REF] Kortunov | Proton-electron mass ratio by highresolution optical spectroscopy of ion ensemble in the resolved-carrier regime[END_REF], approaching or even exceeding the precision of theoretical predictions [START_REF] Korobov | Fundamental transitions and ionization energies of the hydrogen molecular ions with few ppt uncertainty[END_REF]. This has allowed to get improved determinations of the proton-to-electron mass ratio, and to perform tests of QED constraining hypothetical "fifth forces" between hadrons [START_REF] Alighanbari | Precise test of quantum electrodynamics and determination of fundamental constants with HD + ions[END_REF][START_REF] Germann | Three-body QED test and fifth-force constraint from vibrations and rotations of HD +[END_REF]. These results, as well as other ongoing projects [START_REF] Zhong | High-precision spectroscopy of hydrogen molecular ions[END_REF][START_REF] Schmidt | Trapping, Cooling, and Photodissociation Analysis of State-Selected H + 2 Ions Produced by (3 + 1) Multiphoton Ionization[END_REF], and perspectives of reaching higher precision [START_REF] Schiller | Simplest Molecules as Candidates for Precise Optical Clocks[END_REF][START_REF] Karr | H + 2 and HD + : Candidates for a molecular clock[END_REF] e.g. using quantum-logic spectroscopy schemes [START_REF] Schmidt | Spectroscopy using quantum logic[END_REF][START_REF] Wolf | Non-destructive state detection for quantum logic spectroscopy of molecular ions[END_REF][START_REF] Chou | Frequency-comb spectroscopy on pure quantum states of a single molecular ion[END_REF] provide strong motivation to improve the theory further.

In our paper [START_REF] Korobov | Fundamental transitions and ionization energies of the hydrogen molecular ions with few ppt uncertainty[END_REF], we calculated the frequencies of fundamental vibrational transitions in the nonrelativistic QED approach, including corrections up to the mα 8 order. Since then, several new advances in the theory of hydrogenlike atoms have been achieved [START_REF] Yerokhin | Theory of the Lamb Shift in Hydrogen and Light Hydrogen-Like Ions[END_REF][START_REF] Szafron | Virtual Delbrück scattering and the Lamb shift in light hydrogenlike atoms[END_REF][START_REF] Karshenboim | The Lamb shift of the 1sstate in hydrogen: Two-loop and three-loop contributions[END_REF], which allows to get improved results for the corresponding correction terms in the hydrogen molecular ions. In addition, we found that several contributions that had been neglected in our previous consideration were of comparable magnitude to the estimated error bar, and thus should be included. Finally, changes in the recommended values of fundamental constants, the nucleus-to-electron mass ratios but also the Rydberg constant, proton and deuteron radii, between the previous (2014) CODATA adjustment [START_REF] Mohr | CODATA recommended values of the fundamental physical constants[END_REF] and the most recent one (2018) [START_REF] Tiesinga | CODATA recommended values of the fundamental physical constants[END_REF] (see Table I) also affect our theoretical predictions and their uncertainties.

The aim of this work is to reanalyse the theory, with particular emphasis on the evaluation of QED corrections at orders mα 7 and mα 8 in the framework of the adiabatic approximation, where we systematically include "vibrational" corrections, i.e. the second-order perturbation terms due to perturbation of the vibrational wavefunction. Improved theoretical rovibrational transition frequencies are given for experimentally relevant transitions, and the impact of these new results on the determination of the proton-to-electron mass ratio is illustrated.

We use atomic units throughout this paper (h = m e = e = 1). Other constants used in calculations are taken from the CODATA18 adjustment [START_REF] Tiesinga | CODATA recommended values of the fundamental physical constants[END_REF] (see in Table I), including the fine structure constant, α = 7.297 352 5693(11)×10 -3 .

II. NUCLEAR SIZE AND POLARIZABILITY CORRECTIONS

In our previous calculations [START_REF] Korobov | Fundamental transitions and ionization energies of the hydrogen molecular ions with few ppt uncertainty[END_REF][START_REF] Korobov | Leading-order relativistic and radiative corrections to the rovibrational spectrum of H + 2 and HD + molecular ions[END_REF], we only included the leading-order nuclear finite-size correction, see Eq. ( 6) in [START_REF] Korobov | Leading-order relativistic and radiative corrections to the rovibrational spectrum of H + 2 and HD + molecular ions[END_REF]. Some higher-order nuclear corrections are not negligible at the current level of theoretical accuracy, in particular the deuteron polarizability [START_REF] Friar | Higher-order nuclear-polarizability corrections in atomic hydrogen[END_REF]. Here, we follow the notations used in [START_REF] Mohr | CODATA recommended values of the fundamental physical constants[END_REF]. According to Eq. (34) of [START_REF] Mohr | CODATA recommended values of the fundamental physical constants[END_REF]), we write the mα 5 deuteron polarizability contribution as

E (5) pol (D) = [-21.37(8)] πδ(r d ) (kHz). ( 1
)
For example, this results in a 0.33 kHz shift for the frequency of the fundamental vibrational transition (L = 0, v = 0) → (0, 1), which is comparable to the overall theoretical uncertainty of 0.5 kHz for this transition (see Table II). Nuclear finite-size corrections at the same order [START_REF] Friar | Higher-order nuclear-size corrections in atomic hydrogen[END_REF] are written as (see Eq. (59) in [START_REF] Mohr | CODATA recommended values of the fundamental physical constants[END_REF])

E (5) fns (D) = -(2R ∞ c) 2π 3 C η R d a 0 3 πδ(r d ) = [-0.57(3)] πδ(r d ) (kHz), (2)
where a 0 is the Bohr radius, R d = 2.12799(74) the rms charge radius of the deuteron, and C η = 2.0 [START_REF] Alighanbari | Precise test of quantum electrodynamics and determination of fundamental constants with HD + ions[END_REF]. The contribution at the next order (mα 6 ), unlike the previous ones, is not a "state-independent" term proportional to the squared value of the wavefunction at the electron-nucleus, and would require an independent calculation. However, this term can be estimated from its value for the hydrogen atom ground state by using the LCAO approximation for the electronic wavefunction:

ψ LCAO (r) = 1 √ 2 ψ 1s (r p )+ψ 1s (r d ) ,
where ψ 1s is the hydrogen ground state wavefunction. Under this approximation, one gets from Eq. (59) of [START_REF] Mohr | CODATA recommended values of the fundamental physical constants[END_REF]:

E (6) fns (D) = -(2R ∞ c) 2π 3 R d a 0 2 (Zα) 2 C θ -ln ZR d a 0 πδ(r d ) = [3.96(2)] πδ(r d ) (kHz). (3)
with C θ = 0.38 [START_REF] Korobov | Fundamental transitions and ionization energies of the hydrogen molecular ions with few ppt uncertainty[END_REF]. Only the uncertainty from C θ is indicated in this equation. Due to the employed LCAO approximation, one may estimate the uncertainty of E

fns (D) as equal to the nonlogarithmic term, which is still much smaller than the overall theoretical uncertainty.

In the proton case, all the corrections considered above are negligibly small at the current level of accuracy.

III. mα 7 AND mα 8 -ORDER CORRECTIONS IN THE ADIABATIC APPROXIMATION

Relativistic and QED corrections at the orders mα 4 to mα 6 have been evaluated in a full three-body approach using precise variational wavefunctions [START_REF] Korobov | Leading-order relativistic and radiative corrections to the rovibrational spectrum of H + 2 and HD + molecular ions[END_REF][START_REF] Korobov | Bethe logarithm for the H + 2 and HD + molecular ions[END_REF], except for the mα 6 relativistic correction [START_REF] Korobov | Fundamental transitions and ionization energies of the hydrogen molecular ions with few ppt uncertainty[END_REF][START_REF] Korobov | Relativistic corrections of order mα 6 to the two-centre problem[END_REF]. For calculation of mα 7 and higher order one-and two-loop corrections we use the Born-Oppenheimer approach, where the states of the molecule are taken in the form

Ψ BO = φ el (r; R)χ BO (R). (4)
The electronic wavefunction obeys the clamped nuclei Schrödinger equation for a bound electron

[H el -E el (R)] φ el = 0, (5)
where

H el = p 2 2m + V + Z 1 Z 2 R , V = - Z 1 r 1 - Z 2 r 2 .
Here H el is the electronic Hamiltonian, Z 1 and Z 2 are the charges of the nuclei, and r 1 , r 2 are the distances from the electron to nuclei 1 and 2, respectively. The wavefunction χ BO (R) describes the relative nuclear motion, and is a solution of

(H vb -E 0 ) χ BO = - ∇ 2 R 2µ N +E el (R)-E 0 χ BO = 0, (6)
where

µ N = M 1 M 2 /(M 1 +M 2 )
is the reduced mass of the nuclei. Instead of the Born-Oppenheimer solution χ BO (R) we use the adiabatic solution χ ad (R), which includes as well the adiabatic corrections

E ad (R) = E el + dr φ el p 2 8µ N + P 2 2µ N - κ 2µ N pP φ el ( 7
)
where p is the electron impulse in the center-of-mass frame, P the relative impulse of the two nuclei, and κ = (M 1 -M 2 )/(M 1 + M 2 ) is the asymmetry parameter. See Ref. [START_REF] Wolniewicz | The vibration-rotational energies of the hydrogen molecular ion HD +[END_REF], or a review by Carrington et al. [START_REF] Carrington | Spectroscopy of the hydrogen molecular ion[END_REF] for more details.

The one-loop self-energy correction of order mα 7 to the energy of a bound electron in the two-center problem (non-recoil limit) has been determined in [START_REF] Czarnecki | Calculation of the One-and Two-Loop Lamb Shift for Arbitrary Excited Hydrogenic States[END_REF][START_REF] Jentschura | Nonrelativistic QED approach to the Lamb shift[END_REF][START_REF] Korobov | Calculation of the relativistic Bethe logarithm in the two-center problem[END_REF][START_REF] Korobov | mα 7 -Order Corrections in the Hydrogen Molecular Ions and Antiprotonic Helium[END_REF][START_REF] Korobov | Theoretical transition frequencies beyond 0.1 ppb accuracy in H + 2 , HD + , and antiprotonic helium[END_REF]. The electronic part of the correction can be calculated using the effective Hamiltonian of Eq. ( 6) from Ref. [START_REF] Korobov | Theoretical transition frequencies beyond 0.1 ppb accuracy in H + 2 , HD + , and antiprotonic helium[END_REF]:

∆E (7) el-SE = χ ad E (7) 1loop-SE (R) χ ad , (8)
numerical data for the E

1loop-SE (R) effective potential may be found in the Supplemental Material to [START_REF] Korobov | Calculation of the relativistic Bethe logarithm in the two-center problem[END_REF]. The one-loop vacuum polarization (Uehling potential) contribution has been considered in [START_REF] Karr | One-loop vacuum polarization at mα 7 and higher orders for three-body molecular system[END_REF]. The adiabatic approximation was also compared with full three-body calculations, confirming that it is accurate to O(m/M ) (where m/M is the electron-to-nucleus mass ratio). The electronic part of the correction can be written as

∆E (7+) el-VP = χ ad E (7+) 1loop-VP (R) χ ad , (9)
where

E (7+)
1loop-VP (R) is given by Eq. ( 16) of [START_REF] Karr | One-loop vacuum polarization at mα 7 and higher orders for three-body molecular system[END_REF]. As shown in [START_REF] Korobov | Fundamental transitions and ionization energies of the hydrogen molecular ions with few ppt uncertainty[END_REF][START_REF] Karr | One-loop vacuum polarization at mα 7 and higher orders for three-body molecular system[END_REF], beyond the above electronic contributions one also needs to include vibrational corrections. The latter are second-order perturbation terms stemming from perturbation of the vibrational wavefunction by the leading relativistic and radiative corrections to the adiabatic potential E ad (R), namely

E (4) BP (R) = α 2 - p 4 8m 3 + πρ 2m 3 + H so R , E (5)
SE (R) = α 3 4 3 ln 1 α 2 -β(R) + 5 6 Z 1 δ(r 1 ) + Z 2 δ(r 2 ) R , E (5)
VP (R) = -α 3 4 15 Z 1 δ(r 1 ) + Z 2 δ(r 2 ) R . (10)
Here ρ = ∇ 2 V /(4π), H so is the electron spin-orbit Hamiltonian (see [START_REF] Korobov | Relativistic corrections of order mα 6 to the two-centre problem[END_REF] for details), and β(R) is the nonrelativistic Bethe logarithm for the bound electron in the two-center problem, whose values as a function of R may be found in the Supplemental Material to Ref. [START_REF] Korobov | Calculation of the relativistic Bethe logarithm in the two-center problem[END_REF] or in [START_REF] Bukowski | Bethe logarithm and Lamb Shift for the Hydrogen Molecular Ion[END_REF]. The mα 7 -order vibrational correction from one-loop self-energy and vacuum polarization is then obtained via the second-order perturbation formalism as

∆E (7) vb = 2 χ ad E (4) BP (R)Q (E 0 -H vb ) -1 Q E (5) SE (R) + E (5) VP (R) χ ad , (11)
here Q is a projection operator onto a subspace orthogonal to χ ad (R) from Eq. ( 5). The difference with respect to our previous calculations is that E

SE (R) in Eq. ( 10) includes the contribution from the electron anomalous magnetic moment, which was missed in Eq. ( 10) of [START_REF] Korobov | Fundamental transitions and ionization energies of the hydrogen molecular ions with few ppt uncertainty[END_REF].

We now turn to the mα 8 -order, starting with the one-loop self-energy correction. The higher-order remainder (mα 8 and above) for this contribution can be estimated from the hydrogen 1S state results by using the LCAO approximation:

∆E (8+) el-1loop = α 5 χ ad G SE (1S) -A 60 Z 3 1 δ(r 1 ) + Z 3 2 δ(r 2 ) χ ad , (12)
where G SE (1S) = -30.29024(2) is the higher order remainder, and A 60 = -30.924149. Both numbers are taken from Table II of [START_REF] Yerokhin | Theory of the Lamb Shift in Hydrogen and Light Hydrogen-Like Ions[END_REF]. This is more accurate than the treatment of Ref. [START_REF] Korobov | Fundamental transitions and ionization energies of the hydrogen molecular ions with few ppt uncertainty[END_REF], where only the mα 8 -order term was estimated. The theoretical uncertainty is estimated as equal to ∆E

(8+) el-1loop -∆E (8log)
el-1loop , where ∆E

el-1loop is the known mα 8 -order logarithmic term (first term of Eq. ( 14) in [START_REF] Korobov | Fundamental transitions and ionization energies of the hydrogen molecular ions with few ppt uncertainty[END_REF]).

The second-order vibrational contribution is expressed

∆E (8) vb-1loop = 2 χ ad E (4) BP (R)Q (E 0 -H vb ) -1 Q E (6) SE (R) + E (6) VP (R) χ ad , (13)
where

E (6) SE (R) = α 4 139 32 -2 ln 2 π Z 2 1 δ(r 1 ) + Z 2 2 δ(r 2 ) R , E (6)
VP (R) = α 4 5 48 π Z 2 1 δ(r 1 ) + Z 2 2 δ(r 2 ) R .
Finally, we consider two-loop mα 8 -order corrections. For hydrogen-like atoms, the two-loop correction at orders mα 8 and higher is generally expressed in the form [START_REF] Yerokhin | Theory of the Lamb Shift in Hydrogen and Light Hydrogen-Like Ions[END_REF][START_REF] Tiesinga | CODATA recommended values of the fundamental physical constants[END_REF]

E (8+) 2loop = (Zα) 6 π 2 n 3 B 63 ln 3 (Zα) -2 +B 62 ln 2 (Zα) -2 + B 61 ln(Zα) -2 + G 2loop (Zα) , (14)
where G 2loop (Zα) is the higher order remainder calculated in [START_REF] Yerokhin | Two-loop self-energy for the ground state of medium-Z hydrogenlike ions[END_REF][START_REF] Yerokhin | Two-loop self-energy in the Lamb shift of the ground and excited states of hydrogenlike ions[END_REF]. We adopt similar notations for hydrogen molecular ions:

E (8+) 2loop = α 6 π χ ad B 63 (R) ln 3 (α -2 )+B 62 (R) ln 2 (α -2 )+B 61 (R) ln(α) -2 +G 2loop (1S) Z 3 1 δ(r 1 )+Z 3 2 δ(r 2 ) χ ad . ( 15
)
Again, the higher-order remainder is estimated using the LCAO approximation. In our calculations, we adopted the value G 2loop (1S) = -94.5(6.6) from [START_REF] Karshenboim | The Lamb shift of the 1sstate in hydrogen: Two-loop and three-loop contributions[END_REF]. This is more accurate than our previous treatment [START_REF] Korobov | Fundamental transitions and ionization energies of the hydrogen molecular ions with few ppt uncertainty[END_REF], where only the mα 8 -order correction was estimated using B 60 (1S) instead of G 2loop (1S). The theoretical uncertainty is estimated as equal to the term proportional to G 2loop (1S), after subtraction of the known term of order mα 9 ln 2 (α) (Eq. ( 26) of [START_REF] Yerokhin | Theory of the Lamb Shift in Hydrogen and Light Hydrogen-Like Ions[END_REF]). Calculation of the B 6k (R) effective potentials for the two-center problem is described in [START_REF] Korobov | Fundamental transitions and ionization energies of the hydrogen molecular ions with few ppt uncertainty[END_REF]. Since then, a new contribution to the B 61 coefficient in hydrogen-like atoms from light-by-light scattering diagrams has been found, yielding a correction B LbL 61 = 0.830 309 for S states [START_REF] Szafron | Virtual Delbrück scattering and the Lamb shift in light hydrogenlike atoms[END_REF]. For the two-center problem we thus add the following term to B 61 (R):

B LbL 61 (R) = B LbL 61 Z 3 1 δ(r 1 ) + Z 3 2 δ(r 2 ) R . (16)
The second-order contribution due to vibrational motion is expressed as

∆E (8) vb-2loop = 2 χ ad E (4) BP (R)Q (E 0 -H vb ) -1 Q E (6) 2loop (R) χ ad , + χ ad E (5) SE (R) + E (5) VP (R) Q (E 0 -H vb ) -1 Q E (5) SE (R) + E (5) VP (R) χ ad , (17)
with

E (6) 2loop (R) = α 4 π 2 [0.538941] π Z 1 δ(r 1 ) + Z 2 δ(r 2 ) R .
Due to the presence of a logarithmic term in E

(5) SE (R), the vibrational correction is enhanced by a factor of ln 2 (α -2 ) and contributes to the B 62 , B 61 and nonlogarithmic terms. It results in a 1.14 kHz shift for the fundamental vibrational transition, thus the neglect of this term in [START_REF] Korobov | Fundamental transitions and ionization energies of the hydrogen molecular ions with few ppt uncertainty[END_REF] was not justified.

The last corrections requiring new consideration are the muonic and hadronic vacuum polarization corrections. The muonic term is (see Eq. ( 14) of [START_REF] Yerokhin | Theory of the Lamb Shift in Hydrogen and Light Hydrogen-Like Ions[END_REF])

E µVP (R) = m e m µ 2 E (5) VP (R) , (18)
and the hadronic term may be written as (Eq. ( 15) of [START_REF] Yerokhin | Theory of the Lamb Shift in Hydrogen and Light Hydrogen-Like Ions[END_REF])

E hadVP (R) = 0.671(15) E µVP (R). ( 19
)
The sum of these two contributions shifts the fundamental transition frequency by 0.25 kHz.

IV. RESULTS

Our results for the frequency of the fundamental transition (L = 0, v = 0) → (0, 1) in HD + are presented in Table II and compared with previous results from Ref. [START_REF] Korobov | Fundamental transitions and ionization energies of the hydrogen molecular ions with few ppt uncertainty[END_REF]. The change in the nonrelativistic transition frequency ν nr is mainly due to those of the nucleus-to-electron mass ratios (mostly µ p ) and Rydberg constant between the 2014 and 2018 CODATA adjustments. The shift in ν α 2 is due to the nuclear corrections; it stems from the CODATA18 values of the proton and deuteron radii, and from the inclusion of higher-order finite-size and polarizability corrections described in Sec. II. The change in ν α 5 comes from the correction of the vibrational contribution [Eq. [START_REF] Wolf | Non-destructive state detection for quantum logic spectroscopy of molecular ions[END_REF]]. Finally, several improvements have been made in the calculation of ν α 6 as detailed in Sec.III: inclusion of vibrational contributions, estimate of the all-order remainder both in one-loop and too-loop corrections, as well as the inclusion of light-by-light scattering diagrams in the two-loop correction. The total shift of the transition frequency is +6.4 kHz, of which +5.4 kHz is due to the shifts in the fundamental constants (+2.7 kHz for ν nr , and +2.7 kHz for the leading-order finite-size correction in ν α 2 ), and +1.0 kHz comes from the new contributions discussed in Secs. II and III.

The theoretical uncertainty is dominated by the one-loop (Eq. ( 12)) and two-loop (term proportional to G 2loop (1S) in Eq. ( 15)) higher-order remainders (see discussion in [START_REF] Korobov | Fundamental transitions and ionization energies of the hydrogen molecular ions with few ppt uncertainty[END_REF]). In the uncertainty from fundamental constants, the largest contribution by far is that of the proton-to-electron mass ratio µ p (1.7 kHz and 1.1 kHz using 2014 and 2018 CODATA values, respectively).

Theoretical frequencies for the rovibrational transitions measured in recent experiments are presented in Table III. They are in very good agreement with experimental results in all cases. For the (0, 0) → (1, 1) transition the combined uncertainty u = 0.86 kHz (u exp = 0.16 kHz and u theor,spin = 0.85 kHz) is given, for the v = (3, 0) → (3, 9) transition, the revised experimental value from [START_REF] Germann | Three-body QED test and fifth-force constraint from vibrations and rotations of HD +[END_REF] is used. Numerical results of calculations for all the contributions considered in Sec. III (with the mα 6 -order relativistic correction as well) for a wide range of rovibrational states are given in the Supplemental Material [34]. Using these improved theoretical predictions, we give in Table IV some updated determinations of the protonto-electron mass ratio from HD + spectroscopy. We follow the least-squares fitting procedure used in the CODATA adjustments and described in Appendix E of [START_REF] Mohr | CODATA recommended values of the fundamental physical constants[END_REF]. The dependence of HD + transition frequencies on fundamental constants is linearized using a first-order Taylor expansion around their starting (CODATA-2018) values; first-derivative coefficients are obtained by the methods described in [START_REF] Schiller | Tests of time independence of the electron and nuclear masses with ultracold molecules[END_REF][START_REF] Karr | Hydrogen molecular ions for improved determination of fundamental constants[END_REF]. Following the CODATA approach, we adjust the individual particle masses rather than mass ratios. In more detail, we use CODATA-2018 values of m d , R ∞ , r d , r p , include as an additional data point the latest measurement of the electron mass [START_REF] Köhler | The electron mass from g-factor measurements on hydrogen-like carbon 12 C 5+[END_REF], and solve for the electron and proton masses. m p /m e and its uncertainty are then deduced from the adjusted values of m p and m e , taking the correlation between them into account. The covariance matrix of the HD + input data (see [START_REF] Mohr | CODATA recommended values of the fundamental physical constants[END_REF] for definition) is built including the three (uncorrelated) sources of uncertainties: experimental, theoretical, and parametric, the latter stemming from uncertainties of fundamental constants that are not directly involved in the adjustment (m d , R ∞ , r d , r p ). When several HD + measurements are combined, the following assumptions are made regarding correlations: (i) uncorrelated experimental uncertainties, (ii) fully correlated theoretical uncertainties, and (iii) correlations between parametric uncertainties are included taking into account the correlation coefficients between fundamental constants available from [39].

Values obtained in this way from each single HD + experiment are given in the first three lines of Table IV, and the combined result from the three measurements in the fourth line. Note that its uncertainty is only slightly reduced with respect to that of individual lines due to strong correlation between them. The latter value, although slightly higher (by 1.4 combined standard deviations), is in good agreement with that obtained from recent highprecision mass spectroscopy measurements of m p [START_REF] Heiße | High-precision mass spectrometer for light ions[END_REF], m d [START_REF] Rau | Penning trap mass measurements of the deuteron and the HD + molecular ion[END_REF], m(HD + ) [START_REF] Rau | Penning trap mass measurements of the deuteron and the HD + molecular ion[END_REF], and m d /m p [START_REF] Fink | Deuteron-to-Proton Mass Ratio from the Cyclotron Frequency Ratio of H + 2 to D + with H + 2 in a Resolved Vibrational State[END_REF] in Penning traps: m p /m e = 1836.152 673 343(60) [START_REF] Köhler | The electron mass from g-factor measurements on hydrogen-like carbon 12 C 5+[END_REF][START_REF] Rau | Penning trap mass measurements of the deuteron and the HD + molecular ion[END_REF].

Finally, the HD + data can be combined with the mass spectrometry measurements and the CODATA-2018 values of R ∞ , r In conclusion, we have presented a revised and improved theory of spin-averaged transition frequencies in hydrogen molecular ions. Further progress in precision now requires calculations of nonlogarithmic mα 8 -order one-and two-loop corrections.

  d , r p to simultaneously determine the three particle masses, m e , m p and m d . The mass ratios m p /m e and m d /m p deduced from this adjustment are shown in the last line of Table IV. Relative uncertainties of 1.8 × 10 -11 (for m p /m e ) and 1.6 × 10 -11 (for m d /m p ) are obtained, improved by factors of 3.3 and 3.5 with respect to CODATA-2018.

TABLE I

 I

	: Reevaluation of the fundamental constants of atomic physics by CODATA in 2018. CODATA14 (resp. CODATA18)
	values are given in the upper (resp. lower) line.

TABLE II :

 II Fundamental transition frequency ν01 for the HD + molecular ion (in kHz). CODATA14 recommended values of fundamental constants were used in[START_REF] Korobov | Fundamental transitions and ionization energies of the hydrogen molecular ions with few ppt uncertainty[END_REF], and the latest CODATA18 values are used in the present work. Nuclear size and polarizability corrections are included in ν α 2 , and "other" corrections correspond to the muonic and hadronic vacuum polarization. Theoretical uncertainties of contributions at each order in α, if not negligible, are indicated within parentheses. In the final value of the transition frequency, the first error is the theoretical uncertainty, and the second one is due to the uncertainty of fundamental constants.

		[4]	this work
	νnr	57 349 439 952.4	57 349 439 955.1
	ν α 2	958 151.7	958 154.6
	ν α 3	-242 126.3	-242 126.3
	ν α 4	-1708.9(1)	-1708.9(1)
	ν α 5	106.4(1)	105.9(1)
	ν α 6	-2.0(5)	-0.8(5)
	other		0.25
	νtot	57 350 154 373.4(0.5)(1.8)	57 350 154 379.8(0.5)(1.3)

TABLE III :

 III Theoretical and experimental spin-averaged transition frequencies (in kHz). CODATA18 values of fundamental constants were used in the calculations. For theoretical values, the first uncertainty is due to yet uncalculated terms and used approximations in theory, while the second uncertainty is due to inaccuracy in the CODATA18 recommended mass values.

	(L, v) → (L , v )	theory	experiment
	(0, 0) → (1, 0)	1 314 925 752.932(19)(61)	1 314 925 752.910(17)
	(0, 0) → (1, 1)	58 605 052 163.9(0.5)(1.3)	58 605 052 164.24(86)
	(3, 0) → (3, 9)	415 264 925 502.8(3.3)(6.7)	415 264 925 501.8(1.3)

TABLE IV :

 IV Determinations of mass ratios using HD + spectroscopy.

	data	mp/me	m d /mp
	(0, 0) → (1, 0)	1836.152 673 480(63)
	(0, 0) → (1, 1)	1836.152 673 40(10)
	(3, 0) → (3, 9)	1836.152 673 457(73)
	HD +	1846.152 673 466(61)
	HD + /Penning	1836.152 673 454(33)	1.999 007 501 243(31)

Acknowledgements. The authors thank J.C.J. Koelemeij for his help with the least-squares adjustments of fundamental constants. V.I.K. acknowledges support of the Russian Foundation for Basic Research under Grant No. 19-02-00058-a.


Nam Hoang

Abdel Lisser Nguyen

Vikram Vikas

Singh

Nguyen

Abdel Lisser
email: abdel.lisser@centralesupelec.fr

Vikas Vikram Singh
email: vikassingh@maths.iitd.ac.in

Random games under elliptically distributed dependent joint chance constraints

Keywords: Chance-constrained game, Elliptical distribution, Nash equilibrium, Archimedean copulas. Mathematics Subject Classification (2000) MSC 90C15, 90C25, 90C59

We study an n-player game with random payoffs and continuous strategy sets. The payoff function of each player is defined by its expected value and the strategy set of each player is defined by a joint chance constraint. The random constraint vectors defining the joint chance constraint are dependent and follow elliptically symmetric distributions.

Introduction

The publication of the seminal book Researches into the Mathematical Principles of the Theory of Wealth by Cournot in 1838 was the trigger for the widespread use of the equilibrium under market conditions [START_REF] Cournot | Researches into the Mathematical Principles of the Theory of Wealth[END_REF]. Later, the saddle point equilibrium for a two player zero-sum game of von Neumann [START_REF] Neumann | Theory of games and economic behavior[END_REF] follows on this work. In 1950, Nash [START_REF] Nash | Non-cooperative games[END_REF] showed that for a finite strategic game there exists an equilibrium point, know as Nash equilibrium, from which there is no incentive for any player to deviate unilaterally. Despite its practical limitation, the general strategic games are extensively studied in the literature [START_REF] Başar | Dynamic Noncooperative Game Theory[END_REF][START_REF] Debreu | A social equilibrium existence theorem[END_REF][START_REF] Fan | Applications of a theorem concerning sets with convex sections[END_REF]. The theory of Nash equilibrium in deterministic setup faces challenges especially when it comes to deal with real applications with random payoffs and strategy sets. The most commonly used tool to deal with random payoffs is the expectation function [START_REF] Ravat | On the characterization of solution sets of smooth and nonsmooth convex stochastic Nash games[END_REF] which is more appropriate for risk neutral cases. The risk averse games are studied by considering an alternative payoff criterion based on risk measure CVaR [START_REF] Kannan | Addressing supply-side risk in uncertain power markets: stochastic Nash models, scalable algorithms and error analysis[END_REF][START_REF] Ravat | On the characterization of solution sets of smooth and nonsmooth convex stochastic Nash games[END_REF] and chance constraint programming [START_REF] Singh | Existence of Nash equilibrium for chanceconstrained games[END_REF][START_REF] Singh | A characterization of Nash equilibrium for the games with random payoffs[END_REF]. In [START_REF] Singh | Existence of Nash equilibrium for chanceconstrained games[END_REF], the authors studied a finite strategic game where the payoff vector of each player is elliptically distributed, and showed the existence of a Nash equilibrium. The equivalence between the Nash equilibrium of chance-constrained games (CCGs for short) in [START_REF] Singh | Existence of Nash equilibrium for chanceconstrained games[END_REF] and the global optimal solutions of a certain mathematical program is stated in [START_REF] Singh | A characterization of Nash equilibrium for the games with random payoffs[END_REF].

In the above referred games, the players payoff functions are random while the strategy sets are deterministic. However, the strategy sets containing chance constraints are often considered in various applications, e.g., risk constraints in portfolio optimization problem [START_REF] Ji | Risk-budgeting multi-portfolio optimization with portfolio and marginal risk constraints[END_REF] and resource constraints in stochastic shortest path problem [START_REF] Cheng | A second-order cone programming approach for linear programs with joint probabilistic constraints[END_REF]. Recently, the games with chance constraint based strategy sets are introduced in the literature [START_REF] Peng | Games with distributionally robust joint chance constraints[END_REF][START_REF] Peng | General sum games with joint chance constraints[END_REF][START_REF] Peng | Chance-constrained games with mixture distributions[END_REF][START_REF] Singh | A second-order cone programming formulation for two player zero-sum games with chance constraints[END_REF][START_REF] Singh | An equivalent mathematical program for games with random constraints[END_REF]. Singh and Lisser [START_REF] Singh | A second-order cone programming formulation for two player zero-sum games with chance constraints[END_REF] considered a 2-player zero-sum game with individual chance constraints and showed that a saddle point equilibrium problem is equivalent to a primaldual pair of second order cone programs when the random constraint vectors follow elliptically symmetric distribution. Singh et al. [START_REF] Singh | An equivalent mathematical program for games with random constraints[END_REF] considered an nplayer general-sum game with individual chance constraints under elliptically symmetric distributions and showed that a Nash equilibrium problem is equivalent to the global optimization of a nonlinear optimization problem. In the wake of these results, Peng et al. [START_REF] Peng | General sum games with joint chance constraints[END_REF] showed the existence of Nash equilibrium for the n-player general-sum games where the strategy profile set of each player is defined by a joint chance constraint, and the random constraint vectors are either independently normally distributed or follow a mixture of elliptical dis-tributions [START_REF] Peng | Chance-constrained games with mixture distributions[END_REF]. When the probability distributions are not completely known and belong to a given distributional uncertainty set, Peng et al. [START_REF] Peng | Games with distributionally robust joint chance constraints[END_REF] formulated the chance constraints of each player as distributional robust joint chance constrained problem. They consider several uncertainty sets, namely a density based uncertainty set and four two-moments based uncertainty sets where one of them has a nonnegative support. They show that there exists a Nash equilibrium of a distributionally robust joint chance constrained game for each uncertainty set.

In [START_REF] Peng | Games with distributionally robust joint chance constraints[END_REF][START_REF] Peng | General sum games with joint chance constraints[END_REF][START_REF] Peng | Chance-constrained games with mixture distributions[END_REF], the authors assume that the random constraint vectors are independently distributed. However, the random variables are usually dependent in real world applications. In order to study the dependence structure of random variable, the concept of copula was introduced by Abe Sklar in 1959 [START_REF] Sklar | Fonctions de repartition an dimensions et leurs marges[END_REF], as a solution to a probability problem stated by Maurice Fréchet in the context of random metric spaces. Copulas are functions used to distinguish the marginal distributions from a given dependent structure based multivariate distribution. Henrion and Strukgarek [START_REF] Henrion | Convexity of chance constraints with independent random variables[END_REF] introduced the notion of log-exp concavity of copula to investigate the convexity of elliptically dependent distributed joint chance constraints. We refer the reader to [START_REF] Nelson | An introduction to copulas[END_REF] for a detailed introduction to the theory of copulas.

In this paper, we extend the results of [START_REF] Peng | General sum games with joint chance constraints[END_REF][START_REF] Peng | Chance-constrained games with mixture distributions[END_REF] to the general case where the payoff function is random and the strategy profile set of each player is defined by elliptically distributed dependent joint chance constraints. We derive a new reformulation of joint chance constraint with dependent random constraint vectors and show that there exists a Nash equilibrium of the game under mild conditions on the payoff functions.

The rest of the paper is organized as follows. Section 2 contains the definition of an n-player CCG. In Section 3, we prove the existence of a Nash equilibrium of the CCG under elliptical distributions. We conclude the paper in Section 4.

The model and preliminary results

Chance-constrained game

We consider an n-player CCG, where H = {1, 2, .., n} is the set of players.

Let S i ⊂ R di be the strategy set of player i which is a non-empty, convex and compact set. For each i ∈ H, S -i denotes the set of strategy vectors of all players j, j = i. A strategy profile x = (x 1 , x 2 , . . . , x n ) ∈ S is represented as (x i , x -i ) where x i denotes the strategy of player i and x -i denotes the vector of strategies of the players other than player i. The strategy set of player i, i ∈ H, is further restricted by the following joint chance constraint

P(V i x i ≤ D i ) ≥ α i , (1)
where

α i ∈ [0, 1] is a given probability level, D i = (D i,1 , ..., D i,Ki ) T ∈ R Ki is a deterministic vector and V i = [V i,1 , ..., V i,Ki ] T is a K i × d i random matrix,
where V i,k denotes the kth row of matrix V i ; T denotes the transposition.

Let J i = {1, 2, . . . , K i } denotes the index set of ith player's constraints. The feasible strategy set of player i is defined as

S i αi = x i ∈ S i | P(V i x i ≤ D i ) ≥ α i .
We assume that for each i ∈ H, S i αi is a non-empty set. Let α = (α i ) i∈H be the confidence level vector and S α = n i=1 S i αi be the set of all feasible strategy profiles. The payoff function of each player is defined using random variables. For each x ∈ S α , the payoff of player i is given by f i (x, ζ) where ζ is an m-dimensional random vector. We use expected value approach to model the payoff function of each player. Therefore, the payoff function of player i is given by

p i (x) = E[f i (x, ζ)], ∀ x ∈ S α . (2)
We assume that the CCG is of complete information, i.e., the payoff function, the strategy set of each player, and the confidence level vector α are known to all the players. Definition 2.1 A strategy profile y * is a Nash equilibrium of the CCG at confidence level vector α if for each i ∈ H

p i (y i * , y -i * ) ≥ p i (x i , y -i * ), ∀ x i ∈ S i αi .
Assumption 1. For each player i, i ∈ H, the following conditions hold.

1.

f i (., x -i , ζ) is a concave function of x i for every (x -i , ζ) ∈ S -i × R m . 2. f i (•) is a continuous function. 3. p i (x) is finite valued for every x ∈ S.
Under Assumption 1, the payoff function of player i defined by ( 2) is a continuous function of x and it is a concave function of x i for every x -i [20].

Basic concepts and known results

In this section, we present some basic definitions and important results which are used in our subsequent analysis.

Definition 2.2 An n-dimensional random vector X follows a spherical distribution if there exists a function Ψ : R → R such that the characteristic function φ X (t) of X is given by

φ X (t) = E(e it T X ) = Ψ (t T t).
The function Ψ is called a characteristic generator of the spherical distribution. Not all elliptical distributions have a probability density function. Whenever it exists, it has the form

f U (z) = c det(Σ) g s (z -µ) T Σ -1 (z -µ) ,
where g s is a nonnegative function called radial density and c > 0 is a normalization factor which makes f U a probability density function.

Definition 2.4 (Definition 2.2, [START_REF] Henrion | Convexity of chance constraints with independent random variables[END_REF]) A real function f : R → R is rdecreasing for some real number r ∈ R, if f is continuous on (0, +∞) and there exists some strictly positive real number t * such that the function t → t r f (t)

is strictly decreasing on (t * , +∞). of an 1-dimensional real-valued random variable X has r-decreasing density if the probability density function of X is r-decreasing for some real number

r ∈ R.
Table 1 presents some 1-dimensional spherical distributions with r-decreasing densities for some values of r and their thresholds t * .

Table 1 List of selected 1-dimensional spherical distributions with r-decreasing density and

their thresholds t * . Distribution Radial density r t * Normal e -1 2 u 2 r > 0 √ r t 1 + 1 ν u 2 -(1+ν)/2 , ν > 0, ν integer 0 < r < ν + 1 rν ν+1-r Laplace e -|u| r > 0 r √ 2 Kotz type u 2(N -1) e -qu 2s , q, s > 0, N > 1 2 r > 2(1 -N ) 2s 2(N -1)

+r 2qs

Pearson type VII

1 + u 2 m -N , m > 0, N > 1 2 0 < r < 2N rm 2N -r Definition 2.6 A function C : [0, 1] K → [0, 1] is a K-dimensional copula if
C is a joint CDF of a K-dimensional random vector, on the unit cube [0, 1] K , whose marginals are uniformly distributed on [0, 1].

Proposition 2.1 (Sklar's Theorem) Given F : R K → [0, 1] is a joint CDF of a K-dimensional random vector and F 1 , ..., F K are the marginal CDFs, respectively. Then, there exists a K-dimensional copula C such that

F (z) = C (F 1 (z 1 ), ..., F K (z K )) .
Moreover, if F i is continuous for any i = 1, ..., K, then C is uniquely given by

C(u) = F F (-1) 1 (u 1 ), ..., F (-1) K (u K ) .
Proposition 2.2 (Fréchet-Hoeffding upper bound) For any K-dimensional

copula C and K-dimensional vector u = (u 1 , ..., u K ) ∈ [0, 1] K , we have C(u) ≤ min k=1,...,K u k .
Definition 2.7 A K-dimensional copula C is strictly Archimedean if there exists a continuous and strictly decreasing function ψ : (0, 1] → [0, +∞), such that ψ(1) = 0, lim t→0 ψ(t) = +∞, and for any K-dimensional vector

u = (u 1 , ..., u K ) ∈ [0, 1] K , we have C(u) = ψ (-1) K i=1 ψ(u i ) .
The function ψ is called a generator of the copula C.

Table 2 presents a selection of some strictly Archimedean copulas with their generators.

Independent - -log(t) Gumbel-Hougaard θ ≥ 1 [-log(t)] θ Frank θ > 0 -log e -θt -1 e -θ -1 Clayton θ > 0 1 θ (t θ -1) Joe θ ≥ 1 -log[1 -(1 -t) θ ] Definition 2.8 A function f : R → R is K-monotonic on an open interval
I ⊆ R for some positive integer K ≥ 2, if the following three conditions hold:

1. f is differentiable up to the order (K -2) on I,

2. The derivatives of f are satisfied by

(-1) k d k dt k f (t) ≥ 0, 0 ≤ k ≤ K -2,
for all t ∈ I,

3. The function t → (-1) K-2 d K-2 dt K-2 f (t)
is nonincreasing and convex on I.

Proposition 2.3 (Theorem 2.2, [START_REF] Mcneil | Multivariate Archimedean copulas, d-monotone functions and l1-norm symmetric distributions[END_REF]) Given ψ : (0, 1] → [0, +∞) is a strictly decreasing function such that ψ(1) = 0 and lim t→0 ψ(t) = +∞. Then, ψ is the generator of a K-dimensional strictly Archimedean copula if and only if the inverse function ψ (-1) is K-monotonic on (0, +∞) and continuous on

[0, +∞). Definition 2.9 A function f : Q → (0, +∞) is r-concave on a convex set Q ⊂ R s for a given r ∈ (-∞, +∞) if for any x, y ∈ Q and α ∈ [0, 1], f (αx + (1 -α)y) ≥ [αf (x) r + (1 -α)f (y) r ] 1 r , when r = 0, f (αx + (1 -α)y) ≥ f (x) α f (y) 1-α , otherwise.

Existence of Nash equilibrium

For each i ∈ H, we assume that the random constraint vector

V i,k ∼ Ellip(µ i,k , Σ i,k , Ψ i,k ), k ∈ J i . Let λ i,k
,min be the smallest eigenvalue of the positive definite matrix

Σ i,k . Define, Si αi = S i αi \ {0}, then for x i ∈ Si αi , let ξ i,k (x i ) = (V i,k ) T x i -(µ i,k ) T x i (x i ) T Σ i,k x i . g i,k (x i ) = D i,k -(µ i,k ) T x i (x i ) T Σ i,k x i . (3)
It is well known that ξ i,k (x i ) follows 1-dimensional spherical distribution with characteristic generator Ψ i,k [START_REF] Fang | Symmetric multivariate and related distributions[END_REF]. Using the above mentioned notations, the constraint (1) can be written as

P ξ i,k (x i ) ≤ g i,k (x i ), k ∈ J i ≥ α i . (4)
By Proposition 2.1, (4) can be equivalently written as

C i x i (F i,1 • g i,1 )(x i ), ..., (F i,Ki • g i,Ki )(x i ) ≥ α i , (5)
where C i x i is the K i -dimensional copula of the random vector (ξ i,k (x i )) Ki k=1 and F i,k is the cumulative distribution function of ξ i,k (x i ); • denotes the function composition.

Assumption 2. There exists a K i -dimensional copula C i such that C i x i = C i for all x i ∈ S i , and C i is a K i -dimensional strictly Archimedean copula with a generator ψ i such that the inverse function ψ (-1) i is 4-monotonic on (0, +∞).

K i ≥ 4.
Under Assumption 2, we can equivalently write (5) as

i [(F i,1 • g i,1 )(x i ), ..., (F i,Ki • g i,Ki )(x i )] ≥ α i . (6)
Proposition 3.1 If x i ∈ Si αi and Assumption 2 holds, the joint chance constraint (1) is equivalent to

         (i) (F i,k • g i,k )(x i ) ≥ ψ (-1) i (y i,k ψ i (α i )), k ∈ J i .
(ii)

k∈J i y i,k = 1, y i,k ≥ 0, k ∈ J i . (7)
Proof. Let x i ∈ Si αi . Under Assumption 2, the joint chance constraint (1) is equivalent to [START_REF] Debreu | A social equilibrium existence theorem[END_REF]. It is enough to show the equivalence between ( 6) and [START_REF] Fan | Applications of a theorem concerning sets with convex sections[END_REF].

Since, C i is strictly Archimedean copula, (6) is equivalent to

k∈J i (ψ i • F i,k • g i,k )(x i ) ≤ ψ i (α i ). (8)
Define K i -dimensional vector y i = (y i,1 , ..., y i,Ki ) ∈ [0, 1] Ki such that

y i,k = (ψ i • F i,k • g i,k )(x i ) j∈J i (ψ i • F i,j • g i,j )(x i ) , k ∈ J i .
From non-increasing property of ψ -1 i , it follows that (x i , y i ) satisfies [START_REF] Fan | Applications of a theorem concerning sets with convex sections[END_REF]. Conversely, we assume (x i , y i ) satisfies [START_REF] Fan | Applications of a theorem concerning sets with convex sections[END_REF]. By adding all the constraints (i) of ( 7)

after applying ψ i (•) on both sides, we can say that x i satisfies (8) which is equivalent to [START_REF] Debreu | A social equilibrium existence theorem[END_REF].

The convexity of the feasible strategy set S i αi plays a very important role in showing the existence of Nash equilibrium. We show that there exists

α * i ∈ [0, 1] such that S i αi is a convex set for all α i ∈ (α * i , 1]. For each i ∈ H, define an index set I (i) = {k ∈ J i | µ i,k = 0} and a set of real numbers {r i,k | k ∈ J i } such that r i,k > 1, if k ∈ I (i) , r i,k = 1, if k / ∈ I (i) .        (9)
Lemma 3.1 Let Assumption 2 holds and for each k ∈ J i , the CDF F i,k has (r i,k + 1)-decreasing density with a threshold t * i,k , where r i,k is defined by [START_REF] Henrion | Convexity of chance constraints with independent random variables[END_REF] and t * i,k refers to Definition 2.4. Then, S i αi is a convex set for all α i ∈ (α * i , 1],

where

α * i := max 1 2 , max k∈I (i) F i,k r i,k + 1 r i,k -1 λ -1 2 i,k,min ||µ i,k || , max k∈J i F i,k (t * i,k ) . ( 10
)
In order to prove Lemma 3.1, we need the three following lemmas.

Lemma 3.2 Let α i ∈ ( 1 2 , 1] and x i ∈ Si αi . Then, D i,k > (µ i,k ) T x i for all k ∈ J i .
Proof. The proof is given in 4.

Lemma 3.3 Let r i,1 , ..., r i,Ki be the real numbers defined by (9) and for each k ∈ I (i) , define

Ω i,k := x i ∈ S i | D i,k -µ T i,k x i > r i,k + 1 r i,k -1 λ -1 2 i,k,min ||µ i,k || (x i ) T Σ i,k x i .
Then,

Conv( Si αi ) ⊂ k∈I (i) Ω i,k , for all α i ∈ (α * * i , 1],
where

α * * i = max 1 2 , max k∈I (i) F i,k r i,k + 1 r i,k -1 λ -1 2 i,k,min ||µ i,k || , ( 11
)
and Conv represents the convex hull. Moreover, for any convex subset Q i,k of

k∈I (i) Ω i,k such that 0 / ∈ Q i,k , g i,k (x i ) is defined and (-r i,k )-concave on Q i,k for all k ∈ J i
Proof. The proof is given in 4.

Lemma 3.4 Let Assumption 2 holds. Then, ψ

(-1) i (y i,k ψ i (α i )) is a convex function of y i,k for all α i ∈ [0, 1].
Proof. The proof is given in 4.

We present the proof of Lemma 3.1 using the results of Lemma 3.2, Lemma 3.3 and Lemma 3.4.

Proof of Lemma 3.1. Let α i ∈ (α * i , 1], λ ∈ [0, 1] and z 1 , z 2 ∈ S i αi . We need to show that λz 1 + (1 -λ)z 2 ∈ S i αi .
Case 1: Let z 1 = 0 or z 2 = 0. Without loss of generality, we assume that z 2 = 0. This gives D i,k ≥ 0 for all k ∈ J i , which in turn implies that

P(V i λz 1 ≤ D i ) ≥ P(V i z 1 ≤ D i ) ≥ α i . Hence, λz 1 + (1 -λ)z 2 ∈ S i αi .
Case 2: Let z 1 = 0, z 2 = 0 and λz 1 + (1 -λ)z 2 = 0. In this case,

z 2 = -λ 1-λ z 1 ∈ Si αi and z 1 ∈ Si αi . It follows from Lemma 3.2 that (µ i,k ) T z 1 > λ -1 λ D i,k , (µ i,k ) T z 1 < D i,k , ∀ k ∈ J i .
This implies that D i,k ≥ 0 for all k ∈ J i . Therefore,

λz 1 + (1 -λ)z 2 = 0 ∈ S i αi .
Case 3: Let z 1 = 0, z 2 = 0 and 0 ∈ Seg(z 1 , z 2 ), where Seg(z 1 , z 2 ) =

{z 1 + l(z 2 -z 1 ), 0 ≤ l ≤ 1}.
Then, the points on the line segment Seg(z 1 , z 2 )

are either belong to Seg(z 1 , 0) or Seg(0, z 2 ). It follows from Case 1 that Seg(z 1 , 0) and Seg(0, z 2 ) are subset of S i αi . Therefore,

λz 1 + (1 -λ)z 2 ∈ S i αi for all λ ∈ [0, 1]. Case 4: Let z 1 = 0, z 2 = 0 such that 0 / ∈ Seg(z 1 , z 2 ). It is clear that Seg(z 1 , z 2 ) ⊂ Conv( Si αi )
. From Lemma 3.3, g i,k (•) is defined and (-r i,k )concave on Seg(z 1 , z 2 ). Therefore,

g i,k (λz 1 + (1 -λ)z 2 ) ≥ λ(g i,k (z 1 )) -r i,k + (1 -λ)(g i,k (z 2 )) -r i,k -1 r i,k . (12)
Since, z 1 ∈ Si αi , from Lemma 3.2 g i,k (z 1 ) > 0 and it follows from (6) that

C i [(F i,1 • g i,1 )(z 1 ), ..., (F i,Ki • g i,Ki )(z 1 )] > α * i . (13)
By using Proposition 2.2 and the definition of α * i from (10), we get

F i,k (g i,k (z 1 )) > α * i ≥ F i,k (t * i,k ). (14)
This implies that

0 < g i,k (z 1 ) -r i,k < (t * i,k ) -r i,k . Similarly, 0 < g i,k (z 2 ) -r i,k < (t * i,k ) -r i,k .
By applying the non-decreasing function F i,k (•) on both side of ( 12), we can write

(F i,k • g i,k )(λz 1 + (1 -λ)z 2 ) ≥ F i,k λ(g i,k (z 1 )) -r i,k + (1 -λ)(g i,k (z 2 )) -r i,k -1 r i,k . (15)
Since, F i,k (•) has (r i,k + 1)-decreasing density, from Lemma 3.1 of [START_REF] Henrion | Convexity of chance constraints with independent random variables[END_REF], the

function t → F i,k t -1 r i,k is concave on (0, (t * i,k ) -r i,k
). Therefore, we can write

F i,k λ(g i,k (z 1 )) -r i,k + (1 -λ)(g i,k (z 2 )) -r i,k -1 r i,k ≥ λ (F i,k • g i,k )(z 1 ) + (1 -λ) (F i,k • g i,k )(z 2 ) . (16)
From ( 15) and ( 16), we have

(F i,k • g i,k )(λz 1 + (1 -λ)z 2 ) ≥ λ (F i,k • g i,k )(z 1 ) + (1 -λ) (F i,k • g i,k )(z 2 ) . ( 17
)
This implies that the composition function ) k∈J i is also a feasible point of [START_REF] Fan | Applications of a theorem concerning sets with convex sections[END_REF]. Again from the equivalence of Si αi and ( 7),

(F i,k • g i,k )(•)
λz 1 +(1-λ)z 2 ∈ Si αi which in turn implies that λz 1 +(1-λ)z 2 ∈ S i αi .
Next, we prove that S i αi is a closed set.

Lemma 3.5 The probability function x i → P(V i x i ≤ D i ) used in the joint chance constraint (1) of player i is continuous on Si , where Si = S i \ {0}.

Proof. We can write P(V i x i ≤ D i ) as follows

P(V i x i ≤ D i ) = E I {V i x i ≤D i } = E k∈J i I {(V i,k ) T x i ≤D i,k } ,
where I A denotes the indicator function of an event A. Given x i ∈ Si and a sequence x i j ∈ Si such that x i j → x i when j → +∞. For each k ∈ J i ,

let A i,k = ω | (V i,k (ω)) T x i > D i,k and B i,k = ω | (V i,k (ω)) T x i < D i,k .
For ω ∈ A i,k , we have I {(V i,k ) T x i ≤D i,k } (ω) = 0. Since, x i j → x i , there exists a positive integer N (ω) such that for all j > N (ω), we have

(V i,k (ω)) T x i j > D i,k .
In other words, I {(V i,k ) T x i j ≤D i,k } (ω) = 0, for all j > N (ω). Hence,

I {(V i,k ) T x i j ≤D i,k } → I {(V i,k ) T x i ≤D i,k } on A i,k , j → +∞.
Similarly,

I {(V i,k ) T x i j ≤D i,k } → I {(V i,k ) T x i ≤D i,k } on B i,k , j → +∞. Define C i,k = ω | (V i,k (ω)) T x i = D i,k . Using the notations from (3), C i,k = ω | ξ i,k (x i )(ω) = g i,k (x i )
. Since, ξ i,k (x i ) follows an 1-dimensional real continuous distribution with a density function, P(ξ i,k (x i ) = c) = 0, for any c ∈ R.

In other words

, C i,k is a negligible set. Note that A i,k ∪ B i,k ∪ C i,k = Ω,
where Ω is the sample space. Therefore, for each k

∈ J i , I {(V i,k ) T x i j ≤D i,k } → I {(V i,k ) T x i ≤D i,k
} almost everywhere which in turn implies that I {V i x i j ≤D i } → I {V i x i ≤D i } almost everywhere. Moreover, I {V i x i j ≤D i } is upper bounded by a positive integrable function I Ω . Then, it follows from dominated convergence theorem that E(I {V i x i j ≤D i } ) → E(I {V i x i ≤D i } ) as j → +∞. Therefore,

P(V i x i ≤ D i ) is a continuous function on Si .
Lemma 3.6 The feasible strategy set S i αi of player i is a closed set.

Proof. Given x i ∈ S i and a sequence x i j ∈ S i αi such that x i j → x i when j → +∞. If x i = 0, the proof follows from Lemma 3.5. If x i = 0, we need to prove that 0 ∈ S i αi . Let 0 / ∈ S i αi and x i j = 0 for all j ∈ N. Then, there exists k * ∈ J i such that D i,k * < 0. For each j ∈ N, we have

P(V i x i j ≤ D i ) ≤ P((V i,k * ) T x i j ≤ D i,k * ).
It follows from the Cauchy-Schwarz inequality that P((

V i,k * ) T x i j ≤ D i,k * ) ≤ P(-||V i,k * || × ||x i j || ≤ D i,k * ) = P(||V i,k * || × ||x i j || ≥ -D i,k * ). Hence, P(V i x i j ≤ D i ) ≤ P(||V i,k * || × ||x i j || ≥ -D i,k * ). = P ||V i,k * || ≥ -D i,k * ||x i j || . (18)
As -D i,k * > 0 and x i j → 0, we deduce that

P ||V i,k * || ≥ -D i,k * ||x i j || → 0. Then, from (18)
, we have

P(V i x i j ≤ D i ) → 0, j → +∞. (19)
However, as x i j ∈ S i αi , we deduce that P(V i x i j ≤ D i ) ≥ α i , for all j ∈ N which contradicts [START_REF] Peng | General sum games with joint chance constraints[END_REF]. Therefore, 0 ∈ S i αi .

The feasible strategy set S i αi is a compact set because, from Lemma 3.6, it is a closed subset of the compact set S i . Finally, we show that there exists a Nash equilibrium of the CCG. 2. For each i ∈ H and k ∈ J i , V i,k ∼ Ellip(µ i,k , Σ i,k , Ψ i,k ), where Σ i,k is a positive definite matrix.

3. For each i ∈ H and k ∈ J i , suppose CDF F i,k (•) has (r i,k + 1)-decreasing density with a threshold t * i,k , where r i,k is defined by (9) and t * i,k refers to Definition 2.4.

Then, there exists a Nash equilibrium of the CCG for any α ∈ (α * 1 , 1] × ... × (α * n , 1], where α * i , i ∈ H, is defined by [START_REF] Henrion | Convexity of chance constraints with dependent random variables: the use of copulae[END_REF].

Proof. Let α ∈ (α * 1 , 1] × ... × (α * n , 1]. Under Assumption 1, the payoff function p i (x i , x -i ) is a concave function of x i , for every x -i ∈ S -i , and a continuous function of x. It follows from Lemma 3.1 that the feasible strategy set S i αi ,

i ∈ H, is a convex set for all α i ∈ (α * i , 1]. For each i ∈ H, S i αi is a compact set.

Then, the existence of a Nash equilibrium of the CCG follows from Theorem 4 of [START_REF] Fan | Applications of a theorem concerning sets with convex sections[END_REF].

Conclusion

In this paper, we studied an n-player non-cooperative CCG where the strategy sets of each player is given by joint chance constraint with dependent random constraint vectors which follow elliptically symmetric distributions. We propose a new reformulation of the joint chance constraints based on the family of Archimedean copulas. We assume that 1-dimensional spherical distribution function in the reformulation has r-decreasing densities for some values of r.

This condition is satisfied by a list of prominent 1-dimensional spherical distributions. Under mild conditions on the payoff functions, we show that there exists a Nash equilibrium of the CCG.

Since, F i,k (•) is a non-decreasing function, from [START_REF] Rüschendorf | Sklar's theorem, and distributional transform[END_REF] we have

D i,k -(µ i,k ) T x i > r i,k + 1 r i,k -1 λ -1 2 i,k,min ||µ i,k || (x i ) T Σ i,k x i . (23)
Therefore, Si α i ⊂ k∈I (i) Ω i,k . For each k ∈ I (i) , Ω i,k is a convex set which implies that Conv( Si α i ) ⊂ k∈I i Ω i,k . We prove the second part of Lemma 3.3 by considering two cases as below:

Case 1: Let k / ∈ I (i) , then µ i,k = 0. From the definition of α * * i , we have α i > 1 2 . From Lemma 3.2, D i,k > 0. In this case, the proof follows directly from Lemma 3 of [START_REF] Cheng | Second-order cone programming approach for elliptically distributed joint probabilistic constraints with dependent rows[END_REF].

Case 2: Let k ∈ I i . It follows from Lemma 2 of [START_REF] Cheng | Second-order cone programming approach for elliptically distributed joint probabilistic constraints with dependent rows[END_REF] that the function

f i,k (x i ) = (x i ) T Σ i,k x i D i,k -(µ i,k ) T x i r i,k
.

Definition 2 . 3

 23 An n-dimensional random vector U follows an elliptical distribution with location parameter µ, positive definite scale matrix Σ and characteristic generator Ψ , i.e., U ∼ Ellip(µ, Σ, Ψ ), if we have the following representation U d = µ + AX, where X follows a spherical distribution with a characteristic generator Ψ , A ∈ R n×n such that AA T = Σ and µ ∈ R n ; d = implies that the both sides have the same distribution.

Definition 2 . 5

 25 The cumulative distribution function (CDF) F : R → [0, 1]

Remark 3 . 1 (

 31 The 4-monotonicity of ψ follows from Proposition 2.3 that 4-monotonicity condition holds if

Theorem 3 . 1

 31 Consider an n-player CCG defined in Section 2.1, where 1. Assumptions 1, 2 hold.

Table 2

 2 Different types of strictly Archimedean copulas.

	Type of copula	Parameter θ	Generator ψ θ (t)

  is a concave function over Seg(z 1 , z 2 ). It follows from Lemma 3.4 that ψ(y i,k ψ i (α i )) is a convex function of y i,k . Because z 1 , z 2 ∈ Siαi and from Proposition 3.1, Si αi and (7) are equivalent, then there exists vectors (y 1 i,k ) k∈J i and (y 2 i,k ) k∈J i such that z 1 , (y 1 i,k ) k∈J i and z 2 , (y 2 i,k ) k∈J i are feasible points of[START_REF] Fan | Applications of a theorem concerning sets with convex sections[END_REF]. Using the fact that (F i,k • g i,k )(•) is a concave function and ψ

	(-1) i	(•) is a convex func-
	tion, we can say that the convex combination of points z 1 , (y 1 i,k ) k∈J i and
	z 2 , (y 2 i,k
		(-1)
		i

Acknowledgements This research was supported by DST/CEFIPRA Project No. IFC/4117/DST-CNRS-5th call/2017-18/2 and CNRS Project No. AR/SB:2018-07-440.

Appendices

Appendix A: Proof of Lemma 3.2 Let x i ∈ Si α i . By applying Proposition 2.2 on (6), we get

Since, F i,k is the CDF of an 1-dimensional real-valued random variable which is symmetric at

Appendix B: Proof of Lemma 3.3

Let k ∈ I (i) and x i ∈ Si α i . By applying Proposition 2.2 on (6), we get

From the definition of α * * i given in [START_REF] Ji | Risk-budgeting multi-portfolio optimization with portfolio and marginal risk constraints[END_REF], we have is defined and a convex function on k∈I (i) Ω i,k . Therefore, for any y, z ∈ Q i,k and λ ∈ [0, 1],

we have

Note that g i,k 24), we can write

Hence, g i,k is defined and

Appendix C: Proof of Lemma 3.4

The second-order differentiation of U (y i,k ) is given by

) is a convex function of y i,k on (0, 1]. The convexity of U on [0, 1] follows from the continuity of U at 0.


domain: spi math

text:
Zain Ul Abdin

Ahmed Rachid

Bond graph modeling of a water-based photovoltaic/thermal (PV/T) collector

Keywords:

This paper deals with photovoltaic-thermal (PV/T) collectors that combine photovoltaic (PV) module and solar thermal (ST) modules to provide heat and electricity simultaneously. The use of a PV/T system not only enhances PV efficiency but also allows to use solar thermal energy for various heating applications. The performance evaluation of the PV/T collectors had been investigated both theoretically and experimentally for the past few years. In this work, we consider a new water-based PV/T configuration by incorporating a tedlar layer and parallel tubes and further investigate its modeling. The aim is to develop a non-linear dynamic model which gives a realistic and reasonable performance of the system. The fundamental equations are determined for the thermal part using a bond graph approach which is a generic and general tool to represent thermal transfers. Simulations are given to highlight the usefulness of the proposed design and its modeling, considering the influence of some external factors (e.g. wind action) and internal geometrical parameters (e.g. insulator thickness). The results indicate that the rise in wind speed can reduce thermal efficiency from 70% to approximately 40%. The results also show that with an increase in insulator thickness, there is a sudden change in outlet temperature but thicknesses of more than 0.14 m is not required as the energy losses become constant. The assessment of the electrical and thermal performance of a PV/T collector is taken into consideration for a sample weather condition of southern France. Further, our results are compared to existing experimental and numerical ones from the literature on alternative PV/Ts design to validate and demonstrate the viability and consistency of our model.

Introduction

Photovoltaic/thermal (PV/T) systems are panels which provide heat and electricity at the same time by combining photovoltaic (PV) and solar-thermal (ST) modules. These systems are a good solution to improve the efficiency of a solar collectors and are used for various purposes like heating and cooling combined to electricity generation [START_REF] Ul | A survey on applications of hybrid pv/t panels[END_REF]. In the past various types of PV/T collectors had been proposed with theoretical and experimental analysis. There are various approaches used to extract the heat from solar collector which includes water-based, airbased, nanofluids, and PCM based PV/T collectors. The priority of the PV/T systems is to generate electricity and the photovoltaic efficiency highly depends on the cell temperature. PV/T systems come up with the cooling so ultimately the cell temperature is reduced and the efficiency is increased [START_REF] Garg | Transient simulation of conventional hybrid photovoltaic/thermal (pv/t) air heating collectors[END_REF]] and by adding the solar thermal module, overall efficiency is increased [START_REF] Ji | Dynamic performance of hybrid photovoltaic/thermal collector wall in hong kong[END_REF]. Theoretically, 80% efficiency is achievable with more than 60% thermal from a PV/T collector [START_REF] Bergene | Model calculations on a flat-plate solar heat collector with integrated solar cells[END_REF]. The first work on PV/T systems was presented by [START_REF] Lw Florschuetz | Extension of the hottel-whillier model to the analysis of combined photovoltaic/thermal flat plate collectors[END_REF] where the main concept of the PV/T collector was presented by assuming a linear correlation between PV cell temperature and its efficiency. Since then, various theoretical and experimental investigations have been presented. We have recently presented a thermal model for an air-based PV/T collector [START_REF] Ul | Modeling, identification and control of photovoltaic/thermal solar panel[END_REF] using air as a working fluid. The results indicate that fluid (air) has low heat transfer due to low heat capacity and low thermal conductivity. The use of air as a working fluid is cheaper and easier but the use of water as a working fluid is more efficient for hot weather conditions. The performance evaluation of a PV/T system using water as a coolant was conducted by [START_REF] He | Hybrid photovoltaic and thermal solar-collector designed for natural circulation of water[END_REF] where it is concluded that water-based PV/T systems are suitable for the warm climatic applications. In addition, the results obtained were encouraging with a daily thermal efficiency of 40%. A comparative study was done by [START_REF] Cs Malvi | Energy balance model of combined photovoltaic solar-thermal system incorporating phase change material[END_REF], where a 1-D energy balance model is presented incorporating PCM. Using the 1-D energy balance model, the performance of the panel is evaluated using PCM. The model is divided into five layers; glass cover, PV cell layer, copper plate, PCM layer and, insulated frame. It is concluded that by increasing the PCM thickness to a certain limit the output increases by 6.5%, and 9.5% when comparing it with PV only system. Thermal performance of solar installations have been studied by [START_REF] Cristofari | Modelling and performance of a copolymer solar water heating collector[END_REF] where technical characteristics are optimized, and the influence of various parameters is presented theoretically. A number of different configurations and multiple thermal models have been reported with the influence of various parameters on the performance of the PV/T systems. [START_REF] Dubey | Thermal modeling of a combined system of photovoltaic thermal (pv/t) solar water heater[END_REF] consider two flat plate collectors which are interconnected in series and connected to a storage tank. 60% area of one collector is integrated by the PV module and remaining with glazing, the second collector is fully covered by glazing. Under this, there is an absorber, followed by metallic plates through which water flows, and finally at the bottom is an insulator layer. The result indicates that by increasing glazing, efficiency increases from 33% to 64%. [START_REF] Chaabane | Thermal performance of an integrated collector storage solar water heater (icsswh) with phase change materials (pcm)[END_REF] proposed two different numerical models and study the effect of PCM on thermal performance by integrating PCM directly in the collector. Based on control a volume finite difference approach, dynamic model with water as a working fluid is developed by [START_REF] Chow | Performance analysis of photovoltaic-thermal collector by explicit dynamic model[END_REF] and the energy balance matrix is obtained. The collector is divided into various sections and each section consists of a single water tube. The PV cell layer is in direct contact with the water tube. There is direct heat transfer from the PV plate to the tube. A performance evaluation of a water-based PV/T collector is done by [START_REF] Tiwari | Performance evaluation of solar pv/t system: an experimental validation[END_REF] and numerical results are then validated by the experimental results. The simulated results predict a daily thermal efficiency of around 58%, which is slightly less than the experimental value of 61.3%. [START_REF] Herman A Zondag | The thermal and electrical yield of a pv-thermal collector[END_REF] presented four different numerical models with three steady-state models and one dynamical model. The model was built on a single pipe winding (twisting like a snake) to measure the performance of the system. The model is full time-dependent quasi 3D and includes second-order differentiations. Absorber plate and the PV laminate are split up in two directions; parallel to the flow and perpendicular to the flow whereas the top layer is separated; parallel to the flow.

In this paper, a bond graph approach is used to model the thermal part of a water-based PV/T collector by incorporating a tedlar layer and parallel tubes. The performance of a PV/T collector is characterized by its water temperature rise, thermal and electrical efficiencies and the net available energy. At the same time, the sensitivity of geometrical parameters (e.g. insulator thickness) and influence of some external factors (e.g. wind action) is analyzed to improve design and performance of a PV/T collector. Further, the assessment of electrical and thermal performance of a PV/T collector is taken into consideration for a sample weather condition of southern France. Finally, current model is verified by comparing the results with several studies in the literature. The paper is organized as follows; a bond graph modeling technique is presented in section 2 to model the thermal part of a water-based PV/T collector. Section 3 provides the simulation results and discussion. It also provides valuable information regarding the influence of various parameters such as thickness, mass flow rate, wind speed, and solar radiation on the outlet fluid temperature. The paper ends with a conclusion and some perspectives.

Modeling

There exist several papers in which the performance of a PV/T is studied and most of the studies are concentrated on the experimental performance of the PV/T system. The focus is to present a non-linear dynamic model of a waterbased PV/T collector which gives a realistic and reasonable performance of the system. In a PV/T panel configuration, cold water is supplied from one end of the panel, and warm water is collected along with the electrical gain from the other end of the panel as shown in Fig. 1(a). Fig. 1(b) shows the thermal model of the panel, composed of several layers; a protective glass layer for PV cells, solar cells layer, protective tedlar layer under PV cells, absorber, metallic tube with bond, water channel, and thermally insulated frame to reduce heat losses. There are four parallel water tubes available through which fluid flows and these tubes are fitted and pressed with the bond to the insulator layer as shown in Fig. 1(b). Water is inserted from a single point, then it is distributed equally to all the tubes and warm water is collected from the other end of the panel as shown in Fig. 1(d). The governing equations of the thermal model were determined using a bond graph technique. Few assumptions have been considered for the thermal model such as heat transfer between the ground and the panel is over sighted and pressure losses are neglected, the temperature is assumed to be homogeneous for each layer. Radiation heat exchange at the bottom and edge losses are neglected because the temperature difference at the bottom is small.

A bond graph is a graphical sketch of a physical system, represented by symbols and lines. It also consists of elements of resistance, compliance, and inertance which are interconnected by junctions, bonds, and flow of energy or power. Effort ( ) and flow ( ) are two power conjugated variables in the bond graph and in physical domains the way of selection of these two variables is different. There can be multiple domains and this system contains two subsystems i.e. thermal and hydraulic subsystem. In a thermal subsystem, the effort variable is represented by temperature whereas flow variable as heat flow ( ̇ ), solar flow ( ) and enthalpy ( ̇ ). In a hydraulic subsystem, the pressure ( ) is an effort variable and mass flow ( ̇ ) is flow variable. Resistance ( ), compliance ( ), and inertance ( ) are classified as a single port and contain no source of power, convert the supplied energy into stored or dissipated energy, and are classified as passive elements. Sources are the active elements that provide power to the system, they are single port and they can be either flow source ( ) or effort source ( ). Junctions represent system topology, can connect two or more bonds and are of two kinds, "0-junction" and a "1-Junction". In 0-Junction, all connected bonds constrained to have the same efforts whereas flows sum up to zero. 1-junction has equality of flows and the algebraic sum of efforts is zero. A bond graph can include multiple domains, when a fluid is flowing and heat transfer is also involved, thermal and hydraulic energies are combined and are represented by a circle around a bond as shown in a bond graph.

For 0-Junction, effort equation [START_REF] Borutzky | Bond graph modelling and simulation of multidisciplinary systems-an introduction[END_REF] is respectively given as: and has flow equation

∑ =1 = ∑ =1 (2)
For 1-Junction, flow equation is respectively given as:

1 = ... = = 1 = ... = (3)
and has effort equation

∑ =1 = ∑ =1 (4)
Each layer is break up into two pieces: one is dissipated by conduction and the other is stored by the layer. 0junction relates to the compliance element whereas 1-junction relates to the resistance element. Relationship between generalized effort and generalized flow corresponding to resistance and compliance elements are written as ( 5) and ( 6):

( ) = 1 ( ), (5)
where

( ) = 1 ∫ ( ) (6)
When two objects with two different temperatures are in contact with each other, heat energy pass between them. In hot materials, molecules move faster and hit the cold object molecule's which are incontact with the hot object. Heat transfer coefficient (ℎ) is the inverse of thermal heat flow and using (5) thermal heat flow for the resistance element between two layers a and b is given by:

̇ = 1 ( -) (7)
Using ( 6), temperature ( ) of each layer with thermal heat capacity ( ) and thermal heat flow ( ̇ ) accumulated by each layer can be written as:

= 1 ∫ ̇ (8)
Thermal heat capacity with density ( ), volume ( ) and specific heat capacity ( ) for each layer is given by:

= (9)
At a rate of in Watts, solar radiation energy denoted by , absorbed by the glass cover is and the energy is absorbed by the PV cell layer and is , where and are absorptance of glass cover and PV cell whereas is the transmittance of the glass cover. There is a change in temperature and in enthalpy in the water tube during heat exchange and the energy gain [START_REF] Bergene | Model calculations on a flat-plate solar heat collector with integrated solar cells[END_REF] can be written as:

̇ = ̇ -̇ = ̇ ( - ) (10)
Glass cover is the top layer of the PV/T collector and the purpose of adding the glass cover is to protect the PV cell layer and reduce the heat loss to ambient. This will eventually improve the thermal efficiency of the system. The equation for glass layer using ( 7) and ( 8) is respectively given by:

= [ -ℎ , ( -) -ℎ , ( - ) -ℎ , ( -)], (11)
where module area is written as;

= ( 12
)
, is the radiative heat flow between sky and glass cover. Where is Stefan Boltzmann constant and radiation heat transfer coefficient between glass layer and sky [START_REF] Hegazy | Comparative study of the performances of four photovoltaic/thermal solar air collectors[END_REF] is given by:

ℎ , = ( 2 -2 )( 2 + 2 ) ( - ) (13)
To take atmosphere into account, sky temperature ( ) [START_REF] Chaabane | Thermal performance of an integrated collector storage solar water heater (icsswh) with phase change materials (pcm)[END_REF] is given by Swinbank (1963) associate's sky temperature with local air temperature. = 0.0552 1.5 (14)

The heat convection takes place at the top and bottom surfaces. Convective heat coefficient strongly depends on the wind speed ( ) and is given by a dimensional equation [Duffie and [START_REF] Beckman | Solar engineering of thermal processes[END_REF].

ℎ , = 5.7 + 3.8 (15)
, is the conduction heat flow and conductive heat transfer coefficient between two parallel layers (a and b) [START_REF] Slimani | A detailed thermal-electrical model of three photovoltaic/thermal (pv/t) hybrid air collectors and photovoltaic (pv) module: Comparative study under algiers climatic conditions[END_REF] is given as;

ℎ , = 1 + , ( 16
)
The equation for PV cell layer using ( 7) and ( 8) is respectively given by:

= [ + ℎ , ( -) -ℎ , ( -)] - (17)
Electric power ( ) produced by PV panel [START_REF] Lw Florschuetz | Extension of the hottel-whillier model to the analysis of combined photovoltaic/thermal flat plate collectors[END_REF] can be respectively written as (18).

= = [1 -( -, )] (18)
The equation for tedlar layer using ( 7) and ( 8) is respectively given by:

= ℎ , ( -) - ℎ , ( -) - ℎ , ( -), (19)
where

= (20) and = ( -) (21)
The equation for absorber using ( 7) and ( 8) is respectively given by:

= ℎ , ( -) - ℎ , ( -) - ℎ , ( -), (22)
where

= (23) and = ( -) (24)
The equation for tube using ( 7) and ( 8) is respectively given by:

= ℎ , ( -) + ℎ , ( -) - ℎ , ( -) - ℎ , ( -), (25)
where

= + 2 (26)
The inside area of the tube ( ) can be written as:

= (27)
and the area of the bond ( ) with width of bond ( )

= ( 28
)
In the heat transfer processes, convective heat coefficients add inversely in series and overall heat transfer coefficient with contact area, and with a bond thickness ( ) and thermal conductivity ( ) is:

1 ℎ , = 1 ℎ + , ( 29
)
where heat transfer coefficient for fully developed laminar flow can be written as ℎ [START_REF] Chow | Performance analysis of photovoltaic-thermal collector by explicit dynamic model[END_REF].

ℎ = 4.364 (30)
The equation for water channel using ( 7), ( 8) and ( 10) is respectively given by:

= ℎ , ( -) + ̇ ( - ) (31)
is the inlet temperature to the channel whereas is the output temperature of the water channel and the mean temperature ( ) is given as:

= ( + ) 2 (32)
The equation for insulator layer using ( 7) and ( 8) is respectively given by:

= ℎ , ( -) + ℎ , ( -) - ℎ , ( - ) (33)
, is the convective heat flow and convective heat transfer coefficient between insulator and surrounding is given as;

1 ℎ , = 1 ℎ , + 2 (34)
Thermal efficiency is defined as ratio of useful energy gain to the incident solar energy over the same period of time. Applying useful energy gain (10), the thermal performance of PV/T collector [START_REF] Nacer | Dynamic coupled electrical and thermal model for pv-t solar energy collectors[END_REF]

is; ℎ = ̇ ( - ) (35)
PV efficiency of the PV/T panel is defined as ratio of electrical gain to the total incident energy and written as;

= (36)
The overall efficiency is the sum of thermal efficiency (35) and electrical efficiency (36) of a PV/T collector and can be calculated as follows:

= + ℎ (37)
Thermal efficiency mainly depends on the solar radiation available whereas electrical efficiency highly depends on the cell temperature. A high cell temperature will decrease the photovoltaic efficiency so it is important to have not very high cell temperature.

Results and discussion

Simulation studies of a PV/T model

Performance of PV/T collector depends on several parameters: ambient temperature, solar radiation, mass flow rate, wind speed, fluid flowing (e.g. water or air) and thickness of various layers. Ambient temperature (

) and solar radiation (

) are presented under a sample of southern France climatic conditions in Fig. 3. The sample is taken during daylight because of strong solar radiation and for the rest of the day; the solar radiation is hardly at all. Both solar radiation and ambient temperature increase with time and reach a maximum value around noon. The ambient temperature reaches to a maximum value of 35 C and solar radiation goes to 870 W/m 2 , respectively, at the solar noon.

The Matlab program is used for the proposed numerical design and climatic parameters. The obtained results are provided in Fig. 456789. The results obtained from the proposed model are then compared with previous studies in the literature and are provided in Table 1 and Fig. 10. Values of absorptivity, transmissivity and emissivity for the glass and cell layers of a PV/T collector are taken as: = 0.06, = 0.8, = 0.93 and = 0.85. The temperature profile for glass cover, PV cell, tedlar, absorber, insulator, metallic tube and fluid is shown in Fig. 4. Temperatures rise with the increase in solar radiation and ambient temperature, and falls with the drop in ambient temperature and solar radiation as illustrated. Ambient temperature and solar radiation vary, depending on the weather conditions. During the day time, the maximum solar radiation available is 870 W/m 2 which occurs at 13:00 (local time) when the ambient temperature is maximum too. The temperature of PV module components is highest around this point with the PV cell temperature goes to approximately 53 C. It is obvious that the PV cell temperature is higher than the other components temperature whatever the value of solar radiation and ambient temperature. The inclussion of glass cover on top not only protects the PV module but causes an increase in PV cell temperature due to which the fluid temperature rises. However, the PV efficiency drops with rise in every 1 C cell temperature. Further, the fluid temperature rises to 49 C as shown in Fig. 4(b) with the inlet water temperature of 25 C, mass flow rate is kept constant at 0.005 kg/sec and insulaton thickness of 0.12 m. It can also be seen that each component in the PV/T collector has a higher temperature than the ambient temperature. A step-change in solar radiation is used in order to know the dynamic response and is given in Fig. 5. The solar radiation is reduced from 200 W/m 2 to 50 W/m 2 and results indicate that the response is not sudden; it takes time to reach the destination point. Moreover, it can be seen that the temperature rise is faster than the temperature drop. An analysis of Fig. 6 illustrates, increase in an insulator thickness leads to an increase in fluid temperature but there is a sudden change in outlet temperature if the thickness is less than 0.08 m. Additionally, for economical reasons, thicknesses of more than 0.14 m is not required as the energy losses becomes constant. Wind speed has a major influence on the performance of the PV/T collector especially during the summer season. This is mainly due to the fact that convective exchange is higher. The results in Fig. 7 show that an increase in wind speed, leading to degradation of thermal efficiency (35). The thermal efficiency drops from 70% to approximately 40% with the wind speed of 1 m/sec to 5 m/sec. It is also noticed that with solar radiation of 500 and 750 W/m 2 and wind speed of 1 m/sec to 5 m/sec, the thermal efficiency varies between 32% and 70%. For cost effectiveness of a PV/T collector, selection of mass flow rate is an important factor for the pump capacity. Fig. 8 illustrates the effect of mass flow rate on the outlet temperature by varying solar radiation. The results indicate that with an increase in mass flow rate, temperature decreases. An interesting observation is that at lower mass flow rates the temperature rise/drop is getting smaller. Further, an increase in solar radiation, the temperature difference between two different mass flow rates increases. The PV efficiency of a collector is respectively given in Fig. 9. Here, the results showed that with an increase in solar radiation and cell temperature (Fig. 4) the PV efficiency decreases, and by the end of the day, it increases with the decrease in solar radiation and cell temperature under the same mass flow rate. The average PV efficiency during the daytime is 10.82% with values ranging between a minimum of 10.28% and a maximum of 11.56%.

Validation and verification of PV/T model

The numerical results in this work are compared with the ones found in the literature [START_REF] Dubey | Thermal modeling of a combined system of photovoltaic thermal (pv/t) solar water heater[END_REF][START_REF] Tiwari | Performance evaluation of solar pv/t system: an experimental validation[END_REF][START_REF] Chow | Performance analysis of photovoltaic-thermal collector by explicit dynamic model[END_REF][START_REF] Herman A Zondag | The thermal and electrical yield of a pv-thermal collector[END_REF] and are summarized in Table 1.

The hybrid PV/T collector has been experimentally studied by [START_REF] Dubey | Thermal modeling of a combined system of photovoltaic thermal (pv/t) solar water heater[END_REF], equipped with a glass cover encased in an aluminum metallic box with 0.1 m glass wool insulation below the absorber to reduce bottom losses. The experimental value of cell temperature goes to 52 C whereas outlet temperature was found at 43 C with the same insulation and weather conditions. However, if the mass flow rate in the current model is increased from 0.005 kg/sec to 0.09 kg/sec the outlet temperature drops close to 43 C reported experimentally.

The numerical results obtained in the present study were also compared with other numerical and experimental results found in [START_REF] Tiwari | Performance evaluation of solar pv/t system: an experimental validation[END_REF]. The reduction in insulator thickness and solar radiation drops the outlet temperature so with insulator thickness of 0.05 m and solar radiation of 812 W/m 2 , the obtained results are in good agreement.

Furthermore, the values computed by the developed model with the same mass flow rate were compared with the values reported in [START_REF] Chow | Performance analysis of photovoltaic-thermal collector by explicit dynamic model[END_REF]. The outlet water temperature values have a small difference due to the difference in solar radiation.

Our data is also compatible with the experimental results in [START_REF] Herman A Zondag | The thermal and electrical yield of a pv-thermal collector[END_REF]. The two parameters due to which the outlet temperature difference is high are insulation thickness which is 0.12 m (resp. an ambient temperature of 34 C) in this study compared to 0.02 m (resp. 20 C in Zondag et al. [2002]). On the other hand, their results are in good agreement with ours when the solar radiation is taken as 800 W/m 2 with ambient temperature 20 C and insulator thickness of 0.02 m. It is worth it to mention that the highest PV efficiency achieved in the current work is close to the maximum value obtained experimentally in other studies in the literature referred in Table 1. The hybrid PV/T system has also been experimentally studied by [START_REF] Huang | Performance evaluation of solar photovoltaic/thermal systems[END_REF] to understand the performance of an integrated photovoltaic and thermal solar system (IPVTS) and the experimental variations of solar radiation and ambient temperature along with PV cell temperature and fluid temperature obtained is shown in Fig. 10. The experimental climatic data for a sunny day of May is shown in Fig. 10(a) and is used in our proposed model to know the hourly variation of PV cell temperature and fluid temperature as illustrated in Fig. 10(b). It should be noted that in numerical simulation, the wind speed was taken with a value of 3 m/s, insulator thickness is 0.1 m, and water mass flow rate is taken as 0.005 kg/s. The results show that the numerical simulations are fairly close to the experimental data with a small difference. This difference is mainly due to the fact that all the parametric values (e.g. thickness of the layers and material used for absorber and insulator etc), geometry, and design used in this study might be different from experimental setup built-in experimental studies. As a general conclusion, the simulation results founded by the numerical model developed in this study and the experimental studies [START_REF] Huang | Performance evaluation of solar photovoltaic/thermal systems[END_REF] validate the viability of the proposed model. It can be observed in both Table 1 and Fig. 10 that there is good compatibility between the experimental results and the simulated results obtained using the above proposed model.

Conclusion

In this paper, we have addressed a complete thermal model of a water-based PV/T collector including a tedlar layer and parallel water tubes. The modeling has been performed using a bond graph technique which can be adapted to other PV/Ts' design and configuration. Simulations have been implemented to highlight the coherence and the usefulness of the proposed model. Its consistency has been validated by comparison to other existing PV/Ts design alternatives.

A sensitivity analysis of the physical and external variables have also been performed and the following conclusions can be drawn: The outlet water temperature rise is faster than the temperature drop with a step change in solar radiation; The insulator thickness is an important design parameter, influencing the performance of a PV/T collector. The higher outlet water temperature is achieved with the insulator thickness of more than 0.08 m but for economical reasons, thicknesses of more than 0.14 m is not required as the energy losses becomes constant.

The PV/T collector shows higher thermal efficiency if the wind speed is low and the thermal efficiency can be enhanced with the increase of water specific mass flow rate. The thermal efficiency drops from 70% to less than 40% with an increase in wind speed from 1 m/sec to 5 m/sec; The thermal energy obtained from the PV/T collector can also be used for various heating applications and for cooling purpose by utilizing the electricity generated by the PV module to run a heat pump.

Further investigations are in progress to design control algorithms in order to optimize the operational performance of the PV/T.

Figure 1 :

 1 Figure 1: (a) Water-based PV/T collector. (b) cross-sectional view of the collector. (c) metallic tube and bond. (d) PV/T collector flow dynamics.

Figure 2 :

 2 Figure 2: Bond graph model of the PV/T collector.

Figure 3 :

 3 Figure 3: Daily evolution of solar radiation and ambient temperature.

Figure 4 :

 4 Figure 4: (a) Daily evolution of temperatures ( , , , , and ). (b) Daily evolution of output temperature ( ).

Figure 5 :

 5 Figure 5: Fluid temperature with a step change in solar radiation.

Figure 6 :

 6 Figure 6: Influence of insulator thickness on the output temperature ( ).

Figure 7 :

 7 Figure 7: Influence of wind speed on thermal efficiency.

Figure 8 :

 8 Figure 8: Influence of mass flow rate on the output temperature ( ).

Figure 9 :

 9 Figure 9: Hourly variation of PV efficiency and solar radiation.

Figure 10 :

 10 Figure 10: (a) Daily evolution of solar radiation and ambient temperature experimentally. (b) Daily evolution of temperatures ( and ) both experimentally (experimental studies) and theoretically (this study).

Table 1

 1 Summary of results for the current work and other studies in literature.

	Reference	( C)	̇ (kg/sec)	(W/m 2 )	( C)	( C)	(m)	(%)
	Current study	49	0.005	870	34	52	0.12	10.28
	Dubey and Tiwari [2008] 43	0.09	Spring (Delhi)		≃ 50	0.1	10.5 -11
	Tiwari and Sodha [2006] ≃ 44	0.005	812	35	≃ 50	0.05	-
	Chow, 2003	46.4	0.005	800	-	51.72	-	≃ 9
	Zondag et al. [2002]	40	0.005	800	20	-	0.02	9.7

= ... = = 1 = ... = (1)


Alberto Della Torre

Milan Sinobad

Barry Luther-Davis

Pan Ma

Stephen Madden

Sukanta Debbarma

Khu Vu

David J Moss

Arnan Mitchell

Jean-Michel Hartmann

Post Processing Dispersion Trimming for On-Chip Mid-Infrared Supercontinuum Generation

come   L'archive ouverte pluridisciplinaire

Post processing dispersion trimming for on-chip mid-infrared supercontinuum generation

Alberto Della Torre *1 , Milan Sinobad 1,2 , Barry Luther-Davis 3 , Pan Ma 3 , Stephen Madden 3 , Sukanta Debbarma 3 , Khu Vu 3 , David J. Moss 4 , Arnan Mitchell 2 , Jean-Michel Hartmann 5 , Jean-Marc Fedeli 5 , Christelle Monat 1 , Christian Grillet On-chip mid-infrared (3-20 µm) supercontinuum (SC) generation has potential applications in many different fields such as bio imaging, environmental sensors and security [START_REF] Soref | [END_REF]. Recently, the wide transparency window from 3 to 15 µm and CMOS compatibility of germanium, along with the prediction of great nonlinear properties [2][3], have attracted a growing interest toward germanium-based platforms and octave spanning SC generation up to 8.5 µm has been demonstrated in a silicon-germanium on silicon waveguide [4]. The spectral properties and the coherence of the generated SC strongly depend on the waveguide's dispersion profile and a careful design of the group velocity dispersion is required. However, the actual dispersion of the waveguide produced by fabrication is often different from the target one, as it is sensitive to fabrication inaccuracies, surface contamination and the presence of defects. Post-process tuning mechanisms are therefore of great interest to adjust or correct a posteriori the waveguide dispersion to match the target value. Several approaches have been proposed for post-process tuning of optical properties in optical fibers [5][6] and photonic crystals [7][8][9][10]. However, the restricted number of design parameters makes post-process dispersion control of waveguides trickier. Here, we show that it is possible to fine tune the dispersion profile a posteriori by adding a chalcogenide cladding layer on top of a highly nonlinear silicon-germanium on silicon waveguide, introducing a simple post processing tool to control the supercontinuum dynamics and its properties. SC generation was first achieved in a 7 cm long 3.75 x 2.7 µm 2 cross-section air clad waveguide. The waveguide, operating in TE single mode, was pumped in anomalous dispersion with ~200 fs pulses at 4 µm and 63 MHz repetition rate. The resulting SC, spanning from 2.63 up to 6.18 µm, shows asymmetric profile and uneven amplitude across the spectrum, both typical of SC generation in the anomalous dispersion regime. A 1.26 µm thick coating layer of chalcogenide Ge11.5As24Se64.5 was then deposited by thermal vapor deposition. The chalcogenide clad waveguide was pumped with pulses with the same parameters as before at 4.15 µm. A narrower and smoother spectrum (with a -30 dB bandwidth extending from 3.1 to 5.5 µm), typical of SC generation in normal dispersion regime, was obtained. The SC generation process was simulated by numerically solving the nonlinear Schrödinger equation using the split-step Fourier method, obtaining a good agreement with experiments. Finally, a numerical analysis of the impact of the chalcogenide thickness on the dispersion was performed, showing that the dispersion properties can be finely tuned with a reasonable chalcogenide thickness resolution of ≈100 nm. These results introduce the deposition of a chalcogenide cladding layer as a simple post processing tool to trim the waveguide's group velocity dispersion depending on the actual structure produced by fabrication.