Computing solubility and thermodynamics properties of H2O2 in water

Hydrogen peroxide plays a key role in many environmental and industrial chemical processes. We performed classical Molecular Dynamics and Continuous Fractional Component Monte Carlo simulations to calculate thermodynamic properties of H2O2 in aqueous solutions. The quality of the available force fields for H2O2 developed by Orabi&English, and by Cordeiro was systematically evaluated. To assess which water force field is suitable for predicting properties of H2O2 in aqueous solutions, four water force fields were used, namely the TIP3P, TIP4P/2005, TIP5P-E, and a modified TIP3P force field. While the computed densities of pure H2O2 in the temperature range of 253-353 K using the force field by Orabi&English are in excellent agreement with experimental results, the densities using the force field by Cordeiro are underestimated by 3%. The TIP4P/2005 force field in combination with the H2O2 force field developed by Orabi&English can predict the densities of H2O2 aqueous solution for the whole range of H2O2 mole fractions in very good agreement with experimental results. The TIP4P/2005 force field in combination with either of the H2O2 force fields can predict the viscosities of H2O2 aqueous solutions for the whole range of H2O2 mole fractions in good agreement with experimental results. The diffusion coefficients for H2O2 and water molecules using the TIP4P/2005 force field with either of the H2O2 force fields are almost constant for the whole range of H2O2 mole fractions. The Cordeiro force field for H2O2 in combination with either of the water force fields can predict the Henry coefficients of H2O2 in water in better agreement with experimental values than the force field by Orabi&English.


I. INTRODUCTION
Hydrogen peroxide, H 2 O 2 , has attracted considerable interest as it plays a key role in the oxidative chemistry of the troposphere. It can be found both in the gas and in the aqueous phase [1,2], and has several industrial [3], environmental [4], and biological [5] applications. The recombination of hydroperoxyl (HO 2 ) radicals is the most important chemical pathway leading to the production of H 2 O 2 in the troposphere [6][7][8]. Subsequently, H 2 O 2 can lead to the acidification of clouds, rain, and fog by oxidizing SO 2 and converting it into H 2 SO 4 (and to a less extent oxidizing NO 2 and converting it into HNO 3 ) [9][10][11][12][13]. H 2 O 2 also serves as a reservoir of HO x radicals that are key oxidants in controlling the self-cleaning of the atmosphere [14][15][16].
H 2 O 2 was first synthesized by Thenard [17] by the reaction of barium peroxide with nitric acids in 1818 and is now considered an important reagent of green chemistry since it decomposes to water and oxygen as the only reaction products. This feature makes H 2 O 2 an envi-temperature (or non-thermal) plasmas [28,29] which allows H 2 O 2 production at ambient temperatures and pressures [30][31][32][33][34]. This enables direct delivery of H 2 O 2 to different substrates; even to heat sensitive substrates such as living tissues. The latter has led to biomedical applications of low temperature plasmas [35]. For such applications, it is important to know which mechanisms determine the uptake of plasma products (e.g., H 2 O 2 ) in the liquid around the cells. For this, information on solubility and thermodynamics properties of plasma products are necessary so that this can be leveraged into macroscopic plasma fluid models [36] to predict the final concentration of plasma products in the liquid phase. The motivation of this work is to provide such data for H 2 O 2 as limited data are available.
Various computational studies have been carried out which shed light on structural properties of H 2 O 2 monomers as well as its clusters, torsional barrier energies, and vibrational-rotational energy levels [53][54][55][56][57] using quantum mechanical approaches. Structure and dynamics of H 2 O 2 in water were also investigated using quantum mechanical methods in Refs. [58][59][60].
In this work, we use force field based Molecular Dynamics (MD) and Continuous Fractional Component Monte Carlo (CFCMC) simulations with the purpose of obtaining solubilities and thermodynamics properties of H 2 O 2 in water, for the first time, in a systematic manner such that the quality of the available force fields for H 2 O 2 are assessed.
Although several force fields are available for H 2 O 2 [61][62][63][64], only a few of them have been parameterized with respect to the interactions between both H 2 O 2 -H 2 O 2 and H 2 O 2 -H 2 O. One is the ABEEM/MM, the atom-bond electronegativity equalization fluctuating charge molecular force field [65,66], which is computationally very expensive due to its complex potential energy functional form [65,66]. A simple additive potential model for H 2 O 2 was proposed by Orabi & English [67] which was parameterized to account for interactions of H 2 O 2 with itself and with water. The model was calibrated with regard to the experimental density and heat of vaporization of pure liquid H 2 O 2 at 0°C, and was able to reproduce the experimental diffusion coefficient at 0°C and the heat capacity at 25°C of liquid H 2 O 2 . With a combination of the modified TIP3P water force field [68], the H 2 O 2 force field could predict the experimental hydration free energies and densities of aqueous [67]. Another force field parametrization is from the work of Cordeiro [69], wherein the bonded interactions were obtained from ab intio quantum calculations [54,70], and the Lennard-Jones parameters and partial charges were modified to reproduce the properties of pure liquid H 2 O 2 and its hydration free energy. This force field was used to study the distribution, mobility and residence times of H 2 O 2 at the interface of water and phospholipid biomembranes. In addition, there is another parametrization in the paper by Vácha et al. [64] in which the behaviour of H 2 O 2 at the air-water interface was investigated. The force field by Vácha et al. [64] is a rigid force field; it only includes electrostatic and van der Waals interactions which were calibrated against the experimental hydration energies of H 2 O 2 .
In this manuscript, we evaluate the quality of the force fields which were developed by Cordeiro [69], and Orabi & English [67] [68,71] as it performs better in calculating the specific heats of water [72], TIP5P-E as it can capture the thermal conductivities of water [72] and TIP4P/2005 [73] as it can predict the densities and self-diffusion coefficients of water with commendable accuracy [73]. The fourth water force field is a modified version of TIP3P (mTIP3P) [68], which was used in the work by Orabi & English [67]. The results are compared with experimental values. Finally, we compute the Henry coefficients of H 2 O 2 in water at 300 K.
The rest of this manuscript is organized as follows. In section II, details of the force fields which were developed by Cordeiro [69], and Orabi & English [67] are provided, and the MD and CFCMC simulations are described. The results are presented and discussed in section III. Finally, concluding remarks are presented in section IV.

A. Force Fields
Both force fields developed by Cordeiro [69], and Orabi & English [67] for H 2 O 2 consider a non-rigid H 2 O 2 molecule, that is, they incorporate bonds, angles and dihedrals information with the van der Waals (vdW) and electrostatic interactions. The total potential energy  [69], and Orabi & English [67]. E bonds , E angles , E dihedrals , E vdW , and E electrostatic represent the stretching, bending, torsional, van der Waals and electrostatic energies, respectively. The definition of parameters is explained in the text (see section II A). The parameters are provided in Tables S2 and  S3 of the Supplementary Information where E bonds , E angles , E dihedrals , E vdW , and E electrostatic are presented in Table I for both force fields. The bonded interaction parameters (E bonds , E angles , and E dihedrals ) listed in Table I have the following definitions: b is the bond distance, θ is the bond angle, ϕ is the dihedral angle, δ is the multiplicity factor, and ψ is the supplementary angle of ϕ. k b , k θ , k ϕ are the force constants of the bond stretching, angle vibration, and dihedral potentials. b 0 and θ 0 represent the equilibrium bond distance and bond angle, respectively. C n with n ranging from 0 to 5 represents the coefficients for the Ryckaert-Bellemans dihedral potential [74]. q represents the atomic partial charges of the electrostatic energy (E electrostatic ) term. A Lennard-Jones (L-J) potential is used for the long-range van der Waals interactions, in which σ represents the distance at which the particle-particle interaction energy is zero, and ϵ represents the depth of the potential well. The mixing rules for the L-J parameters for two dissimilar non-bonded atoms are given by Lorentz-Berthelot [75] [ for the force field by Cordeiro. The values of these parameters are provided (using the GROMACS convention) in Tables S1 and S2 of the Supporting Information (SI) for both force fields. The cutoff radius for Lennard-Jones and Coulombic interactions was set to 9Å. The Particle-Mesh-Ewald [76,77] method was used to treat long-range electrostatic interactions. Long-range tail corrections were applied to both energies and pressures [78]. We use three different rigid water force fields in this study, namely TIP3P, [68,71]TIP4P/2005, [73] and TIP5P-E [79,80]. We also use a modified TIP3P water force field (mTIP3P) [68] which was used in the work by Orabi & English [67]. In the remainder of this manuscript, the force field developed by Cordeiro [69] is referred to as "Cordeiro" and the force field developed by Orabi & English [67] is referred to as "Orabi". fractions of H 2 O 2 in the range from 0 to 1.0 were performed using the GROningen MAchine for Chemical Simulations (GROMACS) version 2022.4 [81][82][83][84][85]. Each system was prepared in a simulation box with an initial length of 27.6Å, containing 500 molecules. A snapshot of a simulation box containing 250 H 2 O 2 molecules and 250 H 2 O molecules is shown in Figure 1.
After energy minimization using the steepest descent algorithm followed by a conjugate gradient algorithm, the MD simulations were run for 100 ps in the constant number of atoms/molecules, volume and temperature (NVT) ensemble. The simulations were then continued in the constant number of atoms/molecules, pressure and temperature (NPT) ensemble for 25 ns. For calculating the viscosities and self-diffusivities, the simulations were continued in the NVT ensemble for another 20 ns. The temperature was kept fixed by Nosé-Hoover thermostat [86]. The Parinello-Rahman barostat [87] with a time constant of 1 ps and compressibility of 4.5 × 10 −5 bar −1 was used to keep the pressure at 1 bar. In all simulations, the Newton's equations of motion were integrated with a leap-frog [88] algorithm with a time step of 2 fs. Periodic boundary conditions were applied in all Cartesian directions. The parallel linear constraint solver (P-LINCS) [89,90] was used to constrain H bonds.

Continuous
Fractional Component Monte Carlo (CFCMC) simulations [92][93][94] using the opensource Brick-CFCMC software [94][95][96] were performed in the isothermal-isobaric (NPT) ensemble. In the CFCMC technique, fractional molecules (compared to normal or "whole" molecules) are introduced whose interactions with the rest of the system are scaled with a continuous coupling parameter λ (λ ∈ [0, 1]). The minimum value of λ (λ = 0) indicates no interactions between the fractional molecule and the rest of the molecules in the system (i.e., fractional molecules act as ideal gas molecules). λ = 1 represents full interactions between the fractional molecules and the other molecules in the system (i.e., the fractional molecule acts as whole molecules). The coupling parameter λ is biased with a weight function (W (λ)) using the Wang-Landau algorithm [97] to improve molecule transfers (insertions/deletions). This ensures a smooth observed probability distribution of λ. We used 100 bins to construct a histogram for the λ values and its probability of occurrence (p(λ)). The Boltzmann average of any property (A) is then computed using [98] The chemical potential of species i is calculated with respect to its ideal gas chemical potential [95] where µ ideal i and µ ex i are the ideal gas and excess chemical potential of the species i, respectively. The excess chemical potential can be related to the Boltzmann sampled probability distribution of λ by the following equation [95] where p(λ=1) and p(λ=0) are the Boltzmann sampled probability distributions of λ at 1 and 0, respectively. k B is the Boltzmann constant, and T is the absolute temperature. The excess chemical potential at infinite dilution (µ ex,∞ ) can be used to determine the Henry volatility coefficient [99] where ρ is the number density of the solvent. This yields the Henry volatility coefficient (K px v ) in units of [Pa]. The Henry coefficient (H cp s ) in units of [mol/m 3 Pa] can be obtained using the following conversion: in which ρ H2O is the density of water, and M H2O is the molar mass of water [101].
CFCMC simulations contained 300 water molecules in a cubic simulation box with initial dimensions of 21Å. A single fractional molecule of H 2 O 2 was introduced to calculate its excess chemical potential. The cut-off radius for the intermolecular L-J and Coulombic interactions was set to 9Å. The Ewald summation [102] method was used for calculating electrostatic interactions. Longrange tail corrections were applied to the L-J potential. Periodic boundary conditions were applied in all directions.
For CFCMC simulations, 1,000 initialization cycles were carried out followed by 5 × 10 6 equilibration cycles and 5 × 10 6 production cycles. One cycle refers to N number of trial moves, where N is the total number of molecules. Trial moves were selected with the following probabilities: 32% translation moves, 22% rotation moves, 1% volume changes, 5% each of bending and torsion moves, 25% λ changes and 10% hybrid moves that combined swap and identity change moves [95]. Three independent simulations were performed for each combination of water force field and H 2 O 2 force field to obtain an average value and the standard deviation for the Henry coefficients.

A. Densities
The densities of anhydrous H 2 O 2 for a temperature range of 253 K to 353 K (in steps of 20 K) for both the Orabi and Cordeiro force fields are plotted in Figure 2. We used the gmx density tool to compute the average density of each system. The experimental values are shown in black circles [42]. The melting point and boiling point of H 2 O 2 are reported as 272.74 K and 423.15 K, respectively [103]. While the Cordeiro force field underestimates the densities of anhydrous H 2 O 2 by about 3% compared to the experimental values, the densities of anhydrous H 2 O 2 using the Orabi force field are in excellent agreement with the experimental values.
Next, we evaluate which water force field is suitable for predicting the densities of H 2 O 2 aqueous solutions. To this end, we modelled systems of H 2 O 2 aqueous solutions for the whole range of H 2 O 2 mole fractions (0 to 1.0) at T = 298 K and 1 bar using the Orabi or Cordeiro force fields in combination with four different water force fields: TIP3P, TIP4P/2005, TIP5P-E, and the modified TIP3P water force field (mTIP3P) which was used in the work by Orabi & English [67]. The choice of temperature at 298 K and pressure at 1 bar was motivated by the availability of experimental results by which we could validate our models. The results as a function of the H 2 O 2 mole fraction are shown in Figure 3. The experimental values [50] are added for comparison.
The densities of pure water (i.e., a mole fraction of zero), using the four water force fields are in good agreement with the reported data at 298 K and 1 bar [73,104].  We conclude that the Orabi force field is a better force field than the Cordeiro force field for predicting the densities of pure H 2 O 2 in the temperature range of 253 -353 K. In addition, the TIP4P/2005 or the mTIP3P force field in combination with the Orabi force field predicts the densities of H 2 O 2 aqueous solutions for the whole range of mole fractions (0 -1.0) in very good agreement with the experimental values.

B. Viscosities
We used the gmx energy tool to compute the viscosities [105]. This tool uses the Navier-Stokes equation in which an external force is applied to the system. This causes a velocity gradient in the system from which the viscosity can be calculated [105]. The viscosities of H 2 O 2 aqueous solutions for various H 2 O 2 mole fractions (0 to 1.0) at T =293 K were computed. Figure S1 shows the values of viscosity using the Orabi (a) and Cordeiro (b) force fields in combination with the TIP3P, mTIP3P and TIP4P/2005 force fields. The results including the TIP5P-E water force field are shown in Figure S1 of the Supplementary Information. The standard deviation is used to estimate error bars. The experimental values at 293 K are included for comparison [43].
The viscosities of pure water (i.e., mole fraction = 0 ) are in good agreement with the computed values using the TIP3P, mTIP3P, TIP4P/2005, and TIP5P-E force fields [106]. The experimental value of the viscosity of pure H 2 O 2 at 293 K is 1.25 mPa s [43]. The computed value is 1.36 mPa s by using the Orabi force field, and is 1.34 mPa s by using the Cordeiro force field.
The Cordeiro [69] force fields in combination with the TIP3P [71], mTIP3P [68] and TIP4P/2005 [73]. The results including the TIP5P-E water force field are shown in Figure S1 of the Supplementary Information. Error bars are estimated based on the standard deviation. The experimental values [43] at 293 K are added for comparison.
estimates. The combination of Cordeiro force field with the mTIP3P or TIP3P or TIP4P/2005 water force fields follows a similar trend. The Cordeiro-mTIP3P and Cordeiro-TIP3P models underestimate the viscosities up to a mole fraction of 0.8, above which it slightly overestimates. The Cordeiro-TIP4P/2005 model, however, underestimates the values of viscosities by 7% up to a mole fraction of 0.5 while it overestimates by ca. 5% for H 2 O 2 mole fractions higher than 0.8. Contrary to the other water force fields, the TIP5P-E water force field in combination with either the Orabi or Cordeiro force fields predicts a relatively high peak in viscosity at the intermediate mole fractions (mole fraction of 0.5), see Figure S1 of the Supplementary Information. This may be due to structural changes which the TIP5P-E water force field induces in the system. This is addressed in section III D using radial distribution functions.
We conclude that the TIP4P/2005 water force field in combination with the Orabi force field or the Cordeiro force field predicts the viscosities of H 2 O 2 aqueous solutions in better agreement with the experimental values.

C. Self-diffusion coefficients
Diffusion coefficients were calculated from the meansquared displacements (MSD), and were corrected for finite-size effects with the Yeh-Hummer equation [107,108] where D and D MD denote the diffusion coefficient calculated with and without the finite-size effects corrections, respectively. k B is the Boltzmann constant, T is the absolute temperature (in K), ξ is a dimensionless number which for a cubic simulation box is equal to 2.837, L is the length of the cubic simulation box, and η is the viscosity of the system. We used the gmx msd tool to obtain the MSD as a function of time. D MD is obtained by fitting the MSD to where d = 3 is the dimension of the system. The selfdiffusion coefficients were calculated from 1 ns to 20 ns NVT trajectories. Figure S2 of the SI shows an example of MSD versus time on a logarithmic scale. Figure  5 shows the self-diffusion coefficients of H 2 O 2 and water in aqueous H 2 O 2 solutions for the whole range of hydrogen peroxide mole fractions (0 to 1.0). The self-diffusion coefficients of pure water (i.e., mole fraction=0) for the four water force fields are in good agreement with the values reported in Ref. [73,79,109]. The TIP5P-E and TIP4P/2005 water force fields predict the value of the self-diffusion coefficient in better agreement with the experimental value (2.3×10 −9 m 2 /s [110]). The self-diffusion coefficients of both the water and the H 2 O 2 molecules decrease monotonically by increasing the mole fraction of H 2 O 2 using the Orabi-TIP3P or Orabi-mTIP3P models. There is a similar trend for the Cordeiro-TIP3P and Cordeiro-mTIP3P models. The TIP4P/2005 water force field in combination with either the Orabi force field or the Cordeiro force field predicts a relatively constant self-diffusion coefficient for both water and H 2 O 2 for the whole range of mole fractions. This is in agreement with a recent experimental study [33] where it was concluded that the self-diffusion coefficients of H 2 O 2 in solutions are insensitive to its concentration. The TIP5P-E water force field in combination with either the Orabi or Cordeiro force fields predicts a minimum at a mole fraction of 0.5 for the self-diffusion coefficients of both water and H 2 O 2 . This is correlated with its very high value of viscosities (see Figure S1 of the Supplementary Information).

D. Radial Distribution Functions
Structural properties of H 2 O 2 aqueous solution at various mole fractions were investigated using the radial distribution functions (RDF). Note that there are 4 atom types in each system: hydrogen of water (H w ), oxygen of water (O w ), hydrogen of H 2 O 2 (H p ), and oxygen of The RDFs for O w and H w in pure water using the mTIP3P, TIP5P-E and TIP4P/2005 water force fields are shown in Figure S3 (a). The results show a first peak at approximately 0.18 nm, and a second peak at approximately 0.32 nm. The peak heights slightly differ between the water force fields with the TIP4P/2005 predicting a higher value followed by the TIP5P-E, and then the mTIP3P predicting a smaller value.
The respective RDFs in H 2 O 2 aqueous solutions using the Orabi force field in combination with the three water force fields for systems with mole fractions of 0.1, 0.5 and 0.9 are also shown in Figure S3 (b, c, d). The RDFs for O w and H w for systems using the Cordeiro force field in combination with the mTIP3P, TIP5P-E and TIP4P/2005 are shown in Figure S3 of the SI. In the system with a mole fraction of 0.1, the position of the peaks is independent of the water force fields. As the mole fraction of H 2 O 2 is increased to 0.5 and 0.9, the position of the peaks does not change in systems where the mTIP3P or TIP4P/2005 force field is used. In the systems where the TIP5P-E water force field is used, however, the RDF changes: the height of the peaks become smaller and an additional structural correlation appears , and x = 0.9 (d) at 298 K and 1 bar using the Orabi [67] force field in combination with the mTIP3P [68], TIP5P-E [79,80], and TIP4P/2005 [73] water force fields, where x is the mole fraction of H2O2. The first peak (at ca. 0.19 nm) was removed to distinguish the differences between the combinations clearly. The RDF for pure H2O2 is plotted in (d) using the Orabi force field and the Cordeiro force field.
in the system with a mole fraction of 0.5. This additional correlation between the water molecules persists till 0.8 nm, whereas for the systems in which the mTIP3P or TIP4P/2005 is used, the structural correlation persists only till 0.6 nm. In the Orabi-TIP5P-E system with a mole fraction 0.9, the first two peaks disappear. A similar trend can be observed for the combinations involving the Cordeiro force field with the mTIP3P, TIP5P and TIP4P/2005 water force fields (see Figure S3 of SI). In the Cordeiro-TIP5P-E combination with a mole-fraction of 0.5, however, the structural correlation between the water molecules is stronger than that of the corresponding Orabi-TIP5P-E combination. This can be seen from a more prominent peak after 0.4 nm compared to the other models.
The RDFs for O p and H p in H 2 O 2 aqueous solution with mole fractions of 0.1, 0.5 and 0.9 using the Orabi force field in combination with the three water force fields are shown in Figure S5 (a, b, c). Similarly, the respective RDFs using the Cordeiro force field in combination with the three water force fields are shown in Figure S4 of SI. The RDF of O p and H p in pure H 2 O 2 using the Orabi and the Cordeiro force fields are also shown in Figure S5 (d). The first peak in the RDF has a large amplitude, therefore we removed it to be able to distinguish the differences between the systems more clearly (see Figure S5 of SI). RDFs of the system at a mole fraction of 0.9 are almost identical using the three different water force fields with a second peak at 0.37 nm. By decreasing the mole fraction of H 2 O 2 to 0.5, and 0.1, the RDFs remain the same for the systems in which the mTIP3P-E or the TIP4P/2005 is used. For the system in which the TIP5P-E water force field is used, however, the RDF changes drastically. This is also the case with the Cordeiro-TIP5P-E model. The RDF for pure H 2 O 2 using the Orabi force field is almost identical as the system using the Cordeiro force field.
The RDFs for O p and H w , O p and O w , and O w and H p using the Orabi force field with the three water force fields (mTIP3P, TIP4P/2005, and TIP5P-E) for solutions with mole fractions of 0.1, 0.5 and 0.9, are shown in Figure S6. Likewise, RDFs with the Cordeiro force field and the water force fields are shown in Figure S6 of SI. In systems where the mTIP3P and TIP4P/2005 water force fields were used, RDFs have the same structure in which the position of the first peak is in good agreement with X-ray measurements on crystals of H 2 O 2 ·2H 2 O [111] and simulation results [67]. On the contrary, the structural properties in systems where the TIP5P-E force field was used, have changed. A comparable effect can be seen in Figure S6 of SI where the Cordeiro-TIP5P-E combination is used.
The number of water molecules in the micro and first solvation shells of H 2 O 2 molecule were obtained by integrating up to the first and second minima of the RDF for O p -O w , respectively. The results are shown in Table  S3 and S4 of the SI. Orabi & English [67] reported 3.0 and 9.4 water molecules in the micro and first solvation shells, respectively, for a single peroxide molecule in 500 water molecules using the mTIP3P water force field. Authors in Ref. [58] report 6.0 water molecules in the first solvation shell of H 2 O 2 using a hybrid quantum-classical simulation. According to our results, the number of water molecules in the first solvation shell at mole fractions of 0.1 and 0.9 are lower for systems where the TIP5P-E water force field is used. At a mole fraction of 0.5, however, the number of water molecules in the first solvation shell slightly increases.
Our results suggest that the addition of the TIP5P-E water force field to the H 2 O 2 Orabi or Cordeiro models disturbs the structural properties of the systems such that there is a stronger interaction between the water molecules and H 2 O 2 . This effect is more pronounced at a mole fraction of 0.5 which is correlated with the prediction of the Orabi-TIP5P-E and Cordeiro-TIP5P-E models for densities, viscosities, and self-diffusion coefficients (see Figure 3, S1 and 5).

E. Hydrogen Bond analysis
The number of hydrogen bonds (calculated as the summation of hydrogen bonds between H 2 O 2 -H 2 O 2 , H 2 O 2water and water -water) per H 2 O 2 molecule for combinations of the Orabi and Cordeiro force fields with the water force fields is shown in Figure 9. We used the geometric criterion for hydrogen bonds proposed in Ref. [112]. The number of hydrogen bonds for both the Orabi and Cordeiro force fields in combination with the TIP4P/2005, TIP3P and mTIP3P water force fields exhibit a steady increase until a mole fraction of 0.9. Existing literature indicates the existence of about 4 hydrogen bonds between hydrogen peroxide and water molecules [58,60]. For pure H 2 O 2 systems, the number of hydrogen bonds sharply decreases to about 5 for both the Orabi and Cordeiro systems. The Orabi-TIP5P-E and Cordeiro-TIP5P-E combinations are however different compared to the others. These combinations indicate a minimum in the number of hydrogen bonds at a mole fraction of 0.4 (around 3 hydrogen bonds for Orabi-TIP5P-E and 2 for Cordeiro-TIP5P-E). This minima for the number of hydrogen bonds at the intermediate concentrations for both these systems are concurrent with the high viscosities and low diffusion seen earlier. Analysis of the RDFs, solvation shells and the number of hydrogen bonds of the TIP5P-E systems indicate a variation in the arrangement of water molecules around a H 2 O 2 molecule. The effects of this structural difference is reflected in the incongruity of the viscosities and selfdiffusion coefficients of the TIP5P-E systems with respect to other water force fields based systems.  [67] and (b) Cordeiro [69] force fields in combination with the TIP3P [71], mTIP3P [68], TIP4P/2005 [73] and TIP5P-E [79,80] water force fields.

F. Henry coefficients
The Henry coefficients were computed for H 2 O 2 in water using the Orabi and Cordeiro force fields in combination with the various water force fields. It should be noted that we have not further considered the TIP5P-E water force field for the solubility calculations as it has been ascertained from the earlier sections that neither the Cordeiro nor Orabi force fields in combination with the TIP5P-E water force field could accurately predict the densities, viscosities and self-diffusion coefficients of H 2 O 2 aqueous systems. The results of the solubility calculations are provided in Table S6. The reported values range from 670 to 1400 [mol/m 3 Pa] [113][114][115][116]. It is evident that the Cordeiro force field in combination with the TIP3P or mTIP3P water force fields predicts the Henry constants that are within the range of the reported experimental values.

IV. CONCLUSIONS
We performed MD and CFCMC simulations to study thermodynamics properties of aqueous solutions of H 2 O 2 . The quality of the available force fields of H 2 O 2 , Cordeiro [69] and Orabi [67] [33]. We studied the structural properties of H 2 O 2 aqueous solutions using radial distribution functions. These results suggest that the use of the TIP5P-E water force field in combination with either the Orabi or Cordeiro force field, predict a stronger interaction between water molecules and H 2 O 2 molecules. Hydrogen bond analysis indicates a steady increase in the number of hydrogen bonds per H 2 O 2 molecule with increasing solute concentration for H 2 O 2 aqueous solutions. The Cordeiro-TIP5P-E and Orabi-TIP5P-E systems exhibited a minimum at the intermediate solute concentrations. This is in line with the deviation in the dynamic properties (viscosities and self-diffusion coefficients) of these systems. Finally, we computed the Henry coefficients of H 2 O 2 in water. The values using the Cordeiro force field in combination with either of the TIP3P or mTIP3P water force fields are within the range of experimental values. The quantitative data presented in this work can be used by macroscopic plasma fluid models to determine the uptake of H 2 O 2 from the gas phase plasma by liquid [36] or to interpret and complement experimental findings [117].

AUTHOR CONTRIBUTIONS
TS carried out the simulations and data analysis. All authors provided critical feedback on the interpretation of data analysis. BB conceived and supervised the project. TS and BB wrote the manuscript in collaboration with all the authors.

SUPPLEMENTARY INFORMATION
Force field parameters for Cordeiro and Orabi & English using the GROMACS convention are provided in Tables S1 and S2, respectively; Figure S1 shows viscosities of H 2 O 2 aqueous solutions for various mole fractions of H 2 O 2 ; Figure S2 shows an example of the MSD versus time on a logarithmic scale; Radial distribution functions are illustrated in Figures S3-S6; The number of water molecules in the micro solvation shell, and the first solvation shell of H 2 O 2 in aqueous solution for various mole fractions of H 2 O 2 is shown in Tables S3 and S4, respectively.

CONFLICTS OF INTEREST
There are no conflicts to declare.