Cerium oxide (ceria, CeO2) is frequently used as a standard in applications such as synchrotron and x-ray free electron lasers for calibrating x-ray wavelengths and offers the potential for understanding the high pressure properties and deformation mechanisms in a wide range of similar face centered cubic (fcc) materials. In this study, the pressure dependence of the strength of ceria was investigated up to 38 GPa using angle dispersive x-ray diffraction in a radial geometry in a diamond anvil cell. In this experiment, the difference in the stress along the axis of compression and perpendicular to the direction of compression can be determined, giving a quantity known as the differential stress. It was found that the differential stress (t), a measure of the lower bound for yield strength, initially increases rapidly from 0.35 ± 0.06 GPa to 2.2 ± 0.4 GPa at pressures of 1.8 and 3.8 GPa, respectively. Above 4 GPa, t increases more slowly to 13.8 ± 2.6 GPa at a pressure of 38 GPa. The changes in the preferred orientation (texture) of CeO2 with pressure were also measured, allowing for the determination of active deformation mechanisms using an elasto-viscoplastic self-consistent model (EVPSC). It was found that as pressure increased, the  direction had a slight preferred orientation along the axis of compression. Our EVPSC model of experimental fiber (cylindrically symmetric) textures and lattice strains were most consistent with dominant slip activity along .
Cerium oxide (CeO2, ceria) is a face centered cubic (fcc) ceramic material that is used in a large variety of applications, including x-ray wavelength standards,1 fuel cells,2 catalysts,3 solar cells,4 biological applications,5 and automotive exhaust treatment.6 The deformation mechanisms of fcc materials are of general interest in materials science, with many technologically important materials such as high entropy alloys,7 perovskites,8 salts,9 and oxides10 all sharing an fcc structure at ambient and high pressures.
The high pressure behavior of CeO2 has been shown to be affected by temperature, sample environment hydrostaticity, and particle size and morphology. At ambient pressure, CeO2 crystallizes in a face centered cubic fluorite (Fm3m) structure. Under high pressures, CeO2 undergoes a phase transition to an orthorhombic structure.11–16 The structure of this orthorhombic phase has been reported to be either Pnam11,12,17 or Pbnm.18 The pressure range at which this phase transition occurs is known to be affected by temperature,12 pressure transmitting medium used,11 and the particle size.13–16 Increasing the temperature has been shown to decrease the pressure at which the phase transformation begins.12 The pressure transmitting medium used affects the amount of deviatoric stress present in the diamond anvil cell (DAC) sample chamber. Deviatoric stress is the difference between the full stress tensor that a material experiences and the hydrostatic component of that stress. Higher amounts of deviatoric stress, such as when using no pressure transmitting medium, reduce the phase transformation pressure to 31–38 GPa.11 Use of a quasi-hydrostatic pressure transmitting medium is known to increase the phase transition pressure with a phase coexistence region of 35.6–73.0 GPa in He,12 and 31–55+ GPa in a 4:1 methanol:ethanol solution.19 Grain size of CeO2 also affects the phase transformation pressure with the limit of a single crystal sample reducing the coexistence region of the fluorite and orthorhombic phases.13 Compressing CeO2 with particle sizes in the nanoscale regime (4–15 nm particle size) had varied results, with some studies reporting an increase in the phase transformation pressure,13,14 and other studies reporting a decrease in the phase transformation pressure.15,16
In addition to studying the phase transformations of CeO2, two studies have reported the relationship between pressure and yield strength of CeO2. The first of these studies determined the strength from the linewidth of powder x-ray diffraction (XRD) peaks as a function of pressure.18 The linewidth of the peak is related to both the microstrain and crystallite size of a material. If the microstrain and the single crystal elastic moduli are known, a lower bound for yield strength can be determined. Using this method, Liu et al.18 reported that the yield strength of CeO2 increases rapidly until 3.3 GPa and then remains almost constant until the maximum pressure studied of 27 GPa. A second study by Lang et al.20 determined the pressure dependence of the strength of CeO2 using gas-gun shock compression up to 25 GPa. A similar plateau in the strength of CeO2 at pressures above 8 GPa was observed in this study.
In the present work, the relationship between yield strength of CeO2 and pressure is determined using uniaxial compression in a DAC coupled with XRD in a radial geometry, i.e., with incident x rays oriented perpendicular to the compression axis of the DAC. Additionally, the texture development of CeO2 under uniaxial compression is determined through full-profile Rietveld refinement of the diffraction patterns. This texture development was further analyzed to determine the active deformation mechanisms using the physics-based micromechanical elasto-viscoplastic self-consistent (EVPSC) model.21
II. EXPERIMENTAL METHODS AND ANALYSIS
A. Experimental details
A BX-90-style DAC,22 modified to enable x-ray diffraction in the radial geometry, was used in this study. The DAC contained two 300 μm flat culet type Ia diamonds, and an x-ray transparent amorphous boron and epoxy gasket (50 μm diameter hole, 50 μm initial thickness) was used to compress the CeO2 sample to high pressures. The CeO2 powder sample used had 99.99% purity and a 14 μm particle size (Alfa Aesar). A Pt flake, loaded with the sample, was used as an internal pressure standard. The DAC was compressed using a remotely controlled gas membrane system23 at the HPCAT (High Pressure Collaborative Access Team, Sector 16) ID-B beamline at the Advanced Photon Source at Argonne National Laboratory.24 XRD was collected in the radial geometry, which is a commonly used method to determine the yield strength of a material under static high pressures.25 The wavelength of the x-ray beam used was 0.4066 Å, the spot size was 3 × 6 μm2, and the sample to detector distance was approximately 205 mm. The XRD images were collected using a Pilatus 1M-F detector with a 3 s accumulation time. 66 diffraction images were collected at pressures between 1.7 and 39 GPa in roughly 0.5 GPa steps.
B. Radial x-ray diffraction analysis
During compression, a sample is subjected to a cylindrically symmetric stress (uniaxial compression) which can be quantified using three different stress tensor elements that are perpendicular to each other (i.e., ). By convention, the stress along the compression axis is represented by and, by symmetry, the stresses in the two directions perpendicular to the direction of compression are equal (i.e., ). The stress a material experiences during compression can be separated into hydrostatic () and deviatoric stress components [the last term on the right hand side of Eq. (1)],26
Here, is the hydrostatic stress component and is the differential stress component.26
When measuring XRD in the radial geometry in a DAC, the powder diffraction rings from the sample are distorted from the perfectly circular Debye–Scherrer rings typically observed in powder diffraction experiments. Specifically, along the compression axis, the lattice spacings are compressed more than those perpendicular to the axis of compression, thereby resulting in elliptical powder diffraction rings. These changes in the lattice spacings around the diffraction rings are modeled using
Here, dm(hkl) is the measured lattice spacing, dp(hkl) is the predicted lattice spacing under hydrostatic compression, is the azimuthal angle around the diffraction pattern in a counterclockwise sense with denoting the direction of compression.27,28Q(hkl) is the lattice strain parameter that is defined by
In Eq. (3), is the shear modulus in the Reuss limit where stress continuity is assumed across the grains of the sample, and is the shear modulus at the Voigt limit where strain is continuous across the grains of the sample.27,28 Empirically, the experimental value of the shear modulus has been found to be well-approximated by the average of and , implying a value of .26 This allows for t to be determined by
where G is the shear modulus.26 The von Mises yield criterion states that a material will yield if the deviatoric stress of a material is greater than or equal to the yield strength.29 If the deviatoric stress applied to the material during compression in a DAC is less than the material's yield strength, the material will deform elastically rather than plastically. As plastic deformation cannot be directly visualized in situ within the DAC, the von Mises yield criterion implies that t ≤ σy, where σy is the yield strength, and that t can be used as a lower bound for the material's yield strength.
Ultrasound time-of-flight measurements, performed in situ at HPCAT beamline 16 BM-B in a large volume Paris-Edinburgh (PE) style hydraulic press, were used to determine the shear sound speed as a function of pressure.30 The shear modulus was determined from the shear sound speed, VS, using
where is the mass density measured by XRD. Details of the ultrasound experiments can be found in the supplementary material. The slope of Vs(P) changes at roughly 3.3 GPa (Fig. S1 in the supplementary material), and the data below this pressure (low pressure regime) were considered separately from the data above 3.3 GPa (high pressure regime). The shear sound speed and its pressure derivative were determined to be Vs = 3390 ms−1, dVs/dP = 32.17 m s−1 GPa−1 in the low pressure regime and Vs = 3304 m s−1, dVs/dP = 57.95 m s−1 GPa−1 in the high pressure regime. This results in a value of G(P = 0 GPa) = 82.7 GPa using the low pressure regime data and G(P = 3.3 GPa) = 89.1 GPa for the high pressure regime [G(P = 0) = 79.0 GPa if the high pressure regime data were extrapolated to 0 GPa]. The evolution of G (and other elastic constants) with pressure is shown in Fig. S2 in the supplementary material. G = 86.2 GPa at 1.7 GPa and increases to 249 GPa at 38 GPa. The values of G used to determine t were calculated using the exact values of density from the XRD and extrapolation of Vs from the ultrasound measurements.
Figure 1(a) shows an example of a diffraction image collected during the experiment. The 2θ and coordinates are highlighted on this image. The diffraction images collected during the experiment were integrated over ψ using a custom Matlab script. This script integrated the diffraction pattern in 72 sectors around , each spanning 5°. The integration process resulted in 72 one-dimensional diffraction patterns. A representative example of the integrated diffraction patterns with intensities represented as colors is shown in spectrum numbers 0–72 in Fig. 1(b).
The Material Analysis Under Diffraction (MAUD)31 software package was used to perform a full-profile Rietveld refinement on the integrated data as described by Wenk et al.32 From the refinement, the hydrostatic lattice parameters of CeO2 and Pt, Q(hkl) of CeO2, and the texture of the sample were determined. The texture of the sample was refined using the entropy maximization William–Imhof–Matthies–Vinel (E-WIMV) iterative process.33 The orientation distribution function (ODF) resolution was set to 7.5° when refining the texture. Due to the radially symmetric stress distribution, fiber symmetry was imposed when refining the texture. An example refinement is shown in spectrum numbers 73–144 in Fig. 1(b).
C. EVPSC modeling
To interpret the texture development of the sample and determine the activity of different slip systems, elasto-viscoplastic self-consistent (EVPSC) modelling was performed.34 The initial texture for the EVPSC model was generated from the texture observed in the experiment at 2.7 GPa. The elastic moduli at 2.7 GPa (see the supplementary material) were also used as initial parameters for the model. The samples were modeled to contain 2000 grains to approximate a polycrystalline material. In this model, each grain was treated as an inclusion in a host medium having anisotropic elasto-plastic properties. During deformation, interactions between inclusions and the medium are approached iteratively until the macroscopic stress distribution is equal to the average of the microscopic stress distribution in each inclusion and throughout the medium. To achieve a total strain of 20%, we ran the model in 100 steps of 0.2% strain.
EVPSC accounts for both the elastic response of uniaxial compression and the relaxation within a crystal lattice due to plastic deformation along rate-dependent slip systems. A velocity gradient (v) of
was used to achieve a good fit to experimental pressures and to represent the quasi-uniaxial compression occurring within a diamond anvil cell. The plastic strain rate of the grain is the sum of the shear rates contributed by each slip system. The plastic strain rate is given by a visco-plastic law, as a power n of the resolved shear stress (RSS),35
In Eq. (7), the strain rate tensor of the plastic component, “p,” is given by , is the threshold value of shear stress along a slip system “s” required to induce slip at the reference shear rate, . The RSS is given by the product of two tensors: the symmetric Schmid factor , and the local stress tensor for an individual inclusion, . As usual, repeated indices are summed over. Variable n is a stress exponent that describes the non-linearity of rate-sensitive deformation. Using the absolute value of the RSS and carrying the sign separately ensures that the equation is valid for both even and odd powers of n. For the simulations performed in this work, a value of was used because this value provided a good match for the relationship between Q(hkl) and pressure and a large n is generally consistent with room temperature deformation. For the simulations performed in this work, a value of was used because this value provided a good match for the relationship between Q(hkl) and pressure and a large n is generally consistent with room temperature deformation.
The critical resolved shear stress (CRSS, ) is a determined value that is affected by both strain hardening and pressure36 and is modeled as
Here, is the initial CRSS at 2.7 GPa and and are the first and second order pressure derivatives, respectively.
Using a self-consistent model such as EVPSC requires linearization of the rate-sensitive equation (7). Several approaches are possible such as those described in Wang et al.34 and Lin et al.36 A tangent homogenization approach was used for this modeling, allowing for more strain heterogeneity between inclusions. This approach was chosen because it most accurately approximates a uniform microscopic stress state.
The grain rotation and resulting orientations were calculated in EVPSC for each inclusion. This can be used to describe the overall texture development of the modeled polycrystal. Initial textures were imported into EVPSC by creating weighted single-crystal orientations in the software BEARTEX37 from the experimental textures output by MAUD as described in Sec. II B and shown in Fig. 1(b). Importing experimental textures yields more realistic texture development when initial crystallite orientations are non-random rather than starting with a statistically random distribution. Simulated textures are exported from EVPSC also as weighted single-crystal orientations and plotted with the MTEX software package for MATLAB.38
Figure 2(a) shows Q(hkl) as a function of pressure for five different values of (hkl) that had strong enough diffraction peaks to be reliably refined throughout the entire pressure range. The markers represent the values of Q(hkl) determined using MAUD to Rietveld refine the experimental data, while the solid lines show the EVPSC model results. At the lowest pressure of 1.8 GPa, the average value of Q(hkl) and the standard deviation are 7 × 10−4 and 1 × 10−4, respectively. Between 1.8 and 3.2 GPa, the five Q(hkl) are very close in value to each other and rise steeply to an average value of 2.6 × 10−3 and standard deviation of 4.41 × 10−4. At higher pressures (above 3.2 GPa), the values of Q(hkl) for the different planes begin to diverge. The order of Q(hkl) from the highest value to the lowest value [specifically, Q(111) ≈ Q(331) > Q(220) > Q(311) > Q(200)] is generally consistent between 3.2 GPa to the highest pressure of 38.1 GPa. Above 3.2 GPa, the average values of Q(hkl) increase with pressure at a slower rate. It should be noted that at some pressures (20, 30, and 36 GPa) there is a temporary drop in Q(hkl) for all hkl values. This decrease in Q results from realigning the x-ray beam and detector with the sample, which moves slightly during the course of compression. Different regions of the sample experience different local stress conditions, leading to discontinuous changes in Q(hkl) and the differential stress.
Figure 2(b) shows the relationship between pressure and differential stress (t). Initially a sharp increase in t is observed between 1.8 and 3.8 GPa, where t increases from 0.35 ± 0.06 GPa to 2.2 ± 0.4 GPa. Above 3.8 GPa, t increases more slowly with pressure until reaching a maximum of 13.8 ± 2.6 GPa at a sample pressure of 38.2 GPa. Similar to the values of Q(hkl), a noticeable drop in t is observed at sample pressures of 20, 30, and 40 GPa. For comparison, t determined using the shear modulus and its pressure derivative from Ref. 39 and the yield strength of CeO2 from two previous studies are also shown in Fig. 2(b).18,20
Figure 3(a) shows a series of inverse pole figures (IPFs) used to visualize the texture development in CeO2 during compression. The IPFs show how different crystallographic directions are oriented with respect to the direction of compression by visualizing texture in each direction as a multiple of a random distribution (m.r.d.). An m.r.d. of 1 indicates that a crystallographic direction is randomly oriented with respect to the axis of compression. M.r.d values greater or less than 1 indicate alignment along the compression axis, and values less than 1 indicate than an orientation is aligned away from the axis of compression. A smoothing factor of 10 was used to generate the IPFs using MTEX.38 As seen in Fig. 3(a), orientations of sample crystallites are near random at a pressure of 2.7 GPa. At 19 GPa, a slight texturing observed for the  and  directions along and away from the compression axis, respectively. During compression, sample texture continues to develop up to 37 GPa, the maximum pressure reached. At this final pressure, a trace amount of the high pressure orthorhombic phase is observed. This phase is too weak to refine in MAUD.
Figures 2(a), 3(b), and 4 show the modeled change in Q(hkl) with pressure, the texture development and changes in the slip system activity with pressure, respectively. From Figs. 2(a) and 3(b), the modeled results for Q(hkl) and texture are seen to be in qualitative agreement with the experimental data. Tables I and II provide the parameters used for the best fit model, and the elastic moduli and their derivatives used in the model, respectively. The values of the elastic moduli (Cij) were obtained by using G [determined by Eq. (5)] and the bulk modulus (B) using the equations
B was determined by
where vl is the longitudinal sound speed. From the ultrasound measurements, vl = 6357 m s−1, dvl/dP = 76.2 m s−1 GPa−1. The relationships between B, G, and Cij with pressure are provided in the supplementary material (Fig. S2 in the supplementary material). The relationships between B, G, and Cij with pressure are provided in the supplementary material (Fig. S2 in the supplementary material). Because 97% of the XRD measurements are in the high pressure regime (above 3.3 GPa), the elastic moduli used in the modeling were calculated from a second order polynomial fit to the experimental data in the high pressure regime. The fits to Cij vs pressure to the experimental data are provided in Table II. The relationship between pressure and strain from the EVPSC modeling is provided in Fig. S3 in the supplementary material.
|Slip system .||(GPa) .||.||(GPa−1) .|
|Slip system .||(GPa) .||.||(GPa−1) .|
|Elastic modulus .||Cij (GPa) .||Cij׳ .||Cij״ (GPa−1) .|
|Elastic modulus .||Cij (GPa) .||Cij׳ .||Cij״ (GPa−1) .|
Texture results are shown in Fig. 3(b) for 2.7, 19, 29, and 38 GPa. These models have similar textures from the experimental data, with the  direction becoming more oriented toward the axis of compression as pressure increased. Note that the texture is slightly more pronounced in the EVPSC modeling than in the experimental data which is likely because EVPSC does not model intergranular interactions, which tend to impede preferred orientation development.
Figure 4 shows the change in the relative slip system activity of the , , and slip systems with pressure, which can be used to determine the most likely deformation mechanisms responsible for texture development. Deformation initially occurs by shearing along the slip system. As pressure increases, the relative activity of the slip system declines to ∼25% at 6 GPa and ∼15% at 39 GPa (strain = 0.202). Over the same pressure range, the activity of the slip system increases from ∼20% to a maximum value of 80% at 14 GPa (strain = 0.10) at which point it is principally responsible for deformation in the sample. Secondary slip along occurs between 5 and 39 GPa. Crossover of the relative activity of the and slip systems corresponds to experimental pressures at which Q(hkl) values diverge [Fig. 2(a)] and the slope of the differential stress, t, changes [Fig. 2(b)].
The relationship between pressure and yield strength determined by Liu et al.18 is shown in blue diamonds in Fig. 2(b). Similar to the data obtained in this study, Liu et al.18 found, using a technique measuring the width of peaks obtained in the axial XRD geometry (i.e., with incident x rays parallel to the DAC compression axis), that the strength of CeO2 increased rapidly until a sample pressure of 3.3 GPa. Above the sample pressure of 3.3 GPa, the results reported by Liu differed from those of the present study. Instead of continuing to increase with pressure, t was reported to remain relatively constant at 1.6 GPa. The difference in these results can be partially attributed to the difference in the elastic moduli used. Liu et al.18 used values of the elastic constants derived from density functional theory. In this case, the shear modulus and its pressure derivative were G0 = C44 = 68.10 GPa and dC44/dP = 1.13.39 In the present work, the shear modulus (Fig. S2 in the supplementary material) determined by the extrapolation of sound speed measurements and the density measured by XRD were used to calculate t from Eq. (4). To enable comparison to the other values of C44 = G and dC44/dP in the literature, a second order polynomial fit to our values of G gives G0 = C44 = 79.0 GPa, dC44/dP = 2.9, and d2C44/dP2 = 0.0398 GPa−1. These values are significantly higher than those used by Liu et al. Determined values of G(P) can vary significantly depending on several factors, including stock material purity and density, poly- or single-crystal material, the development of preferred orientation of the polycrystal during compression, and the technique used to determine G, i.e., theoretical vs experimental (and which experimental method is used). The values of G(P) determined by ultrasonic measurements in the Paris–Edinburgh press used the same CeO2 stock material as the XRD measurements and were performed at MHz frequencies, where effects of dispersion are minimal. To facilitate comparison of determined values of strength between studies, values for t were also calculated using the shear moduli reported by Liu and Frost, who used published density functional theory (DFT) results and Brillouin spectroscopy, respectively.
Values of t using the G(P) from Liu are shown in Fig. 2(b) in red circles. The pressure dependence of the shear modulus of CeO2 obtained by Frost et al.17 using Brillouin spectroscopy was fit with a linear equation to give the values of G0 = 63 GPa and dG/dP = 1.04. Because the sound velocity was found to be insensitive to grain size, both the microcrystalline and nanocrystalline data from this study were used in this fit. Values of t derived from these values of G are shown by green stars. Though the DFT and Brillouin spectroscopy values of shear modulus result in values of t that are lower than the values obtained in the present work, they are still higher than the values obtained by Liu et al.,18 suggesting that the remaining differences may be due to other factors such as the sample not achieving plastic flow, or differences in grain size or morphology.
Lang et al.20 measured the strength of CeO2 under shock compression in a two stage gas gun. In Lang et al., the strength was determined using both the Tresca and Griffith yield criterion. The Tresca yield criterion is used for when a material yields in a ductile manner and the Griffith yield criterion is for when a material yields in a brittle manner.40 The behavior of strength vs pressure was found to be similar to the behavior reported by Liu et al. with the strength plateauing above a sample pressure of 8 GPa at 2.6 GPa for the Tresca yield strength and 1.04 GPa for the Griffith yield strength. The Tresca yield strength is in agreement with the 2.7 GPa measured at a sample pressure of 4.3 GPa in this experiment. The difference in yield strength might be due to the higher temperatures reached in dynamic compression compared with the samples remaining at a constant 298 K in the quasi-static experiments reported in this work.
Slip system activities in the best-fit EVPSC model suggest that the slip system is the dominant deformation mechanism in CeO2 up to 37 GPa at 300 K with secondary contributions from the and slip systems. Other face centered cubic materials exhibit similar behavior such as Si where slip along was found to be the dominant deformation mechanism.41,42
The differential stress (t) of CeO2 was found to initially increase quickly from 0.35 ± 0.06 GPa to 2.2 ± 0.4 GPa between 1.8 and 3.8 GPa. At higher pressures, t was found to increase at a slower rate up to 13.8 ± 2.6 GPa at 38 GPa. Changes in texture of CeO2 with pressure were observed, where the  direction increasingly aligned with the axis of compression, and the  direction progressively aligned perpendicular to the axis of compression. This change in texture is consistent with the activation of the slip system, which is a commonly activated slip system in face centered cubic materials.
See the supplementary material for the details of the ultrasound measurements used to determine the shear modulus as a function of pressure.
This work was supported by the Laboratory Directed Research and Development program of Los Alamos National Laboratory (LANL) and the G. T. Seaborg Institute. LANL is operated by the Triad National Security, LLC for the DOE-NNSA under Contract No. 89233218CNA000001. Portions of this work were performed at HPCAT (Sector 16), Advanced Photon Source (APS), Argonne National Laboratory. HPCAT operations are supported by DOE-NNSA's Office of Experimental Sciences. The Advanced Photon Source is a U.S. Department of Energy (DOE) Office of Science User Facility operated by Argonne National Laboratory under Contract No. DE-AC02-06CH11357. L.M. acknowledges support from the NSF (No. EAR-1654687) and U.S. Department of Energy National Nuclear Security Administration through the Chicago-DOE Alliance Center (No. DE-NA0003975). We thank Mungo Frost for sharing his CeO2 shear modulus data. We would like to thank Bethany Chidester for helpful technical discussions and the two anonymous reviewers whose feedback improved the manuscript.
Conflict of Interest
The authors have no conflicts to disclose.
Larissa Q. Huston: Conceptualization (equal); Data curation (equal); Formal analysis (lead); Writing – original draft (lead); Writing – review & editing (lead). Samantha C. Couper: Data curation (equal); Formal analysis (equal); Writing – original draft (equal); Writing – review & editing (equal). Matthew Jacobsen: Data curation (supporting); Formal analysis (supporting); Writing – original draft (supporting). Eric K. Moss: Data curation (equal). Lowell Miyagi: Data curation (supporting); Formal analysis (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Jesse S. Smith: Data curation (supporting); Resources (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Nenad Velisavljevic: Data curation (supporting); Formal analysis (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Blake T. Sturtevant: Conceptualization (equal); Data curation (equal); Formal analysis (supporting); Funding acquisition (lead); Writing – original draft (supporting); Writing – review & editing (supporting).
The data that support the findings of this study are available from the corresponding author upon reasonable request.