Tailoring the Mechanical Properties of Metaluminous Aluminosilicate Glasses by Phosphate Incorporation

The characterization of aluminosilicate glasses is highly relevant in geosciences and for engineering applications such as reinforcement fibers or touchscreen covers. The incorporation of phosphate as a third network-forming species into these glasses offers unique opportunities for fine-tuning glass properties via changes in glass structure and polymerization. In this work, we studied melt-quenched aluminosilicate glasses within the system SiO2-Al2O3-Na2O-P2O5 with 50–70 mol% SiO2 and up to 7.5 mol% P2O5. All glasses were metaluminous (Al:Na = 1) in order to maximize the degree of polymerization. Increasing the phosphate content at the expense of NaAlO2 led to reduced glass polymerization and density, resulting in a decrease in elastic moduli and hardness and an increase in strain-rate sensitivity. When increasing the silica content by substituting SiO4 for AlO4 tetrahedra, network polymerization remained mostly unchanged, as confirmed by nearly constant hardness. Densification upon indentation was analyzed by Raman spectroscopy and finite element analysis. We find that the elastic properties and hardness of metaluminous phospho-aluminosilicate glasses are governed by changes in density and network polymerization. Other mechanical properties underlie more complex changes in glass structure.


INTRODUCTION
The mechanical properties of aluminosilicate glasses have been valued by humans since ancient times, as obsidian, a natural aluminosilicate glass vitrified from lava, found use in the fabrication of tools (Cann and Renfrew, 1964;Ericson et al., 1975;Bellot-Gurlet et al., 2004). In the early days of fracture mechanics, the high intrinsic strength of aluminosilicate glasses was demonstrated (Griffith, 1920), paving the way for their use as reinforcement fibers in composites (Botev et al., 1999;Dibenedetto, 2001;Nkurunziza et al., 2005;Sathishkumar et al., 2014). More recently, damage resistant aluminosilicate glasses have been evolving (Bechgaard et al., 2016;Zeng et al., 2016). Understanding the underlying structural origin of their response to mechanical load provides tools for tailoring glass properties (Wondraczek et al., 2011). As an ultimate objective, introducing microscopic ductility as a route to dissipate mechanical energy would enable glasses with unprecedented practical strength (Wondraczek, 2019).
Metaluminous aluminosilicate glasses are highly polymerized as aluminum is present in the form of network-forming AlO 4 tetrahedra, charge-balanced by alkali metal cations such as sodium (Al:Na = 1). While the high degree of polymerization already yields favorable mechanical properties, small amounts of additional glass components may allow for fine-tuning.
The use of P 2 O 5 as a third network-forming oxide besides SiO 2 and Al 2 O 3 is such an option. Phosphate forms tetrahedral structural units similar to those of SiO 4 and AlO 4 . However, phosphate tetrahedra are connected to the surrounding network structure by less than four bridging oxygen bonds (Dupree et al., 1988;Kosinski et al., 1988;Gan and Hess, 1992;Toplis and Schaller, 1998) and may thus reduce glass polymerization. On the other hand, a preferential association of phosphate with Al 3+ and Na + (Ryerson and Hess, 1980;Gan and Hess, 1992;Mysen, 1998) may additionally affect glass structure and properties.
Up to now, few studies have been reported on the effects of phosphate on the mechanical properties of aluminosilicate glasses (Lönnroth and Yue, 2009;Zeng et al., 2016;Tarragó et al., 2018a,b). These studies concerned only with an excess of alkaline ions over aluminum rather than highly polymerized metaluminous glasses. Furthermore, while all of these studies investigated glass hardness, none of them reported on the elastic properties.
The aim of this work was therefore to provide a comprehensive mechanical characterization of phosphatecontaining metaluminous aluminosilicate glasses and to identify the structural origin of these data. For this purpose, the model glass system of SiO 2 -Al 2 O 3 -Na 2 O-P 2 O 5 was investigated. Three glass series were compared: one series with increasing P 2 O 5 content from 0 to 7.5 mol% and constant 60 mol% SiO 2 and two series with increasing SiO 2 content from 50 to 70 mol% and constant P 2 O 5 content of either 0 or 7.5 mol%. The study focused on elastic properties and hardness, complemented by investigations of density, densification upon indentation, strain-rate sensitivity and crack resistance.
transferring the melt inside the crucible to an annealing furnace preheated to temperatures 20 K above the respective glass transition temperature, T g , to reduce residual stresses. At 1 h at this temperature, the crucibles were slowly cooled down to room temperature by switching off the furnace. Glass monoliths were retrieved by breaking the alumina crucible walls: glasses were not cast from platinum crucibles and quenched because of high melt viscosity. From the monoliths, coplanar glass samples were cut and polished on both sides. As leaching of alumina from the crucibles could not be excluded, longer melting times were avoided and sections close to the crucible walls were omitted during sample cutting. The protocol described here allowed for comparable synthesis conditions for all glasses except Si50P0. For this glass, owing to its high Na 2 O content of 25 mol%, crystallization could only be prevented by reducing the annealing temperature to T g -40 K. However, tests with various annealing temperatures showed that even such reduced cooling temperature had no significant effect on glass properties (results not shown).

Compositional Analysis
Actual chemical compositions of all glasses were analyzed by energy-dispersive X-ray spectroscopy (EDX). For each glass the composition was calculated as mean value of four EDX measurements. Errors represent standard deviations. Using different specimens, three of these four measurements were performed with the same scanning electron microscope (Phenom pro X, Phenomworld, Eindhoven, NL). The fourth EDX measurement was performed with a different instrument (Neon 60, Zeiss, Jena, Germany, combined with X-Max X-Ray Detector 80 mm 2 , Oxford Instruments, High Wycombe, UK) at a different laboratory.

Density and Related Properties
Glass density, ρ, was determined as the arithmetic mean out of ten measurements on a single specimen, using a helium pycnometer (AccuPyc 1330, Micromeritics GmbH, Mönchengladbach, Germany). Powder samples with a particle size of 125-250 µm were used to minimize the influence of bubbles. The experimental error is ±0.01 g/cm 3 . Molar volume, V m , was calculated from density and molar weight, M i , of each glass component weighted by its nominal molar fraction, f i .
All calculations presented here were based on nominal glass compositions, justified by a good agreement with the analyzed compositions (Table 1). Atomic packing density, C g , was estimated using the following equation (Makishima and Mackenzie, 1973): Here, V i = 4/3πN xr 3 A + yr 3 B represents the molar volume occupied by an oxide A x B y . N is Avogadro's number and r A as well as r B are the ionic radii of the corresponding cations and anions, respectively. Effective ionic radii were employed (Shannon, 1976), considering Si 4+ (r Si = 26 pm), Al 3+ (r Al = 39 pm), and P 5+ (r P = 17 pm) species in tetrahedral coordination, Na 2+ (r Na = 102 pm) in octahedral coordination and O 2− (r O = 135 pm) in two-fold coordination. It is noted that the radius of oxygen ions may vary depending on the coordinated cations and may be determined by the cation-oxygen bond lengths (Duffy, 1990;Zeidler et al., 2014). However, in the absence of quantitative data on the amount of bridging and non-bridging oxygen species and the connectivity between different polyhedral species, a constant ionic radius of oxygen was chosen for the sake of simplicity. Consequences of this assumption will be addressed in the discussion.
Oxygen packing density, ρ ox , was calculated as the molar amount of oxygen per unit volume of glass (Ray, 1974;Fredholm et al., 2010).
N ox,i is the stoichiometric number of oxygen atoms provided by each glass component, i.e., 2 for SiO 2 , 5 for P 2 O 5 , 3 for Al 2 O 3 and 1 for Na 2 O. Errors of molar volume, packing density and oxygen packing density were estimated by error propagation calculation.

Ultrasonic Echometry and Brillouin Spectroscopy
The longitudinal and transversal sound velocities, v L and v T , were measured by ultrasonic echometry and Brillouin spectroscopy. On that basis, elastic properties were calculated, including shear modulus, G, Young's modulus, E, bulk modulus, K, and Poisson's ratio, ν, under the assumption of isotropic symmetry in the glasses (Rouxel, 2007). Experimental uncertainties were estimated by error propagation.
Ultrasonic echometry was performed on coplanar, opticallypolished samples of ∼2 mm thickness using an echometer (Echometer 1077, Karl Deutsch GmbH & Co. KG, Wuppertal, Germany) equipped with piezoelectric transducers operating at frequencies of 8-12 MHz. Values of v T and v T were derived from sound propagation times, which were recorded with an accuracy of ±1 ns, and from the exact thickness of the 2 mm thick coplanar, optically-polished specimen, as measured with an accuracy of ±2 µm using a micrometer screw. Brillouin spectroscopy was performed at room temperature on coplanar, optically-polished samples of 200-300 µm thickness using a custom-built optical setup as described previously (Veber et al., 2018). Brillouin shifts were averaged between the Stokes and anti-Stokes signals to account for instrumental non-zero position of the Rayleigh line. Experiments were performed in backscattering geometry (through a microscope, 50x magnification, 0.42 numerical aperture) as well as platelet geometry (45 • incidence angle).
In backscattering geometry, v L was determined as v L = f·λ/(2n) (Rabia et al., 2016), where λ (488 nm) is the wavelength of the continuous wave sapphire laser used, and f is the Brillouin shift of the longitudinal acoustic mode, which was determined with an accuracy of ±0.012 GHz. The glasses' refractive index, n, at 488 nm was deduced by linear interpolation between refractive indices at 480 and 546.1 nm as determined by Pulfrich refractometry (PR2, Carl Zeiss, Jena, Germany) with an experimental uncertainty of ±0.001.
In platelet geometry, v L and v T were determined from the Brillouin shifts of the longitudinal and transversal acoustic modes, f L and f T , as v L/T = f L/T · λ/( √ 2) (Whitfield et al., 1976). Although the platelet geometry delivered both v L and v T , the v L value from backscattering geometry is normally more precise because of sample surfaces not being perfectly coplanar and because of precision on the 45 • angle in platelet geometry. The values of v L from backscattering geometry were therefore preferred and further used to correct v T from platelet geometry (v T,corrected = v T × v L,backscatter /v L,platelet ).

Calculation of Elastic Properties
In addition to experimental results, we employed the model of Makishima and Mackenzie (1973) to predict the compositional dependence of Young's modulus. Following this approach, Young's modulus of a glass is determined from the atomic packing density and volume density of bonding energy, <U 0 /V m >.
Note that <U 0 /V m > was calculated by means of the actual glass density rather than the density of the individual components in their crystalline state (Limbach et al., 2017). The average bond energy, U 0 = f i H ai , was calculated from molar dissociation enthalpies, H ai , of the glass components, weighted by their nominal molar fractions, f i . Values of H ai were determined from molar heats of formation, H f , of each oxide component A x B y in its crystalline state and the constituent atoms in gaseous state (Inaba et al., 1999).
Values of H f were taken from Lide (2006), except for H f (P 2 O 5 ) = −306.4 kJ/mol, which was calculated from the dissociation energy of P 2 O 5 , H ai /V m = 62.8 kJ/cm 3 , available from literature (Inaba et al., 1999) and using Equations (1) and (9). The error of the calculated modulus was estimated by error propagation.

Nanoindentation
Instrumented indentation experiments were carried out using a nanoindenter (Nano Indenter G200, Schaefer Technology GmbH, Langen, Germany) equipped with a three-sided Berkovich diamond tip (Synton-MDP, Port, Switzerland) and operating in continuous stiffness measurement (CSM) mode (applying an oscillation of amplitude h = 2 nm and frequency f = 45 Hz). Values of Young's modulus and hardness, H, were determined from indentations with a depth limit, h, of 2 µm created at constant strain-rate,ε, of 0.05 s −1 , following the method proposed by Oliver and Pharr (1992). For statistical relevance, the values of E and H were averaged over at least 15 individual indentation experiments.
Time-dependent indentation deformation was analyzed in a nanoindentation strain-rate jump test, as described in detail elsewhere (Maier et al., 2011;Limbach et al., 2014). This test initially comprises a first segment, where the Berkovich indenter tip penetrates the glass surface to a depth of 500 nm at a constant strain-rate of 0.05 s −1 . Afterwards, with further increasing indenter displacement, the strain-rate is changed in intervals of 250 nm while monitoring the accompanied changes in hardness using the CSM equipment of the nanoindentation setup ( h = 5 nm, f = 45 Hz). In this study, a total number of at least eight strain-rate jump tests were evaluated on each sample with strainrates of 0.05, 0.007, and 0.001 s −1 (in descending order). The strain-rate sensitivity, m, was calculated from the slope of the logarithmic plot of the hardness, lnH, over the indentation strainrate, lnε i , where,ε i is equal toε/2 for materials with a depthindependent hardness value (Maier et al., 2011(Maier et al., , 2013Limbach et al., 2014;Zehnder et al., 2017).
Error bars show standard deviations from averaging over all evaluated indentation experiments.

Vickers Hardness
Vickers hardness, H v , was determined using a Vickers microhardness tester (Duramin-1, Struers GmbH, Willich, Germany). On each sample, 25 indentations with a load, P, of 981 mN were created. Loading duration and dwell time were set to 15 and 10 s, respectively. H v was calculated from the load in N and the length of the projected indentation diagonals, d, in µm according to H v = 1.8544(P/d 2 ) (Tabor, 1970). Error bars show standard deviations.

Crack Resistance
Crack resistance was analyzed by Vickers indentation. Indents with stepwise increasing loads from 981 mN up to 19.6 N were created in air under ambient conditions using a constant loading duration of 15 s and dwell time of 10 s. After unloading, the number of median-radial cracks occurring at the corners (illustrated by inset in Figure 1) of the residual Vickers hardness imprints was counted. On each sample, 25 indentations were performed per load applied and the probability of crack initiation, PCI, was calculated from the number of corners with cracks divided by the total number of corners. The load dependence of PCI is exemplarily shown for glass Si60P7.5 in Figure 1. The experimental data follows the typical sigmoidal trend (orange dashed line) as reported previously (Limbach et al., 2015b). The load value of the fitted sigmoidal curve at 50% PCI was taken as crack resistance, CR, i.e. the load at which on average two cracks were formed. A detailed description of the procedure used for assessing the experimental uncertainty of this method is given in the Discussion.

Raman Spectroscopy
To characterize glass densification, the low frequency region of Raman spectra was analyzed using the method presented by Deschamps et al. (2011Deschamps et al. ( , 2013. After indentation, the low frequency Raman region (blue line in Figure 2), including its center of gravity, COG (dash-dotted gray line in Figure 2), was shifted (illustrated by arrow in Figure 2) toward higher wavenumbers. This COG shift upon indentation was used to quantify densification.
The low frequency envelopes of Raman spectra measured before (black spectrum in Figure 2) and after indentation (orange FIGURE 1 | Probability of crack initiation (PCI) as a function of the applied Vickers indentation load for glass Si60P7.5. The inset shows the typical crack pattern observed after Vickers indentation (load 19.6 N). The orange dashed line illustrates the sigmoidal fit used to determine crack resistance (defined as the load at PCI = 50%). Connecting the end points of the PCI error bars creates an approximated error envelope. The intersections of this error envelope with the 50% PCI line were used to estimate the experimental uncertainty in crack resistance.
FIGURE 2 | Raman spectra of a glass before (black) and after (orange) densification by Vickers indentation. The blue spectrum shows the low frequency region after subtraction of a baseline (dashed black line) and was used to calculate that region's center of gravity (COG, gray dash-dotted line). The densification shifted intensities in the low frequency region. This was quantified by determining the shift of the COG (arrow, not to scale). spectrum in Figure 2) were integrated to obtain their respective COG using the following equation.
Here, the COG is the wavenumber which divides the integrated region into two parts of equal area. I(ω) is the Raman intensity and ω the Raman shift in cm −1 . The lower and upper integration limits are low and up, respectively. The limits were chosen at the positions of the intensity minima which confined the low frequency Raman region. Following the procedure by Deschamps et al. (2011), a linear baseline as illustrated in Figure 2 was subtracted from the low frequency Raman region prior to COG calculation. Furthermore, the use of this baseline was compared to another post-processing routine involving a linear baseline fitted to the signal-free high frequency end of the Raman spectra to remove luminescence, followed by a data reduction by the Long equation (Long, 1977(Long, , 2002. Here, ω 0 is the wavenumber of the excitation laser,h is the reduced Planck constant, k B the Boltzmann constant and T the absolute temperature. The Long reduction required a change of the lower integration limit because it shifted the corresponding intensity minimum to the lowest measured frequency. This enlarged the integrated range for datasets that underwent the Long reduction. However, the impact of the changed integration limits on the COG calculation was considered low as Raman intensities in the affected wavenumber regime dropped significantly during Long reduction. Raman spectra were measured in the centers of each of five residual Vickers hardness imprints (load 981 mN) per sample, as well as five reference spectra of the pristine glasses. Obtained COG shifts were averaged, error bars show standard deviations. Raman measurements were carried out using a Raman microscope (invia Raman Microscope, Renishaw, Wotton-under-Edge, UK). Measurements were performed at room temperature using the 488 nm line of an argon ion laser and a microscope objective lens with 50× magnification and a 0.50 numerical aperture. Spectra were measured in the range from 150 to 1,550 cm −1 . The acquisition time per spectrum was 225 s, summarized over several scans.

Finite Element Analysis
Finite Element Analysis (FEA) was performed using the software package ABAQUS, similar to earlier studies (Bruns et al., , 2020. The indentation process was modeled in a two-dimensional (2D) axisymmetric way using the Berkovich equivalent cone, which exhibits an opening angle of 70.3 • . The indenter/material contact was assumed to be frictionless. All material properties were presumed to represent rate insensitive room temperature values. The series with constant 60 mol% SiO 2 content (Si60Px) was selected as model system, whose corresponding elastic properties (isotropy assumed) can be found in the Results section.
The yield surface of this series was mapped with Drucker-Prager-Cap plasticity according to Bruns et al. (2020). The yield strength under hydrostatic pressure relied on densification hardening curves usually recorded in diamond anvil cell (DAC) experiments. Since corresponding data are not available for the system studied here, the densification capacity was estimated according to Rouxel's model predictions based on Poisson's ratio (Rouxel et al., 2008). For ν = 0.22 (constant value found for the entire Si60Px series) a densification capacity of 8.6 ± 2.5% can be estimated. The experimental densification hardening curve of fused silica (Rouxel et al., 2008;Sonneville et al., 2012;Deschamps et al., 2013) was then compressed to the saturation level determined earlier (8.6%, Figure 3A) and serves as input for cap hardening along hydrostatic axis.
The yield strength under pure shear was fitted based on experimental nanoindentation load-displacement curves using the extrema of the Si60Px series with 0 and 7.5 mol% P 2 O 5 . Since the densification prediction does not exhibit the sensitivity to account for those chemical variations, only one representative yield surface for this series was constructed ( Figure 3B).

Compositional Analysis
Deviations from the nominal compositions are listed in Table 1. Almost all values measured are similar to the nominal compositions within error limit. It can be noticed however that SiO 2 and P 2 O 5 content were systematically lower and Al 2 O 3 and Na 2 O content higher than nominal composition. The overall good agreement justifies the utilization of the nominal compositions for the theoretical considerations on packing efficiency or bond energy density.

Density and Related Properties
Density (ρ) and atomic packing density (C g ) decreased, while molar volume (V m ) and oxygen packing density (ρ ox ) increased for glasses with increasing P 2 O 5 content and constant SiO 2 content (Figures 4A,C,E,G). For increasing SiO 2 content and constant P 2 O 5 content of either 0 or 7.5 mol% density and atomic packing density decreased, and oxygen packing density increased ( Figures 4B,F,H). For phosphate-free glasses the decrease of packing density was within the error limit at lower silica contents (Figure 4F). At 50 mol% SiO 2 , the incorporation of P 2 O 5 did not notably change packing density. The molar volume only slightly decreased for increasing SiO 2 content and remained almost constant for glasses with 7.5 mol% P 2 O 5 ( Figure 4D).

Elastic Properties
Compositional variations in elastic properties, including shear (G), bulk (K) and Young's moduli (E) as well as Poisson's ratio (ν), as obtained from the different experimental methods are illustrated in Figures 5-7. Trends described in the following apply to all methods employed.
For glasses with constant 60 mol% SiO 2 and increasing phosphate content (Figure 5), the values of G, K and E decreased while ν remained constant within error limits.
For increasing silica content in phosphate-free glasses ( Figure 6) as well as glasses with 7.5 mol% P 2 O 5 (Figure 7), the trends were less clear. G remained roughly constant while K and E decreased with increasing silica content, regardless of the phosphate content. The decrease in E was weak, but K showed a decrease of a similar magnitude as the one observed with increasing phosphate content. Poisson's ratio appeared to decrease with increasing silica content.
Elastic properties obtained by ultrasonic echometry and Brillouin spectroscopy were in good agreement (Figures 5-7), reflecting the fact that both methods rely on the propagation of sound waves. The only exception was glass Si60P5 (5 mol% P 2 O 5 , Figure 5), where slightly higher values of G, K and E were probed by Brillouin spectroscopy as compared to ultrasonic echometry.
The trends for E obtained by nanoindentation (green triangles in Figures 5-7) agreed well with the trends determined by ultrasonic echometry and Brillouin spectroscopy. Here again, the sample with 5 mol% P 2 O 5 was slightly off the trend, indicating an outlier. Values for E obtained by nanoindentation were systematically lower by up to 5 GPa than those obtained by  the other two methods. This small discrepancy results from differences in testing conditions and may originate from sink-in effects or a piling-up of material at the periphery of the indenter (Limbach et al., 2015a(Limbach et al., , 2017. Young's modulus as calculated by the model of Makishima-Mackenzie (E M+M ) was relatively close to the values measured by echometry (Figure 8); however, the model systematically underestimated the experimentally determined values by about 2-8 GPa.

Hardness
Compositional trends in the hardness as obtained from microand nanoindentation agreed well, yielding decreasing hardness values with increasing phosphate content (Figures 9A,C). With increasing silica content (Figures 9B,D), hardness remained roughly constant for both methods, regardless of the phosphate content. However, Vickers hardness (H v ) showed a slight decrease for increasing silica content in phosphate-free glasses. The hardness values measured by nanoindentation (H) were systematically higher than the corresponding Vickers hardness by ∼1 GPa. This may be attributed to differences between the two methods: Vickers microindentation probes plastic deformation only, while nanoindentation probes the combined elastic-plastic response of a material, yielding higher hardness values for materials with a pronounced elastic contribution to the indentation deformation (Oliver and Pharr, 1992).

Crack Resistance
Crack resistance did not show any trend for increasing phosphate content ( Figure 10A). Data scattered between values of 4 and 8 N. By contrast, two different trends were found for the glasses with 0 and 7.5 mol% P 2 O 5 , respectively, when increasing the silica content ( Figure 10B). The crack resistance of glasses with 0 mol% P 2 O 5 slightly increased for increasing silica content. For the glasses with 7.5 mol% P 2 O 5 , crack resistance decreased markedly from 50 to 55.5 mol% SiO 2 and then remained constant within error limit for further increasing silica content.

Densification From Raman Analysis
Densification upon indentation resulted in a COG shift by 12-24 cm −1 of the low frequency Raman region in the considered spectral range of 200-650 cm −1 (exemplified in Figure 2, black spectrum before indentation, orange spectrum after indentation). For the two post-processing routes used, i.e., Long data reduction and baseline subtraction (filled and open symbols in Figures 10C,D), the COG shifts upon indentation showed little variation (<2 cm −1 ). For increasing phosphate content ( Figure 10C) the COG shift increased to a maximum value and then slightly decreased toward 7.5 mol% P 2 O 5 , independent of the data processing. For increasing silica content (Figure 10D), the COG shift increased for both glass series. This increase appeared to be linear within the error range for the series with 7.5 mol% P 2 O 5 (Figure 10D, orange circles), while it was non-linear for the series with 0 mol% P 2 O 5 (Figure 10D, black squares).
As an additional observation, an increase in phosphate content did not significantly affect the COG shift for the glasses with 50 mol% SiO 2 (Figure 10D), regardless of the data processing.

Densification From FEA
Previous studies indicate that densification saturation cannot be achieved using indentation testing (Bruns et al., 2020). The predicted saturation densification of 8.6% (see section Finite Element Analysis in the Materials and Methods section) will likely not be achieved in the indentation center. An experimental determination of densification based on the COG shift upon indentation requires correlation factors, usually determined in DAC experiments accompanied by Raman spectroscopy. Those factors are rarely available for glasses, except for fused silica (Sonneville et al., 2012;Deschamps et al., 2013), and it remains unclear if they are applicable to other glasses, such as the aluminosilicate glasses investigated here. In the absence of such data an alternative strategy needs to be proposed.
FEA offers the opportunity to fill this gap and to provide information on the densification field for a glass from the Si60Px series having a Poisson's ratio of 0.22 (Figure 11). Densification is exported from ABAQUS using the parameter PEQC4, a measure for the total volumetric inelastic strain. In this manner, a maximum densification of 5.88 ± 1.39% can be observed in the indentation center.

Strain-Rate Sensitivity
For increasing phosphate content ( Figure 10E) the strainrate sensitivity slightly increased. For increasing silica content ( Figure 10F) the strain-rate sensitivity appears to increase first up to a maximum value in the range between 50 and 60 mol% FIGURE 7 | Elastic properties obtained by Brillouin spectroscopy (black squares), ultrasonic echometry (blue circles) and nanoindentation (green triangles, only E) for glasses with 7.5 mol% P 2 O 5 and varying silica content: (from top to bottom) shear modulus (G), bulk modulus (K), Young's modulus (E) and Poisson's ratio (ν). SiO 2 , followed by a slight decrease with further increasing SiO 2 content for both glass series.

Compositional Analysis
The observed slight increase in Al 2 O 3 content and decrease in SiO 2 content compared to the nominal composition (Table 1) can be attributed to Al 2 O 3 crucible corrosion (Barz et al., 1996;Smith et al., 2014;Palles et al., 2016). The amount of Al 2 O 3 incorporation found here is less than for (silico-)phosphate and tellurite glasses, which have been found to incorporate about 3-9 mol% Al 2 O 3 during melting (Palles et al., 2016;Tagiara et al., 2017;Griebenow et al., 2018;Sawangboon et al., 2020). This difference may be explained by the high melt viscosity of the present glasses.
The lower P 2 O 5 content may originate from P 2 O 5 evaporation during melting (Griebenow et al., 2018;Sawangboon et al., 2020), but P 2 O 5 losses are small considering the high melting temperature. This confirms a previous finding that P 2 O 5 is not likely to leave metaluminous melts (London et al., 1993).

Density and Related Properties
The decreasing density and atomic packing density as well as the increasing molar volume with increasing phosphate content (Figure 4) may be explained by the fact that phosphate was incorporated by substituting P 2 O 5 for NaAlO 2 , i.e., substituting two phosphate tetrahedra for one aluminate tetrahedron chargebalanced by sodium. We note that all aluminate tetrahedra are supposed to have four bridging oxygen atoms (BO) because of the metaluminous composition (Al:Na = 1). Phosphate tetrahedra, on the other hand, may bear fewer BO. Phosphate groups have been reported to partially assume Q 4 P state by forming Al-O-P bonds (Ryerson and Hess, 1980;Gan and Hess, 1992;Kirkpatrick and Brow, 1995). Such bonds between tetrahedral aluminate and phosphate groups have been shown by a bond valence approach to be stable in the absence of an alkaline charge-compensator. In this approach, the Al 3+ provides 0.75 valence units (cation oxidation number / coordination number) to the bridging oxygen, while the P 5+ provides 1.25 valence units. Therefore, the nominal −2 charge of the bridging oxygen is compensated (Brow et al., 1993;Bertmer et al., 2000). Al-O-P bond formation was supported also by heteronuclear dipolar solid state NMR experiments . Formation of these bonds should liberate Na + ions to create new non-bridging oxygen atoms (NBO) rather than charge-balancing aluminate tetrahedra. As a consequence, the progressive substitution of P 2 O 5 for NaAlO 2 introduces more tetrahedral units per mole and simultaneously reduces the average number of BO per tetrahedron. An increase in phosphate content therefore involves a decreasing network polymerization. This depolymerisation was supported by recent results derived from a solid-state 31 P MAS NMR spectroscopic analysis . A glass network with fewer cross-links and a higher number of tetrahedra per mole is likely to be less compact. The observed decrease in density and atomic packing density as well as the increase in molar volume may be attributed to this decreased network polymerization, but also to the decreasing amount of interstitial sodium ions with increasing P 2 O 5 content and, with regard to the molar volume, to a higher number of tetrahedral species per mole. A similar decrease in density has been observed in basaltic glasses (Tarragó et al., 2018b).
For increasing silica content, silica tetrahedra (SiO 4 ) were substituted for aluminate tetrahedra charge-balanced by sodium (Na + [AlO 4 ] − ). As both tetrahedral units bear four BO each, the overall network polymerization is assumed to remain unchanged. This appears to correlate with the nearly constant molar volume, but the latter may also be explained by the constant number of tetrahedral groups, regardless of their degree of polymerization. The decrease in density and atomic packing density for increasing silica content is caused by the requirement for charge-balancing Na + ions in the vicinity of [AlO 4 ] − tetrahedra. Considering the similar ionic radii of Si 4+ (r Si = 26 pm) and Al 3+ (r Al = 39 pm) (Shannon, 1976), their corresponding tetrahedra will also be of comparable size, but a network with more aluminate tetrahedra is likely to occupy more inter-tetrahedral space for the accommodation of the charge-balancing sodium ions. Therefore, both density (ρ) and atomic packing density (C g ) changed, even though the molar volume remained roughly constant.
The invariance of C g within the error limits for all glasses with 50 mol% SiO 2 and for the phosphate-free glass with 60 mol% SiO 2 , even though ρ changed, is attributed to the C g model assumption of densely packed spheres. This model is strongly affected by the selected ionic radii and does not take into account directional bonds and NBO which may be present.
The trend observed for oxygen packing density (ρ ox ) was opposite to that for density (Figure 4). This is counter intuitive for oxide glasses with oxygen as the most abundant element. The concept of oxygen packing density, i.e., the amount of oxygen atoms per defined volume of glass, has been introduced by Ray (1974), who emphasized its usefulness for weakly polymerized alkali silicate glasses with alkaline ions acting as network modifier. The glasses studied here, however, are highly polymerized, and sodium ions take a charge-balancing role rather than a modifying one. Equivalent to C g , the calculation of ρ ox is based on the assumption of densely packed spheres and does not take the directional, covalent character of the oxygen bridging bonding into account. The dependence of oxygen packing density on stoichiometric numbers appears to be a disadvantage in the presence of three different network-forming oxides, all of them with different stoichiometric numbers of oxygen. Oxygen FIGURE 11 | FEA densification field for the system with constant 60 mol% SiO 2 , beneath the residual impression of a Berkovich equivalent conical indentation. As the indenter shape is self-similar, the shape of the densification field is independent from indentation depth. Therefore, no scale is required for the axes X,Y,Z. Densification is exported from ABAQUS using the PEQC4 parameter.
packing density is most useful for glass series with constant numbers of oxygen atoms per formula unit (Tylkowski and Brauer, 2013;Mandlule et al., 2014). Oxygen packing density has been found to depend on the coordination number of the cations (Zeidler et al., 2014;Salmon and Zeidler, 2015). However, recent solid-state MAS NMR spectroscopic investigations indicated that Si, Al and P remained in fourfold coordination in all studied glasses (Grammes et al., in preparation;Nizamutdinova et al., 2020).
Results for packing density and oxygen packing density may additionally be affected by the assumption of a constant ionic radius for oxygen (r O = 135 pm). The ionic radius of oxygen is known to vary depending on the adjacent cationic species (Duffy, 1990;Zeidler et al., 2014). Therefore, values found here could be erroneous. In general, a lower oxygen radius should enable a higher packing density and oxygen packing density. However, once the oxygen is surrounded by cations from more than one species, determination of an exact ionic radius becomes difficult if no quantitative data are available from bond-length measurements. We note that the calculation of oxygen packing density according to Equation (3) does not involve the ionic radius of oxygen [in contrast to the equation used by Zeidler et al. (2014) which requires an a-priori knowledge of the exact bondlengths] while the packing density calculation via Equation (2) does require the oxygen radius. This may also be a reason for the opposing trends between those two properties. As stated above, the trend in packing density appears more valid as it correlates with the measured density. Therefore, the oxygen packing density will not be taken into account in further analyses here.

Elastic Properties
The elastic properties determined for phosphate-free metaluminous glasses agreed well with available literature data on similar compositions. For Si50P0 and Si60P0 glasses, values of Young's modulus ranging from 72 to 74 GPa have been reported (here: 71-75 GPa for Brillouin spectroscopy and ultrasonic echometry) as well as Poisson's ratios from 0.200 to 0.228 (here: 0.213-0.233) (Tiegel et al., 2015;Weigel et al., 2016). To our knowledge, elastic properties have not yet been reported on comparable phospho-aluminosilicate glass compositions. Only peralkaline basaltic glasses with negligible phosphate contents of 0.2 mol% have been investigated (Lönnroth and Yue, 2009). The decrease in E, G and K for increasing phosphate content (Figure 5) may indicate a decrease in network polymerization. This trend may also be related to a decrease in atomic packing density upon phosphate incorporation, as suggested by Makishima and Mackenzie (1973). According to Rouxel, the packing density shows a non-linear correlation with elastic properties for a wide range of glass types (Rouxel, 2006(Rouxel, , 2007. For increasing silica content here, polymerization should remain unchanged, at least in the phosphate-free glasses. For the phosphate-containing glasses, interactions between Na, P and Al may slightly influence polymerization. If little to no change in polymerization is assumed, the elastic properties should remain largely constant. This assumption is well confirmed by the results for G and possibly for E (Figures 6, 7). Results show a clear difference in the trends for G and E with increasing silica content, showing minor changes only, and much more pronounced changes for increasing phosphate content. Despite constant network polymerization with increasing silica content, a clear decrease of K is observed in addition to the very minor decreases of G and E. These decreasing moduli may possibly be explained by the decrease in packing density (Figure 4). The decrease in packing density with increasing silica content indicates an increase in free volume. A glass with more free volume will be more compliant toward elastic deformation as indicated by the trend observed for K (and, possibly, G and E). To conclude, both network polymerization and atomic packing density appear to contribute to modulus changes, but only the change in packing density may explain all measured trends.
Poisson's ratio remained practically constant for increasing P 2 O 5 content ( Figure 5) and decreasing polymerization, even though Rouxel (2006) indicated an increase of ν for glasses with lower degree of cross-linking. However, his data showed large scattering for highly polymerized glasses. Rouxel's observations also include a decrease in Poisson's ratio with decreasing packing density. The results shown here do not confirm this observation, however, as the packing density showed a clear decrease with increasing phosphate content, while Poisson's ratio remained constant. By contrast, the slight decrease of ν for increasing SiO 2 content (Figures 6, 7) matches the relationship found by Rouxel (2006). There is actually a simple mathematical explanation for the constant values observed for ν, as it was calculated by dividing E by G. As E and G showed comparable changes with composition, no change in ν is observed. One additional observation may be noted, though. Rouxel's observation of a correlation between ν and packing density (Rouxel, 2006) was confirmed here when substituting SiO 4 for Na + [AlO 4 ] − , i.e., when the overall amount of network-forming units did not change. With increasing phosphate content, by contrast, this correlation was not observed, most likely because the overall amount of network-forming units increased, as NaAlO 2 was replaced by P 2 O 5 . However, to really test whether the overall ratio of network formers to network modifiers correlates with ν, one would need to test over a much wider range of glass compositions.
Young's moduli calculated using the Makishima-Mackenzie model showed a comparable trend to the experimental values (Figure 8), although the calculated values were lower than the experimental ones. The relatively good agreement in trend highlights that, despite its simplicity, the model provides a useful estimate of elastic properties. This is in agreement with an earlier study of the model on simple silicate, borate and phosphate glasses (Inaba et al., 1999). The deviation between experimental and model values may be caused by various factors.
First of all, the Makishima-Mackenzie model assumes that all elements in the glass have the same coordination as in their corresponding oxide components. This is not true for aluminum which should be present mainly in fourfold coordination in the studied metaluminous glasses, while the Makishima-Mackenzie calculation is based on the dissociation energy of Al 2 O 3 , in which aluminum is exclusively sixfold coordinated. The cation packing fraction of aluminum in fourfold and sixfold coordination is similar, yet the charge of Al 3+ is shared among fewer oxygens in fourfold coordination (She et al., 2019). This is in agreement with early work by Sun (1947) suggesting that the dissociation energy of Al-O bonds is higher when Al is in fourfold coordination compared to sixfold coordination. Therefore, an adjustment of the model taking into account Al coordination may potentially yield larger calculated values of E, leading to better agreement with experimental data.
Secondly, the Makishima-Mackenzie model presents a meanfield approximation. Using quantitative data on the work of deformation and its comparison to predicted bulk modulus, it has recently been shown for a wide range of glassy materials that this is overly simplified (Sawamura and Wondraczek, 2018). Instead, the mechanical response of glasses is strongly dependent on their respective structural heterogeneity (Benzine et al., 2018).
Another shortcoming of the Makishima-Mackenzie model is its implicit assumption that all bonds are bridging oxygen bonds. Neither non-bridging oxygen nor terminal oxygen are taken into account. Inaba et al. (1999) found this to affect the model values for phosphate glasses and successfully adjusted the model to take non-bridging oxygen atoms into account. However, in the present study the model underestimated the experimental values. An adjustment of the model to include non-bridging oxygens connected to phosphate groups would further decrease the calculated moduli, leading to even higher deviation from experimental values.
A further factor influencing the moduli may be bonds between mixed network-forming units, such as Al-O-P bonds, which are not taken into account by the Makishima-Mackenzie model either. A comparison between the molar heats of formation of AlPO 4 (−1733.8 kJ/mol, only Al-O-P bonds) and Al 2 O 3 (−1675.7 kJ/mol, only Al-O-Al bonds) (Lide, 2006) suggests that mixed bonds are stronger than non-mixed ones. This may partially explain why experimental moduli are higher than the calculated ones. An adjustment of the model to take mixed bonds into account would require detailed knowledge of the glass structure, however, to determine the relative amount of each type of bond.
A minor contribution to the deviation between model and experimental data may also result from actual glass composition slightly deviating from the nominal values used in the model. Furthermore, ionic radii used for the calculation of packing density were found to have a major influence on the results of the Makishima-Mackenzie model. For example, the two different sets of "effective ionic radii" and "crystal radii" provided by Shannon (1976) led to differences of more than 10 GPa in the calculated moduli. For this reason, the effective ionic radii used here were chosen as they yielded the best agreement with experimental data.

Hardness
The decrease in hardness with increasing phosphate content can be explained by lower network polymerization, caused by P 2 O 5 incorporation. The constant hardness for increasing silica content may be explained by constant network polymerization upon substitution of SiO 4 for Na + [AlO 4 ] − tetrahedra. Lönnroth and Yue (2009) have noted a weak correlation between hardness and density in peralkaline aluminosilicate glasses. Results here showed a correlation of hardness with packing density (Makishima and Mackenzie, 1973) with increasing P 2 O 5 content but not with increasing SiO 2 content. Instead, the constant hardness with increasing silica content could be related to the constant molar volume (Figures 4D,   9B,D). The decreasing hardness with increasing P 2 O 5 content (Figures 9A,C) was accompanied by increasing molar volume, indicating a relationship between large molar volume and low hardness for the present system. A similar decrease in hardness with increasing phosphate content up to more than 8 mol% has been reported for basaltic glasses (Tarragó et al., 2018b). Zeng et al. (2016) found hardness values between 5 and 5.5 GPa for peralkaline glasses with 5 mol% P 2 O 5 . This range agrees with the Vickers hardness of most phosphate-containing glasses presented here. The hardness values found for the phosphate-free glasses (6 to 6.3 GPa for microindentation) are in good agreement with a value of about 6 GPa reported by Tiegel et al. (2015) for a glass of the Si60P0 composition. Other literature hardness values of metaluminous glasses with various alkaline constituents ranged from 5.6 to 7.4 GPa (Smedskjaer et al., 2013;Tiegel et al., 2015;Bechgaard et al., 2016).

Crack Resistance
To estimate the error of the crack resistance, a method was developed here based on the PCI error bars. In Figure 1, PCI error bars represent the 95% confidence interval after evaluation of 25 indentations (up to 100 potential median-radial cracks) per load. The error of the crack resistance has to be found in the load dimension (horizontal axis). Therefore, the end points of the PCI error bars (blue triangles) were connected by linear segments to create an estimated error envelope of the PCI (red lines in Figure 1). This error envelope intersects the 50% PCI line at two loads above and below the crack resistance, generating two horizontal error bars (not shown) of different length (because of logarithmic axis). The larger error bar out of the two was chosen as the estimated error of the crack resistance.
The slight increase in crack resistance with increasing silica content in phosphate-free glasses ( Figure 10B) confirms previous findings available in the literature (Gross et al., 2009;Hermansen et al., 2013). Bechgaard et al. (2016) have recently reported crack resistances between 1 and 3 N for a glass of Si50P0 composition. This agrees with the value found here (3.1 ± 0.5 N). For metaluminous calcium aluminosilicate glasses with SiO 2 content increasing from 60 to 80 mol%, Gross et al. found an increase in crack resistance with values comparable to those obtained in this study (Gross et al., 2009). However, for >80 mol% SiO 2 , this trend is reverted as cracking changed to silicalike anomalous behavior with different crack systems becoming active. Up to now, there is no literature data available on the crack resistance of phospho-aluminosilicate glasses.
Higher crack resistances have previously been linked to higher silica content in the case of soda-lime silicate glasses (Hermansen et al., 2013), indicating the influence of network polymerization. For phosphate-free sodium aluminosilicate glasses, on the other hand, crack resistance was reported to decrease when increasing the Al/Na ratio from 0 to 1 (Bechgaard et al., 2016), even though this means approaching metaluminous composition and thus increasing polymerisation.
The range of data scattering observed for increasing phosphate content (Figure 10A) is similar to the range spanned by the crack resistance trends for increasing silica content ( Figure 10B). Therefore, the apparent trends in Figure 10B may be coincidental and will not be discussed further.
To increase reliability of crack resistance measurements, an even higher number of indentations would be needed to improve statistics. Also, atmospheric humidity may affect the measurements and should be avoided by measuring under protective atmosphere or in vacuum.

Densification
Changes in the low frequency Raman region have been correlated with bond angle changes and thus with density changes (Galeener, 1979;Rouxel et al., 2008;Deschamps et al., 2011;Cornet et al., 2017). In case of densification, this region shifts toward higher wavenumbers; however, the theoretical basis for this correlation has only been strictly defined for one-component glasses such as vitreous silica (Sen and Thorpe, 1977;Galeener, 1979). The phospho-aluminosilicate glasses studied here carry three different network formers, all of them possessing different Raman scattering cross-sections. This may influence the spectral shifts to a currently unpredictable degree. Even so, in the past the analysis of densification with Raman spectroscopy has not only been exercised on vitreous silica (Kailer et al., 1999;Hehlen, 2010;Cornet et al., 2017), but also on soda-lime silicate glasses (Deschamps et al., 2011;Kassir-Bodon et al., 2012), borosilicate glasses (Winterstein-Beckmann et al., 2014) and aluminosilicate glasses (Poe et al., 2001;Bechgaard et al., 2016). All of these studies correlated changes of the low frequency Raman region with densification. Even the studies on borosilicate and aluminosilicate glasses have not tackled the issue of having more than one network former. Instead, a phenomenological approach has been chosen for these glass systems. As long as silica remained the main network former, the shape of the low frequency Raman envelope was close enough to that of vitreous silica to assume similar behavior. The same assumption will be considered valid for the phospho-aluminosilicate glasses studied here, because of their high degree of polymerisation and the clearly observable changes of the low frequency Raman region upon indentation.
Changes in the medium frequency region [Raman band at 700 cm −1 in Figure 2, assigned to [AlO 4 ] − stretching (Iwamoto et al., 1978;McMillan and Piriou, 1982;Kamitsos et al., 1994)] were not taken into account here. This band would be affected by a change of Al coordination, but there was no reason to assume such a coordination change in the present metaluminous glasses. The formation of higher-coordinated aluminum upon compression is likely for peralkaline compositions but not for metaluminous ones (Hochella and Brown, 1985;Bechgaard et al., 2016). Instead, a weak intensity increase in the medium frequency region upon densification is assumed to be caused by partial overlap with the low frequency bands.
When analyzing densification by Raman spectroscopy, a comparison between the baseline method by Deschamps et al. (2011) and the Long data reduction seemed worthwhile since the use of the Long data reduction is popular in Raman analysis of glasses, yet changes markedly the shape of the low frequency Raman region (Seifert et al., 1982;Neuville and Mysen, 1996;Mysen, 1999;Mysen and Toplis, 2007 Thus, the comparison allowed to test whether the Long data reduction affected the densification analysis. Concluding from the good agreement between the two tested post-processing routes (open and filled symbols in Figures 10C,D), different ways of determining the COG shift may still lead to comparable results regarding densification. Both tested post-processing routes may be regarded as equally applicable.
Raman spectroscopy allowed to track differences in densification between the different compositions, based on different COG shifts. The increasing COG shift upon phosphate incorporation in glasses with 60 mol% SiO 2 (Figure 10C) indicates higher densification upon indentation. This may be explained by the lower density of the phosphate-containing glasses (Figure 4). When equally loaded, a less densely packed glass should be more readily densified than a densely packed one. The maximum trend of the COG shift with increasing P 2 O 5 content ( Figure 10C) may be related to increasing glass depolymerisation upon increasing P 2 O 5 content. As shown by Bechgaard et al. (2016), the degree of polymerisation may affect the deformation mechanism of aluminosilicate glasses.
The increasing COG shift with increasing silica content ( Figure 10D) may be explained by the lower packing density of the high-silica glasses. The apparent non-linearity of this trend may potentially be caused by underlying changes of the deformation mechanism. The slightly higher values of Poisson's ratio for the glasses with lower silica content (Figures 6, 7) support a lower degree of densification and a higher degree of shearing for the low-silica glasses, following Sellappan et al. (2013).
The transformation of COG shifts into absolute values of densification is established for vitreous silica (Deschamps et al., 2013) and soda-lime silicate glass (Deschamps et al., 2011) only. For more complex glass systems such as the given phosphoaluminosilicate glass series no direct correlation is available in the literature. In the present study the degree of densification was estimated using FEA for the 60 mol% SiO 2 glass series and a maximum indentation densification of 5.88 ± 1.39% was observed (Figure 11). The FEA input is strongly linked to Poisson's ratio, which is almost constant around 0.22 for the 60 mol% SiO 2 glass series ( Figure 5). Therefore, the FEA approach does not exhibit the precision to reproduce the compositional influences as they have been found for the COG shift in the given test series. Both approaches should be considered together. The sensitive COG shifts provide insight into how compositional changes lead to relative differences in densification upon indentation. FEA provides a rough idea of the absolute level of indentation densification on which those changes take place.
Raman signal is gathered not only from the sample surface but also from below, where densification is lower (Figure 11) (Bruns et al., 2020). COG shifts therefore correspond to a densification averaged over signal depth. Compositional trends of COG shift should not be affected by this because indent size, load, and Raman depth of field were roughly equal for all samples.

Strain-Rate Sensitivity
All strain-rate sensitivities found here fit well within the relationship between strain-rate sensitivity and Poisson's ratio for glasses shown by Limbach et al. (2014). The strain-rate sensitivities ranged from 0.012 to 0.018 for values of ν between 0.20 and 0.24, agreeing with similar values for soda-lime silicate glasses (Limbach et al., 2014(Limbach et al., , 2017. These strain-rate sensitivities are rather low, yet common over many different glass systems (Limbach et al., 2014). The slightly increasing strain-rate sensitivity for increasing phosphate content ( Figure 10E) may potentially be explained by the decreasing packing density and network polymerisation upon phosphate incorporation. Both trends should account for more free volume within the glass. This additional free volume together with the lower degree of cross-linking may facilitate local structural rearrangements when subjected to sudden changes in strain-rate. Nevertheless, a complete structural explanation of the increased strain-rate sensitivity has to be more complex. For increasing silica content, the results show no correlation between strain-rate sensitivity and density. Here, the non-linearity of strain-rate sensitivity correlates with no other property. A tentative structural hint may be obtained from metaluminous minerals of stoichiometry comparable to that of the phosphate-free glasses. Glass Si50P0 has the same stoichiometric composition as the mineral nepheline (NaAlSiO 4 ). The composition of glass Si60P0 lies between nepheline and jadeite (NaAlSi 2 O 6 ), the latter corresponding to 66.67 mol% SiO 2 . The composition of glass Si70P0 lies between jadeite and albite (NaAlSi 3 O 8 ), with albite corresponding to 75 mol% SiO 2 . The crystal structures of nepheline, jadeite and albite are hexagonal, monoclinic and triclinic, respectively, see e.g., the RRUFF mineral database, entries R040025.1, R050220 and R040068.1 (Lafuente et al., 2015). Potentially, the maximum in strain-rate sensitivity upon increasing silica content may correspond to an underlying change within the glasses' mediumrange structure which bears similarities to the change of the corresponding crystal structure. However, this is a mere hypothesis and requires further study.

CONCLUSIONS
This work has reported the synthesis and study of mechanical properties of glasses in the metaluminous (Al:Na = 1) system SiO 2 -Al 2 O 3 -Na 2 O-P 2 O 5. The incorporation of phosphate into such aluminosilicate glasses by substituting P 2 O 5 for NaAlO 2 resulted in decreasing packing density. A simultaneous decrease in hardness and elastic moduli and increase in strain-rate sensitivity could be attributed to this packing density decrease, and to depolymerisation of the glass network.
Regarding glass densification upon indentation, glasses with lower initial density showed higher shifts of the correlated Raman spectral region, indicating higher densification. The results were non-linear over glass composition, which was tentatively attributed to a change in deformation mechanism.
Variation of the silica content by substituting SiO 4 for Na + [AlO 4 ] − tetrahedra was assumed to maintain nearly constant network polymerisation. This was supported by constant hardness. Simultaneous changes in elastic moduli correlated with packing density, while changes in Poisson's ratio were small. Only in case of increasing SiO 2 content they could be correlated with packing density changes.
Taken together, the findings of this first comprehensive mechanical analysis of metaluminous phospho-aluminosilicate glasses provides a basis for tailoring mechanical properties by simply varying the P 2 O 5 or SiO 2 content. Non-linearities found in densification and strain-rate sensitivity hint toward more complex changes in glass structure, opening up possibilities for further study of the interplay between the three network formers.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
TG, DB, LvW, and DL designed the glass series and the experiments. TG synthesized the samples, performed density, microindentation and Raman spectroscopy experiments, contributed to Makishima-Mackenzie calculations, analyzed and interpreted the results, and wrote the manuscript. RL performed nanoindentation and ultrasonic echometry measurements and contributed to Makishima-Mackenzie calculations, interpretation of mechanical properties, and writing of the manuscript. SB performed the FEA analysis and contributed to the interpretation of glass densification and to writing the corresponding part of the manuscript. LvW, DL, EK, KD, LW, and DB contributed to the interpretation of the results and writing of the manuscript.

ACKNOWLEDGMENTS
We would like to thank Alexander Veber (formerly Friedrich Alexander University Erlangen-Nürnberg) for Brillouin measurements, Johannes Buchheim (Institute of Technical and Environmental Chemistry, Friedrich Schiller University Jena) as well as Detlef Schelle (Institute of Applied Physics, Friedrich Schiller University Jena) for EDX measurements, Timothy Oshiobughie Imogore (Friedrich Schiller University Jena) for implementing the center of gravity calculation into MATLAB, Tadiyos Tarekegn Alemayehu (formerly Friedrich Schiller University Jena) for assistance with crack resistance measurements and two reviewers for their valuable comments. The majority of data shown here, with exception of the Makishima-Mackenzie calculations and the FEA simulation of densification, was published in TG's Ph.D. thesis (Grammes, 2020).