Cooling rate effects on the structure of 45S5 bioglass: Insights Cooling rate effects on the structure of 45S5 bioglass: Insights from experiments and simulations from experiments and simulations

Due to its ability to bond with living tissues upon dissolution, 45S5 bioglass and related compositions materials are extensively used for the replacement, regeneration


Introduction
Since its introduction by Hench, 45S5 bioglass remains one of the most used bioactive glass compositions [1].45S5 bioglass, henceforth referred to as bioglass, is constituted by the oxides of calcium, sodium, phosphorus, and silicon and has the capacity to form strong interfacial bonds with hard tissues in vivo through dissolution of ions when they are hydrated by body fluids [1].Bioactive glass compositions have been widely used for a variety of applications ranging from scaffolds for tissue engineering, [2] bone grafts, [3] periodontal applications, [4,5] and even as a coating for bioinert implants.[6,7] Despite the wide usage of bioactive glasses, some open questions regarding the atomic-level structure and other structural features of bioglass remain, which largely arise from some inconsistencies between previously reported simulation and experiment results.This is exemplified by the fact that in bioglass, while experimental studies suggest that phosphate species exist majorly as isolated tetrahedra, [8,9,[26][27][28][29][30] atomistic simulations suggest the existence of phosphate tetrahedra with higher degrees of crosslinking [9][10][11][12]15,19,20].A detailed understanding of the structure of bioglass is essential to tailor the dissolution processes governing bioactivity.[26,27,31] Bioglass samples, like other glasses, are out-of-equilibrium systems formed by the rapid quenching of an isochemical equilibrium liquid.[32,33] Accordingly, the structure and properties of bioglass strongly depend on the thermal and pressure history they have experienced during melt-quenching or subsequent processing.The structure of bioglasses is typically characterized using techniques such as infrared spectroscopy, Raman spectroscopy, magic angle spinning nuclear magnetic resonance (MAS-NMR) spectroscopy, [26][27][28] or using computational techniques such as molecular dynamics (MD) or density function theory (DFT) simulations.[8,9,19] MD simulations are widely used as they provide some insights into the atomic details that are otherwise invisible to conventional experiments.In particular, classical MD simulations are usually preferred over DFT simulations due to the latter's limitation to small sizes (typically < 500 atoms) and timescales (typically < 10 picosecond).Although classical MD simulations are also restricted to short timescales (typically a few nanoseconds), this timescale can be four to five orders of magnitude higher than that of DFT simulations.In turn, these short timescales (both for MD and DFT) lead to unrealistically high cooling rates in simulations, such as 10 12 K/s, which is almost 10 orders of magnitude higher than typical experimental cooling rates (1-to-100 K/s).Therefore, understanding the effect of the cooling rate on the structure of bioglass provides the key for establishing a realistic, unified description of the atomic structure of bioglass.
Several earlier studies have been conducted to study the effect of cooling rates on model glasses, [34][35][36] vitreous silica, [37][38][39] sodium silicate glasses [40,41], and 45S5 bioglass [10,15] structures.The high cooling rate used in simulations typically results in more structural defects.[40,42,43] Further, the cooling rate is found to have a notable impact on some of the macroscopic properties such as density, thermal expansion, and elastic modulus [40,44]-wherein the extent of the dependence on the cooling rate depends on the glass composition.Moreover, earlier studies based on MD simulations have suggested that, although the short range order is usually weakly affected, the medium range order can present notable changes with varying cooling rates.[40] This is important as, in the case of bioglass, the medium range order is suggested to play a major role in controlling the bioactivity.[27] This raises some concerns on the validity of the structure-property relations generated by MD simulations.
Extensive studies has have been conducted in the past to analyze the structure of 45S5 bioglass .In particular, the phosphorus speciation has been debated when reconciling experimental and simulation studies.Experimental studies suggest that the majority of the phosphate species present in the bioglass exist as isolated orthophosphate units, although non-negligible amount of QP 1 species have been observed [16,20,[23][24][25][26].In some studies, 31 P MAS-NMR spectroscopy revealed a broad resonance at around 9 ppm corresponding to isolated orthophosphate anions (PO4 3-) [26,28,30].These isolated orthophosphate anions associate with metal cations, thereby reducing the network-modifying role of Ca/Na cations in the silicate network.[28] However, multiple molecular dynamics (MD) [9][10][11]13,15] and ab-initio simulation [12,20,22] studies have confirmed the existence of non-orthophosphate groups.The MD studies have been conducted using rigid-ion [15,45] and core-shell potentials [10,11], with the core-shell potential leading to lower non-QP 0 species in comparison to rigid-ion one.While there has been some cooling rate studies using MD simulations on the Q n speciation using both rigid-ion [15,45] and core-shell models [10], no clear picture of the nature of P speciation has emerged [10,11].To reconcile the experimental and simulation results, a later study suggested that, although it is possible that P may exist as Q 1 or Q 2 species, the levels of these species must be significantly lower than those suggested by MD simulations [26].
Herein, using MD simulations, we investigate the structure of bioglass samples prepared using cooling rates ranging over five orders of magnitudes.In line with previous results, [10,40] we find that the cooling rate strongly affects the medium-range order structure, whereas the short-range order remains fairly unaffected.Interestingly, we observe of an increasing propensity phosphate group to exist as isolated orthophosphate ions upon decreasing cooling rate.Finally, by extrapolating the structural features obtained from MD simulations such as Q n distribution to experimental cooling rates, we demonstrate that MD simulations can be reconciled with our NMR measurements of Q n speciation to provide a consistent structure for bioactive glasses.

Methodology Simulation details (i)
Glass simulation 45S5 bioglass samples, with the nominal molar composition (CaO)26.9.(Na2O)24.4.(SiO2)46.1.(P2O5)2.6,were prepared using classical MD simulations by applying the traditional "melt-quench" method.[40,46,47] All simulations were performed using the open-source LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) package.[48] In order to study the effect of the cooling rate, the glass samples were prepared with cooling rates varying over five orders of magnitudes, from 100 K/ps to 0.01 K/ps.The interactions between the atoms were simulated using the well-established Teter potential, [8,40] which has a Buckingham-like form.This potential has been extensively validated against experimental structure data for a wide-range of silicate and phosphate minerals and glasses, [40,44,[49][50][51] including 45S5 bioglass (see Supplementary Information and Ref. [8]).A repulsive term, V(r) = B/r n + Dr, was adopted for r smaller than r0 (variable for each of the atom pairs), to avoid the "Buckingham catastrophe" at high temperature.Here, r0 is defined as the value of r when the third derivative of potential energy approaches zero.At r0, B and n were chosen to make potential energy and its derivatives continuous.The long range Coulombic interactions were calculated with the particle-particle particlemesh (PPPM) solver algorithm with an accuracy of 10 -5 .Initial configurations were built by randomly placing 5673 atoms in a cubical box such that there is no unrealistic overlap.The system was then equilibrated for 1 ns at 3000 K and zero pressure in the isothermal isobaric (NPT) ensemble.This was to ensure that the system reaches an equilibrium state, without any memory of the initial structure obtained by the random placement of atoms.After equilibration, glasses were formed by gradually cooling the systems from 3000 to 300 K linearly at five different cooling rates, namely, 100, 10, 1, 0.1, and 0.01 K/ps under zero pressure in the NPT ensemble.Note that continuous cooling, with no intermittent equilibration, was carried out from 3000 to 300 K using the Nose-Hoover thermostat [52].Finally, all the glass structures were relaxed for 100 ps at 300 K and zero pressure in the NPT ensemble.For statistical averaging, the final glass structure obtained was equilibrated in an NVT ensemble for 100 ps.All the structural and thermodynamic data of the final glass structure were collected during this run.All the properties referring to the "glassy state" are obtained by averaging over 100 frames extracted from this run.

(ii)
Pair distribution function The pair distribution function (PDF) represents the ratio of the local density of atoms with respect to the global density of the atoms in terms of distance from an atom.In particular, the neutron pair distribution function is given by:   () = where   is the fraction of  atoms ( = Na, Ca, Si, P, or O),   is the neutron scattering length of the species, and   are the partial PDFs.Using the partial PDF, the coordination number of each of the species can be obtained by computing the number of neighbors within the first coordination shell of the respective atoms.The cutoff distance was obtained from the first minimum of the respective partial PDFs.
(iii) Structure factor The partial structure factors were calculated from the Fourier transformation of the partial pair distribution functions   () as: where is the scattering vector,  0 is the average atom number density, and  is the maximum value of the integration in real space, which is set to half of the size of one side of the simulation cell.The total structure factor was calculated as: where   and   are the fractions of atoms and   and   are neutron scattering lengths, for elements i and j, respectively.The partial pair distribution functions and structure factors were plotted by taking statistical averages over 100 frames at 300 K.The position of the first sharp diffraction peak (FSDP) and the full width at half maxima (FWHM) of the structure factors were calculated by fitting the FSDP with a Lorentzian function for all the cooling rates.
(iv) Q n distribution The Q n notation quantifies the degree of connectivity of the tetrahedra present in oxide glass networks.In this Q n notation, n represents the number of oxygens that form a bridge between two tetrahedranamely, bridging oxygen (BO)-around a network-forming cation.Thus, Q 4 represents a tetrahedron connected to four other tetrahedra through BO atoms.For 45S5 bioglass, the tetrahedra can be formed by SiO4 and PO4 species.Here, the Q n distribution was computed through an in-house code, which counts the number of BO associated with each Si/P atoms.

(v)
Ring distribution In glassy networks, a ring represents a closed path made of interatomic bonds.In particular, we focused on hetero-polar primitive rings that represent the shortest path made of Si-O and P-O bonds to return to the starting atom without retracing any route.The open-source RINGS package was used to compute the ring size distribution.[53] Note that in oxide glasses, the ring size is represented by the number of network formers in the ring, in this case, Si or P. Thus, a ring size of m contains 2m atoms, namely, with m O atoms and m Si or P atoms.

Experimental details (i)
Glass preparation The 45S5 glass was prepared using the melt-quenching method by mixing a batch of the following reagent grade materials: Na2CO3 (Ph Eur.≥ 99.5%, Sigma-Aldrich), CaCO3 (Reag.Ph Eur, Merck), Ca(H2PO4)2.H2O (Budenheim Ibérica S.L.U., ≥95%), and SiO2 (Sigma-Aldrich, purum p.a.).This batch, with its stoichiometric amounts, was melted in a Pt-Rh crucible at 1350 °C for 1 h in an electrically heated furnace.The melt was cast onto brass plates and the transparent glass was annealed for 30 minutes at its glass transition temperature (Tg) of 523°C.X-ray diffraction analysis shows no signs of crystallization.

(ii)
NMR experiments 29 Si MAS NMR was conducted at 11.7 T (99.27 MHz resonance frequency) using an Agilent DD2 spectrometer and 5 mm MAS NMR probe.The powdered bioglass was packed into a 5 mm outer diameter zirconia MAS NMR rotor, providing sample spinning of 5 kHz. 29Si MAS NMR data were collected using radio-frequency pulse widths of 4.8 s, corresponding to a tip angle of /2, and with a 10 s recycle delay, 3000 acquisitions were collected.The 29 Si T1 value was estimated to be on the order of 1-2 s, resulting from the incorporation of iron (as batch material impurity) in the glass sample.A 29 Si MAS NMR spectrum collected with a much longer recycle delay (60s) was identical, confirming adequate relaxation time for the NMR experiments. 31P MAS NMR data were collected at 16.4 T using an Agilent DD2 spectrometer, 3.2 mm MAS NMR probe (22 kHz sample spinning) and a resonance frequency of 283.27 MHz.600 acquisitions were collected with a /2 pulse width of 3 s and a recycle delay of 5 s, optimized for the short 31 P T1.MAS NMR data were processed without apodization and shift referenced to tetramethylsilane and 85% phosphoric acid solution at 0.0 ppm for 29 Si and 31 P, respectively.

(iii)
Deconvolution of NMR results Spectral deconvolution was made using DMFit [54].The phosphate site (Q n ) concentrations were obtained from the simulations of the isotropic peaks and all associated spinning sidebands of the 31 P MAS spectra.The 31 P quantifications were used to constrain the silica sites, assuming charge neutrality in the glasses.The silica site concentrations were obtained from Gaussian simulations of the isotropic peaks and all of their associated spinning sidebands.To demonstrate the effect of the cooling rate on the features of the glass transition, we are first focusing on the evolution of enthalpy (H) and density of the glass upon quenching computed by MD simulations.Figure 1(a) shows the evolution of the enthalpy with respect to temperature for five different cooling rates, namely, 0.01, 0.1, 1, 10, and 100 K/ps.Upon quenching, we observe that the enthalpy decreases continuously and monotonically with temperature.In particular, we observe that there is no sudden change in the enthalpy of the system, confirming the absence of any first-order phase transition.Further, the enthalpy exhibits a smooth and continuous change in slope at a particular temperature, namely the "fictive temperature" of the glass [55].Note that the fictive temperature decreases with decreasing cooling rate, in agreement with the fundamental nature of the glass transition [40,56].The inset of Fig. 1(a) shows the enthalpy at 300 K of the glasses obtained from different cooling rates.We observe that this enthalpy value decreases with decreasing cooling rate.This suggests that the glasses prepared with lower cooling rates are energetically more stable in comparison to those prepared at faster cooling rates.This can be attributed to the fact that, upon slower cooling rates, the system can sample a larger region of the energy landscape, thereby allowing it to attain lower, more stable energy states [56][57][58][59].

Results and Discussion
Figure 1(b) shows the evolution of density with respect to temperature computed by MD simulations for the five cooling rates considered.Similar to the case of the enthalpy, we observe a continuous and monotonic increase in the density of the glass with decreasing temperature-devoid of any first-order phase transition.Although a change in slope corresponding to the glass transition is also observed for the density, it does not occur at the same temperature as that of the enthalpy for each cooling rate.This is in agreement with earlier studies suggesting that the density and enthalpy relaxations in glassy systems at lower temperature may be decoupled [60,61].The inset of Fig. 1(b) shows the cooling rate dependence of density of the glass at 300 K, which increases monotonically from 2.60 g/cm 3 for 100 K/ps to 2.65 g/cm 3 0.01 K/ps, which suggest an enhanced degree of packing of the atoms in the system with slower cooling rates.Further, the simulated density of the glass tends towards the experimental value of 2.70 g/cm 3 with decreasing cooling rates [8].Overall, these results confirm that our simulations offer a realistic description of the typical features of the glass transition.we observe minor peaks corresponding to second neighbor and third neighbor interactions.At larger distances, the PDF gradually converges to unity, confirming the absence of any long-range order in the system.Further, we observe that the position of the peaks of the total PDF varies only slightly, if any, with the cooling rate.This suggests that the first neighbor environments in the atomic structure are not significantly affected by the cooling rates.To further confirm this, we plot the X-O partial PDFs (with X = Si, P, Na, and Ca).In Fig. 3, we observe that for the X-O partial PDFs, neither the position nor the width of the peaks, are affected significantly by the cooling rate.
Next, we analyze the effect of cooling rate on the homo-nuclear interactions.Figure 4 shows the X-X partial PDFs with X = Si, O, Na, and Ca.We observe that the X-X interactions peak generally become sharper with decreasing cooling rate.This suggests that the overall degree of order among increases with decreasing cooling rate.Although this effect is pronounced in the Si-Si partial PDF (see Fig. 4(a)), it is less pronounced in the other pairs (see Fig. 4(b)-(d)).The Na-Na and Ca-Ca partial PDFs exhibit broader peaks as seen in Figures 4(c) and (d), respectively.This is due to the increased mobility of these atoms, thanks to the non-directional ionic bonds they form with oxygen.Other partial PDFs also exhibit similar behavior wherein the first peak becomes increasingly sharper corresponding to a decreasing cooling rate (see Figs. S2-S5 of the Supplementary Material) suggesting an increased local order.This is in agreement with the earlier observation that the structure attains a lower energy stable state as the cooling rate decreases (see Fig. 1(a)).the average bond angle decreases slightly with decreasing cooling rate up to 1 K/ps.Below 1 K/ps the value exhibits little change.In the case of O-Si-O, we observe that the average intra-tetrahedral bond angle is not affected by the cooling rate.This suggests that the overall tetrahedral structure in the studied bioglass is not significantly affected or distorted by high cooling rates.
Figure 6 shows the inter-tetrahedral BADs.Similar to the intra-tetrahedral BADs, we observe that the Si-O-Si BAD exhibits little variation with different cooling rates.Interestingly, we observe that the peak intensity of BAD for Si-O-P decreases and the distribution becomes broader with decreasing cooling rates.This suggests an increased degree of angular disorder within the Si-O-P bridges with decreasing cooling rate.

(iii) Network connectivity
To analyze the effect of the cooling rate on the silicate network connectivity, we quantify the topology of the oxygen atoms.In oxide glasses, the network formers such as Si or P atoms form tetrahedra, which are connected by some BO atoms through X-O-Y bridges (X or Y = Si or P), thereby increasing the connectivity (or rigidity) of the silicate network.In contrast, network modifiers such as Na or Ca atoms tend to depolymerize the network by breaking the X-O-Y bridges leading to the formation of X-O-Na/Ca bonds.These oxygen atoms are termed as non-bridging oxygen (NBO) atoms.Note that the addition of an Na2O species results in the formation of two NBOs through the formation of two Si-O -Na + bonds.Similarly, the addition of a CaO species result in the formation of two NBOs through the formation of a Si-O -Ca 2+ -O --Si [32].
To this extent, we compute the fraction of Si-O-P and Si-O-Si linkages and compare them with the predictions from a random network model.[62] According to this model, each network former (that is, Si or P) can pick its neighbor with equal probability from the available atoms, irrespective of the atom type.Thus, the fractions of Si-O-Si, Si-O-P, P-O-P, angles are given by where   and   are the numbers of Si and P atoms, respectively.
Figure 7 shows the fraction of the Si-O-Si and Si-O-P linkages predicted from the random model compared with the results from MD simulations for different cooling rates.We observe that, overall, the glasses exhibit a chemically-ordered structure that deviates from the predictions of the random modelwherein the computed fractions of Si-O-P and Si-O-Si linkages are lower and higher than those obtained from a random distribution, respectively.Note that the P-O-P linkages are absent in structures generated from MD simulations, in agreement with experimental observations.As the cooling rate decreases, the deviation from the random model increases and the glass exhibits an enhanced chemically-ordered structure.Further, the fraction of Si-O-P linkages decrease (and Si-O-Si linkages increase) with decreasing cooling rate.Overall, these results suggest that the phosphate species exhibit a tendency to form isolated orthophosphate tetrahedra in bioglass-the manifestation of which is closely linked to the cooling rate used to prepare the glass.To confirm the increased affinity of P toward NBOs, we compute the Si-P partial PDF (see Fig. 8).Upon decreasing cooling rate, the intensity of the Si-P peak diminishes, and the peak becomes notably less intense.Note that this is in contrast to the other partial PDFs (Figs. 3 and 4), which exhibit sharper peaks with decreasing cooling rates.This confirms that the tendency of P atoms to form isolated orthophosphate species, which is magnified upon decreasing cooling rate, i.e., when the glass reaches a more stable state.In other words, for a system obtained with reasonably low cooling rates, such as the experimental glasses, there is a tendency to avoid the formation of X-O-P bridges.Instead, the phosphate tetrahedra tends to remain as isolated [PO4] 3-species with four NBOs.This observation is consistent with previous experimental studies, which report a low fraction of Q 1 P species.[26,27] (iv) To study the connectivity of the tetrahedral network further, we analyze the Q n distribution (see Methodology section) both computationally and from solid-state NMR measurements on a bioglass sample.The fitting parameters used herein for both the 31 P and 29 Si MAS NMR data, as well as the obtained Q n distributions, are presented in Tables 1 and 2.   29 Si and 31 P MAS NMR with deconvolution into respective Q n species.From 29 Si MAS NMR, we observe that the fraction of QSi 2 units is approximately 76%, while that of QSi 1 and QSi 3 are approximately 5% and 17%, respectively (Table 2).The present experimental results from 29 Si MAS-NMR are in agreement with previous studies by Yu et al. [29] and Pedone et al. [20].These results suggest that the silicate network primarily consists of QSi 2 tetrahedra, with some degree of cross-linking provided by QSi 3 species.From 31 P MAS-NMR, we observe that P atoms exists only as Q 0 and Q 1 units, wherein the fraction of Q 0 atoms is around 97% (Table 1).The residuals associated with the fit is provided in the supplementary material.This confirms that the P atoms prefer to exist mostly as isolated Q 0 orthophosphate ions in the bioglass system, although non-trivial amounts (3%) of Q 1 may exist.This observation is in agreement with Yu et al. [29], where 4.1% of QP 1 were observed.Further, the chemical shifts for QP 0 and QP 1 observed by Yu et al. are 8.7 and 0.4 ppm, respectively, which are comparable to the chemical shifts of 8.9 and 2.0 ppm, respectively, observed in the present study.It should be noted that QP 0 and QP 1 exhibit different chemical shift anisotropy (CSA).Since, CSA is highly sensitive to symmetry, QP 0 is expected to exhibit a smaller CSA.Based on an accurate fit 31 P MAS NMR spectrum utilizing all spinning sidebands, we estimate the CSA of QP 0 and QP 1 as 56 and 124 ppm, respectively.Overall, the present work conclusively provides the evidence for the existence of QP 1 species, although in minor fractions of around 3%.Now, we compare the experimental results with the Q n distributions obtained from MD simulations.Figures 10(a) and (b) shows the Q n distributions of the tetrahedral species in the system formed by Si and P atoms, respectively.Note that, in addition to the experimental and simulation results from the present study, experimental results from Yu et al. [29] and Pedone et al. [20] and simulation results from Pedone et al. [20] are included for comparison.In the case of Si, MD simulations results suggest that a large fraction (around 40-45%) of the glass network comprises of QSi 2 tetrahedra, at the expense of QSi 3 or QSi 4 tetrahedra.Further, we observe that the fraction of QSi 0 , QSi 1 , and QSi 4 units decreases with decreasing cooling rate.Now, the simulated values are extrapolated to experimental cooling rates by fitting a line in the least squared sense (see Fig. 10(a)).The extrapolated values are then compared to experimental results as shown in Fig. 10(a).We observe that the trend in Q n distribution of Si is wellpredicted by the MD simulations.In particular, the fraction of Q 3 units exhibit an excellent match when extrapolated.However, Q 2 units are slightly underpredicted while the Q 1 units are overpredicted by the extrapolation from MD simulations.Overall, we observe that the trend predicted by MD simulations exhibits a good match with experiments for the QSi n units of Si atoms.
In the case of the P atoms, we observe that the structure exhibits only Q 0 , Q 1 , and Q 2 units, even at the highest cooling rates (see Fig. 10(b)).Upon decreasing cooling rate, Q 2 units decrease quickly and converge to zero, even within the range of cooling rates accessible by MD simulations.Further, the fractions Q 0 and Q 1 units monotonically increase and decrease with decreasing cooling rate, respectively.This confirms that, as the cooling rate decreases, P atoms increasingly tend to form some isolated [PO4] 3- tetrahedral units.These values are then compared to the experimental results obtained using 31 P MAS-NMR by extrapolating the MD results in a least squared sense.Interestingly, we observe that the extrapolation of MD simulation suggests the existence of Q 0 species only, not accounting for the minor QP 1 species suggested by the experiments.However, this could be due to the low amount of P atoms present in the bioglass system, making an accurate determination of the Q 1 units challenging.Extending MD simulations to lower cooling rates or increasing the size of the system may provide deeper insights into the possible existence of QP 1 units.Overall, we observe that the effect of the cooling rate on the Q n distributions can be reasonably captured when MD simulations are extrapolated to experimental timescales, and furthermore, the experimentally-determined network structure of bioglass is consistent with previous studies of bioglasses.(v) Medium-range order While the PDF is useful for analyzing the structure correlations at smaller distances (≤ 3 Å), the structure factor represents the correlations at larger distances (between 3 Å to 10 Å).In particular, the first sharp diffraction peak (FSDP) captures the extent of correlation in the medium-range structure of glasses.[42,[63][64][65][66] The full width at half-maximum (FWHM) of the FSDP is inversely linked to a mediumrange coherence length L as, L = 7.7/FWHM, [42,[63][64][65][66] which captures the medium-range correlation distance among fluctuations of atomic density.Figure 12(a) shows the neutron structure factors computed by MD simulations with different cooling rates.We observe that the structure factor at larger values of wave vector exhibits little differences.However, FSDP becomes sharper with increasing peak intensity for decreasing cooling rates.This indicates that, in turn, the coherence length in the mediumrange order increases upon decreasing cooling rate.To confirm this, we plot the FWHM and coherence length in Fig. 12(b).We observe that, upon decreasing cooling rate, FWHM decreases and the coherence length increases (see Fig. 12(b)), which denotes an increased degree of order within the medium-range structure of the glass.This echoes the fact that, due to the tendency P atoms to exist as isolated species, the network presents an increasingly chemically-ordered structure upon decreasing cooling rate-in agreement with the experimental and simulation data on Q n as well.Finally, we analyze the effect of the cooling rate on the ring size distribution within the network computed by MD simulations (see Methodology section).As shown in Fig. 13, we find that, in contrast to the short-range order, the ring distribution is largely affected by the cooling rate.First, we observe that there are no two-membered rings, which confirms the absence of any edge-sharing tetrahedra.Further, there are no small three-membered rings for all but the highest cooling rate of 100 K/ps.This arises from the fact that, at high cooling rate, atoms may get frozen in unstable states due to the limited time available for structural rearrangements.With decreasing cooling rate, we observe that larger rings are preferred with a maximum fraction observed for five-membered rings, that is, rings having five Si atoms.In particular, a sharp peak at a ring size of 5 is observed for the cooling rate of 0.01 K/ps, which shows that almost 50% of the rings have a size of five Si atoms.This is in contrast with glassy silica, wherein a ring size of six is preferred.[42,43,67] Finally, we observe that all the rings are formed by Si and O atoms, thereby confirming that the P atoms do not participate in the ring formation.This is consistent with the fact that P atoms mostly form isolated Q 0 and terminating Q 1 units and, hence, are unable to participate in these ring structures.Overall, we conclude that the medium range structure of the bioglass is significantly affected by the cooling rate.The present results suggest a holistic picture reconciling earlier experimental and MD-simulation-based studies.While the short-range structure of bioglass hardly depends on the cooling rate, the medium range structure is significantly affected.In particular, our results confirm that P atoms exhibit an increased affinity to NBOs, thereby leading to the gradual disappearance of Si-O-P bridges at lower cooling rates.In other words, lower cooling rate results in an increase in the concentration of isolated orthophosphate anions, which associate with the Ca and Na cations present in the structure.As such, the extrapolation of MD results toward lower experimental cooling rates would result in: (i) a significant decrease in the fraction of Si-O-P bridges, (ii) the absence of any QP 2 species, as observed in our NMR measurements, and (iii) an increased propensity for P atoms to exist as isolated orthophosphate ions.Indeed, it is well known that the high cooling rates used in MD simulations tend to overestimate the degree of randomness in computed Q n distributions.Overall, the present study emphasizes the importance of using a cooling rate study in MD simulations covering different orders of magnitude to reconcile experimental and simulation structural studies.
From a more general perspective, these results can have significant ramifications in the understanding of the structure and bioactivity of bioglass.It is well known that bioactivity is closely related to the leaching of ions from the glass during the dissolution process.[68] This occurs through the breakage of interatomic bonds present in the glass structure.As such, the absence of a large number of Si-O-P bridges, which were otherwise postulated by previous MD studies, suggests that isolated phosphate ions may exhibit an increased mobility and, hence, can easily leach out from the glass.This is further confirmed by the NMR experiments wherein approximately 97% of P atoms are shown to exist as Q 0 units.This may in turn contribute to the bioactive behavior of the glass.However, confirming this hypothesis would require future experimental studies on the degree of congruency in the leaching rates of Si and P species in bioglass.

Conclusions
The present study reveals the importance of the cooling rate in governing the structure of bioactive glasses.In particular, we show that the medium-range order structure is largely controlled by the thermal history of the glass.Interestingly, we observe an increased affinity of P atoms toward NBOs leading to a significant decrease in the fraction of Si-O-P bridges upon lower cooling rates.This observation, along with the propensity for P atoms to form isolated tetrahedra, confirms the earlier experimental findings that phosphorus atoms primarily exist as isolated orthophosphate ions in bioglass, although minor quantities of QP 1 may also exist.More importantly, we show that, through a detailed analysis of cooling rate simulations covering several order magnitudes, the experimental results can be reasonably extrapolated using MD simulations.This, in turn, can provide a handshake between MD simulation and experiments, which is otherwise challenging due to the extremely high cooling rates involved in MD simulations.

Figure 1
Figure 1 (a).Computed values of enthalpy (H) as a function of temperature during cooling under different rates by MD simulations.The inset shows H at 300 K as a function of the cooling rate.The line is a guide for the eye.(b) Computed values of density as a function of temperature during cooling under different rates by MD simulations.The inset shows density at 300 K as a function of the cooling rate.The line is a guide for the eye.

Figure 2 .Figure 3 .
Figure 2. Total pair distribution functions in the glassy state for the studied cooling rates computed by MD simulations.

Figure 4 .
Figure 4. Homo-nuclear partial pair distribution functions in the glassy state, computed by MD simulations for the studied cooling rates, of: (a) Si-Si, (b) O-O, (c) Ca-Ca, and (d) Na-Na.

Figures 5 (
Figures 5(a) and (b) show the intra-tetrahedral O-P-O and O-Si-O bond angle distributions (BADs),respectively.We observe that the average intra-tetrahedral angles (i.e., around 109°) are not significantly affected by the cooling rate, although the peaks become slightly sharper for lower cooling rates.This suggests an increased degree of angular order upon lower cooling rate.In the case of O-P-O,

Figure 5 .Figure 6 .
Figure 5. Intra-tetrahedral (a) O-P-O bond and (b) O-Si-O bond angle distribution (BAD) in the glassy state computed by MD simulations for selected cooling rates.

Figure 7 .
Figure 7. Fraction of Si-O-Si and Si-O-P linkages predicted from the random model compared with the results from MD simulations for different cooling rates.

Figure 8 .
Figure 8. Si-P partial pair distribution functions in the glassy state for selected cooling rates computed by MD simulations.
Figures 9 (a) and (b) show the29 Si and31 P MAS NMR with deconvolution into respective Q n species.From29 Si MAS NMR, we observe that the fraction of QSi 2 units is approximately 76%, while that of QSi 1 and QSi 3 are approximately 5% and 17%, respectively (Table2).The present experimental results from29 Si MAS-NMR are in agreement with previous studies by Yu et al.[29] and Pedone et al.[20].These results suggest that the silicate network primarily consists of QSi 2 tetrahedra, with some degree of cross-linking provided by QSi 3 species.From 31 P MAS-NMR, we observe that P atoms exists only as Q 0 and Q 1 units, wherein the fraction of Q 0 atoms is around 97% (Table1).The residuals associated with the fit is provided in the supplementary material.This confirms that the P atoms prefer to exist mostly as isolated Q 0 orthophosphate ions in the bioglass system, although non-trivial amounts (3%) of Q 1 may exist.This observation is in agreement with Yu et al.[29], where 4.1% of QP 1 were observed.Further, the chemical shifts for QP 0 and QP 1 observed by Yu et al. are 8.7 and 0.4 ppm, respectively, which are comparable to the chemical shifts of 8.9 and 2.0 ppm, respectively, observed in the present study.It should be noted that QP 0 and QP 1 exhibit different chemical shift anisotropy (CSA).Since, CSA is highly sensitive to symmetry, QP 0 is expected to exhibit a smaller CSA.Based on an accurate fit 31 P MAS NMR spectrum utilizing all spinning sidebands, we estimate the CSA of QP 0 and QP 1 as 56 and 124 ppm, respectively.Overall, the present work conclusively provides the evidence for the existence of QP 1 species, although in minor fractions of around 3%.

Figure 9 .
Figure 9. (a)29 Si MAS NMR spectrum of bioglass with deconvolution into Q 1 , Q 2 and Q 3 silicate groups and (b)31 P MAS NMR spectrum of bioglass with deconvolution into Q 0 and Q 1 phosphate groups.

Figure 10 .
Figure 10.Q n distribution of (a) Si and (b) P tetrahedra in 45S5 Bioglass computed by MD simulations for five different cooling rates.Experimental results obtained by MAS-NMR and from Yu et al.[29] and Pedone et al.[20] and simulation results from Pedone et al.[20] are added for comparison.The lines are fitted using the MD simulation data in a least-squared sense.

Figure 12 .
Figure 12.(a) Computed total neutron structure factors in the glassy state for selected cooling rates.(b) Full width at half maxima (FWHM) of the total structure factor as a function of cooling rate (left yaxis) and the corresponding coherence length (right y-axis).The line is a guide for the eye.

Figure 13 .
Figure 13.Ring distribution for the bioglass at various cooling rates computed by MD simulations.

Table 1 .
Parameters used for deconvolution of 31 P MAS NMR spectra of 45S5 bioglass.The QP n populations (atom %) include intensity from the spinning sidebands.

Table 2 .
Parameters used for deconvolution of29Si MAS NMR spectra of 45S5 bioglass.The QSi n populations (atom%) also include all spinning sideband intensity.