Impact-sensitivity predictions based on the vibrational up-pumping model show a strong polymorph dependency for RDX and highlight that one of the high-pressure forms, which forms during shock-wave experiments, is appreciably more susceptible to mechanical initiation. The origin of the predicted impact sensitivity variation can be attributed to vibrational mode hardening by pressure and to differences in the molecular conformation of RDX in the four polymorphs studied. These polymorphs present different distributions of molecular vibrations within their respective up-pumping windows, which leads to their varying ability to up-pump and trap the vibrational energy that arises from mechanical insult.
INTRODUCTION
The ease of initiating energetic materials (EMs, explosives, propelants and gas generators) by impact and shock are essential parameters from both an application and safety standpoint.1 Impact sensitivity (IS) is typically measured by a drop-hammer apparatus, where a known mass is dropped from a variable height until the minimum energy threshold required to induce initiation is established. However, despite the existence of well-established testing protocols, it is not uncommon for variable results to be reported by different laboratories, with temperature, humidity, sample purity, crystallinity, particle size, and operator experience all known to affect the outcome of the binary “go/no-go” call.2 In many EM research labs, RDX (formally hexahydro-1,3,5-trinitro-1,3,5-triazine, see Fig. 1) has informally been adopted as a laboratory internal standard for IS measurements, allowing the sensitivities of novel materials to be ranked in comparison to this widely characterized EM. However, there is substantial variability in the reported IS values of RDX itself.2
Structures of RDX showing differences in (a)–(c) molecular conformations and (d)–(g) crystal packing arrangements. (a) AAA, (b) AAE, and (c) AAI conformations, denoting A-axial, E-equatorial and I-intermediate positions of the nitro groups. (d) α-, (e) β-, (f) ε-, and (g) γ-polymorphic forms.
Structures of RDX showing differences in (a)–(c) molecular conformations and (d)–(g) crystal packing arrangements. (a) AAA, (b) AAE, and (c) AAI conformations, denoting A-axial, E-equatorial and I-intermediate positions of the nitro groups. (d) α-, (e) β-, (f) ε-, and (g) γ-polymorphic forms.
The issues associated with IS measurements have prompted a response from the modeling community to generate reliable structure/prediction models. This has included investigating correlations between IS and bond dissociation energies,3 electron density topologies,4,5 crystal void space and compressibility,6 and electronic band gaps.7,8 However, all of these methods tend to offer a more qualitative rationale for impact sensitivity and are generally restricted to chemically similar energetic molecules. In contrast, vibrational up-pumping has emerged as a reliable tool capable of successfully ranking a broad range of EMs according to their experimental IS values.8–18
In essence, the up-pumping model describes the ease with which the vibrational modes of a crystalline material [represented by the phonon density of states, g(ω)] can channel the energy from a mechanical impact event through its low energy external lattice modes to reach the localized molecular vibrations. A general schematic of the vibrational up-pumping process is shown in Fig. 2. The initial mechanical impact induces a compressive wave in the crystalline solid, depositing energy in the three acoustic (translational) motions of the crystal. From here the energy quickly equilibrates across the continuum (or “bath”) of external phonon modes, q, with an upper bound Ωmax that typically falls around 200 ± 50 cm−1.11 The vibrationally excited phonon bath modes subsequently combine and conduit (up-pump) excess energy into high frequency internal (localized) molecular vibrations, Q.
The up-pumping process defines permitted phonon–phonon scattering pathways, in accordance with a general vibrational Hamiltonian as follows:
where low-energy external lattice modes (Hq) and the higher-energy internal molecular modes (HQ) are treated as unique, separable entities. The third term of Eq. (1) refers to the conversion of energy between external and internal modes and can be expanded to first order by a Fermi golden rule formalism as follows:
where is the cubic anharmonic coupling constant that describes the strength of interaction between three vibrational modes and ensures energy conservation. While the values of are system specific, high-resolution Raman spectroscopy and computational modeling work have found that the magnitude of is largely conserved for organic molecular EMs, including for HMX, RDX, TATB, and PETN.19–21 The explicit calculation of this anharmonic term is, therefore, omitted in this work. The remainder of Eq. (2) is the two-phonon density of states, Ω(2), which describes the scattering (combination) of two phonons (q1, q2) to create a third, higher energy phonon (Q). The permitted pathways include the scattering of two identical phonon bath modes (akin to the generation of an overtone state in vibrational spectroscopy), two non-identical phonon bath modes (akin to the creation of a combination mode vibrational state), or one phonon bath mode combining with a doorway vibrational mode, QD, that resides between 1 and 2 Ωmax. These QD modes are typically low-lying angle bend/torsional motions, often associated with e.g., C–NO2 functional groups. The upper limit accessible by two phonon scattering events is, therefore, 3Ωmax, and this defines the upper boundary of the vibrational up-pumping window used in this work. The molecular vibrations contained within the window induce distortions of the weakest bonds in the energetic molecule. Thus, the up-pumping process localizes the initial mechanical energy through the external vibrations into the local modes, resulting in the activation of trigger linkages, leading to bond-breaking and initiation.22–24 In our work, we take the projection of Ω(2) onto g(ω) over 1–3 Ωmax as a metric to describe how efficiently the crystal lattice of the EM can trap mechanical energy in its molecular vibrations.
With the up-pumping model having already demonstrated success for numerous well-known EMs,8–18 we now turn our attention to RDX. This has been the subject of a number of previous reports. McNesby applied an up-pumping model using Raman spectroscopy data that successfully ranked the measured impact sensitivities of a number of EMs, including RDX.14 Chauduri et al. applied an up-pumping model to study the energy transfer at a hot–spot interface in RDX.25 Chung et al. demonstrated that lattice phonons dominate the thermal energy transfer to modes with significant nitrogen–nitrogen (NN) character in RDX,26 and later quantified the three-phonon scattering rates and mode-to-mode scattering rates for energy transfer from phonons expressed across the Brillouin zone into the NN molecular vibrations.27 Our own interest relates to the polymorphic nature of RDX, and in this, regard this study was also motivated by our recent combined experimental and theoretical report that another well-known EM, FOX-7, was predicted to show a polymorphic dependency for IS.28 Note in this work experimental validation proved impossible to obtain, as the act of measuring the IS for γ-FOX-7 induced a phase transition to α-FOX-7. Chaudhuri et al. have also predicted using an overtone-based up-pumping model that IS varies for polymorphs of both CL-20 and HMX.18
The two crystalline phases (known as the α-29 and β-forms30) of RDX have been structurally characterized at ambient conditions. A third polymorph, the ε-form,31 can be produced under high-pressure (5.2 GPa) and high-temperature (450 K) conditions. This form can be recovered under ambient pressure and low temperature conditions (<230 K), where it has been structurally characterized. A fourth phase, the γ-phase, has been obtained by compressing a single crystal of α-RDX to 3.9 GPa at ambient temperature.32 A fifth polymorph, the δ-phase, was identified in a high-pressure Raman spectroscopy study, but its crystal structure has not yet been determined.33 Crystal-packing arrangements of the four structurally characterized polymorphs are shown in Fig. 1.
As well as exhibiting different crystal-packing arrangements, the polymorphs of RDX also show different molecular conformations. While the triazine ring always maintains a chair conformation, the three pendant nitro groups vary in their orientations (see Fig. 1). In α-RDX, all molecules have two nitro groups in the axial (A) position and one in the equatorial (E) position, to give the AAE conformation. For the β- and ε-forms, all molecules adopt the AAA conformation, whereas the asymmetric unit of γ-RDX contains one molecule in the AAA conformation and a second in the AAI conformation. The AAI conformation can be regarded as intermediate between AAE and AAA.
In this paper, we report on the application of the vibrational up-pumping model to predict the IS for four polymorphs of RDX. The model provides a convenient platform to discuss structure/property relationships, permitting the effects of molecular conformation, crystal packing, and hydrostatic compression on IS to be probed in some detail for this important EM. Finally, we show how the predicted IS values for the RDX polymorphs rank alongside other EMs previously reported by the vibrational up-pumping model.
COMPUTATIONAL MODELING
Input geometries for the four polymorphs of RDX were obtained from the Cambridge structural database34 (codes: CTMTNA11,35 CTMTNA04,30 CTMTNA07,31 and CTMTNA0232) and subjected to full optimization using CASTEP version 17.21 (PBE functional with dispersion correction TS for the α- and β-polymorphs and G06 for the ε- and γ-polymorphs, norm-conserving pseudopotentials coupled to a plane wave basis set expressed at 950 eV, Brillouin zone sampling of 0.05 Å−1 and an FFT grid set to 2.0).36–38 Input file preparation was aided using seek-path.39 Geometry optimization convergence criteria were as follows: atomic forces <5.0 × 10−3 eV/Å, change in energy per atom <2.0 × 10−6 eV/atom, unit-cell stress <5.0 × 10−3 GPa and maximum atomic displacement 5 × 10−4 Å. For the γ-form, an external pressure of 3.9 GPa was applied during optimization. In addition, to elucidate whether any differences in the calculated phonon spectra for the high-pressure phase could be attributed directly to the external pressure (i.e., mode hardening), optimizations were repeated for γ- and α-RDX in the absence/presence of a 3.9 GPa external pressure, respectively.
Following geometry optimization (unit-cell parameters quoted in the supplementary material; the resulting deviations in unit cell volumes are of the order 4%–7%, indicating that reliable optimized structures have been obtained),40 phonon calculations were carried out at the Brillouin zone gamma-point only using density functional perturbation theory (DFPT).41 The acoustic sum rule was applied analytically. All structures returned all positive vibrational frequencies with the exception of β-RDX, which initially showed five low-lying imaginary frequencies (at −18, −17.5, −15.6, −6.1, and −2.6 cm−1). Tightening the fine FFT grid (from 2.0 to 4.0) reduced this to one (at −9 cm−1). As this is an experimentally characterized polymorph, it is likely that the imaginary frequency has arisen due to further numerical instabilities in the geometry optimization process, rather than an indication of a metastable state. However, as the predicted IS metric increased by less than 3% upon removal of four of the five imaginary frequencies, it is unlikely that further pursuing the calculation quality, by further tightening the convergence criteria, would result in any appreciable differences in the predicted impact sensitivity value.
All resulting g(ω) plots were presented with a Gaussian smearing width of 5 cm−1. Assignment of Ωmax was made in each case by tracking the displacement of the center of mass (CoM) for all molecules in a unit cell for each eigenvector, backed up with mode visualization using Jmol.42 When the CoM displacements fell below 10% of the maximum value, the phonons were considered to be more localized in nature than delocalized, thus marking the transition from external lattice mode to molecular-based vibrations. These plots are available in the supplementary material. The shock temperature, Tshock, i.e., the superheated phonon quasi-temperature adopted to simulate the mechanical impact event, is obtained from the ratio of the bulk heat capacity to the phonon bath heat capacity, Ctot/Cph (also shown in the supplementary material). Based on previous work a Ctot/Cph ratio of 5.00 was set to 3278 K (calculated from the adiabatic compression of a model organic crystal).16
To successfully capture the IS of RDX, our up-pumping model has evolved to include an additional scattering pathway. Our previous model considered: (1) the scattering of two external phonons into the doorway region (i.e., 1–2 Ωmax), followed by (2) the scattering of an external phonon with a doorway mode into the region 2–3 Ωmax. Although successful for many systems, this model severely underestimated the predicted IS of RDX. We have resolved this issue by now allowing the scattering in the pathway (2) to also up-pump density onto the doorway region, without sacrificing the predictive power across the previously reported EM dataset. Further discussion is provided in the supplementary material.
RESULTS AND DISCUSSION
The calculated g(ω) for all optimized structures is shown in Fig. 3(a). The same computational approach has previously been shown to provide excellent agreement with the experimental inelastic neutron scattering (INS) spectrum of α-RDX.43 From g(ω), we obtained the phonon scattering (up-pumping) pathways, Ω(2), also shown in Fig. 3(a) using the parameters listed in Table I. To ensure that this process generates up-pumping values that can be compared across different crystal structures, it is important that the data are normalized. In this work, g(ω) was normalized to the number of modes in the phonon bath prior to generating Ω(2). This allows us to compare different phonon spectra on an absolute scale. In addition, the integral of the overlap between Ω(2) and g(ω) within the up-pumping window is divided by the number of molecules in the unit cell. The latter normalization step is to account for the localization of phonon energy after up-pumping. Further details are available in the supplementary material.
Simulated vibrational spectra for RDX polymorphs. (a) g(ω) (gray) of the RDX polymorphs with Ω(2) overlayed (red), vertical blue lines represent sequential multiples of Ωmax, with the up-pumping window defined by the 1–3 Ωmax limits. (b) Analogous data for α-RDX compressed at 3.9 GPa. (c) Partial g(ω) showing the contributions to the vibrational modes from conformers AAA and AAI in γ-RDX.
Simulated vibrational spectra for RDX polymorphs. (a) g(ω) (gray) of the RDX polymorphs with Ω(2) overlayed (red), vertical blue lines represent sequential multiples of Ωmax, with the up-pumping window defined by the 1–3 Ωmax limits. (b) Analogous data for α-RDX compressed at 3.9 GPa. (c) Partial g(ω) showing the contributions to the vibrational modes from conformers AAA and AAI in γ-RDX.
Parameters used to calculate Ω(2) for polymorphs of RDX. Ωmax denotes the top of the phonon bath, Z is the number of molecules in each unit cell, Y is the number of amalgamated vibrations per molecule in the phonon bath regions, and Tshock is the temperature adopted for the phonon bath modes. The up-pumped density is the projection of Ω(2) onto g(ω) per molecule in the up-pumping window 1–3 Ωmax.
Polymorph . | ΩMAX . | Z . | Z(6 + Y) . | Y . | Tshock/K . | Up-pumped density/×103 A.U. . |
---|---|---|---|---|---|---|
α-RDX | 164 | 8 | 96 | 6 | 3265 | 48.1 |
β-RDX | 164 | 8 | 96 | 6 | 3390 | 26.7 |
ϵ-RDX | 164 | 4 | 48 | 6 | 3395 | 15.9 |
γ-RDXa | 268 | 8 | 112 | 8 | 2820 | 176.2 |
α-RDXa | 257 | 8 | 112 | 8 | 2835 | 67.0 |
Polymorph . | ΩMAX . | Z . | Z(6 + Y) . | Y . | Tshock/K . | Up-pumped density/×103 A.U. . |
---|---|---|---|---|---|---|
α-RDX | 164 | 8 | 96 | 6 | 3265 | 48.1 |
β-RDX | 164 | 8 | 96 | 6 | 3390 | 26.7 |
ϵ-RDX | 164 | 4 | 48 | 6 | 3395 | 15.9 |
γ-RDXa | 268 | 8 | 112 | 8 | 2820 | 176.2 |
α-RDXa | 257 | 8 | 112 | 8 | 2835 | 67.0 |
Optimized with 3.9 GPa external pressure.
It is immediately apparent that the β- and ε-forms have very similar g(ω) [and hence Ω(2)], which arises due to the same molecular conformation (AAA) being present in their respective crystallographic unit cells. This also suggests that the intermolecular interactions contained in these two polymorphs are weak. This agrees with the lower lattice energy previously reported for β-RDX compared to α-RDX43 and the observation that ε-RDX readily transforms to α-RDX on heating above 230 K.31 α-RDX has a notably different up-pumping window, which reflects the different molecular conformation (AAE) in this polymorphic form. Overall γ-RDX is the biggest outlier, with a significantly higher Ωmax and a more densely populated up-pumping window that contains molecular vibrations for two molecular conformations (AAA and AAI).
These differences in g(ω) and Ω(2) are carried through to the predicted up-pumping intensities (Table I). Our calculations suggest that the IS values for the β- and ε-forms are very similar, with α-RDX being higher, and γ-RDX being even higher. These predictions, therefore, suggest that the molecular conformation of RDX is critically important: maintaining the AAA conformation and varying the crystal packing observed for the β- and ε-polymorphs has little effect while switching to the AAE conformation presents more vibrational states to boost and trap Ω(2). The same observation holds for γ-RDX, where the two molecular conformations (AAA and AAI) display even more vibrational modes in the up-pumping window [Fig. 3(a)].
The next question that arises from the up-pumping model is why the phonon bath for γ-RDX extends so high compared to the other three polymorphs. The animation of the vibrational modes leading up to Ωmax (∼100–164 cm−1 for α-, β- and ε−RDX; 100–268 cm−1 for γ-RDX) shows the expected amalgamated mode (Y) behavior, where molecular deformation modes combine with the lattice motions (the six translational/rotational degrees of freedom, per molecule). For the ambient-pressure phases, these eigenvectors (six per molecule) involve various combinations of –NO2 twisting motions, while the first clusters of vibrational modes that fall within the up-pumping windows describe various combinations of N–NO2 out-of-plane motions. For the high-pressure γ-phase, both of these molecular motions are contained in the phonon bath, resulting in a greater number of amalgamated modes (eight per molecule) for γ-RDX compared to the ambient-pressure phases (Table I). Moreover, the extent of lattice mode behavior is significantly more pronounced for γ-RDX, as demonstrated by tracking the changes in the CoM for each of the eigenvectors (see the supplementary material).
Thus, the question now becomes: why do the lattice and molecular vibrations separate less readily for γ-RDX? The most likely explanation is that this polymorph was subjected to external pressure during geometry optimization. It has been noted that pressures of up to 4 GPa can be expected to produce a blue shift in molecular-based vibrational frequencies on the order of 5–20 cm−1,43 which is confirmed in the observed changes in mode distributions for the polymorphs of RDX noted above. It is also well known that external phonon frequencies tend to harden substantially with pressure, as noted by the Debye temperature dependence on pressure.44 In an attempt to quantify directly the effect of the external pressure on predicted IS values within the up-pumping model, the geometry optimization for γ-RDX was repeated in the absence of the external pressure. Unfortunately, this caused the unit cell to expand by over 30%, resulting in an unrealistic representation of the high-pressure phase. As an alternative, α-RDX was re-optimized under 3.9 GPa external pressure conditions to observe how this changes g(ω), Ω(2), and the predicted up-pumped density [Fig. 3(b)]. The unit-cell volume of compressed α-RDX decreased to levels comparable with γ-RDX (see supplementary material). On analysis of the resulting g(ω), a considerable increase in the value of Ωmax is immediately apparent, placing Ωmax at 257 cm−1, in close alignment with γ-RDX (Table I). The amalgamated mode count for compressed α-RDX also matches that for γ-RDX (Table I), while the remaining external modes of vibration have experienced mode hardening of ∼5–10 cm−1. The corresponding predicted up-pumping metric derived from the projection of Ω(2) onto g(ω) shown in Fig. 3(b) has increased compared to the non-compressed form (Table I). Thus, the effect of applying external pressure to α-RDX is significant, and by extension will likely also play a similar role on the high-pressure γ-form.
While the external pressure undoubtedly plays a role in boosting the predicted IS for γ-RDX, other factors must be at work to account for the near three-fold difference still observed between compressed α-RDX and γ-RDX. To account for this we now turn to the doorway region, which is the first half of the up-pumping window (i.e., 1–2 Ωmax). This is a particularly important region because, in the up-pumping model, the doorway modes both contribute to and capture the up-pumped energy. Given that the γ-form contains two molecular conformations, partial g(ω) plots were constructed to show the contribution from both conformations [see Fig. 3(c)]. Crucially this analysis shows that peaks at 330–360 cm−1, which are attributed to the AAI conformation alone, help to populate the doorway region. These eigenvectors are best characterized (visually) as N–NO2 bond stretching and ring deformation modes. The analogous vibrations for the AAA conformation fall at 360–390 cm−1, where they occur alongside other AAI ring deformation modes. Similar eigenvectors are observed for the AAA conformer at 350–360 cm−1 in β-RDX and at 360–370 cm−1 in ε-RDX, while those associated with the AAE conformer in α-RDX appear at 340–350 cm−1. Thus, these modes fall outside the doorway regions for all phases except for the γ-form.
It is, therefore, apparent that two factors are responsible for the prediction of increased IS for the γ-form of RDX: (i) direct compression, which leads to shorter and stronger intermolecular interactions that influence the structure of the g(ω) bath region, and (ii) changes in molecular conformation that alter the distribution of molecular vibrations that fall in the up-pumping window.
As the predicted IS values for the different polymorphs of RDX are relative values, it is possible to rank them alongside those of other EMs that have been investigated using the same up-pumping model (see the supplementary material). This is presented in Fig. 4, from which a clear relationship is drawn – the more sensitive a material is to impact (i.e., the lower the mechanical stimulus needed to initiate the EM), the higher the calculated up-pumped density. Taking an experimental IS value for RDX of 13 J,2 and assuming this corresponds to α-RDX, our prediction sits close to the curve. Note that, as experimental values are not known for the other polymorphs, we have plotted all corresponding data points at the same value on the x-axis. For completeness we also include the predicted variability for the FOX-7 polymorphs according to our updated up-pumping model (see supplementary material); this represents an update from our earlier publications.16,28
Experimental IS values vs vibrational up-pumped densities for a range of EMs (see the supplementary material), alongside predictions for the polymorphs of RDX, shown in red. Filled symbols correspond to experimentally measured data points, unfilled to predicted values only.
Experimental IS values vs vibrational up-pumped densities for a range of EMs (see the supplementary material), alongside predictions for the polymorphs of RDX, shown in red. Filled symbols correspond to experimentally measured data points, unfilled to predicted values only.
The differences in predicted impact sensitivities shown in Fig. 4 are quite stark. β-RDX sits below α-RDX, but arguably still close to the curve, despite the presence of one imaginary frequency in its g(ω) which likely means its predicted IS is slightly underestimated. ε-RDX sits lower still, with a predicted sensitivity closer to α-FOX-7. γ-RDX sits well above it, approaching a value close to that of the highly sensitive ε-CL-20. This raises important questions considering the handling of RDX under shock-wave conditions. At 3–5.5 GPa shock loading, Patterson et al. reported in situ Raman spectra that corresponded to the α → γ phase transition.45 Similarly, molecular dynamics simulations have also demonstrated that the α → γ transformation can occur readily under shock loading.46 Our simulations, therefore, raise the intriguing possibility that if the α → γ phase transition could be suppressed, then the shock sensitivity of RDX could be improved. In principle, this might be achieved by doping α-RDX with an additive that increases the α → γ transition pressure, thereby suppressing the formation of the more sensitive γ-form. There are precedents that pressure-induced polymorphism can be tuned through crystal engineering strategies. For example, doping the well-known EM ammonium nitrate with a group 1 nitrate suppresses an undesirable temperature-induced phase transition that is responsible for the deterioration of its mechanical properties.47 Doping has permitted high-pressure polymorphs to be stabilized at ambient conditions,48 and phase transitions have been pushed to higher pressures.49
CONCLUSIONS
Through the application of the vibrational up-pumping model, this work predicts that the impact sensitivity of RDX should show a strong polymorphic dependency. The high-pressure γ-form, which is accessible under shock-loading conditions, is predicted to be more prone to mechanochemical initiation. This is due to mode hardening through the presence of external pressure, and to the presence of two molecular conformations in the crystallographic unit cell, AAA and AAI, which boosts the number of molecular vibrations in the up-pumping window to both contribute to and trap the two phonon density of states energy, Ω(2). This suggests that if the α → γ phase transition could be suppressed, the shock sensitivity of RDX could be reduced, thereby enhancing safety under operational conditions. The β- and ε-forms, which both comprise the AAA molecular conformation of RDX, are predicted to be less efficient at trapping Ω(2) compared to the AAE conformation present in a-RDX. Consequently, vibrational up-pumping predicts that β- and ε-RDX should be less prone to mechanochemical initiation than α-RDX. This work further highlights the power of the vibrational up-pumping method to give new insights into important performance and safety metrics for energetic materials.
SUPPLEMENTARY MATERIAL
See the supplementary material for further information on the vibrational up-pumping model, geometry optimization, and phonon processing data.
ACKNOWLEDGMENTS
We thank Dr. Ranko Vrcelj (Center for Defense Chemistry, Cranfield University, Defense Academy of the United Kingdom) and Dr. Malcolm J. Burns (Los Alamos National Laboratory) for helpful discussions. This material is based upon work supported by the Air Force Office of Scientific Research under Award No. FA8655-20-1-7000. We are grateful for the computational support from the United Kingdom Materials and Molecular Modeling Hub, which is partly funded by EPSRC (Grant Nos. EP/PO20194 and EP/TO22213), for which access was obtained via the UKCP consortium and funded by EPSRC Grant No. ref EP/PO22561/1. I.L.C. thanks the University of Edinburgh, School of Chemistry for the award of a Ph.D. scholarship.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Imogen L. Christopher: Investigation (equal); Writing – original draft (equal); Writing – review & editing (supporting). Colin R. Pulham: Conceptualization (supporting); Validation (equal); Writing – review & editing (supporting). Adam A. L. Michalchuk: Conceptualization (equal); Formal analysis (supporting); Writing – review & editing (equal). Carole A. Morrison: Conceptualization (lead); Formal analysis (equal); Funding acquisition (lead); Investigation (lead); Methodology (equal); Project administration (lead); Writing – original draft (lead); Writing – review & editing (equal).
DATA AVAILABILITY
The data that support the findings of this study are available within the article and its supplementary material.