domain: phys sde sdu

text:
Julien Joliat

Sylvain Picaud
email: sylvain.picaud@univ-fcomte.fr

Pál Jedlovszky

Competitive Adsorption of Trace Gases on Ice at Tropospheric Temperatures: A Grand Canonical Monte Carlo Simulation Study

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Gas/ice interfaces are ubiquitous in atmospheric environments, where they influence the local chemistry by acting as trace gas scavengers or as catalysts for heterogeneous surface reactions that may strongly differ from gas phase processes. [1][2][3][START_REF] Mcneill | Organics in Environmental Ices: Sources, Chemistry, and Impacts[END_REF] Thus, many laboratory studies have been devoted to the characterization of the trapping properties of various atmospherically relevant molecules, ranging from weakly (e.g., small organic compounds, such as formaldehyde) to very strongly (e.g., inorganic acids, such as HNO 3 ) interacting species on ice, at tropospheric temperatures. 2,3 Meanwhile, some of these experimental data have nicely been complemented by results coming from numerical simulations, based on the description of the gas/ice interactions at the molecular scale. [START_REF] Picaud | Molecular-Scale Simulations of Organic Compounds on Ice: Application to Atmospheric and Interstellar Sciences[END_REF] In these framework, Monte Carlo simulations performed in the grand canonical ensemble (GCMC) [START_REF] Allen | Computer Simulation of Liquids[END_REF] have proven to be one of the most suitable modeling approaches, [START_REF] Picaud | Molecular-Scale Simulations of Organic Compounds on Ice: Application to Atmospheric and Interstellar Sciences[END_REF][START_REF] Jedlovszky | Determination of the Adsorption Isotherms of Methanol on the Surface of Ice. An Experimental and Grand Canonical Monte Carlo Simulation Study[END_REF] because GCMC directly provides the adsorption isotherm, being one of the observables most commonly recorded in the experiments. Thus, maximum surface coverage as well as adsorption energy can be provided on both experimental and theoretical sides for comparison, which allows, for instance, assessment of the interaction potential models used in the simulations.

Laboratory studies of atmospherically relevant adsorption processes have most often been conducted with a single species, and so has also been done by the modeling approaches.

However, taking into account the large variety of organic species emitted in the atmosphere, it is also of fundamental importance to characterize how the competition for occupying similar adsorption sites may change the trapping efficiency of the ice surface when different species are simultaneously present. [START_REF] Sokolov | Competitive Adsorption of Atmospheric Trace Gases Onto Ice at 228 K: HNO 3 /HCl, 1-Butanol/Acetic acid and 1-Butanol/HCl[END_REF] Despite the expected importance of such co-adsorption processes, it is disappointing to note that, as far as we know, only two, relatively old such laboratory studies have been published so far. [START_REF] Sokolov | Competitive Adsorption of Atmospheric Trace Gases Onto Ice at 228 K: HNO 3 /HCl, 1-Butanol/Acetic acid and 1-Butanol/HCl[END_REF][START_REF] Hessberg | Interaction of Formic and Acetic Acid with Ice Surfaces Between 187 and 227 K. Investigation of Single Species-and Competitive Adsorption[END_REF] This could, nevertheless, be explained by the corresponding experimental constraints. However, the situation is even worse from the theoretical side, as no study has provided a detailed, molecular-level understanding of the simultaneous adsorption of several compounds on ice.

Again, the GCMC method is an ideal tool to tackle such a challenge, as it has proven to give a realistic description of the selective trapping of small volatile species in multi-guest clathrates hydrates, in various contexts relevant for astrophysical situations. [START_REF] Lectez | A 32-70 K Formation Temperature Range for the Ice Grains Agglomerated by Comet 67 P/Churyumov-Gerasimenko[END_REF][START_REF] Patt | A Grand Canonical Monte Carlo Study of the N 2 , CO, and Mixed N 2 -CO Clathrate Hydrates[END_REF][START_REF] Petuya | Molecular Selectivity of CO-N 2 Mixed Hydrates: Raman Spectroscopy and GCMC Studies[END_REF][START_REF] Patt | Molecular Selectivity of CH 4 -C 2 H 6 Mixed Hydrates: A GCMC Study[END_REF] Thus, based on the experience we have acquired by using the GCMC method to characterize, on one hand, the adsorption of one single organic species on ice [START_REF] Picaud | Molecular-Scale Simulations of Organic Compounds on Ice: Application to Atmospheric and Interstellar Sciences[END_REF] and, on the other hand, the competition between different species to occupy the cages of clathrate hydrates, [START_REF] Lectez | A 32-70 K Formation Temperature Range for the Ice Grains Agglomerated by Comet 67 P/Churyumov-Gerasimenko[END_REF][START_REF] Patt | A Grand Canonical Monte Carlo Study of the N 2 , CO, and Mixed N 2 -CO Clathrate Hydrates[END_REF][START_REF] Petuya | Molecular Selectivity of CO-N 2 Mixed Hydrates: Raman Spectroscopy and GCMC Studies[END_REF][START_REF] Patt | Molecular Selectivity of CH 4 -C 2 H 6 Mixed Hydrates: A GCMC Study[END_REF] here we present the results of the first molecular scale simulation study of a competitive adsorption process on ice under tropospheric conditions.

As a first application, we have chosen to investigate the case of mixtures of 1-butanol and acetic acid molecules, i.e., two small, polar, oxygenated organics that are abundant in the atmosphere. [START_REF] Cavalli | FTIR Kinetic, Product, and Modeling Study of the OH-Initiated Oxidation of 1-Butanol in Air[END_REF][START_REF] Andersen | Atmospheric Chemistry of i-Butanol[END_REF][START_REF] Khan | Investigating the Tropospheric Chemistry of Acetic Acid Using the Global 3-D Chemistry Transport Model, STOCHEM-CRI[END_REF] The choice of these adsorbates is also dictated by the fact that their co-adsorption on ice has previously been investigated in coated-wall flow tube experiments, [START_REF] Sokolov | Competitive Adsorption of Atmospheric Trace Gases Onto Ice at 228 K: HNO 3 /HCl, 1-Butanol/Acetic acid and 1-Butanol/HCl[END_REF] opening thus the possibility of comparing the present results to the experimental data.

Computational Details

Water molecules have been described by the TIP4P/Ice model, [START_REF] Abascal | A Potential Model for the Study of Ices and Amorphous Water: TIP4P/Ice[END_REF][START_REF] Conde | High Precision Determination of the Melting Points of Water TIP4P/2005 and Water TIP4P/Ice Models by the Direct Coexistence Technique[END_REF] whereas the 1-butanol and the acetic acid molecules have been modeled with the flexible AUA4 potential. [START_REF] Ferrando | Transferable Force Field for Alcohols and Polyalcohols[END_REF][START_REF] Ferrando | Prediction of Phase Equilibrium and Hydration Free Energy of Carboxylic Acids by Monte Carlo Simulations[END_REF] Note that previous simulation studies have shown the relevance of this combination of interaction potential models to accurately describe the behavior of the 1-butanol molecules on ice. [START_REF] Joliat | Adsorption of C2-C5 Alcohols on Ice: A Grand Canonical Monte Carlo Simulation Study[END_REF] Both the TIP4P/Ice and the AUA4 models are defined as sums of pairwise dispersion-repulsion and electrostatic terms, represented by Lennard-Jones (LJ) and Coulomb interactions. Heteronuclear interactions have been calculated in the same way, using the cross LJ parameters determined according to the standard Lorentz-Berthelot combination rules. [START_REF] Frenkel | Understanding Molecular Simulation: From Algorithms to Applications[END_REF] Analytical tail correction [START_REF] Allen | Computer Simulation of Liquids[END_REF] has been applied for the LJ contributions to the interaction potential, while the long-range part of the electrostatic interaction has been accounted for by means of the Ewald summation method. [START_REF] Frenkel | Understanding Molecular Simulation: From Algorithms to Applications[END_REF] Ice grains have been simulated by a slab made of 2880 water molecules, arranged in 18 solid layers stacked along the 0001 crystallographic axis (defining the z axis of the system). This ice slab has been placed in a rectangular simulation box of the dimensions of L x = 35.926 Å, L y = 38.891 Å, and L z = 100 Å, in such a way that the slab exhibits two interfaces in contact with the gas phase, consisting of the 1-butanol and acetic acid molecules, along the z axis. Periodic boundary conditions have been applied along all the three axes and a cut-off radius equal to half the smallest size of the simulation box has been applied for the calculation of the interactions. Thus, two molecules pertaining to the two opposite interfaces have not interacted with each other directly. The simulations have been performed using the GIBBS software, which is a general Monte Carlo simulation code, developed jointly by three French Institutions (Université Paris Sud, CNRS and IFPEN), [START_REF] Ungerer | Applications of Molecular Simulation in the Oil and Gas Industry: Monte Carlo Methods[END_REF] and which has been used here under a Licence freely provided.

To characterize the adsorption properties of a gas phase containing a mixture of 1-butanol and acetic acid molecules, GCMC runs have been performed, allowing the number of the adsorbate molecules to fluctuate in the simulation box with the value of the total chemical potential µ. For each µ value considered, the number of molecules of each species that are in contact with the ice surface has been determined. This number of adsorbed molecules as a function of µ corresponds, by definition, to the adsorption isotherm for a given composition of the gas phase and a given temperature. For comparisons with experimental data, the results of such simulations have, however, been expressed in terms of pressure instead of chemical potential, using a conversion procedure consistent with the interaction potential models, for the selected compositions of the gas phase, as in our previous works [START_REF] Patt | Molecular Selectivity of CH 4 -C 2 H 6 Mixed Hydrates: A GCMC Study[END_REF][START_REF] Joliat | Adsorption of C2-C5 Alcohols on Ice: A Grand Canonical Monte Carlo Simulation Study[END_REF][START_REF] Joliat | Adsorption of Organic Compounds at the Surface of Enceladus' Ice Grains. A Grand Canonical Monte Carlo Simulation Study[END_REF] (see also the Supporting Information).

The GCMC simulations have started with a period of equilibration of typically 2-4×10 8 Monte Carlo steps. Once the energy of the system and the number of molecules in the basic box started to fluctuate around their equilibrium values, a production run of 2 × 10 8 Monte Carlo steps has been performed. The 1-butanol and acetic acid molecules have been flexible, while the water molecules have been kept rigid, in accordance with the parametrization of the TIP4P/ice model. [START_REF] Abascal | A Potential Model for the Study of Ices and Amorphous Water: TIP4P/Ice[END_REF][START_REF] Conde | High Precision Determination of the Melting Points of Water TIP4P/2005 and Water TIP4P/Ice Models by the Direct Coexistence Technique[END_REF] Thus, in the simulations, the water molecules have only been subject to translational and rotational moves, performed with equal probabilities, while for the 1-butanol and acetic acid molecules, insertion (20 %), deletion (20 %), translation (15 %), rotation (15 %), configurational-biased regrowth (i.e., change of the internal molecular configuration, 15%) and identity swap (i.e., one molecule of a given type is replaced by a molecule of the other type, at the same position, 15%) moves have been attempted.

Results and Discussion

Our previous works showed that the simulated adsorption isotherms of small alcohols on ice compare very well with experimental data, provided that a small (i.e., about 5 K) temperature shift is considered in the calculations. [START_REF] Jedlovszky | Determination of the Adsorption Isotherms of Methanol on the Surface of Ice. An Experimental and Grand Canonical Monte Carlo Simulation Study[END_REF][START_REF] Joliat | Adsorption of C2-C5 Alcohols on Ice: A Grand Canonical Monte Carlo Simulation Study[END_REF] For this reason, the present simulations have not only been performed at the temperature of the experiments of 228 K, but also at 233 K. However, only the results obtained at this slightly higher temperature are detailed here, because they indeed show an overall better agreement with the experimental data, in particular, when looking at the low pressure part of the adsorption isotherms, as discussed below. Note that we have not studied the adsorption process at higher temperatures, the main goal being, in this first simulation study of the competitive adsorption process, to stay as close to the experimental conditions as possible, for the purposes of validating our GCMC approach.

Adsorption isotherms

The simulated adsorption isotherms at 233 K are shown in Figure 1 for four various compositions of the gas mixture, corresponding to the gas phase 1-butanol:acetic acid molar ratios of 100:0, 62:38, 25:75, and 0:100. These compositions have been chosen in accordance with the experimental conditions, considering both species separately, and also two mixtures, charac-terized by gas partial pressure ratio of P 1-butanol = (1.65 ± 0.15) × P acetic acid (corresponding to the molar ratio of 62:38) and P acetic acid = 3.1 × P 1-butanol (i.e., molar ratio of 25:75). [START_REF] Sokolov | Competitive Adsorption of Atmospheric Trace Gases Onto Ice at 228 K: HNO 3 /HCl, 1-Butanol/Acetic acid and 1-Butanol/HCl[END_REF] Note that the experimental data at 228 K, 8 extracted from the original publication using the WebPlotDigitizer software, 25 are also reported on Figure 1, for comparisons.

At first, we discuss the adsorption of 1-butanol and acetic acid in the respective neat systems. As seen from Figure 1 (top row), the adsorption behavior of the two molecules is rather similar to each other in the low pressure range, where a rapid increase of their surface concentration, Γ, corresponding to the building up of the adsorption layer at the surface, is observed. It should, however, be noted that the adsorption process starts at much lower pressures for acetic acid than for 1-butanol and, as a consequence, in the low pressure range, up to about 5 × 10 -3 Pa, the amount of acetic acid molecules trapped at the ice surface ice is considerably larger than that of the alcohol molecules. The increase of the pressure above 2 × 10 -2 Pa only leads to a small variation of Γ for 1-butanol, which is a strong indication of the formation of a stable monolayer. [START_REF] Joliat | Adsorption of C2-C5 Alcohols on Ice: A Grand Canonical Monte Carlo Simulation Study[END_REF] On the other hand, it leads to a continuous increase of Γ for acetic acid due to the adsorption of an increasing amount of molecules at the ice surface up to the completion of the layer, occurring just before the condensation occurs (this condensation is indicated by the sudden vertical jump of the isotherm). Further, in the intermediate pressure range of about 10 -2 -1 Pa, the number of adsorbed alcohol molecules is always larger than that of acetic acid.

The simulated isotherms, calculated at 233 K, agree nicely with the experimental ones, measured at 228 K, [START_REF] Sokolov | Adsorption to Ice of n-Alchohols (Ethanol to 1-Hexanol), Acetic Acid and Hexanal[END_REF] especially in the low pressure range and for 1-butanol (see Figure 1a). At larger coverage, the agreement is less satisfactory for these alcohol molecules, although the large dispersion of the experimental data makes the comparison with the simulation results not so easy. As previously discussed, [START_REF] Joliat | Adsorption of C2-C5 Alcohols on Ice: A Grand Canonical Monte Carlo Simulation Study[END_REF] this slight discrepancy could come from an underestimation of the lateral interactions between alcohol molecules with the AUA4 potential model, especially when long aliphatic chains are involved. Also, our simulations seem to somewhat overestimate the amount of adsorbed acetic acid molecules in the whole pressure range (Fig-ure 1b), a feature that could come from a slight overestimation of the interaction between the COOH group and the ice surface, when using the AUA4 model. Nevertheless, it is very satisfying to see that the subtle balance between adsorbate-ice and adsorbate-adsorbate interactions, which governs the adsorption process, can be reproduced quite well by using a simple combination of empirical models originally parametrized for pure compounds only.

Properties of the adsorbed molecules

As in our previous work, [START_REF] Joliat | Adsorption of C2-C5 Alcohols on Ice: A Grand Canonical Monte Carlo Simulation Study[END_REF] the analysis of these results has been performed in terms of number density profile, adsorption energy, and orientation of the adsorbed molecules (some of the corresponding curves are provided as Supplementary Information, in addition to those given below). This analysis allows us to draw the following conclusions on the individual adsorption of the two species at 233 K.

At the low surface coverage corresponding to the pressure value of 0.003 Pa, the distribution of the molecular C-C axis orientations (P (cosθ 1 )) for the 1-butanol molecules (see Figure 2a) is characterized by a broad peak, the maximum of which corresponds to the value of θ 1 ∼ 105 • . This indicates that the 1-butanol molecules are adsorbed with their molecular axis being slightly (by about 25 • ) tilted with respect to the ice surface, the OH group being closer to the water molecules than the methyl group. This configuration is very similar to that obtained at 228 K. [START_REF] Joliat | Adsorption of C2-C5 Alcohols on Ice: A Grand Canonical Monte Carlo Simulation Study[END_REF] The distribution of the interaction energy between one adsorbed 1-butanol molecule and the ice surface, U alc-w , exhibits only one, broad, peak, the maximum of which corresponds to about -67 kJ.mol -1 (Figure 2c), in a nice agreement with the experimental value of the adsorption energy of -67.8 ± 3.8 kJ.mol -1 , measured at very low coverage. [START_REF] Sokolov | Adsorption to Ice of n-Alchohols (Ethanol to 1-Hexanol), Acetic Acid and Hexanal[END_REF] Meanwhile, the distribution of the interaction energy of one adsorbed 1-butanol molecule with the rest of the adsorption layer, U alc-alc , also exhibits one single peak, the maximum of which corresponding to a small value close to zero energy, indicating that the lateral interaction between the adsorbed alcohol molecules is very weak in this low pressure range (Figure 2e). By contrast, at pressures corresponding to saturation coverage of the ice surface, most of the adsorbed molecules tend to move their apolar tail away from, and pointing their polar head towards the ice surface, their molecular axis now being almost perpendicular to the surface, as indicated by the main peak of the distribution P (cosθ 1 ) around the value of cosθ 1 = -0.9 (Figure 2b). For simplicity of the following discussions, this orientation will be referred as orientation Alc A . In this new orientation, the U alc-w and U alc-alc distributions of the interaction energies have their maximum around -44 kJ.mol -1 (Figure 2d) and -57 kJ.mol -1 (Figure 2f), respectively, indicating a strong increase of the lateral interactions between the adsorbed 1-butanol molecules at the expense of their interaction with ice. It is, however, worth noting that a few adsorbed 1-butanol molecules may also stay almost flat on the surface at this high coverage, as shown by the small peak around cosθ 1 = 0 in Figure 2b (below, this orientation will be referred to as orientation Alc B ). As an illustration, snapshots of the 1-butanol/ice system are given on Figure 3 at low (Figure 3a) and saturation (Figure 3b) coverages of the ice surface, in which two molecules exhibiting orientations Alc A and Alc B are indicated.

When considering the acetic acid molecules, the analysis of the simulation results in the low pressure range, when only a few molecules are adsorbed at the ice surface, indicates that all these molecules adopt almost the same alignment, as indicated by the large peak in P (cosθ 1 ) between -0.85 and -0.65 (Figure 4a). In this configuration (referred to as orientation Ac α ), the C-C axis of the acetic acid molecule is tilted by about 30-50 • from the surface normal, the carboxyl group pointing to the water molecules. The corresponding distribution of the interaction energy between one adsorbed acetic acid molecule and the ice phase, U ac-w , is characterized by one single peak, around -73 kJ.mol -1 (Figure 4c). In this situation, the lateral interaction between adsorbed acetic acid molecules is vanishingly small (Figure 4e), indicating that the adsorbed molecules stay isolated from each other.

Note that the obtained value of U ac-w agrees well with the experimental value of -73.1 ± 11.7 kJ.mol -1 . [START_REF] Sokolov | Adsorption to Ice of n-Alchohols (Ethanol to 1-Hexanol), Acetic Acid and Hexanal[END_REF] These conclusions remain valid at larger pressures, in the region where the adsorption isotherm slope changes, as indicated by the analysis of the simulation data performed at the pressure of 0.02 Pa. In particular, there is still only one single peak in the P (cosθ 1 ) and P (U ac-w ) distributions, which can be unambiguously related to configuration Ac α .

On the contrary, the analysis of the results obtained at higher pressures, just below the point of condensation, indicates the preference for two different orientations, as evidenced from the two peaks in the P (cosθ 1 ) distribution of the molecular C-C axis orientations (Figure 4b). It should also be noted that, consistently with this picture, the corresponding energy distributions exhibit two separate peaks (Figures 4d and4f). The first of these orientations is, in fact, very similar to Ac α , observed also at low coverage, while in the second orientation (referred to as orientation Ac β ) the molecular C-C axis lies more or less parallel to the ice surface. This new orientation, Ac β , can thus be associated with the second peak of the U ac-w distribution, occurring at -38 kJ.mol -1 , while the first peak at -68 kJ.mol -1 is still associated with Ac α . The distribution of the lateral interactions, U ac-ac , has its maxima at -22 kJ.mol -1 (orientation Ac α ) and -45 kJ.mol -1 (orientation Ac β ), indicating that these two orientations differ not only in the interaction with the water molecules, but also in that with the surrounding adsorbed molecules. These results can thus be related to the existence of two preferred adsorption configurations at the ice surface for the acetic acid molecules, as previously evidenced for formic acid molecules. [START_REF] Jedlovszky | Adsorption Isotherm of Formic Acid on the Surface of Ice, as Seen From Experiments and Grand Canonical Monte Carlo Simulation[END_REF] The first one, at which acetic acid molecules adopt orientations Ac α , is preferentially evidenced at low surface coverage, i.e., at the low pressure range (Figure 5a). At higher coverages, additional adsorption configuration of a different type is also adopted by the acetic acid molecules, characterized by the orientation Ac β (Figure 5b). When the gas phase contains a mixture of 1-butanol and acetic acid molecules, the corresponding simulated adsorption isotherms (blue and green curves in Figure 1) overlap with those of the non-mixed species in the very low pressure range, indicating that the adsorbed molecules of one species are not significantly influenced by those of the other species, when the total number of adsorbed molecules remains quite weak (low coverage of the ice surface). This is confirmed, for instance, by the analysis of the orientational distributions P (cosθ 1 ) (Figure 6a and 6b) which do not show any strong difference at different gas phase compositions.

The situation is clearly different when the pressure increases. Indeed, the adsorption isotherms for the mixtures all show a plateau, for both molecules at higher pressures (blue and green curves in Figure 1), indicating that, before the condensation occurs, one molecular adsorbed layer is stabilized at the surface of ice in a broad range of pressures. [START_REF] Picaud | Molecular-Scale Simulations of Organic Compounds on Ice: Application to Atmospheric and Interstellar Sciences[END_REF] Moreover, although this monolayer contains co-adsorbed 1-butanol and acetic acid molecules, it is characterized by a significantly lower number of both species than in the non-mixed situation at similar pressures. This shows that the 1-butanol and acetic acid molecules indeed compete for being adsorbed at the surface.

It is worth noting that the simulated behavior of the mixtures agrees, in general, well with the experimental results (grey symbols on Figures 1e and1f). This agreement is even quantitative when the gas phase contains 25 % of 1-butanol, as evidenced in Figure 1e. On the other hand, the agreement is only qualitative when the gas mixture contains 38 % of acetic acid, as here the simulation results overestimate the amount of adsorbed acetic acid molecules, as already observed also in the non-mixed case. In any case, when the two species can simultaneously be trapped at the ice surface, the maximum number of both molecules that can be adsorbed at the surface is considerably less than in the corresponding non-mixed situations, in accordance with the experimental observations.

To understand the situation more deeply, the distributions of molecular orientations at the ice surface have thus been analyzed at saturation coverage (corresponding to the pressure of 1 Pa), just before the occurrence of the condensation. The corresponding curves are shown in the bottom row of Figure 6, for the mixed (blue and green curves) and non-mixed (red curves) systems (given for reference). These distributions clearly show that, within the mixed adsorbed layer, preference for orientation Ac β nearly vanishes and, as a consequence, acetic acid molecules preferentially adopt orientation Ac α . The situation is somewhat more balanced for 1-butanol, as both orientations Alc A and Alc B are still evidenced when the alcohol molecules are co-adsorbed with acetic acid. However, it is clear that the higher the concentration of 1-butanol is at the ice surface, the weaker is their preference for orientation Alc B . In addition, the surface distribution of the XY positions of the carbon atom to which the alcoholic or acidic functional group is attached has been calculated for both species in 10,000 equilibrium sample configurations, and superimposed in the same graph, shown in Figure 7, for both gas mixtures considered. First, as evidenced by the various spots on these spatial distributions, the entire surface is covered by the adsorbed molecules, a feature which could be related to the plateau observed in the corresponding adsorption isotherms (Figures 1c andd). Further, 1-butanol and acetic acid molecules clearly prefer to occupy well-defined adsorption positions on the ice surface. It it also seen that the two species are not fully mixed at the ice surface, instead they tend to form well-separated self-aggregates of the same species. Snapshots taken from the simulations are given in Figure 8, to illustrate the orientational configurations adopted by the 1-butanol and acetic acid molecules in the mixed adsorption layers at the ice surface, at low and saturation coverages.

Selectivity of the adsorption process

The spatial distributions shown in Figure 7 also suggest that the composition of the adsorption layer is different from that of the gas phase, due to the competition in the adsorption process between the 1-butanol and acetic acid molecules. This difference can be more quantitatively characterized by calculating the selectivity of the ice surface for 1-butanol with respect to acetic acid, defined as [START_REF] Petuya | Molecular Selectivity of CO-N 2 Mixed Hydrates: Raman Spectroscopy and GCMC Studies[END_REF][START_REF] Guo | Competitive Adsorption of Methanol-Acetone on Surface Functionalization (-COOH,-OH,-NH 2 , and-SO 3 H): Grand Canonical Monte Carlo and Density Functional Theory Simulations[END_REF]

α C4H10O/C2H4O2 = N C4H10O /N C2H4O2 y C4H10O /y C2H4O2 (1)
where N j represents the molar fraction of species j (j = C 4 H 10 O or C 2 H 4 O 2 ) in the adsorption layer, while y j is the corresponding molar fraction in the gas phase. Thus, following the definition given by Eq. 1, values of α C4H10O/C2H4O2 larger than 1 indicate that the relative proportion of 1-butanol with respect to acetic acid is larger at the ice surface than in the gas phase, i.e., the adsorption process is selective for 1-butanol. By contrast, α C4H10O/C2H4O2

values smaller than 1 indicate that the adsorption process is selective for acetic acid.

As shown in Figure 9, the adsorption process remains always selective in favor of the acetic acid molecules (i.e., α C4H10O/C2H4O2 < 1) when 1-butanol molecules are in majority (i.e., 62 %) in the gas phase (green symbols). Nevertheless, it should be noted that the fraction of 1-butanol at the surface increases with pressure up to saturation, as indicated by the increase of α C4H10O/C2H4O2 up to a more or less constant value. This behavior is even more pronounced for the other mixture considered here (25 % of 1-butanol). Indeed, in this case, the value of α C4H10O/C2H4O2 remains lower than 1 only up to a threshold pressure, above which the adsorption process, initially being selective for acetic acid, becomes progressively more favorable to the alcohol molecules, as indicated by the selectivity values larger than 1 at higher pressures. Note that the threshold pressure at which the selectivity becomes larger than 1 roughly corresponds to the pressure at which the adsorption isotherms start to exhibit a plateau (i.e., when a more or less saturated mixed monolayer is formed at the surface).

Conclusions

Summarizing, the results of this first GCMC simulations devoted to the characterization of competitive adsorption clearly show that the trapping of acetic acid molecules is favored at the surface of ice with respect to that of 1-butanol at low pressures, irrespective of the gas phase composition, in a clear accordance with earlier experimental findings. [START_REF] Sokolov | Competitive Adsorption of Atmospheric Trace Gases Onto Ice at 228 K: HNO 3 /HCl, 1-Butanol/Acetic acid and 1-Butanol/HCl[END_REF] By contrast, at higher pressures, more 1-butanol molecules can be adsorbed on the expense of the adsorption of the acetic acid molecules, preventing the acetic acid molecules from adopting orientation Ac β . Moreover, the adsorption process leads to a gradual increase of the 1-butanol mole fraction in the adsorption layer with respect to that in the mixed gas phase, as shown by the increasing value of the ice selectivity for 1-butanol. Similar conclusions are obtained at the two temperatures considered, i.e., 228 K and 233 K. The present results clearly point out that interactions between co-adsorbing gases can play a significant role in their trapping on atmospheric ice particles, and modify their relative molar fractions, especially in the saturated portion of the adsorption isotherm, as already observed in the experiments.

Further, they underline that competitive adsorption effects cannot be disregarded when considering real atmospheric situations. The present study also shows that the GCMC method is a powerful and versatile tool for characterizing this competition, provided that accurate interaction potential models are available.

Figure 1: Average number of 1-butanol (circles) and acetic acid (triangles) molecules adsorbed on ice as a function of the pressure (note that the errors bars are smaller than the symbols), as calculated from the present GCMC simulations at 233 K, for various proportions of these two species in the gas phase. Top row: Red curves represent the results obtained for the neat (non mixed) systems. Middle row: Results obtained for the gas phase 1-butanol:acetic acid molar ratios of 25:75 (blue curves) and of 62:38 (green curves). Bottom row : Comparison between results obtained for the non-mixed and the mixed systems, for 1-butanol (left) and acetic acid (right) molecules. Lines connecting the points are only guides to the eye. Comparisons with the experimental data 8 are shown as grey symbols. Note that for the mixtures, these experimental data are available for only one species in each case, i.e., 1butanol for the gas phase 1-butanol:acetic acid molar ratios of 25:75, and acetic acid for the gas phase 1-butanol:acetic acid molar ratios of 62:38. For reference, the outer tail of the number density profile of the water centers of mass is also given (dashed blue line). In addition, snapshots issued from the simulations at saturation coverages are given, to illustrate the adsorption geometry when the molecules form a monolayer at the ice surface. In these snapshots, red, white, and black (gray) circles represent oxygen, hydrogen and carbon (bottom carbon) atoms, respectively. Note that hydrogen atoms of the aliphatic chains are not explicitly included in the united atom AUA4 potential model and, thus, they are not shown.

Figure 2 :

 2 Figure2: Top row : Distributions of the molecular orientations of the adsorbed 1-butanol molecules for the neat (non mixed) system. Insets show the definition of the θ 1 angle and, also, examples of the preferential orientations (Alc A and Alc B ) adopted by the 1-butanol molecules at the ice surface (see text for more explanation). In these representations of the molecules, red, white, and black (gray) circles represent oxygen, hydrogen and carbon atoms, respectively. Note that hydrogen atoms of the aliphatic chains are not shown for clarity. Middle and bottom rows : Distribution of the interaction energy of an adsorbed 1-butanol molecule with the ice phase P (U ads-w ) and with the other adsorbed molecules P (U ads-ads ). All these analyses have been performed at 233 K and at pressure values of 0.003 Pa (low coverage, left hand side of the figure) and 1 Pa (saturation coverage, right hand side).

Figure 3 :

 3 Figure 3: Snapshots taken from the simulations at (a) low and (b) saturation coveragesfor the 1-butanol molecules adsorbed on ice, at 233 K. Two molecules exhibiting the typical orientations Alc A and Alc B adopted by the 1-butanol species at the saturation coverage of the ice surface are indicated in blue, as an illustration. In these snapshots of the molecules, red, white, and black (gray) circles represent oxygen, hydrogen and carbon atoms, respectively. Note that hydrogen atoms of the aliphatic chains are not shown for clarity. The ice crystal has been drawn as a continuous grey material.

Figure 4 :

 4 Figure4: Top row : Distributions of the molecular orientations of the adsorbed acetic acid molecules for the neat (non mixed) system. Insets show the definition of the θ 1 angle and, also, examples of the preferential orientations (Ac α and Ac β ) adopted by the acetic acid molecules at the ice surface (see text for more explanation). In these representations of the molecules, red, white, and black (gray) circles represent oxygen, hydrogen and carbon atoms, respectively. Note that hydrogen atoms of the aliphatic chains are not shown for clarity. Middle and bottom rows : Distribution of the interaction energy of an adsorbed acetic acid molecule with the ice phase P (U ads-w ) and with the other adsorbed molecules P (U ads-ads ). All these analyses have been performed at 233 K and at pressure values of 0.0005 Pa (low coverage, left hand side of the figure) and 1 Pa (saturation coverage, right hand side).

Figure 5 :

 5 Figure 5: Snapshots taken from the simulations at (a) low and (b) saturation coverages for the acetic acid molecules adsorbed on ice, at 233 K. Two molecules exhibiting the typical orientations Ac α and Ac β adopted by the acetic acid species at the saturation coverage of the ice surface are indicated in blue, as an illustration. In these snapshots of the molecules, red, white, and black (gray) circles represent oxygen, hydrogen and carbon atoms, respectively. Note that hydrogen atoms of the aliphatic chains are not shown for clarity. The ice crystal has been drawn as a continuous grey material.

Figure 6 :

 6 Figure 6: Distributions of the molecular orientations of the adsorbed 1-butanol (left) and acetic acid (right) molecules, for the neat (i.e., non mixed) systems (red curves) and for 1-butanol:acetic acid gas phase molar ratios of 62:38 (green curves) and 25:75 (blue curves).Top and bottom rows correspond to low and saturation coverages of the ice surface by the adsorbate, respectively. Insets show the definition of the θ 1 angle. The preferential orientations for the 1-butanol (Alc A and Alc B ), and for the acetic acid (Ac α and Ac β ) molecules are also indicated (see text for more explanation). In these representations of the molecules, red, white, and black (gray) circles represent oxygen, hydrogen and carbon atoms, respectively. Note that hydrogen atoms of the aliphatic chains are not shown for clarity.

Figure 7 :

 7 Figure 7: Surface distribution of the XY positions of the carbon atom of the 1-butanol (purple dots) and acetic acid (blue dots) molecules to which the functional group is attached, at saturation coverages of the ice surface, and for the two compositions of the mixed gas phase investigated here.

Figure 8 :

 8 Figure8: Snapshots taken from the simulations of mixtures of 1-butanol and acetic acid molecules adsorbed on ice, at 233 K and at low (left hand side) and saturation (right hand side) coverages. Typical orientations adopted by the adsorbed species at the ice surface (see text) are indicated in blue, as an illustration. In these snapshots of the molecules, red, white, and black (gray) circles represent oxygen, hydrogen and carbon atoms, respectively. Note that hydrogen atoms of the aliphatic chains are not shown for clarity. The ice crystal has been drawn as a continuous grey material.

Figure 9 :Figure S1 :

 9S1 Figure9: Selectivity of the ice surface as a function of the total pressure, as defined by Eq. 1. Results have been calculated at 228 K (open symbols) and 233 K (full symbols), for the gas phase 1-butanol:acetic acid molar ratio of 62:38 (green symbols) and 25:75 (blue symbols). The limit between acetic acid and 1-butanol selective ice surfaces, corresponding to the selectivity value of 1, is indicated by the horizontal dashed line.

Figure S2 :Figure S3 :

 S2S3 FigureS2: Distribution of the molecular angle θ 2 formed between the CO vector and the normal z to the ice surface (as illustrated by the inserts), for 1-butanol (left) and acetic acid (right) molecules adsorbed on ice at low (top row) and saturation (bottom row) coverages of the ice surface. The results for the non-mixed systems are given in red, whereas blue and green curves represent the results obtained for the 25:75 and 62:38 mixed systems, respectively.

Acknowledgements

The Region Bourgogne Franche-Comté is gratefully acknowledged for its financial support through the grant PASCOA 2020-0052. Computations have been performed using resources of the Mésocentre de calcul de Franche-Comté. PJ acknowledges financial support from the Hungarian NKFIH Foundation under Project No. 134596. The license to use the GIBBS software has been granted under the contractual letter ref. IFPEN n°2020-0744.

Conflicts of interest

There are no conflicts to declare.

Supporting Information

Includes information on the NPT simulations performed for the conversion between chemical potentials and pressures (Widom procedure, see also Ref. [START_REF] Patt | Molecular Selectivity of CH 4 -C 2 H 6 Mixed Hydrates: A GCMC Study[END_REF][START_REF] Joliat | Adsorption of C2-C5 Alcohols on Ice: A Grand Canonical Monte Carlo Simulation Study[END_REF][START_REF] Joliat | Adsorption of Organic Compounds at the Surface of Enceladus' Ice Grains. A Grand Canonical Monte Carlo Simulation Study[END_REF] ), and additional results on the structural and orientational analyses performed to better understand the behavior of the molecules adsorbed at the surface of ice.

Supporting information

Competitive Adsorption of Trace Gases on Ice at Tropospheric Temperatures : A Grand Canonical Monte Carlo Simulation Study

Julien Joliat † , Sylvain Picaud † , Pál Jedlovszky ‡ † Institut UTINAM -UMR 6213, CNRS / Université de Franche-Comté, 25000 Besançon, France ‡ Department of Chemistry, Eszterházy Károly Catholic University, Leányka u. 6, H-3300 Eger, Hungary

S1. NPT simulation details

To determine a consistent and accurate relation between pressure and chemical potential, we have performed a series of Monte Carlo simulations in the (N,P,T) ensemble, over the same ranges of pressure and temperature than in the GCMC simulations, in a cubic box (of side 3.10 4 Å) containing 300 molecules (with varying proportions of each species depending on the simulated gas phase mixture). In these simulations, the molecules have been subjected not only to translation and rotation MC moves but also to changes in their internal configuration. In addition, identity exchange move have been considered in the case of mixtures, allowing exchange of positions of two randomly selected molecules. These four moves have been performed with equal probabilities of 24% each. The system as a whole has also been subjected to the volume change move, with a 4% probability. Once the systems have reached equilibrium after a total of 10 8 Monte Carlo steps, the Widom insertion test phase took place over an additional 10 8 Monte Carlo steps, with probabilities of Monte Carlo movements during this phase set at 15% for translation and for rotation, 15% for internal configuration change, 15% for identity exchange, 1% for volume change, and 39% for insertion tests. Note that the interactions between the molecules have been calculated using the same models than in the GCMC simulations, with cut-offs equal to half the size of the simulation box. Then, the conversion between chemical potential and pressure has been achieved as in our previous works. [1-3]


Valérie Berthé
email: berthe@irif.fr

Herman Goulet-Ouellet

⋆⋆

On Substitutions Preserving their Return Sets

Keywords: Substitutions, Return words, Sturmian substitutions

We consider the question of whether or not a given primitive substitution preserves its sets of return words-or return sets for short. More precisely, we study the property asking that the image of the return set to a word equals the return set to the image of that word. We show that, for bifix encodings (where images of letters form a bifix code), this property holds for all but finitely many words. On the other hand, we also show that every conjugacy class of Sturmian substitutions contains a member for which the property fails infinitely often. Various applications and examples of these results are presented, including a description of the subgroups generated by the return sets in the shift of the Thue-Morse substitution. Up to conjugacy, these subgroups can be sorted into strictly decreasing chains of isomorphic subgroups weaving together a simple pattern. This is in stark contrast with the Sturmian case, and more generally with the dendric case (including in particular the Arnoux-Rauzy case), where it is known that all return sets generate the free group over the underlying alphabet.

Since the seminal work of Durand at the end of the 90s-including a generalization of Cobham's theorem [START_REF] Durand | A generalization of Cobham's theorem[END_REF], a characterization of sequences defined by primitive substitutions [START_REF] Durand | A characterization of substitutive sequences using return words[END_REF], and with Host and Skau, a description of the dimension groups of substitution dynamical systems via sequences of Kakutani-Rohlin partitions [START_REF] Durand | Substitutional dynamical systems, Bratteli diagrams and dimension groups[END_REF]-return words have proved to be a useful tool in combinatorics on words and symbolic dynamics. They have been used, among other things, to characterize Sturmian words [START_REF] Vuillon | A characterization of Sturmian words by return words[END_REF], study maximal bifix codes [START_REF] Berstel | Bifix codes and Sturmian words[END_REF][START_REF] Berthé | Maximal bifix decoding[END_REF], or else study the Schützenberger group, a topological/algebraic invariant of minimal shift spaces which takes the form of a projective profinite group [START_REF] Almeida | Profinite groups associated with weakly primitive substitutions[END_REF][START_REF] Almeida | Presentations of Schützenberger groups of minimal subshifts[END_REF][START_REF] Almeida | A geometric interpretation of the Schützenberger group of a minimal subshift[END_REF].

A return set is a set of (first) return words to a given word in a given set of words, for instance the language of a substitution. The Return Theorem, a striking result by a group of authors which includes the first author of this paper, states that in a dendric shift space (a common generalization of episturmian shifts and codings of interval exchanges), every return set is a basis for the free group over the alphabet of the shift [START_REF] Berthé | Acyclic, connected and tree sets[END_REF]. This kind of stable behavior for subgroups generated by return sets was observed in the more general class of suffix-connected shift spaces, introduced by the second author [START_REF] Goulet-Ouellet | Suffix-connected languages[END_REF]. The behavior of return sets in a given shift space has also been linked, by Almeida and Costa, to the Schützenberger group; in particular, global stability of the subgroups generated by the return sets entails freeness of the Schützenberger group [START_REF] Almeida | A geometric interpretation of the Schützenberger group of a minimal subshift[END_REF]. Together with results from [START_REF] Costa | Conjugacy invariants of subshifts: An approach from profinite semigroup theory[END_REF][START_REF] Goulet-Ouellet | Pronilpotent quotients associated with primitive substitutions[END_REF], this can be used to show that dendric shift spaces are not topologically conjugate to shift spaces of primitive uniform substitutions-though the same conclusion can also be reached using other means [START_REF] Berthé | Rigidity and substitutive dendric words[END_REF].

In addition to the algebraic regularity from the Return Theorem, let us recall that Durand proved in [START_REF] Durand | A characterization of substitutive sequences using return words[END_REF] a combinatorial regularity property for return sets of substitutive shift spaces. His result uses derived sequences, which are obtained by recoding an infinite word with respect to some return set. Durand showed that having a finite number of derived sequences precisely characterizes substitutive sequences. Coming from a dynamical point of view, this might set some expectation that return sets in substitutive shift spaces should also show some algebraic regularity. However, looking at concrete examples quickly reveals that, outside of the well-charted dendric case, the subgroups generated by return sets can exhibit a more complicated behavior-and in fact, little is known about how these subgroups behave in general.

In a roundabout way, this was the initial motivation for the work presented in this paper, which started as a study of the subgroups generated by the return sets in the Thue-Morse shift. In this case the subgroups, far from all being equal, form instead several decreasing chains of subgroups, crossing and weaving into each other. Moreover, these chains can be completely described-details are given in Proposition 1 and in § 5-thanks to a simple preservation property which sparked our curiosity: with finite exceptions, the Thue-Morse substitution maps its return sets to other return sets. Our first main result (Theorem 1) states that this property holds for every primitive aperiodic bifix substitution (and by bifix, we mean that the letter images form a bifix code). Our second main result (Theorem 2) states that, on the other hand, it fails for primitive Sturmian substitutions, at least up to conjugacy.

Preliminaries

Throughout this paper, A is a finite alphabet, A * is the set of words over A equipped with concatenation, and ε is the empty word. For a subset B ⊆ A * , let B * denote the submonoid of A * generated by B. We use u • v as a shorthand for the pair (u, v), in particular when u and v are in A * . We say that s

• t refines u • v if u is a suffix of s and v is a prefix of t, i.e. s • t = s ′ u • vt ′ for some s ′ , t ′ ∈ A * ; we then write u • v ⪯ s • t.
Let A Z be the set of two-sided infinite word over A. Recall that a shift space is a subset X ⊆ A Z invariant under the shift map (x n ) n∈Z → (x n+1 ) n∈Z and its inverse, and closed for the product topology of A Z (A being equipped with the discrete topology). The language of X, denoted L(X), is the set of all finite words occurring as factors (i.e. consecutive blocks) in the elements x ∈ X. A left extension of u ∈ L(X) is a word x ∈ A * such that xu ∈ L(X). Likewise, y is a right extension of u if uy ∈ L(X). We say that u is left special if it admits at least two left extensions of length 1, and right special if those are instead right extensions. We say that u is bispecial if it is both left and right special.

Given a pair u • v with uv ∈ L(X), a word r ∈ A * is a return word (sometimes called first return word ) to u • v in X if urv ∈ L(X) ∩ uvA * ∩ A * uv and contains exactly two occurrences of uv. We denote by R u•v the set of return words to u • v in X (the shift space is easily inferred from context). Those are called the return sets of X. Note that this definition differs from the one-perhaps more common in the literature-where return words are one-sided ; in our notation, this would correspond to return sets of the form R ε•v or R u•ε . Concrete instances of return sets are given, for instance, in Example 3. The lemma below collects some standard properties of return sets that will be useful later.

Lemma 1 ([2,14,15]). Let X be a shift space and u • v be a pair with uv ∈ L(X).

(i) The set R u•v generates freely the submonoid R * u•v of A * , i.e. it is a code. (ii) The submonoid R * u•v contains {w ∈ L(X) : uwv ∈ L(X) ∩ uvA * ∩ A * uv}. (iii) R u•v = u -1 R ε•uv u = vR uv•ε v -1 .
The item (iii) above says that we can pass from two-sided to one-sided return sets and back simply by conjugating. This however should not be construed as saying that two-sided return sets are useless: there are circumstances in which the two-sided version plays an important role, as in [START_REF] Almeida | Presentations of Schützenberger groups of minimal subshifts[END_REF][START_REF] Almeida | A geometric interpretation of the Schützenberger group of a minimal subshift[END_REF][START_REF] Durand | Substitutional dynamical systems, Bratteli diagrams and dimension groups[END_REF]. In the present paper, it is helpful in handling some aspects of the proof of our first main result (Theorem 1).

Let σ : A * → A * be a primitive substitution. We denote by X σ the two-sided shift space associated with σ (as per the standard definition; see e.g. [28, § 5] or [17, § 1.4]) and by L(σ) the language of X σ . By primitivity, all the return sets in X σ are finite. We say that σ is aperiodic if the shift space X σ does not consist of a single finite orbit of the shift map. We say that σ is an encoding if σ is injective on A * . In particular, encodings preserve codes [START_REF] Berstel | Codes and Automata[END_REF]Corollary 2.1.6]. We say that σ is a bifix encoding, or simply bifix, if for all distinct letters a ̸ = b ∈ A, σ(a) is neither prefix nor suffix of σ(b). Every bifix substitution is an encoding [7, Proposition 2.1.9]. Finally, we say that σ is left proper (or right proper ) if there is a letter a 0 ∈ A such that σ(a) ∈ a 0 A * (or σ(a) ∈ A * a 0 ) for all a ∈ A.

A Sturmian substitution is a binary substitution, here defined on the alphabet {0, 1}, which preserves Sturmian sequences-aperiodic sequences in {0, 1} N of minimal factor complexity (see [START_REF] Lothaire | Algebraic combinatorics on words[END_REF]). By [START_REF] Berstel | A characterization of Sturmian morphisms[END_REF]Theorem 5], this is the same as being weakly Sturmian (the image of some Sturmian sequence is Sturmian); in the primitive case, this is also the same as generating a Sturmian shift space. In particular, primitive Sturmian substitutions are aperiodic.

In the next example, we recall two important substitutions which will reappear several times in this paper.

Example 1. Consider the following primitive binary substitutions:

τ : 0 → 01, 1 → 10 and ϕ : 0 → 01, 1 → 0.
The former is known as the Thue-Morse substitution, and it is bifix and aperiodic; it is neither left nor right proper. The latter, the Fibonacci substitution, is Sturmian, hence aperiodic, and left proper; it is neither bifix nor right proper.

Let F (A) be the free group over A. It consists, we recall, of the words over the alphabet A ∪ A -1 (where A -1 = {a -1 : a ∈ A} is a disjoint copy of A) which are irreducible under the rewriting rules aa -1 → ε for all a ∈ A. There is a natural embedding A * → F (A) which allows us to view words as elements of F (A) and substitutions as endomorphisms of F (A)-and we do so whenever convenient. The endomorphisms of F (A) that are also substitutions are called positive; the positive automorphisms of F ({0, 1}) are precisely the Sturmian substitutions [START_REF] Mignosi | Morphismes sturmiens et règles de Rauzy[END_REF][START_REF] Wen | Local isomorphisms of the invertible substitutions[END_REF].

For a given subset B ⊆ F (A), we let ⟨B⟩ denote the subgroup of F (A) generated by B. Similar to codes in A * , we call B free if the inclusion B → F (A) extends to an isomorphism F (B) ∼ = ⟨B⟩. We also say that B is a basis of the subgroup ⟨B⟩. In contrast with free monoids, every subgroup of the free group admits a basis (this is the celebrated Nielsen-Schreier theorem; see [START_REF] Lyndon | Combinatorial Group Theory[END_REF]Proposition 2.11]), and in fact infinitely many bases whenever it is not cyclic. We will be interested in the subgroups generated by the return sets of a shift space, dubbed the return groups of the shift space.

Example 2. If X ⊆ A Z is a minimal dendric shift space, then all return groups are equal to F (A), and moreover every return set viewed as a subset of F (A) is free; this is the aforementioned Return Theorem [START_REF] Berthé | Acyclic, connected and tree sets[END_REF]Theorem 4.5]. This applies in particular to Sturmian shift spaces, as these are all dendric. More generally, if X is suffix-connected, then it has only finitely many return groups which all belong to the same conjugacy class [19, Theorem 1.1].

Preservation of Return Words

Let σ be a primitive substitution over A. Given a pair u • v with uv ∈ L(σ), we consider the following preservation property:

σ(R u•v ) = R σ(u)•σ(v) . (P)
Let us state our two main results.

Theorem 1. Let σ be a primitive aperiodic bifix substitution. There exists a constant K > 0 such that the property (P) holds for every pair u•v with uv ∈ L(X) and |uv| ≥ K.

Theorem 2. Let σ be a primitive Sturmian substitution. There is a substitution in the conjugacy class of σ for which the property (P) fails for infinitely many pairs u • v with uv ∈ L(σ).

The proofs are given later, in § 3 and § 4 respectively. The first main result will be used to describe the return groups in the Thue-Morse shift.

Proposition 1. Up to conjugacy, all return groups in the shift space of the Thue-Morse substitution are equal to one of the following subgroups of F ({0, 1}):

⟨τ n (0), τ n (1)⟩, ⟨τ n (1), τ n (00), τ n (0110)⟩, ⟨τ n (0), τ n (11), τ n (1001)⟩,
with n ≥ 0. Moreover, the given generating sets are bases of these subgroups.

The proof of Proposition 1 will be given in § 5. The remainder of this section explores various aspects of the property (P), starting with an example.

Example 3. For the Thue-Morse substitution τ : 0 → 01, 1 → 10, the pair 0 • ε fails the property (P) since

R 0•ε = {0, 10, 110} and R 01•ε = {01, 001, 101, 1001}.
On the other hand, R 0•1 = {10, 100, 110, 1100} and

τ (R 0•1 ) = {1001, 100101, 101001, 10100101} = R 01•10 .
In fact, every pair u • v with uv ∈ L(τ ) and |uv| ≥ 2 satisfies the property (P). This is a consequence of the proof of Theorem 1 detailed in Proposition 4.

Next is a simple lemma that will be useful for the proof of Theorem 1.

Lemma 2. Let σ be a primitive substitution which is also an encoding. A pair u • v with uv ∈ L(σ) satisfies the property

(P) if and only if R σ(u)•σ(v) ⊆ σ(R * u•v ). Proof. The inclusion σ(R u•v ) ⊆ R * σ(u)•σ(v) always holds by Lemma 1 (ii). Since σ(R u•v ) and R σ(u)•σ(v)
(R u•v ) = Card R σ(u)•σ(v) .
The lemma below shows that the property (P) depends only on the word uv, rather than on the pair u • v. Lemma 3. Let σ be a primitive substitution. A pair u • v with uv ∈ L(σ) satisfies the property (P) if and only if the pair uv • ε does.

Proof. This is straightforward by Lemma 1 (iii); for instance, assuming that u • v satisfies the property (P), we get

σ(R uv•ε ) = σ(v) -1 σ(R u•v )σ(v) = σ(v) -1 R σ(u)•σ(v) σ(v) = R σ(uv)•ε . ⊓ ⊔
The next proposition gives a sufficient condition for a pair of words to fail the property (P). It is used to prove Theorem 2.

Proposition 2. Let σ be a primitive substitution which is also left proper. Let u • v be a pair with uv ∈ L(σ). If σ(uv) is right special, then u • v fails the property (P).

Proof. Thanks to the previous lemma, we may assume that v = ε. Since σ is left proper, there is a letter a

0 ∈ A such that σ(R u•ε ) ⊆ a 0 A * . As σ(u) is right special, it has a right extension b ∈ A with b ̸ = a 0 .
Then, there is a return word r ∈ R σ(u)•ε such that r ∈ bA * , and clearly r /

∈ σ(R u•v ). ⊓ ⊔ Remark 2.
In Proposition 2, the assumption that σ is left proper can be replaced by the local condition that uv is not right special.

Of course, similar statements hold with the left special property for right proper substitutions or under the local assumption that uv is not left special. We thus deduce the following result as a straightforward consequence of Theorem 1; to the best of our knowledge, this is new.

Corollary 1. Let σ be a primitive aperiodic bifix substitution. There is a constant K > 0 such that, for every u ∈ L(σ) with |u| ≥ K, u is left special, right special or bispecial whenever σ(u) is. If σ is left proper (or right proper), then Im(σ) ∩ L(σ) contains only finitely many right special words (or left special words).

Finally, we give an example illustrating failure of the property (P); it is a special case of Theorem 2.

Example 4. Consider the Fibonacci substitution ϕ : 0 → 01, 1 → 0. Let u ∈ L(ϕ) be right special and consider the word v = u1. We claim that ϕ(v) is right special. Indeed, we have u0 and v0 ∈ L(ϕ), hence ϕ(u0) = ϕ(v)1 and ϕ(v1) = ϕ(v)0 both belong to L(ϕ) as well. Since ϕ is left proper, it follows from Proposition 2 that ϕ fails the property (P) for all such words v. Take for instance u = 0010. Then, v = 00101 and σ(v) = 01010010 have for return sets

R v•ε = {00100101, 00101} and R σ(v)•ε = {01010010, 10010},
confirming that v • ε fails the property (P). Remark 3. Example 3 shows that the property (P) may fail for reasons other than Proposition 2 and Remark 2, as τ is neither left nor right proper and 0 is bispecial in X τ .

The Bifix Case

This section presents the proof of Theorem 1, which relies on the following notions and subsequent lemma. Let σ : A * → A * be a primitive substitution. A parse of u • v (under σ) is a pair x • y such that xy ∈ L(σ) and σ(x) • σ(y) ⪰ u • v; we call parsable a pair that admits a parse. On the other hand, an interpretation of u ∈ L(σ) is a triple (s, w, t) such that w ∈ L(σ) and σ(w) = sut. A pair u • v is called synchronizing if, for every interpretation (s, w, t) of uv, there exists a pair

x ′ • y ′ such that w = x ′ y ′ and σ(x ′ ) • σ(y ′ ) = su • vt.
The next lemma is a reformulation of Mossé's celebrated recognizability theorem [START_REF] Mossé | Puissance de mots et reconnaissabilité des points fixes d'une substitution[END_REF] which is due to Kyriakoglou [START_REF] Kyriakoglou | Iterated morphisms, combinatorics on words and symbolic dynamical systems[END_REF]Proposition 3.3.20] (see also [START_REF] Durand | Dimension Groups and Dynamical Systems[END_REF]Proposition 1.4.38]).

u ′ q w v ′ r σ(u) σ(v) s t σ σ

Lemma 4 ([22]

). Let σ be a primitive aperiodic substitution. There exists a constant L > 0 such that every parsable pair u • v with |u| and |v| ≥ L is synchronizing.

Proof (Theorem 1). Let L > 0 be the constant given by the previous lemma.

σ(u)•σ(v) ⊆ σ(R * u•v ) holds. Fix a return word r ∈ R σ(u)•σ(v) ; we wish to show that r ∈ σ(R * u•v
). Thanks to Lemma 3, we may assume that |u| and |v| are both ≥ ⌊K/2⌋.

We start by showing that r admits a preimage under σ. Let (s, w, t) be an interpretation of σ(u)rσ(v). Note that σ(u) • σ(v) is parsable with |σ(u)| and |σ(v)| ≥ L, therefore it is synchronizing. Since moreover σ(u)r

• σ(v) ⪰ σ(u) • σ(v), there is a factorization w = w ′ v ′ such that σ(w ′ ) • σ(v ′ ) = sσ(u)r • σ(v)t.
(u) • r ⪰ σ(u) • σ(p)
, there exists a factorization w ′ = u ′ q such that σ(u ′ ) = sσ(u) and σ(q) = r (see Fig. 1).

We finish by showing that q ∈ R * u•v . For this purpose, recall that bifix codes generate submonoids that are biunitary [7, Proposition 2.2.7]. In the case at hand, it means that the two following properties hold:

x, xy ∈ Im(σ) =⇒ y ∈ Im(σ) and x, yx ∈ Im(σ) =⇒ y ∈ Im(σ).

Applying these properties to σ(u ′ ) = sσ(u) and σ(v ′ ) = σ(v)t respectively yields words s ′ , t ′ such that u ′ = s ′ u and v ′ = vt ′ . In particular, uqv ∈ L(σ). Since r ∈ R σ(u)•σ(v) , there are words x, y that σ(u)rσ(v) = σ(uv)x = yσ(uv). Applying the above properties once again, this time to σ(uqv) = σ(uv)x = yσ(uv), produces words x ′ , y ′ such that uqv = uvx ′ = y ′ uv. By Lemma 1 (ii), q ∈ R * u•v .

⊓ ⊔

For instance, in the case of the Thue-Morse substitution, the constant K from the proof above equals 2 (see Proposition 4). Using a computability result of Durand and Leroy, we also deduce the following. Proposition 3. The constant K of Theorem 1 is computable.

Proof. The constant L of Lemma 4 is related in a straightforward way to the recognizability constant (see the proof of [START_REF] Kyriakoglou | Iterated morphisms, combinatorics on words and symbolic dynamical systems[END_REF]Proposition 3.3.20]), which is computable by a result of Durand and Leroy [START_REF] Durand | The constant of recognizability is computable for primitive morphisms[END_REF]. As for the constant M , note that it can be set equal to any m ∈ N such that min{|r| : r ∈ R w•ε } ≥ |σ|⌈L/⟨σ⟩⌉ for all w in the finite set L(σ) ∩ A m . Since the language of a primitive substitution and its return sets are computable, this property is decidable for each natural number, and therefore the constant M is computable.

⊓ ⊔

The Sturmian Case

In preparation for the proof of Theorem 2, let us set up some notation and recall a handful of facts about Sturmian substitutions. Given a pair u • v of words in {0, 1} * , we let [u • v] be the substitution 0 → u, 1 → v. We say that two substitutions σ and ρ defined on the same alphabet A are rotationally conjugate, or simply conjugate, if there is a word w ∈ A * such that either σ(a)w = wρ(a) for all a ∈ A, or wσ(a) = ρ(a)w for all a ∈ A. Then, we write respectively σ = w -1 ρw or σ = wρw -1 . Consider the set of standard pairs: it is the smallest subset of {0,

1} * × {0, 1} * which contains {0•1, 1•0} and which is closed under u•v → u•uv and u•v → vu•v. Substitutions of the form [u•v] with u•v
σ 0 = σ, σ 1 = a -1 0 σ 0 a 0 , σ 2 = a -1 1 σ 1 a 1 , • • • σ n = a -1 n-1 σ n-1 a n-1 ,
with n = |x| + |y| -1, and where σ i (0) and σ i (1) both start with the letter a i for

0 ≤ i < n. Let σ i = [x i • y i ].
Let E be the automorphism of {0, 1} * exchanging 0 and 1 and consider the substitutions

σ E i = E • σ i • E.
They also form a Sturmian conjugacy class where L(σ E i ) = E(L(σ i )). Thus, up to replacing σ by σ E if needed, we may assume that x = ys for some word s, which is non-empty since σ is not periodic. Next, note that x n and y n start with distinct letters by [START_REF] Séébold | On the conjugation of standard morphisms[END_REF]Lemma 5]; in particular, y n is not a prefix of x n . Since y 0 is a prefix of x 0 , we may let j be the largest index satisfying 0 ≤ j < n such that y j is a prefix of x j but y j+1 is not a prefix of x j+1 . Write x j = y j s j and let b j be the first letter of s j ; the choice of j implies that a j ̸ = b j .

Next, observe that there exists ℓ ≥ 0 such that 01 ℓ 0 and 01 ℓ+1 0 ∈ L(σ) with these being the only two factors of the form 01 k 0 in L(σ). (Note that ℓ may be equal to 0, i.e. the letter 1 does not need to be the most frequent letter.) In particular, 1 ℓ is right special and admits 0 and 10 as right extensions. Let u be a right special factor in X σ . Since the right special factors in X σ are all suffixes of one another, either u = 1 k with k < ℓ or 1 ℓ is a suffix of u. Either way, it follows that u10 = v0 ∈ L(σ). Since σ j is left proper, by Proposition 2, we are done if we can show that σ j (v) is right special. But note that σ j (u0) = σ j (u)y j s j = σ j (v)s j ∈ L(σ) and σ j (v0) = σ j (v)x j ∈ L(σ).

Hence, σ j (v)a j and σ j (v)b j ∈ L(σ). Since a j ̸ = b j , this shows that σ j (v) is indeed right special. ⊓ ⊔

The next example illustrates the construction in the above proof.

σ 0 = [1101 • 110], σ 1 = [1011 • 101], σ 2 = [0111 • 011], σ 3 = [1110 • 110], σ 4 = [1101 • 101], σ 5 = [1011 • 011].
Here, σ 3 is the first member of the conjugacy class whose second component is not a prefix of the first, so j = 2. Moreover, 0110 and 01110 ∈ L(σ), so ℓ = 2.

According to the argument above, if u is right special, then v = u1 is such that σ 2 (v) is also right special. And indeed, σ 2 (v0) = σ 2 (v)011 and σ 2 (u01) = σ 2 (v)1011 both belong to L(σ). Since σ 2 is left proper, we can apply Proposition 2 to conclude that v fails the property (P) with respect to σ 2 .

For a concrete example, take u = 11. In this case, v = 111, σ 2 (v) = 011011011, and the fact that v fails the property (P) is apparent from the equalities

R v•ε = {0110110111, 0110110110111}, R σ2(v)•ε = {011, 1011011011}.

Return Words of Thue-Morse Substitution

In this section, we give precise formulas for the return sets in the shift of the Thue-Morse substitution τ : 0 → 01, 1 → 10 and proceed to give the proof of Proposition 1. Our starting point is the following simple proposition, which is a consequence of the proof of Theorem 2. Proposition 4. Every pair u • v with uv ∈ L(τ ) and |uv| ≥ 2 satisfies the property (P) with respect to τ .

Proof. First, we claim that the constant L from Lemma 4 can be set to L = 2. To establish this, it suffices to show that every parsable pair u • v with |u| = |v| = 2 is synchronizing for τ . There are four such pairs, namely, 01 • 10, 10 • 01, 01 • 01 and 10 • 10. That the first two pairs are synchronizing is obvious: otherwise, 11 and 00 would be images of letters, and this is not the case. If, say, 01 • 01 was not synchronizing, then there would be a pair x • y such that 010 • 1 ⪯ x • y, xy ∈ L(τ ) and x, y ∈ Im(τ ). This would then imply 1010•10 ⪯ x•y, and therefore 101010 ∈ L(τ ), a contradiction. The fact that 10 • 10 is synchronizing is proved similarly. This proves the claim. Finally, since

L = |τ | = ⟨τ ⟩ = 2, we find that R(|τ |⌈L/⟨τ ⟩⌉) = R(2) = 2, so K = max(2, 2) = 2. ⊓ ⊔ uv u • v Ru•v
wn wn- As noted in Balková et al. [4, § 3.1], every return set is in fact determined by a bispecial pair-a pair u • v such that uv is bispecial. To see why, fix a minimal aperiodic shift space X ⊆ A Z and consider a pair u • v with uv ∈ L(X). If u • v is not left special, then we may replace it by au • v where a ∈ A is the unique left extension of uv of length 1 in L(X), and this does not change the return set. A similar phenomenon occurs with the unique right extension b ∈ A of uv of length 1 when uv is not right special. We can therefore consider the minimal refinement s • t ⪰ u • v which is bispecial; let us call this the bispecial closure of u • v. What Balková et al. observed is stated (with an extra maximality property) in the next lemma.

Lemma 5. Let X ⊆ A Z be a minimal aperiodic shift space and u • v be a pair with uv ∈ L(X). The bispecial closure of u • v is the greatest refinement s

• t ⪰ u • v with st ∈ L(X) such that R s•t = R u•v .
Next, we recall the classification of the bispecial factors in the Thue-Morse shift given by de Luca and Mione in 1994 [START_REF] De Luca | On bispecial factors of the Thue-Morse word[END_REF]. Let us mention in passing that the idea behind de Luca and Mione's result has been significantly generalized by Klouda, who gave an algorithm for computing similar classifications for any primitive aperiodic substitution [START_REF] Klouda | Bispecial factors in circular non-pushy D0L languages[END_REF].

Proposition 5 ([24]). Every bispecial factor of length ≥ 2 in X τ is equal to one of the following words, for some n ≥ 0:

w n = τ n (01), w n = τ n (10), z n = τ n (010), z n = τ n (101).
Let also w -1 = 0 and w -1 = 1. Note that w n = τ (w n-1 ), z n = τ (z n-1 ), and so on. Note moreover the following factorizations for n ≥ 0:

w n = w n-1 w n-1 , z n = w n w n-1 , w n = w n-1 w n-1 , z n = w n w n-1 .
Thanks to Proposition 4, computing the return sets to the bispecial factors is as simple as computing them for w 0 , z 0 , w 0 and z 0 , and taking direct images by τ n . Table 1 gives the return sets to each of the pairs given by the factorizations above.

Next, we proceed to describe the return groups as stated in Proposition 1, again noting that we only need to consider those that are determined by bispecial pairs. For this purpose, it is useful to have in mind the following lemma. A quick and easy way to establish this is to show that {01, 10} is a free subset of F ({0, 1}), which is indeed equivalent to τ being injective. This can be done using Stallings' algorithm, a powerful algorithm that we also use in the proof of Proposition 1 below. This algorithm, given a finitely generated subgroup of a free group, produces an automaton, known as the Stallings automaton of the subgroup, which allows to easily test for membership. It can also be used to test whether a subset is a basis of the subgroup, and has numerous other applications beyond (calculate subgroup indexes, test for conjugacy, compute intersections, etc.). See the paper by Kapovich and Myasnikov [START_REF] Kapovich | Stallings foldings and subgroups of free groups[END_REF] for a detailed description of Stallings' algorithm and many of its applications; see also the paper by Touikan [START_REF] Touikan | A fast algorithm for Stallings' folding process[END_REF] for an efficient implementation of the algorithm.

Proof (Proposition 1). Consider the following return groups in X τ :

W n = ⟨R wn-1•wn-1 ⟩, Z n = ⟨R wn•wn-1 ⟩, W n = ⟨R wn-1•wn-1 ⟩, Z n = ⟨R wn•wn-1 ⟩.
By Lemma 5 and Proposition 5, every return set R u•v with uv ∈ L(τ ) and |uv| ≥ 2 is conjugate to one of these. It remains to show that these subgroups have the bases stated in the proposition. First, we claim that, for all n ≥ 0, {τ n (0), τ n (1)} forms a basis of W n and {τ n (1), τ n (00), τ n (0110)} forms a basis of Z n . This can be checked directly in the case n = 0 using Stallings' algorithm (the relevant Stallings automata are found in Fig. 2) and it then holds for all n ≥ 0 thanks to Proposition 5 and Lemma 6.

To finish, let E be the automorphism exchanging 0 and 1, viewed as an automorphism of F ({0, 1}). On the one hand, it is clear that W n = E(W n ) = W n , which shows that the return groups W n are redundant. On the other hand, E(Z n ) = Z n and thus Z n have the required basis for all n ≥ 0. This concludes the proof. Up to conjugacy, every return group in the Thue-Morse shift occurs in one (and only one) of these sequences. It is not hard to see that they are all strictly decreasing with respect to inclusion. The next proposition clarifies how these sequences are related to one another. The proof is a straightforward application of Stallings algorithm; all relevant Stallings automata are found in Fig. 2. Proposition 6. For all n ≥ 0,

W n = ⟨Z n ∪ Z n ⟩ > Z n ∩ Z n > W n+2 .
The Hasse diagram of the poset formed by the return groups W n = W n , Z n and Z n ordered by inclusion is depicted in Fig. 3.

Conclusion

The statement of Theorem 2 is somewhat unsatisfactory, in the sense that we are left wondering what happens with the other Sturmian substitutions within the conjugacy class. Based on examples, we expect all Sturmian substitutions to satisfy the conclusion of Theorem 2. It is not clear at the moment how generalizations, like dendric substitutions, might behave.

On the other hand, Theorem 1 shows that being bifix is a sufficient condition for having the property (P) for all but finitely many pairs, but we are not sure whether or not it is necessary. We would like to know if other natural classes of primitive substitutions satisfy the property (P) for all but finitely many pairs (clearly, by Theorem 2, such classes cannot contain the class of primitive Sturmian substitutions). It might also be interesting to investigate the property (P) for non-primitive substitutions.

Finally, we wonder whether or not the sophisticated pattern formed by the return groups of the Thue-Morse shift is common among primitive aperiodic bifix substitutions. While classifications of bispecial factors can always be obtained in the primitive case thanks to the aforementioned algorithm of Klouda, the return preservation property might not suffice to determine all return groups, since the general form of Klouda's algorithm involves taking more than direct images-it uses what Klouda calls f B -images, which in the case of the Thue-Morse substitution reduces to the direct image. In this regard, the Thue-Morse substitution seems to be quite special.

Fig. 1 .

 1 Fig. 1. An illustration of the proof of Theorem 1.

  Let us write |σ| = max{|σ(a)| : a ∈ A} and ⟨σ⟩ = min{|σ(a)| : a ∈ A}. Let M > 0 be such that min{|r| : r ∈ R u•v } ≥ |σ|⌈L/⟨σ⟩⌉ for all pairs u • v with uv ∈ L(σ) and |uv| ≥ M . Such a constant M exists by primitivity [14, Lemma 3.2]. Finally, let K = max M, 2⌈L/⟨σ⟩⌉).Fix a pair u • v with uv ∈ L(σ) and |uv| ≥ K. By Lemma 2, it suffices to show that the inclusion R

  Let p be the prefix of v of length ⌈L/⟨σ⟩⌉. Then |σ(p)| ≥ L, and since |r| ≥ |σ|⌈L/⟨σ⟩⌉, σ(p) is a prefix of r. The pair σ(u) • σ(p) is synchronizing, being parsable with |σ(u)| and |σ(p)| ≥ L. Since σ

  standard are called standard substitutions. Every standard (binary) substitution is Sturmian [5, Proposition 2.6]. If σ is a Sturmian substitution, then its conjugacy class contains |σ(01)| -1 substitutions, all of which are Sturmian and one of which is standard [29, Propositions 9 and 10]. Moreover, all the primitive Sturmian substitutions in the same conjugacy class share the same shift space [29, Lemma 8]. Proof (Theorem 2). Let σ = [x • y] be a primitive standard substitution whose conjugacy class is given by

Example 5 .

 5 Let σ = [1101 • 110]. It is a primitive standard substitution with conjugacy class

Fig. 2 .Lemma 6 .

 26 Fig. 2. Stallings automata of the subgroups W0, W 0, Z0, Z0 and Z0 ∩ Z0.

Fig. 3 .

 3 Fig. 3. Hasse diagram of return groups of Xτ ordered by inclusion.

  are both codes (by Lemma 1 (i)), they are equal if and only if they generate the same submonoid (by [7, Corollary 2.2.4]).

⊓ ⊔ Remark 1. In the previous lemma, we can replace the assumption that σ is an encoding by the local condition Card σ

Table 1 .

 1 1 • wn-1 {τ n (10), τ n (110), τ n (100), τ n (1100)} wn wn-1 • wn-1 {τ n (01), τ n (001), τ n (011), τ n (0011)} Return sets in Xτ .

	zn	wn • wn-1	{τ n (001), τ n (01101), τ n (0011001), τ n (011001101)}
	zn	wn • wn-1	{τ n (110), τ n (10010), τ n (1100110), τ n (100110010)}

⋆ This work was supported by the Agence Nationale de la Recherche through the project "Codys" (ANR-18-CE40-0007). ⋆⋆ Corresponding author.


Fernando Terroso-Sáenz
email: fernando.terroso@upct.es

Andres Muñoz
email: andres.munoz@uca.es

Philippe Roose
email: philippe.roose@univ-pau.fr

War & Music: The impact of the Ukrainian War on the Music Listening Behaviour in Eastern Europe*

Keywords: Music-listening behaviour, Spotify, timeseries, warfare, Ukrainian War

According to many analysts, February 2022 is dated as the exacerbation of the armed conflict between Ukraine and Russia. Of course, this shocking event has an impact on mood of people. In this context, we studied the music listening habits of people in Eastern European countries. We applied timeseries and statistical analysis of mood-relate features (danceability, valence and energy) to a Spotify dataset in order to define new paradigms of emotion-aware systems that can be applied to detect the population's state of mind. The results show that the trend in the type of music changes notably towards songs with less genre diversity and reflecting less positivity.

I. INTRODUCTION

Due to music digitization, songs can now be defined in terms of quantitative features describing certain sound-and moodrelated properties such as their danceability or valence. Since these features can be analyzed from multiple points of view, it has encouraged an interesting course of action to capture people's emotions through their music listening habits [START_REF] Anderson | just the way you are": Linking music listening on spotify and personality[END_REF]. In turn, this has allowed to develop and assess novel emotionaware systems like music recommender systems [START_REF] Dai | Personalized knowledge-aware recommendation with collaborative and attentive graph convolutional networks[END_REF], [START_REF] Holtz | The engagement-diversity connection: Evidence from a field experiment on spotify[END_REF].

Apart from that, on 24th February 2022, the conflict between Russia and Ukraine was exacerbated with a large military deployment, which has turned into a long-term warfare causing major social and economic disturbances, especially in Eastern Europe [START_REF] Kurapov | Toward an understanding of the russian-ukrainian war impact on university students and personnel[END_REF], [START_REF] Lo | The russo-ukrainian war and financial markets: The role of dependence on russian commodities[END_REF]. Such major conflicts inevitably have a substantial impact on the mood of the population.

Considering the above, this work focuses on analyzing whether the emotional impact of the Ukrainian war is reflected in the music listening habits of the citizens of some of the most directly affected Eastern European countries (Estonia, Latvia, Lithuania, Belarus, Ukraine, Romania and Bulgaria) through Spotify. To do so, we apply different techniques from timeseries and statistical analysis to a range of mood-relate features of the most streamed songs in such countries. In utility terms, the findings of this work could be useful for defining new paradigms of emotion-aware systems capable of responding to users' distress in severe situations.

Financial support for this research has been provided under grant PID2020-112827GB-I00 funded by MCIN/AEI/10.13039/501100011033. The remainder of the paper is structured as follows. Section 2 gives an overview about existing trends for music streamingdata analysis. Then, section 3 states the data and tools used to measure the music listening patterns from Spotify. The extracted patterns are described and evaluated in section 4. Lastly, section 5 summarizes the main findings of this work.

II. RELATED WORK

A broad line of research in the music streaming data mining field have evaluated the impact of external and disruptive events on the listening habits of the users of music streaming services. For example, some works have assessed the impact of the COVID-19 pandemic on the music consumption. As a matter of fact, authors in [START_REF] Hennessy | Music and mood regulation during the early stages of the covid-19 pandemic[END_REF] showed a remarkable positive correlation between the levels of distress and tension of a set of volunteers and their willingness to listen to music on a regular basis during the first months of the pandemic. Besides, the study stated in [START_REF] Kalustian | evacuate the dancefloor": Exploring and classifying spotify music listening before and during the covid-19 pandemic in dach countries[END_REF] developed a classifier able to label the listening behaviour of Germany, Austria and Switzerland in the non-pandemic and pandemic periods to demonstrate the behavioural change that occurred during these months. Finally, the work in [START_REF] Terroso-Saenz | Evolution of global music trends: An exploratory and predictive approach based on spotify data[END_REF] also performed a correlation analysis between the evolution of the pandemic and the features of most popular songs in several countries showing quite different results depending on the geographical area.

Unlike the previous works, we intend to analyze the actual impact of wars on the listening behaviour of a population in terms of the mood-relate features of the most popular tracks among that population. It is true that previous works have compared the mood-related listening habits of a large population with respect to the COVID-19 pandemic, which can also be considered an extreme and longstanding event. However, we believe that a pandemic and a state of war are very different in nature and in the manner they are perceived and handled by the civil population. To the best of the authors' knowledge, this is the first attempt to compare the evolution of war and the music listening behavior of the directly affected population.

III. EXTRACTION OF THE FEATURES TO CHARACTERIZE

THE MUSIC LISTENING PATTERNS In this section we put forward the general approach followed to model the music-listening behaviour of the target countries.

To do so, we first describe the Spotify dataset used in the study. Next, we state the features extracted from such a dataset that will be analyzed later on to study the effect of the Ukrainian war.

A. Dataset description

For our analysis, we took as input the Top 200 daily rankings released by Spotify, which contain the 200 most played songs in 69 different countries per day 1 . At this point it is important to note that we have opted to center our study on Eastern Europe countries due to their geographical proximity to the Ukrainian conflict. In particular, we focused on Estonia, Latvia, Lithuania, Romania Bulgaria, Belarus and Ukraine itself as target countries. Note that there are other surrounding countries that might also be interesting for this study, such as Moldavia, but the Spotify platform did not provide access to their rank data. Besides, Russia has not been included as Spotify stopped to release its daily lists on March, 3rd 2022.

Using an ad-hoc crawler, daily rankings were extracted for each single country for a 2-year period from 2020-07-17 to 2022-10-11. As a result, we composed a raw dataset of 1,310,738 records comprising 14,816 unique songs. Each record is represented as ⟨s id , l, pos, co, d⟩, indicating that the song with identifier s id was listened to l times in country co occupying the pos position in the co ranking on date d.

B. Mood-related Timeseries

For each song in the dataset, we extracted its associated features provided by the Spotify Developer Platform (SDP) through its public API 2 . For this study we focused on three mood-relate features (according to the Spotify terminology) of a song, namely

• Danceability (D). This feature indicates how suitable a song is for dancing. It takes values from 0 (not danceable) to 1 (highly danceable). • Valence (V). It indicates how positive and upbeat a song is. Again, it ranges from 0 to 1. Hence, songs with a high valence feel cheerful and happy whereas low-valence songs feel sad and depressing. • Energy (E). This attribute measures the intensity and activity of a song. A very energetic song will likely be fast, loud and noisy whereas a low-energetic song will be much more quiet and calm such a Johann Sebastian Bach's prelude score. Again, it ranges from 0 to 1. Next, for each target country it was composed 3 different timeseries, one for each mood feature. Each series took the form for the last day of the covered period). By means of these series, it is possible to measure the evolution through time of the mood-relate features of the most popular songs in each country. The resulting timeseries are shown in Fig. 1. As we can see from this figure, the danceability of Romania was slightly higher than the other countries from the beginning of 2022. Moreover, similar patterns are observed for the valence and energy in Bulgaria.

F c f = ⟨f c 1 → f c 2 → .. → f c L-1 → f c L ⟩

IV. RESULTS

In this section, we evaluate how the timeseries generated in the previous section have been affected by the Ukrainian War.

A. Homogeneity Test of the Mood Timeseries

The first experiment was to perform a homogeneity test on each of the timeseries of the target countries. This test is a statistical method that checks whether two (or more) datasets originate from the same statistical distribution. In a timeseries, it allows to detect one (or more) breakpoint in the series. This breakpoint occurs when the data set changes its distribution.

Table I shows the breakpoints for each country and feature detected by the test for a p-value of 0.05. As can be seen, a single breakpoint was identified for each pair of countries and features. Hence, some interesting findings can be extracted.

Most notably, Ukraine's danceability and energy timeseries changed significantly during April and May 2021. It is important to note that during these months Russia initiated a major military deployment near the Ukrainian border according to some analysts. Regarding the valence feature, its break point occurred in November 2020. During that time, Russia also conducted a second reinforcement along the border 3 . Thus, these breakpoints reflect a shift in the features of the most streamed songs in Ukraine, with a slight drop in their danceability and positivity (e.g, the mean of the danceability feature changed from 0.73 to 0.70 according to Table I) and an increment of their energy.

Comparing these results with those of other target countries, Table I shows that most of them had a slight increment in the energy of the songs a few months after the breakpoint date of Ukraine (2021-05-23). Regarding the danceability and valence features, they changed in all countries except Belarus by the end of 2021 or the first months of 2022 as long as the likelihood of a war increased. In the case of danceability, this involved a decrease in its mean value, which is also consistent with the behaviour observed in Ukraine. For the valence feature, its behaviour was more heterogeneous: it increased in certain countries (Latvia, Bulgaria, Lithuania and Belarus) but declined in others (Romania, Bulgaria and Estonia). Moreover, most of these valence breakpoints occurred in the second half of 2021, during which the Ukrainian conflict was exacerbated. For the sake of completeness, Fig. 2 shows the breakpoints of each mood timeseries for Ukraine along with the mean values after and before the detected change.

B. Clustering of Countries based on the Mood Timeseries

In the previous section, the timeseries of all the target countries were analyzed in an isolated manner. In order to study the relationships among these countries, we developed a second experiment to group the target countries based on the similarity of their mood timeseries.

First, each timeseries F c f was divided into 2 different sub-sequences, namely F c f,prewar that comprised the feature series before the war outrage (from 2021-07-15 to 2022-02-23) and F c f,war that comprised the part of the timeseries covering the war period (from 2022-02-24 to 2022-11-10). Then, we applied the Hierarchical Clustering algorithm using the Dynamic Time Wrapping (DTW) method as a distance metric [START_REF] Giorgino | Computing and visualizing dynamic time warping alignments in R: The dtw package[END_REF] to each type of sub-sequence for all the target countries. DTW is a well-known distance metric for timeseries that exploits the shape of the input timeseries. Fig. 3 shows the resulting clusters of this procedure. Regarding the danceability feature, Figs. 3a andb show that the conflict caused a large aggregation of the target countries. In particular, a large cluster comprising Estonia, Latvia, Lithuania, Ukraine and Bulgaria was created after the conflict began. A similar pattern was observed in the case of the energy feature as Ukraine joined the cluster that already comprised Latvia, Lithuania and Romania during the state of war (Figs. 3e andf). Last, the behaviour of the valence feature was slightly different from the previous ones as, in this case, Ukraine remained in an isolated cluster before and after the war began. This is consistent with the fluctuation and previous differences observed in this feature with respect the danceability and energy patterns.

As a result, this clustering analysis reveals that, overall, the mood timeseries of the target countries has become less diverse after the war. It shows large clusters comprising almost all the target countries, at least for the energy and danceability features.

C. Variation in the Number of Daily Streams

The last experiment studied the impact of the war with respect to the daily number of song stream in the Top-200 rankings of the target countries. As observed in Fig. 4, all countries suffered a sharp decline on the same date as the beginning of the war and in the following days. For example, Ukraine fell from 2,122,435 total streams in its Top-200 rank on February 23rd, 2022 to 1,442,228 (about 32% less) the day after Russia started the conflict. This may reflect the obvious impact on the daily lives of Ukrainians during that particular day. Then, the number of streams steadily increased since that date.

V. CONCLUSIONS AND FUTURE WORK

In this article we study the effect that an extreme event such as the recent war between Ukraine and Russia may have on the state of mind of the Eastern European population most affected by the conflict. To this end, we have analyzed the evolution of the music listened to by the citizens of these countries on the Spotify platform in periods before and after the start of the war. The results show that whilst the flow of songs between the two countries at war has not been altered, a change in the type of music listened to in Ukraine has been detected, specifically related to the danceability and energy during the first months of the conflict. In addition, we found that the type of music listened to in these countries has clearly homogenized since the beginning of the conflict. Finally, we detected that a clear drop in the number of songs played just at the beginning of the war, this number has been gradually increasing again.

  where f c i indicates the mean value of feature f of the songs in the Top-200 list of the country c at the the i-th day of study (L stands

Fig. 1 .

 1 Fig. 1. Timeseries of 3 mood features (danceability, valence and energy) considering all the target countries.

  (a) Breakpoint of the F U KR D timeseries. (b) Breakpoint of the F U KR V timeseries. (c) Breakpoint of the F U KR E timeseries.

Fig. 2 .

 2 Fig. 2. Results of the the homogeneity test made on the 3 mood-related timeseries of Ukraine. Each breakpoint is depicted as a vertical dashed red line, whereas the mean of the corresponding mood feature of each half is represented by yellow and green horizontal lines, respectively.

  (a) Cluster of the countries before the war outrage based on the danceability feature. (b) Cluster of the countries after the war outrage based on the danceability feature. (c) Cluster of the countries before the war outrage based on the valence feature. (d) Cluster of the countries after the war outrage based on the valence feature. (e) Cluster of the countries before the war outrage based on the energy feature. (f) Cluster of the countries after the war outrage based on the energy feature.

Fig. 3 .

 3 Fig. 3. Clustering of the target countries based on the DTW algorithm applied on their mood-related timeseries.

Fig. 4 .

 4 Fig. 4. Evolution in the total number of daily song stream per country. The vertical dashed line indicates the beginning of the war (February 24th, 2022).

https://charts.spotify.com/charts/overview/global

https://developer.spotify.com/discover/

https://taskandpurpose.com/news/russian-military-equipment-whitemarkings/

CODE AVAILABILITY

The source code of the project is available at https://github. com/fterroso/war and music