Yield strength of CeO₂ measured from static compression in a radial diamond anvil cell

Cerium oxide (CeO₂, ceria) is a face centered cubic (fcc) ceramic material that is used in a large variety of applications, including x-ray wavelength standards,¹ fuel cells,² catalysts,³ solar cells,⁴ biological applications,⁵ and automotive exhaust treatment.⁶ The deformation mechanisms of fcc materials are of general interest in materials science, with many technologically important materials such as high entropy alloys,⁷ perovskites,⁸ salts,⁹ and oxides¹⁰ all sharing an fcc structure at ambient and high pressures.

The high pressure behavior of CeO₂ has been shown to be affected by temperature, sample environment hydrostaticity, and particle size and morphology. At ambient pressure, CeO₂ crystallizes in a face centered cubic fluorite (Fm3m) structure. Under high pressures, CeO₂ undergoes a phase transition to an orthorhombic structure.^11–16 The structure of this orthorhombic phase has been reported to be either Pnam^11,12,17 or Pbnm.¹⁸ The pressure range at which this phase transition occurs is known to be affected by temperature,¹² pressure transmitting medium used,¹¹ and the particle size.^13–16 Increasing the temperature has been shown to decrease the pressure at which the phase transformation begins.¹² 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.¹¹ 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,¹² and 31–55+ GPa in a 4:1 methanol:ethanol solution.¹⁹ Grain size of CeO₂ also affects the phase transformation pressure with the limit of a single crystal sample reducing the coexistence region of the fluorite and orthorhombic phases.¹³ Compressing CeO₂ 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 CeO₂, two studies have reported the relationship between pressure and yield strength of CeO₂. The first of these studies determined the strength from the linewidth of powder x-ray diffraction (XRD) peaks as a function of pressure.¹⁸ 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.¹⁸ reported that the yield strength of CeO₂ 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.²⁰ determined the pressure dependence of the strength of CeO₂ using gas-gun shock compression up to 25 GPa. A similar plateau in the strength of CeO₂ at pressures above 8 GPa was observed in this study.

In the present work, the relationship between yield strength of CeO₂ 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 CeO₂ 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.²¹ 

A BX-90-style DAC,²² 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 CeO₂ sample to high pressures. The CeO₂ 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 system²³ at the HPCAT (High Pressure Collaborative Access Team, Sector 16) ID-B beamline at the Advanced Photon Source at Argonne National Laboratory.²⁴ 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.²⁵ The wavelength of the x-ray beam used was 0.4066 Å, the spot size was 3 × 6 μm², 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.

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., $σ1,σ2andσ3$⁠). By convention, the stress along the compression axis is represented by $σ3$ and, by symmetry, the stresses in the two directions perpendicular to the direction of compression are equal (i.e., $σ1=σ2$⁠). The stress a material experiences during compression can be separated into hydrostatic (⁠$σp$⁠) and deviatoric stress components [the last term on the right hand side of Eq. (1)],²⁶ 

Here, $σp$ is the hydrostatic stress component and $t=σ3−σ1$ is the differential stress component.²⁶ 

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, d_m(hkl) is the measured lattice spacing, d_p(hkl) is the predicted lattice spacing under hydrostatic compression, $ψ$ is the azimuthal angle around the diffraction pattern in a counterclockwise sense with $ψ=0°$ denoting the direction of compression.^27,28Q(hkl) is the lattice strain parameter that is defined by

In Eq. (3), $GR$ is the shear modulus in the Reuss limit where stress continuity is assumed across the grains of the sample, and $GV$ 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 $GR$ and $GV$⁠, implying a value of $α=1/2$⁠.²⁶ This allows for t to be determined by

where G is the shear modulus.²⁶ 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.²⁹ 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.³⁰ The shear modulus was determined from the shear sound speed, V_S, 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 V_s = 3390 ms⁻¹, dV_s/dP = 32.17 m s⁻¹ GPa⁻¹ in the low pressure regime and V_s = 3304 m s⁻¹, dV_s/dP = 57.95 m s⁻¹ GPa⁻¹ 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 V_s 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)³¹ software package was used to perform a full-profile Rietveld refinement on the integrated data as described by Wenk et al.³² From the refinement, the hydrostatic lattice parameters of CeO₂ and Pt, Q(hkl) of CeO₂, 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.³³ 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).

To interpret the texture development of the sample and determine the activity of different slip systems, elasto-viscoplastic self-consistent (EVPSC) modelling was performed.³⁴ 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),³⁵ 

In Eq. (7), the strain rate tensor of the plastic component, “p,” is given by $ε˙ijp$⁠, $τs$ is the threshold value of shear stress along a slip system “s” required to induce slip at the reference shear rate, $γ0˙$⁠. The RSS is given by the product of two tensors: the symmetric Schmid factor $mkls$⁠, and the local stress tensor for an individual inclusion, $σkl$⁠. 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 $n=10$ 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 $n=10$ 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, $τs$⁠) is a determined value that is affected by both strain hardening and pressure³⁶ and is modeled as

Here, $τ0s$ is the initial CRSS at 2.7 GPa and $dτdP$ and $d2τd2P$ 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.³⁴ and Lin et al.³⁶ 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 BEARTEX³⁷ 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.³⁸ 

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⁻⁴ and 1 × 10⁻⁴, 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⁻³ and standard deviation of 4.41 × 10⁻⁴. 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 CeO₂ 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 CeO₂ 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.³⁸ 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 [001] and [011] 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 (C_ij) were obtained by using G [determined by Eq. (5)] and the bulk modulus (B) using the equations

B was determined by

where v_l is the longitudinal sound speed. From the ultrasound measurements, v_l = 6357 m s⁻¹, dv_l/dP = 76.2 m s⁻¹ GPa⁻¹. The relationships between B, G, and C_ij with pressure are provided in the supplementary material (Fig. S2 in the supplementary material). The relationships between B, G, and C_ij 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 C_ij 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.

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 [011] 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 ${111}⟨11¯0⟩$⁠, ${110}⟨11¯0⟩$⁠, and ${100}⟨011⟩$ 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 ${110}⟨11¯0⟩$ slip system. As pressure increases, the relative activity of the slip system ${110}⟨11¯0⟩$ declines to ∼25% at 6 GPa and ∼15% at 39 GPa (strain = 0.202). Over the same pressure range, the activity of the ${111}⟨11¯0⟩$ 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 ${100}⟨011⟩$ occurs between 5 and 39 GPa. Crossover of the relative activity of the ${111}⟨11¯0⟩$ and ${110}⟨11¯0⟩$ 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.¹⁸ is shown in blue diamonds in Fig. 2(b). Similar to the data obtained in this study, Liu et al.¹⁸ 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 CeO₂ 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.¹⁸ used values of the elastic constants derived from density functional theory. In this case, the shear modulus and its pressure derivative were G₀ = C₄₄ = 68.10 GPa and dC­₄₄/dP = 1.13.³⁹ 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 C₄₄ = G and dC₄₄/dP in the literature, a second order polynomial fit to our values of G gives G₀ = C₄₄ = 79.0 GPa, dC₄₄/dP = 2.9, and d²C₄₄/dP² = 0.0398 GPa⁻¹. 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 CeO₂ 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 CeO₂ obtained by Frost et al.¹⁷ using Brillouin spectroscopy was fit with a linear equation to give the values of G₀ = 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.,¹⁸ 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.²⁰ measured the strength of CeO₂ 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.⁴⁰ 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${111}⟨11¯0⟩$ slip system is the dominant deformation mechanism in CeO₂ up to 37 GPa at 300 K with secondary contributions from the ${110}⟨11¯0⟩$ and ${100}⟨011⟩$ slip systems. Other face centered cubic materials exhibit similar behavior such as Si where slip along ${111}⟨11¯0⟩$ was found to be the dominant deformation mechanism.^41,42

The differential stress (t) of CeO₂ 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 CeO₂ with pressure were observed, where the [011] direction increasingly aligned with the axis of compression, and the [001] direction progressively aligned perpendicular to the axis of compression. This change in texture is consistent with the activation of the ${111}⟨11¯0⟩$ 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 CeO₂ shear modulus data. We would like to thank Bethany Chidester for helpful technical discussions and the two anonymous reviewers whose feedback improved the manuscript.

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.