Molecular dynamics simulation of surface phenomena due to high electronic excitation ion irradiation in amorphous silica

We studied by means of an atomistic model based on molecular dynamics the thermal evolution of surface atoms in amorphous silica under high electronic excitation produced by irradiation with swift heavy ions. The model was validated with the total and differential yields measured in sputtering experiments with different ions and ion energies showing a very good quantitative prediction capability. Three mechanisms are behind the evolution of the surface region: (1) an ejection mechanism of atoms and clusters with kinetic energy exceeding their binding energy to the sample surface, which explains the experimentally observed angular distributions of emitted atoms, and the correlation of the total sputtering yield with the electronic stopping power and the incidence angle. (2) A collective mechanism of the atoms in the ion track originated by the initial atom motion outwards the track region subsequently followed by the return to the resulting low-density region in the track center. The collective mechanism describes the energy dissipation of bulk atoms and the changes in density, residual stress, defect formation and optical properties. (3) A flow mechanism resulting from the accumulation and subsequent evolution of surface atoms unable to escape. This mechanism is responsible for the crater rim formation. Evolution of the molecular dynamics simulation of amorphous silica for 12 keV/nm electronic stopping power with a hot cylinder of radius 3.0 nm and an incidence angle of 20º Evolution of the molecular dynamics simulation of amorphous silica for 12 keV/nm electronic stopping power with a hot cylinder of radius 3.0 nm and an incidence angle of 20º


Introduction
High electronic excitation induced in dielectric materials, such as silica and quartz, by ions with energies exceeding 0.1 MeV/amu has gained interest for a number of new applications, e.g., waveguide manufacturing, etching techniques and fine tuning of physical properties [1][2][3][4][5][6][7]. High electronic excitation also plays a role in a number of important phenomena in new irradiation facilities, e.g., nuclear fusion research facilities, large accelerators and intense laser facilities [8][9][10][11][12] and it is the basis for studies of interstellar objects wandering in space for millions of years before impacting on Earth [13,14]. Ions in this regime transmit their kinetic energy to the electronic system of the irradiated materials through ion-electron collisions that produce high energy electrons (δ-electrons). Subsequently, a myriad of events takes place in a sub-ps timescale: a e-mail: fra.sperez@gmail.com (corresponding author) impact ionization, carrier scattering, Auger recombination and exciton formation [15][16][17][18]. Energy stored in the electronic system is eventually transmitted to the atomic system at longer times (from tens to hundreds of ps) [15]. We can consider that the energy transfer to the lattice and the subsequent evolution takes place in a much longer timescale than the initial electronic evolution. This evolution eventually leads to electronic sputtering, defect formation and changes in the physicochemical properties of the material [10,12,15,19] and it is the object of our study. These final changes can be characterized by several experimental techniques [10][11][12][19][20][21][22][23] and described with different theoretical models [23][24][25][26][27][28][29][30][31][32]. Unfortunately, no model can accurately describe all the processes due to the disparate timescales involved. Thus, a description of the electronic evolution can be done through ab initio models, appropriate for sub-ps times and small systems (thousands of atoms) [33][34][35]. With a lesser degree of accuracy, kinetic models can follow the electronic evolution for longer timescales and larger systems [36]. On the other hand, atomistic models based on molecular dynamics are useful to describe hot lattice evolution and therefore, effects such as amorphization, defect production, modification of properties or emission of surface atoms [37][38][39][40][41][42][43].
We discussed ion irradiation-induced high electronic excitation effects in bulk silica in a previous paper [44] with the aid of an atomistic model that provides quantitative results on track formation and explains the observed density changes, stress generation, defect production and, consequently, the modification of the optical properties. The model was successfully validated with experimental results.
Surface effects on amorphous silica were experimentally characterized by atom emission distributions [19-22, 37-40, 45, 46]. These works analyzed the influence of the electronic stopping power, ion charge, incidence angle, polar and azimuthal emission angle, etc. A thermal emission was inferred from the results, as opposed to the jet emission observed in other materials, like LiF [19-22, 37-40, 45, 46]. In addition, the final crater structure was observed at different fluences [47,48] and it was related to the density variations and residual stress [5]. However, the mechanisms responsible for these surface effects are not well described. Theoretical models were developed to this aim, but they provided only a partial explanation [19,25,32,46].
In this work, we study silica surface emission induced by irradiation with ions of high electronic stopping power (S e > 4 keV/nm). We use an atomistic model to quantify the atom emission induced by electronic sputtering and to study the resulting crater formation. The model was validated with experimental results obtained with different ions over a broad range of ion energies. We ascribe the observed surface evolution to a thermal process consisting of three mechanisms: (1) an ejection mechanism characterized by a dissipation of energy based on atom emission; (2) a flow mechanism produced by the surface atom energy dissipation and (3) a collective mechanism of energy dissipation in bulk silica by radial displacements.

Swift ion-induced atom emission experiments on silica
In this section, we briefly summarized the available experimental data on surface emission upon swift ion irradiation. The total atom emission yield is directly measured by means of in situ ERDA (elastic recoil detection analysis) experiments during irradiation [49]. Experiments to determine differential emission yield are based on a catcher to capture the emitted atoms [19-23, 46, 50]. A schematic drawing of a catcher for this kind of experiments is shown on Fig. 1a. The catchers are analyzed ex situ with different techniques, e.g., ion beam techniques as ERDA, RBS (Rutherford backscattering spectrometry) or NRA (nuclear reaction analysis).
A few experimental campaigns were devoted to the study of emission from silica surface upon swift ion irradiation. For modeling purposes, the most interesting ones were performed with the ion charge state in equilibrium, [21-23, 49, 50] because in this situation the electronic stopping power (S e ) can be estimated with the SRIM code [51]. The results show that the emitted oxygen to silicon ratio is 2:1, i.e., the stoichiometric relation of the material is preserved [23]. In the case of glancing incidence (α < 30º), one can find experiments carried out with ions of stopping powers in a broad range (up to S e ∼ 20 keV/nm) [22,23,49]. There are fewer reported experiments for higher incidence angles due to the significant difficulties associated with those experiments. However, some experiments are available, even at normal incidence (α = 90°), which requires, a special annular detector located a few mm over the irradiation spot [21,50]. As shown in Fig. 2, the total emission yield becomes important only above a certain threshold value for electronic stopping power. The relation between the total emission yield (Y ) and the stopping power follows the expression: where E is a fitting parameter and n ∼3, typically, according to experiments [24]. The incidence angle plays a very important role, as shown in Fig. 2. The analysis of series obtained with constant stopping power (ion energy) and varying incidence angle, indeed, show that the lower the incidence angle, the higher the emission that is measured. In this case, we can describe the dependency through the following expression: where Y 90 is the total yield in normal incidence experiments and d ∼1.65 according to Ref. [23]. We fitted the experimental values obtained with α = 19°-20°with the aid of Eq. (1). The result is the continuous line in Fig. 2. Then, by means of Eq. (2), we obtained a curve for Y 90 (purple dashed line). The few experimental points obtained for α = 90°fairly follow the curve, as expected.
Once we obtained a curve for Y 90 , we can find the yield for any incidence angle by means of Eq. (2), as shown for α = 25°-40°(dashed curves). This way, we were able to predict the total yield for any incidence angle using the existing experimental values. A few experimental series were carried out to measure the differential emission yield dY/dΩ as a function of the emission angle θ [19,23]. In all these cases, a catcher perpendicular to the incidence plane was used, schematically shown in Fig. 1a. In this configuration, a symmetric emission is observed, as expected due to the symmetry of the arrangement. In Fig. 3b, the symmetric emission is shown for silica experiments irradiated with 190 MeV Au ions (S e = 17.1 keV/nm) at an incidence angle α = 19°. The differential yield in this  [19,23]; α is the ion incidence angle and θ the emission angle. b Schematic 3D representation of the simulation box showing the typical emission cone characterized by a polar emission angle (ϕ) and an azimuthal emission angle (φ). c Top view of the simulation box. The incidence angle (α) of the ion beam is arbitrary. Thus, the generated hot cylinder is inclined. The ellipse is the intersection of the inclined hot cylinder and the sample surface. A thermal bath at 300 K is set at the 4 lateral sides. d Side view of the simulation box. The height is 21.6 nm, which includes the silica atoms in the center (10 nm high) and two regions of 5.8 nm for evolution of the emitted atoms. The atoms that reach the control surfaces located 5.8 nm over the sample top and back surfaces are removed from the simulation configuration can be described by the expression: where A and b are fitting parameters and the exponent b ≥ 1.

Atomistic model
Surface effects produced by ion irradiation-induced high electronic excitation on amorphous silica were simulated with an atomistic model based on molecular dynamics (MD). The atomistic model was described in our previous paper [44] to study the bulk properties of silica under ion irradiation. Some changes were included in this work to account for surface effects. The simulations were carried out with the parallel computing code MDCASK. We used the Feuston-Garofalini [52] inter-atomic potential modified by the Ziegler-Biersack-Littmark (ZBL) potential to represent close range interactions [53]. The silica samples were rectangular prisms containing more than 1.1 × 10 6 atoms (60 × 30 × 10 nm 3 ) in a simulation box larger than the sample to provide free space for atom emission. The simulation times reached 100 ps with a time step of 0.5 fs, found to be appropriate to avoid convergence problems [44,53]. Simulations were performed on the supercomputer CESVIMA-MAGERIT-2 (UPM) using 512 processors, which provided a computation time per run of ∼20 h. The initial silica sample was made from a cristobalite structure. It was heated to 7000 K during 25 ps. Then, it was cooled down to 300 K by a repetitive process consisting of reducing the temperature in steps of 1000 K, keeping each temperature constant for 25 ps. This process was performed with periodic boundary conditions (PBC) as described in the previous paper [44] and elsewhere [54,55]. Finally, the temperature was maintained at 300 K for 25 ps with non-PBC in Z direction, in order to obtain two free surfaces. The box center coincided with the origin of the coordinate system. Such a thin box is appropriate to simulate thin samples irradiated with ions impinging on the front surface parallel to the plane XY. In such a configuration, both the front and back surfaces parallel to the plane XY are emission surfaces. Since we are interested in surface effects, using a thin box is a good way to reduce the number of atoms and, therefore, the computation time. In addition, the emission from both surfaces is equivalent due to the symmetry of the system. Thus, combining the emission results from both surfaces is a way to improve the statistics. The sample length along the Z direction (normal to the irradiation surface) was 10 nm, chosen to be as small as possible but in such a way that the emission from both surfaces is not affected by each other as confirmed by studying the initial and final positions of every atom. We assume that the energy is transferred from the incoming ion to electrons in silica in a very short time. Therefore, we ignore the details of the process. We simply define a cylindrical region (hot cylinder) around the ion trajectory and increase the total kinetic energy of the atoms in the hot cylinder according to the electronic stopping power (S e ) of the ion. The hot cylinder radius (a = 3.0 nm) is related to the deposited energy density, which, in turn, is related to the specific energy of the ion [44]. The hot cylinder axis coincides with the irradiation direction, and it is defined by the incidence angle (α), as indicated in Fig. 1. In order to simulate the energy dissipation from the ion track to the surrounding massive material, we set a thermal bath reset to 300 K every 10 steps, i.e., 5 fs, by rescaling the velocity of the atoms in a layer of 0.716 nm at the four lateral surfaces [44]. The ion arrival was simulated with PBC for lateral surfaces and non-PBC otherwise. In order to study atom emission, we let emitted atoms and atom clusters move freely in a 5.8 nm thick empty layer around the (front and back) emission surfaces (Fig. 1). The sputtered atoms are characterized by their polar (ϕ) and azimuthal (φ) angles, see Fig. 1b. Once single atoms (or the first atom of a cluster) reach a control surface (5.8 nm above the front surface or below the back surface), we consider that the atom or cluster has escaped the sample and will keep on moving in a rectilinear trajectory. Therefore, these atoms or clusters are removed from the simulation. After the simulation, we can calculate the total yield as the total number of atoms reaching the control surfaces from the front and back surfaces divided by two (thus, every simulation is equivalent to irradiation with two ions). To obtain the differential emission, we follow the trajectories of the emitted particles approaching the control surfaces and calculate the trajectory coordinates in the spherical coordinate system defined in Fig. 1b by means of least squares minimization using a singular value decomposition algorithm. Note that the trajectories of free atoms are straight lines; however, atoms forming part of clusters follow a helical orbit around the straight trajectory of the center of mass (Fig. 4). This trajectory is the relevant one because it reproduces the experimental trajectory of emitted clusters and it can be well estimated by the fitting procedure. The problem of these fits is that, due to in-flight collisions, sudden trajectory changes take place. The goodness of the fit makes possible to identify and filter these events. In this paper, we excluded the trajectories with a coefficient of determination R 2 < 0.93. This value is a compromise to leave out trajectories that underwent in-flight collisions but avoiding to exclude atoms in clusters following orbits too far from the linear trajectory of the center of mass (which translates to a lower R 2 ). Following this procedure, we assign a straight trajectory to every atom (if not excluded because R 2 < 0.93). In case of single atom emission, the fit provides the straight atom trajectory and in case of cluster emission, we assign the trajectory obtained from the fit (good approach to the trajectory of the center of mass). Once we have the trajectory of every individual atom and every cluster, we can easily estimate the number of particles emitted per solid angle unit in any given direction. In particular, to obtain the differential emission (Fig. 3), we calculate the number of particles emitted in cones (solid angle ∼0.1 sr) oriented toward a catcher parallel or perpendicular to the incidence plane, as in the experiments. A schematic representation of a catcher perpendicular to the incidence plane with an emission cone forming an angle θ with the normal direction to the irradiation surface is shown in Fig. 1a.

Results and discussion
We carried out simulations as described in the previous section at incidence angles α = 20°, 25°, 40°and 90°with an initial hot cylinder of radius a = 3.0 nm and stopping powers from 4 to 20 keV/nm. The total yield obtained in each case by counting the number of atoms that reach the control surfaces is plotted in Fig. 2 (open symbols). The vast majority of the experimental data were obtained at incidence angles α = 19°-20°( fit represented by the continuous line). The agreement of our results to these experimental data is very good, especially considering the broad range of stopping powers. Fewer data (and with poorer quality) are available for higher incidence angles. However, the dashed lines obtained according to the procedure described in the previous section are a good estimate of the experimental yield. Our simulations also reproduce fairly well the expected trends shown in Fig. 2, i.e., our results on total yield follow the power law dependence of the stopping power according to Eq. (1) and the inverse sine law of the incidence angle according to Eq. (2). Moreover, our simulations also predict the oxygen to silicon ratio, ∼2, i.e., following the material stoichiometry. A significant fraction of oxygen atoms escapes as single atoms (around 1/4 of them at 12 keV/nm increasing up to around 1/3 at 20 keV/nm for an incidence angle of 20º). The rest escapes in clusters of different sizes dominated by clusters of one silicon atom and two or three oxygen atoms. For an incidence angle of 20º and stopping power of 12 keV/nm, nearly 30% of the clusters are of this type, while for a stopping power of 20 keV/nm the fraction of these clusters reduces to 20% and a higher fraction of larger clusters is observed. Increasing the incidence angle leads to a higher fraction of single oxygen atoms and prevents the formation of large clusters, in this sense, the trend is similar to reducing the stopping power. The fact that the simulations predict the Si/O ratio observed in experiments constitutes an additional validation of our model. Recall that the only free parameter in our simulations was the initial radius of the hot cylinder, set to a = 3.0 nm for all the simulations. Considering the complicated underlying physical mechanisms and the broad range of experimental conditions, it is quite impressive that our atomistic model can provide quantitative agreement with the experiments and enables the use of the model to draw a picture about the underlying mechanisms.
In addition to the total yield, we obtained from our simulations the differential yield as described in the previous section. We plotted on Fig. 3 all the differential yields obtained for two catcher orientations with respect to the incidence plane. The dashed lines are the best fit to Eq. (3). The exponent b is larger than 1, indicating a preferential emission in the direction perpendicular to the irradiation surface but not so high as to consider jet emission like in other materials, mainly alkali halides [19-22, 37-40, 45, 46]. When the catcher is oriented perpendicularly to the incidence plane, as in Fig. 1a, one expects a symmetrical differential yield distribution as observed. On the other hand, when the catcher is contained in the incidence plane, one cannot expect a symmetrical distribution in agreement with our simulations. The good agreement between our simulations at 17 keV/nm (green symbols) and the only available experimental dataset (190 MeV Au ions, S e = 17.1 keV/nm) is remarkable. Indeed, the continuous line is the fit of the experimental data to Eq. (3) and coincides fairly well with the fit of the simulation data. We consider this an additional validation of our model that makes possible to extract conclusions about the underlying mechanisms as described next.
As mentioned above, our simulations show that only oxygen atoms are emitted as single atoms, whereas silicon atoms are always emitted as part of a cluster. The most common type of cluster is the one formed by one silicon atom and two or three oxygen atoms. This can be explained by considering that oxygen atoms have fewer bonds (two bonds per atom) than silicon atoms (four bonds per atom), easing their emission as single atoms. According to our model, the fraction of single atoms increases with increase in stopping power, but is nearly independent of the incidence angle.
We present in Fig. 5, by means of snapshots, two examples of MD simulations to follow the surface evolution upon irradiation. These examples correspond to normal incidence, α = 90°, (left) and glancing incidence, α = 20°, (right). The snapshots correspond to a thin central slab (thickness 2.5 nm), i.e., the region most affected by the ion passage. By means of the snapshots, we can qualitatively understand the irradiation effects. Initially, 1 ps after depositing the energy in the hot cylinder, massive atom displacement takes place outwards the hot cylinder. Bond rupture occurs to the extent that we observe single oxygen atoms escaping from both the front and back surfaces. Subsequently, cluster emission takes place. The pressure exerted by the hot cylinder atoms outwards produces surface deformation and compression around the track. At 3 ps, the collective atom displacement generates a crater. One might be tempted to think that the crater has a cylindrical shape, but this is not the case because deep atoms can hardly reach the surface level due to lattice reconstruction. This way, the crater shape is roughly a double cone with the vertex in the box center. Note that the snapshots might be misleading in this point because, for illustrative purposes, they correspond to a very thin slab, which may look like a cylinder at a first glance. The process continues until the atom energy is not high enough to allow them to escape from the sample. The result is a well-defined crater with a very evident crater rim. The crater rim undergoes a continuous deformation as a consequence of the exerted pressure. In addition, it acts as a nucleation center for atoms escaping with glancing trajectories. The continuous trapping of these atoms increases the height of the rim. In the case of glancing incidence (right panels), the rim is very asymmetric due to the higher fraction of atoms and clusters escaping toward the + X direction than toward the − X direction. At 12 ps, the process is nearly over and only lattice relaxation is observed further on. At 50 ps, the temperature is so low that no appreciable changes occur any longer.
We can identify three mechanisms contributing to the surface evolution: (1) Ejection mechanism The observed ejection mechanism accounts for atom displacement toward the surface followed by atom emission from the surface. This Fig. 5. Molecular dynamics simulations with a hot cylinder of radius a = 3.0 nm for silica irradiated with ions of stopping power 12 keV/nm. The incidence angle α = 90°(left) or 20°(right). The side views correspond to a 2.5 nm-thick slab centered in the plane XZ for different times. The top panels are a 3D view obtained for times longer 50 ps to observe the permanent crater formed upon irradiation. Colors indicate height along Z axis from 4 nm (red) to 8 nm (yellow). Red arrows indicate the ion incidence direction is possible when the kinetic energy of atoms at the surface exceeds the binding energy. Based on our results, we can conclude that the ejection mechanism is constituted by two contributions. Firstly, the energy transmitted by the ion is enough to produce bond rupture releasing single oxygen atoms. These atoms are emitted almost perpendicular to the irradiation surface, as can be observed in Fig. 3 and for this reason b is > 1 in Eq. (3). The second contribution to the emission of atoms is a consequence of the first one. The high excitation generated leads to a secondary rupture of oxygen bonds, but these atoms are affected by a dominant radial displacement. Consequently, a broader differential sputtering yield is generated.
(2) Collective mechanism Atoms out of the surface region and inside the hot cylinder move outwards the hot cylinder dissipating their energy to the surrounding atoms. The collective motion in the radial direction produces a mechanism originated by the dissipation of the thermal energy introduced in the hot cylinder. This mechanism is responsible for mass transport, residual stress and defect generation, as explained for bulk silica in our previous paper [44].
(3) Flow mechanism The majority of the surface atoms are not emitted due to the continuous energy loss before reaching the sample surface. Part of the surface atoms do not have enough energy to produce any bond rupture. These atoms are displaced following a heat transfer process. They dissipate their energy in the axial direction, to the surface (this component is negligible in bulk process), and in the radial direction, to the bulk. The contribution of both components and the displacements in both directions produces a spiral movement, which is responsible for the energy dissipation. Thus, many atoms near the surface cannot escape and are affected by a spiral displacement that we can observe in Fig. 5.
Amorphous silica irradiated with a high electronic excitation exhibits permanent surface effects in the form of craters. The craters are characterized by a volume loss due to atom emission (ejection mechanism), a crater rim (flow mechanism) and a residual stress (collective mechanism). The final structure of the crater is classified by the electronic stopping power of the incident ion. We distinguish two different threshold values for the stopping power: to produce permanent changes in silica (damage threshold) [22,46] and to produce electronic sputtering (sputtering threshold), the first one always being lower than the second one. The damage threshold can be obtained from bulk irradiations and it is between 2 and 5 keV/nm depending on the specific ion energy (velocity effect, discussed in our previous paper [44]). The sputtering threshold is harder to obtain because it is influenced by more parameters, including incidence angle, velocity effect, surface preparation, preexisting defects, etc. An ion with electronic stopping power above the damage threshold value, but under the sputtering threshold, is able to produce a surface depression caused by the residual stress [5], but a crater rim is not formed. An ion with an electronic stopping power higher than the sputtering threshold causes a crater depression due to the volume loss generated upon atom emission (ejection mechanism), and the residual stress generated in bulk silica by energy dissipation (collective mechanism). On the other hand, surface atoms with enough energy to produce a partial rupture of their bonds are emitted and relocated around the depression, creating the crater rim. When electronic stopping power of the ion is enough to produce ablation, a hole is produced due to the large number of atoms that are emitted. The final surface structure is mainly influenced by the ejection mechanism.

Conclusions
Surface effects on amorphous silica irradiated with ions in the regime of high electronic excitation were studied with atomistic simulations based on molecular dynamics. We analyzed the evolution of atoms in the surface region upon swift ion irradiation with the aid of a simple model that provides good results as demonstrated through quantitative validation to the total and differential emission yields obtained by means of sputtering experiments. The fact that our model is able to quantitatively predict experimental results in a broad range of stopping powers indicates that the surface evolution is dominated by thermal effects. The predictive capability of our model led us to identify three mechanisms responsible for the final fate of the involved atoms: (1) an ejection mechanism originated by the gain of thermal energy followed by bond rupture and emission of single atoms and clusters; (2) a collective mechanism originated by the initial atom motion in radial direction outwards and then inwards the track region exchanging energy with the surrounding atoms and producing permanent modifications in density, defect concentration and optical properties; (3) a flow mechanism affecting surface atoms that do not have enough energy to produce bond rupture and responsible for the crater rim formation.
The methodology developed in this work can be applied to other materials that experience surface evolution upon irradiation dominated by the gain of thermal energy. Thus, it is expected that good quantitative predictions can be obtained on the important phenomenon of electronic sputtering, provided that there are appropriate inter-atomic potentials. On the other hand, for materials with a surface emission dominated by athermal effects, a more sophisticated approach than that described in this paper is expectedly needed to reach quantitative predictions and a good understanding of the involved phenomena.
The results presented in this paper are relevant in different fields. For example, the model is able to predict morphological changes at the nanometer scale (crater and crater rim formation). Further work complemented with AFM experiments of irradiated samples can provide useful information on the crater formation mechanisms as well as quantitative predictions on the crater geometry, of interest for applications based on the fabrication of nano-pores [56]. Another application of this work is related to the study of interstellar dust. In this regard, this methodology can be applied, for instance, to the analysis of iron ion irradiation of interstellar origin. These ions can reach energies of 100 MeV and larger (equivalent to a stopping power of ∼7 keV/nm in the silica surface) [57]. According to our results, this has a significant impact on the irradiated material because every impinging ion leads to an emission yield of around 1000 atoms in grazing incidence. This effect accumulated over millions of years could have a significant influence on the interstellar dust evolution [58]. The role of energetic particles in space weathering is also important for debris disks [59] and for material in our solar system [60]. Additionally, preferential ejection of isolated O would impact dust charging and astrochemistry evolution [61].
Acknowledgements This work was funded by the projects Radiafus-5 (PID2019-105325RB-C32) of Spanish Ministry of Science, Technofusion (S2018/EMT-4437) of Madrid Regional Government and Eurofusion (EH150531176). The authors acknowledge the computer resources and technical assistance provided by the Centro de Supercomputación y Visualización de Madrid (CeSViMa) CESVIMA-MAGERIT. AP acknowledges the support of FONDECYT under grants 3190123. EMB thanks support from grant ANPCyT PICTO-UUMM-2019-00048. JK was supported by the Beatriz Galindo Program (BEA-GAL18/00130) from the Ministerio de Educación y Formación Profesional of Spain Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: All the data used in this work are offered by the authors of the article to all those who request it.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.