Surface reactivity and cation non-stoichiometry in BaZr 1 xYxO 3 ( x = 0 – 0 . 2 ) exposed to CO 2 at elevated temperature

Surface formation of BaCO3 mainly proceeds by transfer of Y to the Ba-site while a BaZrO3 formula unit is consumed.


Introduction
A II B IV O 3 perovskites have been widely studied as materials for solid state electrochemical devices. [1][2][3] Alkali earths such as strontium and barium on the A-site implies the possibility of reaction with CO 2 to form stable carbonates at low to intermediate temperatures in electrochemical systems for the conversion of hydrocarbons. 1,3-10 The CO 2 -sensitivity of the materials can originate from thermodynamic phase instability or surface reactivity, and minor compositional modications can affect the chemical stability. 11 Consequently, traditional solid-state reaction synthesis can be insufficient to avoid compositional heterogeneities that exhibit poorer resistance to CO 2 11 and powder processing can therefore also have a signicant impact on stability as well as sinterability. [12][13][14] The high proton conductivity of Y-doped BaZrO 3 (BZY) has made it the state-of-the-art solid-state electrolyte in proton ceramic fuel cells (PCFCs) 3,15 and electrochemical membrane reactors. 16,17 BZY has a cubic crystal structure, 18 and BZY containing 10 mol% yttria (BZY10) has become one of the most studied compositions. 19 Ce-containing compositions, BaZr 1ÀxÀy Ce x Y y O 3 (BZCY), possess higher proton conductivity and improved sintering properties. 20 However, the chemical stability towards CO 2 -, H 2 O-, or H 2 S-containing atmospheres is reduced with increasing Ce-content. 21 The presence of carbonates can have a signicant effect on the transport properties of proton conductors, 22,23 as well as on the catalytic activity in the presence of hydrocarbons. 9 Nevertheless, a recent study by Duan et al. showed excellent coking resistance and sulfur tolerance of BaZrO 3 -based PCFCs. 24 Structural and chemical degradation has been reported even in yttria stabilized zirconia (YSZ) in carbon-rich fuel gases such as CH 4 , CO or syngas, inducing coking/ graphitization and (oxy)carbide formation [25][26][27][28] that led to irreversible changes in the transport properties. 27,28 Thus, a thorough understanding of potential degradation processes at the surface of BZY materials is highly desirable, and could lead to valuable insights into the surface chemistry and structure of this important proton-conducting electrolyte.
While bulk BZY has been shown to be chemically stable in 1 atm CO 2 above about 600 C in accordance with thermodynamic considerations, 29,30 the surface of the material may still exhibit some reactivity. Computational studies have recently shown that CO 2 exhibits strong chemisorption on BaZrO 3 surfaces by formation of carbonate species and that a carbonate overlayer of BaCO 3 can be thermodynamically preferred even in 400 ppm CO 2 above 100 C. 31 The reactivity between BZY and pure CO 2 has been demonstrated to result in the formation of BaCO 3 and deterioration of the mechanical properties at 550-750 C. 32 It was determined that the reaction with CO 2 resulted in Ba-deciency in the exposed surface rather than formation of ZrO 2 crystallites, 32 and the reaction could thereby reasonably be described according to where v 00 Ba and v O denote Ba-and O-vacancies, respectively, in Kröger-Vink notation. 33 Similar degradation has also been reported for Sr-based perovskites. 5,34 Alternative reaction pathways can be hypothesized based on further consideration of the phase stability and cation defect chemistry of BZY. In this respect, the coexistence of a Y-rich composition has been reported within single grains of Ba-decient BZY. [35][36][37] Han et al. showed by synchrotron X-ray diffraction that excess Y can occupy the A-site, i.e., Y Ba . 18 Another possible reaction mechanism between BZY and CO 2 towards BaCO 3 and Ba-decient BZY can therefore be expressed by relocation of Y from the Zr to the Ba-site and consumption of one BaZrO 3 formula unit according to In this study, we present a thorough investigation of the reactivity of dense BaZr 1Àx Y x O 3Àd (x ¼ 0-0.2) ceramics with CO 2 as the gas relevant for anodic fuel conversion electrochemical devices. The BZY samples were exposed to 1 atm CO 2 at 650 C for 10-1000 h, and the extent of BaCO 3 formation and changes in the BZY surface microstructure and chemical composition were investigated by using several techniques. Reactions (1) and (2) both result in Ba-decient BZ(Y) according to the change in A/B-site ratios, Ba/Zr and Ba/(Zr + Y), respectively. However, the reactions may be differentiated by the presence of the Y Ba species as well as the level of deterioration of the BZY grains due to the consumption of a BaZrO 3 formula unit in reaction (2). Furthermore, density functional theory (DFT) calculations were performed to evaluate relevant point defect energetics and to compare the thermodynamic feasibility of reactions (1) and (2).

Sample preparation
Dense BaZr 1Àx Y x O 3Àd (x ¼ 0-0.2) ceramics were prepared by the procedure described in two recent contributions. 32,38 Powders made by spray pyrolysis (CerPoTech, Norway) were calcined at 950 C for 12 h, ball milled and pressed into pellets by uniaxial pressing followed by cold isostatic pressing at 200 MPa. The pellets were sintered at 1600 C for 10 h in a powder bed consisting of 90 wt% BZY and 10 wt% BaCO 3 . The pellets were polished with SiC papers and diamond suspensions down to 1 mm in order to obtain surfaces with low roughness. The density of the materials was measured by the Archimedes method using isopropanol at room temperature. The polished disc samples were broken into several pieces and annealed under 1 atm of dry owing CO 2 at 650 C for total exposure times of 10, 20, 100 and 1000 h. The 10 and 20 h annealing experiments were performed in an alumina tube furnace, and the longer annealings of 100-1000 h were performed in a ProboStat measurement cell (NORECS, Norway) with a sample holder of alumina and platinum wires. The samples are identied by their yttrium content (BZ, BZY10 and BZY20 for 0, 10 and 20% Y doping, respectively) and their exposure time to CO 2 at 650 C.
X-ray powder diffraction (XRD) of the calcined powders and sintered materials was performed using a Bruker D8 Advance DaVinci diffractometer using CuK a1 radiation. Grazing incidence XRD (GIXRD) was performed using CuK a1 and parallel beam optics on pristine samples and samples exposed to CO 2 to characterize the formation of BaCO 3 on the surface. Rietveld renement of the XRD patterns for powders was carried out using TOPAS V4.1 soware using a cubic structure model (Pm 3m) for all the BZ(Y) materials. The microstructure of the materials was investigated with Hitachi S3400N and Hitachi FEG Zeiss Ultra scanning electron microscopes (SEM) on gold coated samples. The chemical composition of phases was investigated by energy dispersive spectroscopy (EDS) with Oxford Instruments AZtec Energy analysis soware.

Transmission electron microscopy
Cross-sectional transmission electron microscopy (TEM) samples were made using a focused ion beam on a dual-beam JEOL SEM. TEM was carried out using a FEI Titan G2 60-300 operated at 300 kV with a high brightness XFEG, probe corrector, Gatan imaging lter for electron energy loss spectroscopy (EELS) analysis and super-X EDS detector. EELS results were acquired with STEM spectral imaging, using dual EELS, acquiring both high loss and low loss, for energy referencing.

Secondary ion mass spectrometry
The distribution of the isotopes was determined by secondary ion mass spectrometry (SIMS) using a Cameca IMS 4f instrument. A primary ion beam of 10 keV Cs + was applied while qualitative optimization of sputtering and mass intensity was performed with a 12.5 keV O 2 + ion beam. The primary ion current varied between 30 and 60 nA, and the secondary ion intensity ranged between 5 and 6 nA. Depending on the secondary ion intensities, an area with a diameter of $30 mm in the center of a sputtered 100 Â 100 mm 2 area was gated for signal detection and analysis in order to avoid negative interference (edge and wall effects). An electron shower with an acceleration of 4.5 kV for $1.2 mA sample current was used to charge compensate the insulating samples. Dynamic transfer optics setting (DTOS) of 60% was used. The secondary molecular ion signals were measured as a function of cycles per second (cps) and sputter time.

X-ray photoelectron spectroscopy
X-ray photoelectron spectroscopy (XPS) was performed on both pristine and CO 2 -exposed BZY samples using a Kratos Axis Ultra DLD XPS system with an incident monochromated Al Ka X-ray beam. The working pressure in the XPS chamber at room temperature was below 5 Â 10 À9 mbar. High resolution scans were recorded for the Ba 3d and 4d, Zr 3d, Y 3d, O 1s and C 1s energy regions, with a step size of 0.1 eV. All of the obtained binding energies (BE) were calibrated based on the C 1s peak for aliphatic bonds from adventitious carbon, set to 284.8 eV BE. The resulting XPS data were analysed using CasaXPS soware. A Gaussian-Lorentzian line shape GL (30) was employed to t the XPS peaks aer background subtraction based on the Shirley algorithm. The peak separation and area ratio of the Y 3d 5/2 and the Y 3d 3/2 peaks were constrained to be 2.05 eV and 3 : 2, respectively.

Computational methods
Y-doped BaZrO 3 was modeled as a 3 Â 3 Â 3 supercell (135 atoms) containing 4 yttrium acceptors and 2 oxygen vacancies for charge neutrality, corresponding to a dopant concentration of approx. 15% (denoted as BZY15). Interactions between Y 0 Zr and v O were thereby taken into account for a specic conguration of the highly doped material. The binding energies between isolated defects were calculated based on the single fully ionized defects in 4 Â 4 Â 4 supercells. The BZY15 cell was constructed with a pair of Y 0 Zr on next neighbor sites and two Y 0 Zr on second next neighbor sites. Based on the calculated association energy between isolated Y 0 Zr and v O of À0.29 eV, the two chargecompensating v O were placed in the most stable neighboring sites to Y 0 Zr . The same procedure was followed for introducing additional barium Schottky defect pairs, or the introduction of two Y Ba and removal of the least stable v O , according to eqn (1) and (2), respectively. The tendency of the defects to segregate to the surface was investigated using 11-layer (0 0 1) BaO-terminated slabs with 3 Â 3 expansion perpendicular to the surface. Defect segregation energies were calculated as the total energy difference between charge neutral cells with the defects residing in the bulk and surface regions, respectively. The entropies of reactions (1) and (2) were estimated from the tabulated entropies of CO 2 and the defect-free solid phases. 39,40 The DFT calculations were performed with VASP 41 based on the projector-augmented wave (PAW) method 42 and the generalized gradient approximation (GGA) functional proposed by Perdew, Burke and Ernzerhof. 43 Geometric optimization of the Y-doped supercell was performed with a 500 eV plane-wave energy cut-off and the subsequent calculations of the Ba-decient cells were performed with xed lattice parameters to represent bulk continuation to the surface. A 2 Â 2 Â 2 Monkhorst-Pack scheme was used for k-point sampling of the supercells. 44 The ionic positions were relaxed until the residual forces were within 0.05 eVÅ À1 , while the self-consistent electronic optimization was performed with an energy convergence criterion of 10 À6 eV. Core level shis for Y 3d in the different congurations and clusters were obtained according to the initial state approximation as implemented in VASP.

Formation of BaCO 3
The extent of BaCO 3 formation on the surface of BZ and BZY samples aer annealing in CO 2 at 650 C for 10 to 1000 h was investigated by GIXRD and SEM. Fig. 1a and b show the GIXRD patterns of the pristine and annealed specimens. For both compositions, reections corresponding to BaCO 3 (Pmcn) appear aer annealing, and the amount of BaCO 3 increases with annealing time. On the other hand, there is no clear evidence of formation of ZrO 2 other than possibly small amounts for the BZY10 sample aer 1000 h. Based on the relative intensities of the main reections, the formation of BaCO 3 was more prominent in the Y-doped composition. As shown in the insets, there is a continuous shi of the BZ and BZY Bragg reections towards smaller lattice parameters with increasing annealing time and BaCO 3 formation. Fig. 1c and d show the surface of BZY10 before and aer exposure to CO 2 at 650 C for 10 h. Crystals of BaCO 3 were formed on the surface by exposure to CO 2 , resulting in a rough surface. The average size of the crystals was about 4 mm. SIMS analyses were performed on the BZ and BZY samples annealed in CO 2 at 650 C for 1000 h. The depth proles of the elements normalized to Zr intensity are given in Fig. 2c and d. According to these proles, the BZ sample exhibits Ba-deciency at the surface, while BZY10 and BZY20 seem to be quite stoichiometric. However, considering Y-segregation to the surface region shown in Fig. 2b, the normalization to Zr does not reveal a consistent A/B-ratio due to the change in all of the elements from the surface towards the bulk.

Cation depth proles
The carbon concentration proles in BZY10 are shown in Fig. 3. The amount and depth of carbon increase as a function of annealing time in CO 2 . Due to the low atomic number, the carbon intensities were quite low compared to the intensities of Ba, Y and Zr.

TEM investigation of the surface microstructure and phases
TEM was used to perform microstructure and phase characterization and in particular to investigate possible formation of 3850 | J. Mater. Chem. A, 2019, 7, 3848-3856 This journal is © The Royal Society of Chemistry 2019 secondary phases at the surface. Fig. 4 shows a surface crosssectional TEM image of BZY20 annealed for 1000 h in CO 2 at 650 C, and electron diffraction patterns from the surface region and 150 nm into the bulk. The grains corresponding to BaCO 3 are clearly visible on the surface. Other secondary phases like ZrO 2 or Y 2 O 3 were not observed.
The bulk (red) diffraction pattern shows the h100i zone axis of the BZY perovskite structure (Pm 3m) with a lattice parameter of 4.1Å. The diffraction pattern from the surface (blue) shows the same crystal structure, albeit with a slightly increased lattice parameter of 4.4Å, in addition to diffraction patterns from BaCO 3 and/or Ba(OH) 2 . The increase in the BZY lattice parameter might suggest a higher Y concentration in the top 150 nm. However, no difference in Y concentration was detected by EDS. One might speculate whether the apparent increase in C content towards the surface observed by EDS proling may reect interstitial carbon species and associated defects that expand the lattice. In this respect, DFT calculations have shown that carbon interstitials can be stabilized as carbonate species in either Ba-or Zr-vacancies in BaZrO 3 . 22 While their equilibrium concentrations were determined to be negligible in the acceptor doped material, the reaction may be facilitated by preexisting cation vacancies as has been shown for dissolution of interstitial Ni into BaZrO 3 from NiO. 45 In Fig. 5a, a magnied high angled annular dark eld (HAADF) STEM image shows that smaller pores of 5-20 nm were present inside the grains near the surface of BZY20 annealed for 1000 h in CO 2 at 650 C. This is also shown in Fig. 5b, in addition to structural dislocations. These nanoscopic voids exhibit high surface energies and cannot be expected to be stable aer sintering of the material, and it is reasonable to assume that they may have been formed during the annealing

XPS and computational analyses of chemical states
XPS was performed to characterize changes in the chemical state of Y in the BZY samples aer reaction with CO 2 to emulate possible formation of Y Ba according to reaction (2). The XPS spectra of Ba 4d and Y 3d are compared for BZY10 and BZY20 annealed in CO 2 for 100 h and 1000 h in Fig. 6. The bulk spectra were taken from the annealed samples (1000 h) fractured under vacuum, and they may have been affected by long-range diffusion processes during the annealing although the main changes in stoichiometry were observed within 10 mm from the surface based on the EDS and SIMS proles. The Ba 4d peak shows that at least two chemical states are present (Fig. 6a). The lower binding energy (BE) peak can be related to BZY while the higher BE peak is typical of BaCO 3 and Ba(OH) 2 . The intensity of the latter peak clearly increases with exposure time to CO 2 , well in agreement with the BaCO 3 formation on the surface of BZY10 samples as detected by GIXRD and EELS (Fig. S2 †). The Zr element is quite stable in CO 2 atmosphere, showing the unchanged position of Zr 3d peaks with increasing exposure time (Fig. S3 †). On the other hand, the Y 3d peak shows a complex structure that changes with increasing time of CO 2 exposure (Fig. 6b and c). Two spin-orbit pairs describing at least two different chemical states were required for a satisfactory tting of the Y 3d energy region. The spin-orbit pair at the highest BE exhibited a core level shi of about 2.5 eV relative to the pair at the lowest BE for BZY10. In the case of BZY20, the core level shi was quite similar for the bulk, while it decreased to about 1.6 eV aer exposure to CO 2 (Fig. 6c). The peak intensities for the chemical state at lower BE decreased with CO 2 exposure time in favor of the chemical state at higher BE (particularly evident for BZY20). The relative Y 3d core level shi of various Y defects and clusters was investigated computationally for comparison with results from XPS. Fig. 7 shows several congurations of Y 0 Zr (Fig. 7b-d) and Y Ba (Fig. 7e-g) with the calculated Y 3d core level shis (D) relative to the fully coordinated Y 0 Zr octahedra (Fig. 7a). The Y 3d core level shi for Y 0 Zr associated with a protonic defect, a nearest neighbor Y 0 Zr and/or an oxygen vacancy was 0.18-0.88 eV to higher binding energy relative to the isolated Y 0 Zr . In comparison, the core level shi for Y Ba was more signicant, 2.21-2.40 eV towards higher binding energy. Similar Y 3d core level shis were obtained for Y 0 Zr and Y Ba defects and clusters on the (0 0 1) BaO-terminated surface, i.e., there was no signicant difference between the core level shis for the bulk and surface.
Despite the rather large difference in size between Y and 12coordinated Ba, Y Ba attained a 6-coordinated conguration by   relaxing towards the edge (Fig. 7f) or face (Fig. 7g)  Zr and Y Ba species corresponds well with the separation of about 1.6-2.5 eV observed by XPS. In the case of dehydrated materials due to dry atmosphere, and a relatively low concentration of oxygen vacancies due to reaction (2), the calculated shi would be closer to 2.4 eV (Fig. 7a and e).

DFT calculations of reaction thermodynamics
DFT calculations were used to calculate and compare the enthalpy of reactions (1) and (2), while the reaction entropies were estimated from the tabulated entropies of CO 2 and the solid phases. The relaxed structures of the computational cells with an optimized conguration of additional defects are shown in Fig. 8. The BZY15 reference cell, shown in Fig. 8a, contains one oxygen vacancy in between two nearest-neighbor Y 0 Zr , i.e., (Y Zr -v O -Y Zr ) Â , and one oxygen vacancy associated with a single Y 0 Zr .
The cell in Fig. 8b contains a Schottky pair where the most favored site for v 00 Ba was adjacent to the (Y Zr -v O -Y Zr ) Â complex, with the oxygen vacancies associated with the remaining fully coordinated Y 0 Zr . The cell in Fig. 8c contains two Y Ba , both of which were found to be most stable adjacent to ( The total energy differences between the cells in Fig. 8b and c with respect to the reference cell in Fig. 8a were used to evaluate the enthalpies of reactions (1) and (2). Table 1 shows the calculated Gibbs energies for the reactions at 650 C based on cells containing different amounts of defects as indicated by the degree of Ba-deciency, x. While the obtained Gibbs energy of reaction (1) was positive for all the considered cells and congurations, it was signicantly lower and close to zero for reaction (2). This indicates that reaction (2) involving formation of BaCO 3 by accommodating Ba-deciency and Yenrichment in the bulk perovskite structureand hence consumption of BaZrO 3can occur to a signicant extent at 650 C in 1 atm CO 2 .
The tendency of yttrium to segregate to the surface was evaluated for BaO-terminated (0 0 1). The calculated segregation energy of Y Ba was À0.33 eV, indicating a considerable stabilization of the defect at the surface, shown in Fig. S4 (ESI †). The segregation energy for an additional Y 0 Zr to the Y Ba at the surface was À0.50 eV, resulting in an overall segregation energy of À0.82 eV for the associated defect pair from the bulk to the surface (see the ESI † for additional details).

Discussion
Observations by several techniques (SEM, TEM, and GIXRD) show that BaCO 3 crystallites form on the surfaces of BZ and BZY aer annealing in 1 atm CO 2 at 650 C. The amount of BaCO 3 increases with annealing time from 10 to 1000 h and it becomes larger for Y-doped compositions. The surface becomes Ba-decient based on SIMS and SEM-EDS analysis of BZ and BZY, respectively (it is difficult to assess Ba-stoichiometry in BZY by SIMS since there is no constant reference due to simultaneous changes in the Y and Zr contents). The Ba-deciency reaches several micrometers into the material, and it is likely that Ba is transported along grain boundaries. 46 Moreover, the corresponding carbon prole from SEM-EDS indicates BaCO 3 formation along grain boundaries and/or dissolution of carbon species into the perovskite. The BZY s-urfaces become enriched in Y as indicated by SIMS and XPS, and the Y-enrichment appears to extend several micrometers into the material, i.e., similar to the corresponding deciency in Ba. XPS shows two distinctly different chemical states of Y, which based on DFT calculations were identied as Y substituted on the A-site, Y Ba , and B-site, Y 0 Zr , respectively. Furthermore, the DFT calculations indicate that these Y-species tend to cluster, both in the bulk and at the surface. The resulting local enrichment in Y becomes structurally similar to Y 2 O 3 and may explain previous observations of Y-rich regions within BZY grains. [35][36][37] The change in stoichiometry near the surface is accompanied by a substantial increase in the lattice parameter of the remaining perovskite, as observed by TEM.
All in all, formation of BaCO 3 on BZY in CO 2 atmosphere at 650 C appears to mainly proceed according to reaction (2), in which Y 0 Zr ends up as a charge compensating donor, Y Ba , while a BaZrO 3 formula unit is consumed. The latter ts with the observation of small intragrain pores by TEM and the general absence of Y 2 O 3 or ZrO 2 secondary phases. In comparison to BZ, where BaCO 3 formation is lower, the Y-dopant thereby provides an extra compositional degree of freedom that facilitates BaCO 3 formation according to reaction (2). Now, the question is to what extent the formation of BaCO 3 by reaction (2) will proceed. Under the conditions used it is possible that it will proceed until complete decomposition to BaCO 3 and for instance Y-saturated ZrO 2 and Zr-saturated Y 2 O 3 . This process will take thousands of hours to complete, while larger BZY grains will probably retard the transport of Ba to the surface and hence also the decomposition and degradation. Another alternative is that it stops at equilibrium when a certain content of Y Ba is reached. The observation that there is a certain Ba-deciency and Y excess quite far into the materialinstead of complete decomposition at the very surfacemay support this alternative. These and other degradation paths could be studied by the reaction between for instance BaCO 3 and YSZ in 1 atm CO 2 .
In any case, while we have shown that BZY suffers from severe BaCO 3 formation in 1 atm CO 2 and probably degradation of its functional and thermomechanical properties, many applications of BZY electrolytes may not be signicantly affected. Lower CO 2 partial pressures or impurity levels of CO 2 should mitigate the degradation processes in hydrogen fuel cells and steam electrolyzers. However, such devices operating at steam pressures within the stability regime of Ba(OH) 2 may exhibit similar degradation phenomena and mechanisms as we have observed here. 47 A recent study on CO 2 and H 2 O coadsorption has shown preferential adsorption of H 2 O on BZY and indicated that the presence of H 2 O prevents the formation of a carbonate overlayer on BZY surfaces. 48 For applications with high CO 2 partial pressures due reforming of hydrocarbons, the presence of steam may therefore limit the BaCO 3 formation and the associated degradation processes. Further degradation studies should therefore take into account both CO 2 and H 2 O. Temperatures lower than 650 C are desirable for many applications of BZY, which worsens the thermodynamics of stability but slows down diffusion, so the effect of temperature may go both ways.

Conclusions
Over a time-scale of 1000 h, BaZr 1Àx Y x O 3Àd (x ¼ 0-0.2) ceramics exhibit severe degradation under 1 atm CO 2 at 650 C due to formation of BaCO 3 on the surfaces and possibly along the grain boundaries, accompanied by Ba-deciency in the subsurface regions. For the Y-doped compositions, the surface and subsurface become enriched in Y, which may be ascribed to a degradation reaction in which Y ends up on the Ba-site while BaZrO 3 formula units are consumed. Accordingly, the degradation proceeds faster with Y as an additional compositional degree of freedom and a more favourable degradation reaction. The local enrichment of Y both in the bulk and on the surface attains a structure similar to Y 2 O 3 . Electrochemical devices based on Y-doped BaZrO 3 electrolytes may be stable towards BaCO 3 formation under lower CO 2 partial pressures and/or in the presence of H 2 O, while further work is necessary to be conclusive.

Conflicts of interest
There are no conicts to declare.  (1) and (2) at 650 C showing the defects formed and consumed in the reactions, and the extent of Ba-deficiency (x) in the defective cells. The reaction enthalpies were obtained from DFT and the reaction entropies were estimated from the tabulated entropies of CO 2  and electrolysers" (FOXCET, 228355) conducted by SINTEF, University of Oslo and The Norwegian University of Science and Technology (NTNU), is gratefully acknowledged in addition to FME NCCS (257579). Computational resources were provided by the Norwegian Metacenter for Computational Science (NOTUR) under the project nn9259k.