domain: sdu sde

text:
A V Shvarov

J Castillo

Schott

A Yu Bychkov
email: bychkov@geol.msu.ru

P Bénézeth

O S Pokrovsky

Yu V Shvarov

A Castillo

J Schott

Experimental determination of calcite solubility and the stability of aqueous Ca-and Na-carbonate and -bicarbonate complexes at 100-160°C and 1-50 bar pCO 2 using in situ pH measurements

Keywords: calcite, pH, high temperature, pCO 2, carbonate, solubility

Experimental determination of calcite solubility and the stability of aqueous Ca-and Na-carbonate and -bicarbonate complexes at 100-160 °C and 1-50 bar pCO2 using in situ pH measurements

INTRODUCTION

The reduction of CO 2 emissions to the atmosphere is one of the main challenges of this century.

To reduce or at least slow down the atmospheric CO 2 concentration, the geologic sequestration of CO 2 seems quite promising due to the large storage capacity and geographic ubiquity of deep geological formations (e.g. [START_REF] Bachu | Aquifer disposal of CO 2 : Hydrodynamic and mineral trapping[END_REF][START_REF] Holloway | Storage of fossil fuel-derived carbon dioxide beneath the surface of the earth[END_REF][START_REF] Oelkers | Geochemical aspects of CO 2 sequestration[END_REF]Broecker 2005;[START_REF] Bénézeth | Dawsonite synthesis and reevaluation of its thermodynamic properties from solubility measurements: Implications for mineral trapping of CO 2[END_REF]Bénézeth et al., , 2009a)). Mineral trapping, which involves the incorporation of CO 2 into a solid phase, for example via the precipitation of carbonate minerals, is a safer way of CO 2 storage compared to solubility trapping in saline aquifers, for instance. Its modelling, in turn, requires knowing the thermodynamic parameters of the main carbonate minerals and especially calcite and the Ca aqueous species present in the fluids resulting from CO 2 injection. Although a large amount of work has been devoted to the calcite-water system at ambient conditions, the thermodynamics of this system, including the stability constants of the aqueous CaCO 3 ° and CaHCO 3 + complexes, remain poorly characterized at elevated temperatures and pCO 2 pressures most pertinent to CO 2 geological sequestration. Moreover, although sodium chloride is considered as the main component of deep geological fluids, the stability of sodium-carbonate and bicarbonate complexes has been mostly investigated at room temperature and their extrapolation to elevated temperatures is not straightforward. Only the recent study of [START_REF] Stefánsson | Carbonic acid ionization and the stability of sodium bicarbonate and carbonate ion pairs to 200 °C -A potentiometric and spectrophotometric study[END_REF] has determined the parameters of these complexes to 200°C.

The solubility of calcite has been thoroughly investigated. Starting from [START_REF] Engel | Sur la solubilité des sels en présence des acides, des bases et de sels[END_REF], many researchers conducted solubility experiments at both room temperatures (e.g., [START_REF] Leather | The system water, calcium carbonate, carbonic acid[END_REF][START_REF] Kendall | The solubility of calcium carbonate in water[END_REF][START_REF] Johnston | The solubility-product constant of calcium and magnesium carbonates[END_REF][START_REF] Frear | Solubility of calcium carbonate in certain aqueous solutions[END_REF][START_REF] Greenwald | The dissociation of calcium and magnesium carbonates and bicarbonates[END_REF][START_REF] Martynova | Dissociation constants of certain scale forming salts[END_REF], 1989;[START_REF] Millero | The solubility of calcite, strontianite, and witherite in NaCl solutions at 25°C[END_REF] and super-ambient temperatures [START_REF] Miller | A portion of the system calcium carbonate-carbon dioxide-water, with geological implications[END_REF]Weyl, 1959;[START_REF] Ellis | The solubility of calcite in carbon dioxide solutions[END_REF][START_REF] Ellis | The solubility of calcite in sodium chloride solutions at high temperatures[END_REF][START_REF] Segnit | The solubility of calcite in aqueous solutions-I. The solubility of calcite in water between 75° and 200° at CO 2 pressures up to 60 atm[END_REF][START_REF] Malinin | The solubility of calcite in homogeneous H 2 O-NaCl-CO 2 systems in the 200-600 degrees C temperature interval[END_REF][START_REF] Plummer | The solubilities of calcite, aragonite and viterite in CO 2 -H 2 O solutions between 0 and 90 o C, and evaluation of aqueous model for the system CaCO 3 -CO 2 -H 2 O. Geochim[END_REF][START_REF] Wolf | Solubility of calcite in different electrolytes at temperatures between 10 and 60°C and at CO 2 partial pressures of about 1 kPa[END_REF]. In contrast, the stability constants of the aqueous Ca hydroxo-carbonate complexes have been mostly investigated at room temperatures (i.e., [START_REF] Hopkins | The Solution Thermochemistry of Polyvalent Electrolytes. I. Calcium Hydroxide[END_REF][START_REF] Nakayama | Calcium activity, complex and ion-pair in saturated CaCO 3[END_REF][START_REF] Lafon | Calcium complexing with carbonate ion in aqueous solutions at 25°C and 1 atmosphere[END_REF][START_REF] Dyrssen | Ionic medium effects in sea water -A comparison of acidity constants of carbonic acid and boric acid in sodium chloride and synthetic sea water[END_REF][START_REF] Reardon | Thermodynamic properties of the ion pairs MgCO 3 o and CaCO 3 o from 10 to 50 o C[END_REF] with only few measurements up to 90°C [START_REF] Plummer | The solubilities of calcite, aragonite and viterite in CO 2 -H 2 O solutions between 0 and 90 o C, and evaluation of aqueous model for the system CaCO 3 -CO 2 -H 2 O. Geochim[END_REF]. To determine the distribution of the Ca aqueous species in equilibrium with calcite one has to know the solution pH value in addition to the total Ca concentration and pCO 2 . pH measurements at temperatures above 100°C are, however, technically challenging which explains why very few studies have been devoted to the temperature dependence of aqueous Ca-carbonate complexes stability. Using various modeling approaches, a number of authors found good agreement between calculated and measured calcite solubilities [START_REF] Gal | Calcium carbonate solubility: a reappraisal of scale formation and inhibition[END_REF][START_REF] Garcia | Prediction of mineral scale formation in geothermal and oilfield operations using the Extended UNIQUAC model Part II. Carbonatescaling minerals[END_REF]). Yet, the mostly frequently used SUPCRT92 software packages that encode the revised-HKF equation of state [START_REF] Johnson | SUPCRT92 -A software package for calculating the standard molal thermodynamic properties of minerals, gases, aqueous species, and reactions from 1 bar to 5000 bar and 0 degrees C to 100 degrees[END_REF] and HCh code [START_REF] Shvarov | HCh: new potentialities for the thermodynamic simulation of geochemical systems offered by Windows[END_REF], yield large disagreements with the experimental result and as such require improvement of their thermodynamic database as shown below. [START_REF] Duan | Coupled phase and aqueous species equilibrium of the H 2 O-CO 2 -NaCl-CaCO 3 system from 0 to 250°C, 1 to 1000 bar with NaCl concentrations up to saturation of halite[END_REF] attempted to calculate thermodynamic parameters and activity coefficients of solution components for description of calcite solubility up to 250°C and 1000 bars. Based on available literature data, they calculated the Gibbs free energy of formation, enthalpy and entropy of CaCO 3 °(aq) and CaHCO 3 + (aq), and suggested Pitzer parameters for activity coefficients.

Potentiometric measurements at elevated temperatures and partial CO 2 pressures allow determination of the solubility products of carbonate minerals as well as the formation constants of aqueous complexes under rigorous control of experimental parameters. There are only a few electrode systems that are capable of providing in situ pH values above 100 °C and pressures higher than saturated water pressures. The most reliable technique known today to measure the pH at hydrothermal conditions is the hydrogen electrode concentration cell [START_REF] Palmer | Aqueous high temperature solubility studies. I. The solubility of boehmite at 150 o C as a function of ionic strength and pH as determined by "in situ" measurements[END_REF]Bénézeth et al., 2009b), which has been used to measure the solubility of various carbonate minerals [START_REF] Bénézeth | Dawsonite synthesis and reevaluation of its thermodynamic properties from solubility measurements: Implications for mineral trapping of CO 2[END_REF](Bénézeth et al., , 2009b[START_REF] Bénézeth | Experimental determination of the solubility product of magnesite at 50 to 200°C[END_REF][START_REF] Bénézeth | Experimental determination of the solubility product of dolomite at 50 to 250°C[END_REF][START_REF] Lindner | On the growth of anhydrous Mg-bearing carbonates -Implications from norsethite growth kinetics[END_REF]. The pH measurements in hydrothermal systems at temperatures above 300°C are possible with an Y 2 O 3 -doped zirconia ceramic electrode [START_REF] Lvov | Progress on yttria-stabilized zirconia sensors for hydrothermal pH measurements[END_REF]. A persistent problem is the high resistance of the zirconium electrode at low temperatures and the design of the reference electrode. The solid-contact glass electrode with special glasses developed to sustain both high pressures and temperature also has not been used widely (i.e., [START_REF] Pokrovski | Experimental determination of the stability constants of NaSO 4and NaB(OH) 4 0 in hydrothermal solutions using a new high-temperature sodium-selective glass electrode -Implications for boron isotopic fractionation[END_REF]Pokrovsky et al., 2009b). The key factor in the successful use of these glass electrodes is the reliability of the reference electrode system. The internal Ag/AgCl electrodes immersed in experimental Cl --containing solution cell without transfer were used for carbonate dissolution kinetics studies at temperatures below 100-150°C and pCO 2 up to 55 bars (Pokrovsky et al., 2009a,b;[START_REF] Golubev | Siderite dissolution kinetics in acidic aqueous solutions from 25 to 100°C and 0 to 50 atm pCO 2[END_REF]. However, the stability of this electrode system at temperatures above 100°C is quite poor due to the dissolution of the reference electrode components in the test solution. The systems with external reference electrode use a liquid junction salt bridge with the flow direction from the reactor to the external media [START_REF] Pokrovski | Experimental determination of the stability constants of NaSO 4and NaB(OH) 4 0 in hydrothermal solutions using a new high-temperature sodium-selective glass electrode -Implications for boron isotopic fractionation[END_REF]). An alternative system of liquid junction electrodes employs the injection of electrolyte in the reactor via a high-pressure pump [START_REF] Reukov | The determination of dissociation constant of HCl at the temperature 350°C and the pressure 200 bars with ceramic electrode[END_REF][START_REF] Tagirov | A potentiometric study of the stability of aqueous yttrium-acetate complexes from 25 to 175 °C and 1-1000 bar[END_REF]. However, the large difference in pressure and temperature between the internal solution and the external reservoir requires a high flux of the electrolyte to the system, which limits the wide use of such devices, especially for long-term solubility or kinetic measurements. Spectrophotometric pH determination using indicators (such as 2-napthol/4-nitrophenol) that are thermally stable has been used recently to determine the first and second ionization constants of carbonic acid and the ion pair formation constants for NaHCO 3 °(aq) and NaCO 3 - (aq) up to 200°C [START_REF] Stefánsson | Carbonic acid ionization and the stability of sodium bicarbonate and carbonate ion pairs to 200 °C -A potentiometric and spectrophotometric study[END_REF] as well as the formation constants for MgHCO 3 + (aq) and MgCO 3(aq) up to 150°C [START_REF] Stefánsson | Potentiometric and spectrophotometric study of the stability of magnesium carbonate and bicarbonate ion pairs to 150°C and aqueous inorganic carbon speciation and magnesite solubility[END_REF]. This technique, however, requires the knowledge of the ionization constants of the indicators as well as their thermal stability.

The present study is devoted to the determination of calcite solubility at 100, 120, 140 and 160°C and 1-50 bar pCO 2 in 10 -3 -0.1 mol•kg -1 NaCl solutions. For this purpose, we used a glass sodium-selective electrode as a reference electrode in cells without liquid junction in solutions having constant sodium activity. These solid-contact, high pressure and high temperature electrodes are commercially available from various Russian companies (Potential®, Econix-Expert®) and are proven to operate in the circumneutral range of pH to temperatures and pressures up to 200-220°C and 300-500 bars, respectively [START_REF] Pokrovski | Experimental determination of the stability constants of NaSO 4and NaB(OH) 4 0 in hydrothermal solutions using a new high-temperature sodium-selective glass electrode -Implications for boron isotopic fractionation[END_REF][START_REF] Zotov | Potentiometric study of the stability of Eu 3+ acetate complexes as a function of pressure (1-1000 bar) at 25°C[END_REF]. The measurements of pH and pNa electrode potentials in the Na-Ca-Cl-CO 2 -H 2 O system using different buffer solutions up to 160°C allowed the determination of the stability constants of the sodium and calcium carbonate and bicarbonate ion pairs.

EXPERIMENTAL METHODS

2.1

In-situ pH measurements using a pH/pNa -Ag/AgCl reference electrode

The solution pH was measured in-situ using a commercial solid contact Li-Sn alloy pH electrode coupled with a home-made Ag/AgCl reference electrode placed in an extensible Teflon container (called "silfon" below) filled with 3.5 M KCl and saturated with AgCl (Fig. 1). The outlet of the silfon had a 45 µm porous filter. Before filling up, the silfon was mechanically expanded such that, during its contraction, a flux of electrolyte from the silfon to the outside solution did not exceed 1 mL per day. This addition of electrolyte increased the salinity of the initial solution by no more than 0.003 mol•kg -1 or 0.3 -3%. After 3-5 days, the flux of electrolyte ceased and new filling of silfon was necessary. This electrode system allowed measuring the potential in the cell:

Sn-Cu, Li-Sn alloy │ H + -selective glass │ test solution ││ 3.5 M KCl -AgCl/Ag

This electrode system had an E° of ≤ -2200 mV and was connected to a high input impedance, high resolution pH meter (Econix-Expert (R) , Russia). The calibration of this system was performed using three buffer solutions: 0.01 mol•kg -1 HCl + 0.1 mol•kg -1 NaCl, 0.1 mol•kg -1 K-phthalate and 0.05 mol•kg -1 Na-tetraborate (Fig. A1 of the Electronic Annex A). The pH values of the calibration solutions were taken from [START_REF] Galster | pH measurement: fundamentals, methods, applications, instrumentation[END_REF]. Fast Nernstian response was observed and electrode potentials exhibited sufficient stability to provide uncertainties of ±0.01 pH units at 25°C and ±0.03 to 0.05 pH units at 150°C. Between two subsequent calibrations, the electrode shift was less than 0.5 mV/day at 25°C, 1 mV/day over 3 days at 60°C, 2-3 mV/day at 100°C and 3-5 mV/day at 150°C. The experimental slope was close to the Nernstian one at each temperature. Since the potential of the Ag/AgCl electrode strongly depends on temperature, there is no isopotential point on this calibration plot (Fig. A1 of the Electronic Annex A). Instead, there is a range of pH values over which lines intersect at different temperatures. The cell potential weakly depends on temperature at pH from 6 to 9, which is convenient for studying carbonate equilibria.

Calibration of the glass sodium-selective electrode

The following liquid-junction electrode system was used to calibrate the electrode:

Sn-Cu, Li-Sn alloy │ Na + -selective glass │ test solution ││ 3.5 M KCl -AgCl/Ag

The calibration was performed in four NaCl solutions (0.001, 0.01, 0.1 and 1.0 mol•kg -1 ) as illustrated in Fig. A2 of the Electronic Annex A. Na + activity in these solutions was calculated using the HCh code [START_REF] Shvarov | HCh: new potentialities for the thermodynamic simulation of geochemical systems offered by Windows[END_REF]. To evaluate the effect of the liquid junction potential in this electrode cell, the silfon electrode was filled with KCl solutions of different concentrations (3, 4 and 5 M). We did not detect a change in cell potential among the different solutions within the precision of our measurements (± 1 mV). Therefore, one can neglect the effect of diffusion and liquid junction potential in this system. The selectivity of the sodium electrode with respect to protons was tested in 0.1 mol•kg -1 NaCl solutions at 4 ≤ pH ≤ 7 and 140°C (Fig. A3 of the Electronic Annex A). It can be seen that at pH < 4.3, the Na electrode potential increases due to competition with H + and thus this electrode is not suitable for pNa measurements at pH < 4.3.

Calibration of the pH electrode in cell without transfer with a Na-electrode as a reference

To calibrate the glass pH-electrode coupled with a Na-selective reference electrode, buffer solutions were prepared from sodium acetate and sodium carbonate having constant Na + (aq) activity of 0.08. The pH values of these solutions were calculated considering the complexation of sodium with acetate and carbonate ions [START_REF] Fournier | Experimental determination of aqueous sodium-acetate dissociation constants at temperatures from 20 to 240°C[END_REF]and this work, respectively). Resulting pH values for each buffer solution are listed in Table A1. For comparison, the pH values measured directly by the pH-glass electrode are also listed in this table. A good agreement can be observed between the measured and calculated data. The differences between E pH and E Na depend only on pH and sodium activity. The direct measurement of these potentials thus provides a new reliable method for pH measurements in solutions having a constant activity of sodium ions in the cell without transfer:

Sn-Cu, Li-Sn alloy │ H + -glass │ test solution │ Na + -glass │ Sn-Cu, Li-Sn alloy.

The calibration curve is presented in Fig. A4 of the Electronic Annex A. The slopes of the calibration lines are close to theoretical Nernst slopes, both reported in Table A1. Using a cell with two glass electrodes provides reliable pH measurements over several days under high temperature and pressure. The maximum operating temperature of the cell is limited by the tin alloy melting point ( 200°C) and the stability of the glass membrane itself. In the circumneutral pH range of our experiments, the glass remained stable for more than 6 days at 140°C, but its corrosion rate increased at 160°C.

Experimental measurement of calcite solubility

Pure Prolabo calcite having a surface area of 0.38 m 2 •g -1 and a grain size of 10-50 µm (Pokrovsky et al., 2009b) was used for the solubility experiments. The powder (~20 grams) was placed in a titanium 450-mL Parr reactor continuously stirred by a magnetically driven Parr stirrer.

Measurements of pH in experiments carried out in pure water were conducted using a cell with transfer:

Sn-Cu, Li-Sn alloy │ H-selective glass │ test solution ││ 3.5 M KCl -AgCl/Ag, whereas, for experiments in 0.1 mol•kg -1 NaCl, pH measurements were performed with a cell without transfer:

Sn-Cu, Li-Sn alloy │ H + -glass │ test solution │ Na + -glass │ Sn-Cu, Li-Sn alloy

The reacting fluid was sampled through a porous 2 µm Ti filter using a manual Ti valve. Experiments were performed at 100, 120, 140 and 160°C and 1 -50 bar pCO 2 . Typical time for achieving equilibrium did not exceed 2 hours as verified by preliminary kinetic experiments (Fig. A7).

Analyses

All input and output solutions were analyzed for total dissolved calcium ([Ca]), alkalinity ([Alk]), and pH as a function of time. Calcium was measured by flame atomic absorption with an uncertainty of ±1 % and a detection limit of 7•10 -8 mol•kg -1 . Alkalinity was determined following a standard HCl titration procedure with an uncertainty of ±2 % and a detection limit of 5•10 -5 mol•kg -1 .

Thermodynamic calculations

We used the HCh software package [START_REF] Shvarov | A general equilibrium criterion for an isobaric-isothermal model of a chemical system[END_REF][START_REF] Shvarov | UNITHERM: A Software Package for Geochemical Equilibrium Modeling. User's Guide[END_REF][START_REF] Shvarov | HCh: new potentialities for the thermodynamic simulation of geochemical systems offered by Windows[END_REF], employing algorithms of free energy minimization. This code uses the Unitherm thermodynamic database containing Gibbs free energies of the simple ions, some polyatomic ions, and aqueous complexes, which are calculated according to the revised Helgeson-Kirkham-Flowers equations of state (HKF) [START_REF] Helgeson | Theoretical prediction of the thermodynamic behaviour of aqueous electrolytes at high pressures and temperatures: Calculation of activity coefficients, osmotic coefficients, and apparent molal and standard and relative partial molal properties to 600°C and 5 kb[END_REF][START_REF] Tanger | Calculation of the thermodynamic and transport properties of aqueous species at high pressures and temperatures: Revised equations of state for the standard partial molal properties of ions and electrolytes[END_REF][START_REF] Shock | Calculation of the thermodynamic properties of aqueous species at high pressures and temperatures[END_REF].

Activity coefficients of all ions were computed using the Helgeson B-dot equation [START_REF] Helgeson | Thermodynamics of hydrothermal systems at elevated temperatures and pressures[END_REF] with parameters from [START_REF] Shvarov | UNITHERM: A Software Package for Geochemical Equilibrium Modeling. User's Guide[END_REF]. The CO 2 fugacity coefficient was calculated following the Peng-Robinson equation of state [START_REF] Shvarov | UNITHERM: A Software Package for Geochemical Equilibrium Modeling. User's Guide[END_REF] which agrees with the values of Spycher and Read (1988). To insure the consistency of calculations in the H 2 CO 3 -H 2 O system, we calculated the apparent dissociation constants of carbonic acid (Q 1 , Q 2 ) in 0.2 m NaCl solution at 50-200°C, using the thermodynamic parameters of Patterson et al. (1982Patterson et al. ( , 1984)). The agreement was within 0.01-0.04 log units, except for the calculated Q 1 at 150 and 200°C where the HKF evaluation was 0.09 and 0.14 log units lower, respectively. The aqueous species and solid component HKF parameters in the CaCO 3 -H 2 O system were taken from the Slop07 database [START_REF] Shock | Calculation of the thermodynamic and transport properties of aqueous species at high pressures and temperatures: Correlation algorithms for ionic species and equation of state predictions to 5 kb and 1000°C[END_REF][START_REF] Shock | Metal-organic complexes in geochemical processes: Estimation of standard partial molal thermodynamic properties of aqueous complexes between metal cations and monovalent organic acid ligands at high pressures and temperatures[END_REF][START_REF] Shock | Inorganic species in geologic fluids: Correlations among standard molal thermodynamic properties of aqueous ions and hydroxide complexes[END_REF][START_REF] Sverjensky | Prediction of the thermodynamic properties of aqueous metal complexes to 1000°C and 5 kb[END_REF] [START_REF] Geopig | [END_REF].

RESULTS

Determination of sodium-carbonate and sodium-bicarbonate stability constants at 100-160°C.

To use the Na + -selective electrode as a reference for studying carbonate systems, the complexation of the sodium ion with carbonate and bicarbonate ions and other solution constituents should be taken into account. There are four major complexation reactions for Na + (aq):

Na + + HCO 3 -= NaHCO 3 °(1)      3 3 3 HCO Na NaHCO NaHCO a a a K Na + + CO 3 2-= NaCO 3 - (2)      2 3 3 3 CO Na NaCO NaCO a a a K Na + + Cl -= NaCl o (3)      Cl Na NaCl NaCl a a a K Na + + OH -= NaOH o (4)      OH Na NaOH NaOH a a a K
where a i and K i stand for the activity and stability constant of the ith aqueous species or complex, respectively. The stability constants of the sodium chloride complex (reaction 3) were determined in a wide range of temperature using conductivity techniques [START_REF] Quist | The electrical conductances of some alkali metal halides in aqueous solutions from 0-800°C and at pressures to 4,000 bars[END_REF]. The stability constants of the NaOH°(aq) complex (reaction 4) were also studied at high temperature by [START_REF] Ho | Conductivity measurements of dilute aqueous LiOH, NaOH, and KOH solutions to high temperatures and pressures using a flow-through cell[END_REF] and [START_REF] Fuangswasdi | Enthalpies of dilution of NaOH, KOH, and HCl and thermodynamic quantities for the formation of these species from their constituent ions in aqueous solution from 300 to 370[END_REF]. In contrast and until recently, sodium-carbonate andbicarbonate complexes (reactions 1 and 2) were studied only at room temperatures, reviewed in the study of [START_REF] Stefánsson | Carbonic acid ionization and the stability of sodium bicarbonate and carbonate ion pairs to 200 °C -A potentiometric and spectrophotometric study[END_REF] who measured the carbonic acid ionization and the stability of sodium bicarbonate and carbonate ion pairs up to 200°C using potentiometric and spectrophotometric approaches.

In the present study, the stability constants were determined using an isothermal cell with transfer operating simultaneously with two electrode systems: Sn-Cu, Li-Sn alloy │ Na + -selective glass │ test solution ││ 3.5 M KCl -AgCl/Ag, Sn-Cu, Li-Sn alloy │ H-selective glass │ test solution ││ 3.5 M KCl -AgCl/Ag Four tests solutions having constant ionic strength of 0.1 and 1.0 were used (0.9 mol•kg -1 NaCl + 0.05 mol•kg -1 Na 2 CO 3 ; 0.6 mol•kg -1 NaCl + 0.2 mol•kg -1 Na 2 CO 3 ; 0.5 mol•kg -1 Na 2 CO 3 ; 0.1 mol•kg -1 NaCl). The electrode cell potential in solutions having different Na 2 CO 3 molalities is illustrated in Fig. A5 of the Electronic Annex A, and the results of the calculations of NaCO 3 -stability constant in the different buffer solutions are listed in Table A2 of the Electronic Annex A. The study of the Na + -bicarbonate ion pair stability was performed in three buffer solutions with a constant CO 2 partial pressure of 10 bars (0.9 mol•kg -1 NaCl + 0.05 mol•kg -1 Na 2 CO 3 ; 0.6 mol•kg -1 NaCl + 0.2 mol•kg -1 Na 2 CO 3 ; 1 mol•kg -1 NaHCO 3 ). The dependence of the electrode cell potential as a function of sodium bicarbonate molality is illustrated in Fig. A6, and the results of the calculation of NaHCO 3 ° stability constants are reported in Table A3 of the Electronic Annex A.

To calculate the NaHCO 3 ° and NaCO 3 -stability constants from our measurements, activity coefficients of all ions were computed using the Helgeson B-dot equation [START_REF] Helgeson | Thermodynamics of hydrothermal systems at elevated temperatures and pressures[END_REF]. All calculations were performed using the HCh computer code [START_REF] Shvarov | HCh: new potentialities for the thermodynamic simulation of geochemical systems offered by Windows[END_REF] and the Gibbs free energies of formation of sodium -carbonate and bicarbonate complexes were derived from the obtained constants. The calculated values of logK NaHCO3° and logK NaCO3-are reported in Table 1, and plotted as a function of temperature in Fig. 2 together with literature data. The values obtained in this study are in excellent agreement with the experimental values reported by [START_REF] Stefánsson | Carbonic acid ionization and the stability of sodium bicarbonate and carbonate ion pairs to 200 °C -A potentiometric and spectrophotometric study[END_REF].

The determination of NaHCO 3 °(aq) and NaCO 3 -(aq) stability constants allowed precise calculation of the pH values of the carbonate buffer solutions and produced a practical (operational) scale of pH measurement using the sodium-selective electrode as a reference in the cell without transfer.

Calcite solubility

The results of calcite solubility measurements in pure water at 120, 140 and 160 °C as a function of pCO 2 are shown in Fig. 3 and listed in Table A4 of the Electronic Annex A. It can be seen that the solubility of calcite in pure water increases with pCO 2 and decreases with temperature, in agreement with the data of [START_REF] Ellis | The solubility of calcite in carbon dioxide solutions[END_REF] as shown in Fig. 3A. Since [START_REF] Ellis | The solubility of calcite in carbon dioxide solutions[END_REF] did not use a uniform temperature step in his experiments, we reported in Fig. 3A the results he obtained at the temperatures closest to those investigated in the present study (120, 140, 160 °C). Prediction of calcite solubility using the Slop07 database [START_REF] Geopig | [END_REF] significantly overestimates the calcite solubility measured in the present study. It should be noted that our measured pH values are the same within uncertainties as those calculated using Slop07 (see Fig. 3B). This suggests that the discrepancy between measured and calculated calcite solubility arises from the values of HKF parameters of Ca aqueous species in the Slop07 database.

Results of calcite solubility in 0.1 mol•kg -1 NaCl are listed in Table A5 of the Electronic Annex A and illustrated in Fig. 4A. As in the case of calcite solubility in pure water, calculations by the HCh software using the parameters in the slop07 database, together with activity coefficients calculated with the B-dot equation, overestimates calcite solubility in 0.1 mol•kg -1 NaCl, especially at T ≤ 140°C.

In addition, calculated pHs with Slop07 database systematically overestimate the experimental pH values at all temperature (Fig. 4B).

Results of calcite solubility measurements in Na 2 CO 3 solutions without CO 2 pressure are illustrated in Fig. 5 and the raw experimental data are listed in Table A6 of the Electronic Annex A.

It can be seen that there is no dependence of calcite solubility on Na 2 CO 3 concentration up to 0.2 mol•kg -1 Na 2 CO 3 .

DISCUSSION: THERMODYNAMIC MODELING

The thermodynamic properties of calcite are well known [START_REF] Bäckström | The thermodynamic properties of calcite and aragonite[END_REF][START_REF] Robie | Thermodynamic properties of minerals and related substances at 298.15K and 1 bar (105 Pa) pressure and at higher temperatures[END_REF][START_REF] Garrels | A chemical model for seawater at 25°C and one atmosphere total pressure[END_REF][START_REF] Plummer | The solubilities of calcite, aragonite and viterite in CO 2 -H 2 O solutions between 0 and 90 o C, and evaluation of aqueous model for the system CaCO 3 -CO 2 -H 2 O. Geochim[END_REF]. The thermodynamic parameters of Ca 2+ (aq) are based on the enthalpy of CaO dissolution [START_REF] Thorvaldson | Studies on the thermochemistry of the compounds 23 occurring in the system CaO-Al 2 O 3 -SiO 2 . I. The heat of solution of calcium oxide in hydrochloric acid[END_REF] and the temperature dependence of calcite solubility product [START_REF] Reardon | Thermodynamic properties of the ion pairs MgCO 3 o and CaCO 3 o from 10 to 50 o C[END_REF]. The HKF parameters for Ca 2+ originate from [START_REF] Shock | Calculation of the thermodynamic and transport properties of aqueous species at high pressures and temperatures: Correlation algorithms for ionic species and equation of state predictions to 5 kb and 1000°C[END_REF].

The first main result of the present study is that the HKF model with the Slop07 database significantly overestimates calcite experimental solubility. This discrepancy increases with pCO 2 and temperature . Therefore, we revised the calcite solubility product and the stability constants of calcium bicarbonate and carbonate complexes based on the new experimental data obtained in this study. From our experiments, we can define the total Ca concentration by taking into account the stability constants for the reactions below (5-9). Ca 2+ + HCO 3 -= CaHCO 3 +

(5)

     3 3 3 HCO Сa СaHCO СaHCO a a a K Ca 2+ + CO 3 2-= CaCO 3 o (6)      2 3 2 3 3 CO Сa СaCO СaCO a a a K
The solubility of calcite can be expressed as:

[Ca] T =[Ca 2+ ]+[CaHCO 3 + ]+[CaCO 3 o ]
where [i] stands for the molailty of the ith aqueous species.

The reactions describing calcite solubility with their corresponding stability constants can be represented in the following way:

CaCO 3 (s)+2H + =Ca 2+ +H 2 O(l)+CO 2 (g) (7) 2 H CO Ca 7 2 2    a P a K CaCO 3 (s)+H + =CaHCO 3 + (8)    H CaHCO 8 3 a a K CaCO 3 (s)=CaCO 3 o (9)   3 CaCO 9 a K
The total concentration of calcium in equilibrium with calcite is then given by: (10)

0 9 1 H 8 CO 2 2 H 7 2 ] Ca [    K a K P a K T     
where  0,  1 , and  2 are the activity coefficients of neutral, singly-charged and divalent species, respectively. The value of K 9 (Eqn. 9) could be directly derived from calcite solubility in alkaline solutions as illustrated in Fig. 5 which presents calcite solubility in Na 2 CO 3 solutions without added carbon dioxide. The logarithm of K 9 value was found to be equal to -5.10, -5.13, and -5.16 at 120, 140 and 160 °C, respectively. However, the precision of pH measurements at elevated temperatures is not sufficient to rigorously resolve Eqn. 10 for K 7 and K 8 . As such, we used a computer program OptimA to determine the thermodynamic parameters in the system. The program OptimA defines simultaneously at different temperatures the Gibbs free energies of a limited number of species in aqueous solutions based on experimental data using the HCh package for Windows [START_REF] Shvarov | A suite of programs, OptimA, OptimB, OptimC, and OptimS compatible with the Unitherm database, for deriving the thermodynamic properties of aqueous species from solubility, potentiometry and spectroscopy measurements[END_REF].

Details on operating the OptimA program and examples of calculations are given in the Electronic

Annex B. To accurately define the Gibbs free energy of the calcium aquo ion and its carbonate and bicarbonate complexes, the program processed simultaneously all data on the solubility of calcite measured in the present work as a function of temperature and pCO 2 in pure water, 0.1 mol•kg -1 NaCl and Na 2 CO 3 solutions listed in Table A4, A5 and A6 of the Electronic Annex A. The program adjusts the Gibbs free energy of formation of Ca 2+ (aq) to provide the best fit of the full set of experimental data. The results of the calculations demonstrate that, in order to reproduce the solubility of calcite above 100° C, one has to take into account the formation of Ca-bicarbonate complexes. The values of the formation constants of the sodium and calcium bicarbonate and carbonate complexes calculated based on our experimental data are listed in Table 1.

The temperature dependence of the calcite solubility product (CaCO 3 (s) = Ca 2+ +CO 3 2-) and the stability constants of CaHCO 3 + (aq) and CaCO 3 °(aq) complexes, reported in Table 1, are depicted in Figs. 6, 7 and 8, respectively. These values are compared to the results of previous literature studies (reported as symbols in these figures) as well as prediction based on values extrapolated from the work of [START_REF] Plummer | The solubilities of calcite, aragonite and viterite in CO 2 -H 2 O solutions between 0 and 90 o C, and evaluation of aqueous model for the system CaCO 3 -CO 2 -H 2 O. Geochim[END_REF], [START_REF] Duan | Coupled phase and aqueous species equilibrium of the H 2 O-CO 2 -NaCl-CaCO 3 system from 0 to 250°C, 1 to 1000 bar with NaCl concentrations up to saturation of halite[END_REF]Li (2008), Facq et al. (2014), [START_REF] Huang | Extended Deep Earth Water Model for predicting major element mantle metasomatism[END_REF] or using the HKF model with the Slop07 database [START_REF] Shock | Calculation of the thermodynamic and transport properties of aqueous species at high pressures and temperatures: Correlation algorithms for ionic species and equation of state predictions to 5 kb and 1000°C[END_REF][START_REF] Shock | Metal-organic complexes in geochemical processes: Estimation of standard partial molal thermodynamic properties of aqueous complexes between metal cations and monovalent organic acid ligands at high pressures and temperatures[END_REF][START_REF] Shock | Inorganic species in geologic fluids: Correlations among standard molal thermodynamic properties of aqueous ions and hydroxide complexes[END_REF][START_REF] Sverjensky | Prediction of the thermodynamic properties of aqueous metal complexes to 1000°C and 5 kb[END_REF]. The HKF parameters reported by [START_REF] Huang | Extended Deep Earth Water Model for predicting major element mantle metasomatism[END_REF] and [START_REF] Facq | In situ Raman study and thermodynamic model of aqueous carbonate speciation in equilibrium with aragonite under subduction zone conditions[END_REF] have been used here to recalculate the values log K calcite , log K CaHCO 3 + and log K CaCO 3 ° using the HCh code and Unitherm database [START_REF] Shvarov | HCh: new potentialities for the thermodynamic simulation of geochemical systems offered by Windows[END_REF] as plotted in Figs. 678. It can be seen that the values calculated from the data of [START_REF] Facq | In situ Raman study and thermodynamic model of aqueous carbonate speciation in equilibrium with aragonite under subduction zone conditions[END_REF] for the solubility product of calcite are almost identical to the predictions made using Slop07.

However, their data are in very good agreement with our experimental values for log K CaHCO 3 + and log K CaCO 3 °, which is not the case for the data generated from the revised data provided by [START_REF] Huang | Extended Deep Earth Water Model for predicting major element mantle metasomatism[END_REF].

Based on the new values of calcite solubility product and CaHCO 3 + and CaCO 3 ° stability constants generated in the present study, we calculated calcite solubility (as total calcium concentration) and equilibrium pH in pure water and in 0.1 mol•kg -1 NaCl as a function of pCO 2 (Figs. 9 and 10, respectively). New recommended log K values (Table 1) allow an accurate description of experimental data at 100-160°C within the uncertainty of measurements. Experimental pH values of 0.1 mol•kg -1 NaCl solutions in equilibrium with different partial pressure of pCO 2 at 120, 140 and 160°C agreed (± 0.1-0.2 pH units) with those calculated using the HCh code with implemented stability constants of NaCO 3 -and NaHCO 3 ° or those taken from [START_REF] Stefánsson | Carbonic acid ionization and the stability of sodium bicarbonate and carbonate ion pairs to 200 °C -A potentiometric and spectrophotometric study[END_REF] as illustrated Fig. 10B.

These results show that the proposed potentiometric set-up with a pH electrode in a cell without transfer and a Na-electrode as a reference can be used in aqueous solutions where the activity of the sodium ion is known or can be calculated. This cell is highly stable could be used, for example, to monitor pH in the ocean, groundwaters, soda lakes, and other natural environments with constant activity of sodium ion.

CONCLUSIONS

The use of a solid-contact pH glass electrode coupled with Ag/AgCl/KCl (sat) liquid-junction reference electrode allows robust and precise measurement of pH in CO 2 -H 2 O systems up to 160°C and 60 bar pCO 2 . Conceptually new in situ pH-measurement system consisting of a pH-glass electrode and a reference Na + electrode (solid contact system without liquid junction) has been also described which allows continuous pH measurements in circumneutral and alkaline solutions with fixed activity of Na + ion. This electrode system was used to measure calcite solubility and to assess calcite solubility product and the stability of CaCO 3 °, CaHCO 3 + and NaCO 3 -complexes at 100, 120, 140 and 160°C and 1 to 50 atm pCO 2 .

It is hoped that the new experimental sets-ups and the results of this study could renew the interest for modeling carbonate equilibria in low-temperature hydrothermal conditions to better understand serpentinization and other natural processes, and predict the consequences of CO 2 sequestration and storage efforts to overcome climate change. These data also could help in the modeling of the subsurface biosphere in continental and marine hydrothermal systems, deep sediments, and basement rocks on the seafloor.

Acknowledgments

The authors thank Sergey Golubev for assistance in analytical measurements, Alexander Zotov for discussing the results. The manuscript was improved thanks to thorough and constructive comments of Everett Shock, Dimitri Sverjensky and two anonymous reviewers. This work was supported by CNRS/INSU and the Russian Science Foundation, grant 19-17-00200.


Ben Faten

Chaabane

Thomas Lopez

Laetitia Pradere

Béatrice Dagens

Aloyse Degiron

Modelling active dipolar media in photonics and optoelectronics with the finite element method

Keywords: finite element method, nonlinear active media, nonlinear nanophotonics, strong coupling, collective states

We introduce a simple and efficient numerical method to model the nonlinear response of active dipolar media in arbitrary photonic environments in the harmonic regime. At variance with existing approaches that require a deep knowledge of time-domain resolution techniques, our model is based on a reformulation in the frequency domain of the rate equations that govern the behavior of active dipolar media. The resulting equations can then be readily solved with popular solvers based on the finite element method. We discuss the validity of the method for a single layer of organic molecules and then illustrate its potential by studying the strong coupling between the molecules and an array of plasmonic nanoparticles. The model correctly predicts the formation of strongly coupled states and also captures more subtle features linked to dipole-dipole interactions.

Introduction

Many research topics in photonics and optoelectronics involve interactions between optical nanostructures and dipolar media (e.g. layers of organic molecules, semiconducting nanocrystals, quantum wells…). This is of course the case for lasers, with current developments focused on topologically protected cavities [START_REF] Bandres | Topological Insulator Laser: Experiments[END_REF] and parity-time symmetry for enhanced modal purity [START_REF] Feng | Single-Mode Laser by Parity-Time Symmetry Breaking[END_REF]. One can also cite the staggering breakthroughs accomplished with optical cavities in the strong coupling regime, such as the possibility of using hybrid light-matter states to enhance the optoelectronic properties of the dipolar medium [START_REF] Orgiu | Conductivity in Organic Semiconductors Hybridized with the Vacuum Field[END_REF] or to engineer the kinetics of chemical reactions in radically new ways [START_REF] Hutchison | Modifying Chemical Landscapes by Coupling to Vacuum Fields[END_REF]. Also worth mentioning are novel avenues for second harmonic generation in nanoscale systems [START_REF] Cox | A Nonlinear Switching Mechanism in Quantum Dot and Metallic Nanoparticle Hybrid Systems[END_REF][START_REF] Gigli | Quasinormal-Mode Non-Hermitian Modeling and Design in Nonlinear Nano-Optics[END_REF] and the emerging potential of collective superradiant states for infrared optoelectronic devices [START_REF] Vasanelli | Ultra-Strong Light-Matter Coupling and Superradiance Using Dense Electron Gases[END_REF].

To model the interactions of a dipolar medium and its surrounding photonic environment, a semiclassical treatment based on the matrix density formalism can be applied in many cases of interest [START_REF] Sukharev | Optics of Exciton-Plasmon Nanomaterials[END_REF]. This general framework makes it possible to evaluate how the different energy levels of the system are populated and, in turn, how these populated states influence the optical properties of the structure under investigation. In mathematical terms, the problem can be formulated as a system of rate equations coupled to Maxwell's equations that must be solved in a self-consistent way. For nanostructures, it is generally not possible to solve these equations analytically due to the complexity of the electromagnetic fields associated with subwavelength features. It is therefore necessary to rely on numerical simulations. Because the rate equations are time-dependent, they are usually solved in the time domain, typically using the Finite-Difference Time-Domain (FDTD) method [START_REF] Sukharev | Optics of Exciton-Plasmon Nanomaterials[END_REF].

Although FDTD has proven highly successful, its implementation is everything but straightforward for the class of problems involving active dipolar media. In fact, there are many pitfalls and subtleties that are not handled by commercial FDTD software and therefore advanced skills are necessary to develop home-built codes [START_REF] Sukharev | Optics of Exciton-Plasmon Nanomaterials[END_REF]. Furthermore, for all the cases that do not deal with transient processes, the photonics community tends to favor another numerical approachthe finite element method (FEM). Contrarily to FDTD, the FEM is usually formulated in the frequency domain. It is used by many groups in the world to solve a large variety of electromagnetic problems such as the scattering of subwavelength nanoparticles [START_REF] Knight | Nanoshells to Nanoeggs to Nanocups: Optical Properties of Reduced Symmetry Core-Shell Nanoparticles Beyond the Quasistatic Limit[END_REF], the mode propagation in optical and plasmonic waveguides [START_REF] Zenin | Directional Coupling in Channel Plasmon Polariton Waveguides[END_REF] or the design of topological metamaterials [START_REF] Yves | Crystalline Metamaterials for Topological Properties at Subwavelength Scales[END_REF]. Although most of the FEM models built by the community deal with passive linear problems, the method has been successfully adapted to active, nonreciprocal and nonlinear systems such as Parity Time symmetric structures [START_REF] Benisty | Implementation of PT Symmetric Devices Using Plasmonics: Principle and Applications[END_REF], nonreciprocal metamaterials [START_REF] Degiron | Efficient Finite Element Resolution of Gyromagnetic and Gyroelectric Nonreciprocal Electromagnetic Problems[END_REF], nonlinear plasmonic waveguides [START_REF] Degiron | Nonlinear Long-Range Plasmonic Waveguides[END_REF] and nonlinear metasurfaces [START_REF] Gigli | Quasinormal-Mode Non-Hermitian Modeling and Design in Nonlinear Nano-Optics[END_REF]. However, none of these developments tackles the important class of problems involving interactions between electromagnetic fields and dipolar media described by rate (or population) equations.

In this article, we address this problem by rewriting the population equations of a dipolar medium in the frequency domain. The expressions are then included in a self-consistent electromagnetic model solved with a commercial finite element solver (Comsol Multiphysics). We first discuss the accuracy and limitations of the model for a simple organic layer surrounded by air. We then illustrate the potential of our approach by examining the strong coupling between this organic layer and an array of plasmonic nanoparticles.

Methods

We consider the case of a two-level dipolar medium and write n1 and n2 the populations in the ground and excited states, respectively. Furthermore, we assume an isotropic distribution of the transition dipole moments and write their norm as µ12. Under these assumptions, the response of the material to a local electric field E ⃗ ⃗ can be expressed using the density matrix method as [START_REF] Sukharev | Optics of Exciton-Plasmon Nanomaterials[END_REF]:

dn 1 dt -γ 21 n 2 = - 1 ℏΩ 0 E ⃗ ⃗ . ∂P ⃗ ⃗ ∂t (1.a) dn 2 dt = - dn 1 dt (1.b) ∂ 2 P ⃗ ⃗ ∂t 2 +( γ 21 + 2γ d ) ∂P ⃗ ⃗ ∂t + Ω 0 2 P ⃗ ⃗ = -κ(n 2 -n 1 )E ⃗ ⃗ , (1.c)
where γ21 is the spontaneous decay rate of the excited state, γd is the pure dephasing rate, ħΩ0 is the transition energy, P ⃗ ⃗ is the macroscopic polarization and κ= 2Ω0µ12 2 /3ħ is a coupling constant.

Note that these expressions take into account the dipole-dipole interactions between individual molecules-a point that will be critical later on.

We look for an equivalent formulation of these equations in the harmonic regime. To this end, we formally express the x, y and z components of the polarization and the electric field as:

p i = 1 2 [|P i | e j[arg(P i )+ωt] + c. c. ] (2.a)
e i = 1 2 [|E i | e j[arg(E i )+ωt] + c. c. ] (2.b)
where Pi and Ei are the i th components of the complex polarization and electric field phasors and i= x, y, z. After some algebraic manipulations, and by noting that in the harmonic regime, the n1 and n2 populations are time-independent on average, we can recast the original system of equations as:

〈n 2 〉 = 1 ℏ Ω 0 γ 21 〈E ⃗ ⃗ . ∂P ⃗ ⃗ ∂t 〉 = -ω 2ℏ Ω 0 γ 21 ∑ |P i | |E i | sin[arg(P i ) -arg(E i )] i (3.a) 〈n 1 〉 = n 0 -〈n 2 〉 (3.b) [Ω 0 2 -ω 2 + jω(γ 21 + 2γ d )]P i = -κ(2〈n 2 〉 -n 0 )E i (3.c)
where < . > denotes time-averaged quantities, i = x, y, z and n0 is the total molecular, atomic or electronic density.

The last step to include equations (3.a)-(3c) in a finite element model in the frequency domain is to express the relative complex permittivity εr of the dipolar medium as a diagonal tensor:

εr,ii = 1 + Pi/ ε0Ei , i=x, y, z (4)

Results and Discussion

Validation of the model

We test the validity of our approach by examining the properties of a molecular layer illuminated by a plane wave under normal incidence. The layer has a thickness t=15 nm and is surrounded by air on both sides. In this example, we consider molecular parameters already considered in the literature [START_REF] Fauché | Plasmonic Opals: Observation of a Collective Molecular Exciton Mode Beyond the Strong Coupling[END_REF], namely, a transition energy ħΩ0 = 2.94 eV, a spontaneous decay rate γ21 = 1 ps -1 , a pure dephasing rate γd = 730 fs -1 and a transition dipole moment µ12=10 D.

We build our model in the commercial software Comsol Multiphysics. We use the predefined environment offered by the RF module to set the electromagnetic part of the problem based on Maxwell's equations. To include the molecular response, we enter eqs 3a-3c manually by creating an additional "Domain ODEs and DAEs" node within the model. The two sets of equations are coupled by ascribing the complex permittivity given by eq 4 to the molecular layer. Because the problem is intrinsically nonlinear, we compute the solution iteratively, using a multifrontal massively parallel sparse direct solver at each step and converging toward the solution with the Newton's method (this is the strategy implemented by default in Comsol).

Figure 1a shows the evolution of the populations <n1> and <n2> as a function of the incident power.

The total molecular density chosen in this first example is n0=10 25 m -3 . In all cases, the molecules are resonantly pumped around the frequency Ω0 of the molecular transition, which is the expected behavior. As the incident power increases, the fundamental level is gradually depleted and it is possible to reach a point where approximately half of the population is in the excited state.

Importantly, our model does not converge for incident powers that exceed this threshold because no population inversion can be observed with a two-level system in the steady state. The reason why the model does not find solutions for <n2> larger than <n1> is encoded in eq 3c. The term 2<n2>-n0 = <n2>-<n1> on the right side of this equation must always be negative, otherwise the molecular susceptibility and permittivity εr are not self-consistent (i.e. the susceptibility and permittivity would exhibit an anomalous dispersion outside the resonance, as can be seen by noting that eq 3c can be recast as the i th component of the electric susceptibility Pi/ε0Ei).

Figures 1b-1c show that the real and imaginary parts of εr have a resonant shape, as expected for a material with resonant transitions around Ω0. Furthermore, the strength of the resonance weakens as the power increases, which is also consistent with the fact that the interactions with the incident field decrease as more and more molecules leave the ground state.

The molecular film considered in this first example is sufficiently thin that all the solved quantities are almost invariant along the propagation direction. Thus, the same problem can also be directly solved in the time domain with eqs 1a-1c because the latter admit tractable solutions (i.e. without using brute-force FDTD) in the limit of vanishingly thin films. The results of the time-domain calculations, performed with the help of a small code written in Octave, are shown as filled symbols on Figures 1a-1c. They are in excellent agreement with the Comsol simulations, validating our FEM model.

We close this section by noting that the other key parameter that influences the response of the molecular layer is the total molecular density n0. This point is illustrated on Figure 1d where the real part of the permittivity is plotted for a fixed incident power P=12.5 MW/m and increasing values of n0. One can see that the Lorentzian-like features sharpen as n0 increases, which is not surprising given that more molecules are involved in the interactions. Interestingly, εr can even reach negative values for the highest molecular concentrations, which is indicative of a metallic behavior. Note that in this case also, the model only finds solutions characterized by <n1> > <n2>.

We remind the reader that this is not an intrinsic limitation of our approach but a consequence of the fact that we consider a two-level system in the steady state. Had we considered a medium with three levels or more, population inversion would have been possible.

Application to Strong Coupling

We now illustrate the potential of our approach by examining the strong coupling between the generic two-level molecular emitters investigated above and an array of resonant plasmonic nanoparticles. We have chosen this example because it has a rich physics that has already been well studied with FDTD simulations [START_REF] Fauché | Plasmonic Opals: Observation of a Collective Molecular Exciton Mode Beyond the Strong Coupling[END_REF][START_REF] Sukharev | Plasmon Transmission Through Excitonic Subwavelength Gaps[END_REF]. In other words, we can use this system to verify that our FEM model captures all the features predicted with simulations performed in the time domain-and also provide additional quantitative insight on this system.

The structure under investigation, shown in Figure 2, differs from the existing literature on two points. First, the plasmonic array has not the same geometry. Second, the coupled organic/plasmonic system is embedded within an inverted organic LED architecture, so that the results presented below may be used to develop electroluminescent devices in the future. The full stack consists of an Al cathode, a TiO2 film to inject the electrons, the plasmonic array, the generic two-level dipolar medium and a transparent cover with a refractive index close to that of ITO. It may be argued that this stack lacks of hole injection layer (HIL) and transport hole layer (HTL), but studies on quantum dot LEDs have shown that these layers and the transparent ITO anode can be optically modeled as a single material [START_REF] Wang | Electroluminescence of Colloidal Quantum Dots in Electrical Contact with Metallic Nanoparticles[END_REF].

Here, we will solely focus on the passive optical properties of the device. Specifically, we model the system following the guidelines exposed previously and examine the reflection properties of the device when the latter is illuminated by an incoming plane wave under normal incidence.

Rather than simulating the full structure, we model a single unit cell and apply periodic boundary conditions on the lateral sides of the computational domain. The periods that will be considered below are always smaller than the wavelength so that no diffraction effect modulates the optical response of the system.

Figure 3a shows the result of the simulations for four molecular densities n0 ranging from 10 23 m - 3 to 5×10 23 m -3 . The reflection spectra exhibit two dips on each side of the molecular transition frequency. Furthermore, the separation between these two minima increases with n0, suggesting that they are the signature of the strong coupling between the molecular system and the plasmonic resonances of the Al nanoparticules. This interpretation is confirmed with the calculations of Figure 3b which shows, for the molecular density n0 = 5×10 23 m -3 , the position of the two minima as a function of the array periodicity P. Also plotted are the energy of the molecular transition ħΩ0 and the dispersion of the plasmonic resonance Epl(P) without the organic layer. A clear anticrossing with a Rabi splitting of 20 meV is observed. Furthermore, the branches have the expected polaritonic behavior, as they can be fitted by the analytical eigenmodes resulting from the hybridization of ħΩ0 and Epl(P) via a coupling term being half of the Rabi splitting [START_REF] Bellessa | Strong Coupling between Surface Plasmons and Excitons in an Organic Semiconductor[END_REF] (continuous lines on Figure 3b).

We have then repeated the simulations for molecular concentrations ranging from 10 24 m -3 to 4×10 24 m -3 . As shown on Figure 3c and Figure 3d, the splitting between the two polaritonic branches further increases and a new minimum emerges near the transition energy of the molecules (2.94 eV). This nearly dispersionless feature has been noted by other authors who studied similar structures in the past [START_REF] Fauché | Plasmonic Opals: Observation of a Collective Molecular Exciton Mode Beyond the Strong Coupling[END_REF][START_REF] Salomon | Strong Light-Molecule Coupling on Plasmonic Arrays of Different Symmetry[END_REF][START_REF] Antosiewicz | Plasmon-Exciton Interactions in a Core-Shell Geometry: From Enhanced Absorption to Strong Coupling[END_REF]. They argued that it corresponds to the emergence of a new collective mode at high molecular concentrations n0. This assertion was supported by a qualitative analytical model and by the fact that the mode emerges when the permittivity of the molecular layer becomes negative. In our case also, we have verified that this mode appears when the real part of εr reaches negative values near the resonance frequency (Figure 1d). As discussed in the previous section, a negative permittivity implies that the medium acts as a metal at optical frequencies, which, by definition, implies a screening of the incoming field-and therefore constructive interferences resulting from the cooperation of all the dipoles of the medium. As a corollary, the results of Figures 3c-3d suggest that dipole-dipole interactions are correctly taken into account by our model-without them, the cooperative behavior between the molecules would not be captured and the third mode would not occur [START_REF] Salomon | Strong Coupling between Molecular Excited States and Surface Plasmon Modes of a Slit Array in a Thin Metal Film[END_REF].

The third mode should not be confused with the signature of the dipoles that remain uncoupled in real experiments [START_REF] Aberra Guebrou | Coherent Emission from a Disordered Organic Semiconductor Induced by Strong Coupling with Surface Plasmons[END_REF] or with the collective subradiant modes resulting from inhomogeneous broadening [START_REF] Houdré | Vacuum-Field Rabi Splitting in the Presence of Inhomogeneous Broadening: Resolution of a Homogeneous Linewidth in an Inhomogeneously Broadened System[END_REF]. To support this claim, we have repeated the calculations after setting the emitter's dephasing time to zero and reducing their decay rate by an order of magnitude, so as to emulate a dipolar medium with vanishing spectral broadening and almost no losses. The results, plotted in the inset of Figure 3d, prove that the third mode still exists in this case-a point that has never been verified to our knowledge.

Conclusion

In conclusion, we have presented a numerical method to include the response of a two-level dipolar medium in electromagnetic finite element simulations in the frequency domain. Our approach has a general validity and can be adapted to other nonlinear problems involving media described by rate/population equations, including those with more than two energy levels. In this regard, it may be of interest for anyone using the finite element method as their primary design tool in photonics, optoelectronics and material sciences. The period is treated as a parameter that will be subsequently varied between 290 nm and 340 nm. The refractive indices of Al and TiO2 are taken from two references of the literature [START_REF] Wang | Electroluminescence of Colloidal Quantum Dots in Electrical Contact with Metallic Nanoparticles[END_REF][START_REF] Rakić | Optical properties of metallic films for vertical-cavity optoelectronic devices[END_REF] while the top layers on top of the molecular film are modeled as a single material with the refractive index of ITO multiplied by 1.25 for the real part [START_REF] Wang | Electroluminescence of Colloidal Quantum Dots in Electrical Contact with Metallic Nanoparticles[END_REF]. Note that we model the Al nanoparticles with a bulk refractive index, but since we deal with nanorods that are separated by hundreds of nanometers, deviations from the bulk (such as nonlocal effects) are weak. . The new collective mode that appears close to the molecular transition is seen as yellow symbols. Inset: result of the simulations when the dephasing rate of the emitters is set to zero and their decay rate reduced by one order of magnitude.

FIGURESFigure 1 .

 1 FIGURES

Figure 2 .

 2 Figure2. Schematic of one unit cell of the periodic structure examined for strong coupling. The Al nanoparticles have a length L = 75 nm and a thickness t = 25 nm. The period is treated as a parameter that will be subsequently varied between 290 nm and 340 nm. The refractive indices of Al and TiO2 are taken from two references of the literature[START_REF] Wang | Electroluminescence of Colloidal Quantum Dots in Electrical Contact with Metallic Nanoparticles[END_REF][START_REF] Rakić | Optical properties of metallic films for vertical-cavity optoelectronic devices[END_REF] while the top layers on top of the molecular film are modeled as a single material with the refractive index of ITO multiplied by 1.25 for the real part[START_REF] Wang | Electroluminescence of Colloidal Quantum Dots in Electrical Contact with Metallic Nanoparticles[END_REF]. Note that we model the Al nanoparticles with a bulk refractive index, but since we deal with nanorods that are separated by hundreds of nanometers, deviations from the bulk (such as nonlocal effects) are weak.

Figure 3 .

 3 Figure 3. (a)Reflection spectra for a structure with a period P = 320 nm. The molecular concentrations n0 are 10 23 m -3 for the light blue curve, 2×10 23 m -3 for the blue curve, 3×10 23 m -3 for the red curve and 5×10 23 m -3 for the gray curve. The vertical dashed line indicates the molecular transition energy ħΩ0 and the black spectrum reveals the uncoupled plasmonic resonance Epl when the simulations are performed without the organic layer. (b) Dispersion relation for n0 = 5×10 23 m -3 . The filled symbols are the output of the FEM simulations for different periods P while the continuous curves result from the diagonalization of the 2×2 Hamiltonian hybridizing ħΩ0 and Epl(P). The uncoupled parameters ħΩ0 and Epl(P) are visualized as gray and black dashed lines, respectfully. (c) Reflection spectra for molecular concentrations n0 = 10 24 m -3 , 2×10 24 m -3 , 3×10 24 m -3 and 4×10 24 m -3 . (d) Dispersion relation for n0 = 4×10 24 m -3 . The new collective mode that appears close to the molecular transition is seen as yellow symbols. Inset: result of the simulations when the dephasing rate of the emitters is set to zero and their decay rate reduced by one order of magnitude.

ACKNOWLEDGMENT

This work was supported by Groupe PSA through Contract CIFRE No 2017/0535 as part of the projects of the OpenLab PhOVeA.

AUTHOR INFORMATION

Corresponding Author

*E-mail: (A.D.) aloyse.degiron@u-paris.fr