Theoretical estimations of equilibrium isotope fractionations between aqueous and crown ether coordinated Gadolinium (III)

Abstract

Gadolinium (Gd), as one of the rare earth elements, plays a significant role in planetary science and studies of early solar system evolution, since its isotopic composition can record both fractionation processes within nebulae and preserve primordial nucleosynthetic heterogeneities inherited from presolar materials. Among the isotopes of gadolinium, ¹⁵⁵Gd and ¹⁵⁷Gd have exceptionally high neutron absorption cross-sections, making the separation of these isotopes highly valuable as excellent neutron absorbers. However, with increasing applications, it has become an emerging contaminant in aqueous environments. This study employs quantum chemical software packages Gaussian09 and DIRAC19 to systematically investigate the coordination behavior of four functionalized crown ether resins (abbreviated as PMADB15C5, PMADB18C6, PMADB21C7, and PMADCH18C6) with Gd³⁺ ions, and evaluates their contributions to Gd isotope fractionations caused by both mass-dependent and the nuclear volume effects in terms of ¹⁶⁰Gd/¹⁵⁵Gd and ¹⁶⁰Gd/¹⁵⁷Gd. Our results indicate that in aqueous solutions, functionalized crown ether resins preferentially enrich light Gd isotopes compared to the [Gd(H₂O)₉]³⁺. Among them, PMADB15C5 exhibits the greatest potential for Gd isotope fractionations and is identified as a more reliable adsorbent for gadolinium isotope purification. Using the Rayleigh fractionation model at standard conditions (298.15 K, 1 atm) and taking the ¹⁶⁰Gd/¹⁵⁵Gd pair as an example, the isotopic fractionations induced by the adsorption process can be limited to within 0.138‰ when PMADB15C5 is used as the adsorbent with yield ≥ 95%. This study elucidates, at the molecular level, the coordination mechanisms between functionalized crown ether resins and Gd³⁺ ions. The theoretical findings on isotopic fractionation provide important insights for designing and optimizing chemical procedures in Gd isotope analysis.

Graphical abstract

Keywords

gadolinium / functionalized crown ether / equilibrium isotope fractionation / density functional theory / nuclear volume effect

1 Introduction

As an important rare earth element, gadolinium and its isotopes are widely used in various fields, including the nuclear industry, medicine (Sears, 1984), and planetary science. However, the widespread use of gadolinium-based materials has led to their increasing accumulation in the environment, with medical wastewater and industrial discharges potentially rendering gadolinium an aquatic contaminant (Blaum et al., 2002; Knappe et al., 2005; Kulaksız and Bau, 2007). Nevertheless, gadolinium can still serve as an effective tracer in karst aquifer systems for investigating surface water–groundwater interactions and the impacts of human activities (Boester and Rüde, 2020). Against the backdrop of rapid development in global emerging industries, the demand for critical metals such as gadolinium continues to grow. However, their limited reserves and uneven geographical distribution pose significant risks to supply chains. Among gadolinium isotopes, ¹⁵⁵Gd and ¹⁵⁷Gd are particularly prominent in multiple applications due to their exceptionally high thermal neutron absorption cross-sections. For instance, ¹⁵⁵Gd and ¹⁵⁷Gd serve as neutron absorbers in nuclear reactors to regulate reaction rates (Moreno et al., 1996; Oettingen and Cetnar, 2014; Dumazert et al., 2018). Variations in the ¹⁵⁶Gd/¹⁵⁵Gd and ¹⁵⁸Gd/¹⁵⁷Gd ratios can reflect thermal neutron capture processes or interactions between cosmic rays and planetary materials (Hidaka et al., 1999, 2000). ¹⁵⁷Gd holds promise for applications in neutron capture therapy (NCT) (Bufalino et al., 2006). Furthermore, since ¹⁶⁰Gd is primarily produced via the r-process, while ¹⁵⁵Gd and ¹⁵⁷Gd receive contributions from both s- and r-processes (Bisterzo et al., 2014). Hu et al. (2021, 2025) proposed that Gd isotopes can record both fractionation processes within the nebula and preserve inherited nucleosynthetic heterogeneities from presolar sources, making them important tracers for early solar system evolution. Additionally, since rare earth elements commonly exist in the + 3 oxidation state and exhibit highly similar geochemical behaviors, Gd occupies a central position in the lanthanide series. Its unique electronic configuration renders the Gd³⁺ ion (4f ⁷) an ideal reference system for studying the chemical behavior of other rare earth elements. For these reasons, theoretical studies on gadolinium isotopes hold significant scientific importance.

As a typical lanthanide element, Gd exhibits pronounced relativistic effects and a large nuclear volume isotope effect (Bigeleisen, 1996; Schauble, 2007). The nuclear volume effect (also known as the nuclear field shift effect) arises from changes in nuclear size and shape due to differences in neutron numbers among isotopes. These changes alter the nuclear charge radius, leading to shifts in the electronic ground state energy and ultimately resulting in mass-independent fractionation (MIF). In previous studies, the classical formulae established by Urey (1947) and Bigeleisen and Mayer (1947) failed to adequately explain the anomalous fractionation observed in uranium and zinc isotopes (Fujii et al., 1989; Nishizawa et al., 1996). Research by Nishizawa et al. (1995) prompted Bigeleisen (1996) to revise the original formula by incorporating the nuclear volume effect (NVE) as one of the key correction terms. Subsequently, Schauble (2007) suggested that pronounced nuclear volume effects should be present in heavy elements such as mercury. Although earlier work proposed that nuclear volume effects might contribute to gadolinium isotope fractionation in cation-exchange columns (Ismail et al., 2000), theoretical investigations on this topic remain scarce. Even in the limited theoretical studies available (e.g., Schauble, 2023), gadolinium isotopes have not been a primary focus. In Boda et al. (2022) study, they employed Turbomole package and density functional theory (DFT) to perform geometric optimizations, frequency analyses, and single-point energy calculations for complexes of functionalized crown ether resins with Gd³⁺ ions at the B3LYP/SVP and B3LYP/TZVP levels. Their work systematically investigated the geometric structures and thermodynamic properties of these systems, and calculated isotopic fractionation factors using the classical Bigeleisen–Mayer (Urey, 1947; Bigeleisen and Mayer, 1947) equation. Additionally, they performed atom-in-molecule calculations at the B3LYP/TZP/ZORA level using the ADF quantum chemistry package, thereby incorporating relativistic effects to characterize the electronic structure of the Gd-O bonds. The results indicated that PMADB18C6 exhibited the highest fractionation factor, favoring the enrichment of lighter isotopes (¹⁵⁵Gd and ¹⁵⁷Gd) relative to [Gd(H₂O)₉]³⁺. However, these calculations considered only for mass-dependent fractionation (MDF) and did not extensively investigate the influence of nuclear volume effects.

Therefore, in the case of simultaneously considering the two key factors of relativistic effects and nuclear volume effects, this study builds upon the work of Boda et al. (2022) and employs the Gaussian09 (Frisch et al. 2009) and DIRAC19 (Gomes et al. 2019) programs to perform quantum chemical calculations, aiming to theoretically evaluate Gd isotope fractionation abilities of four functionalized crown ether resins during chromatographic separation. Particular attention is given to the role of the nuclear volume effect, and the isotope analysis is through the Rayleigh fractionation model. This approach systematically elucidates the fractionation mechanisms of gadolinium isotopes during separation processes from a theoretical perspective.

2 Method

2.1 Theories of equilibrium isotope fractionation

Isotope effects originate intrinsically from mass-dependent variations in molecular vibrational frequencies that subsequently influence the physicochemical properties of compounds — a quintessential quantum mechanical phenomenon. In this study, the mass-dependent fractionation factor is derived based on thermodynamic equilibrium fractionation theory (Urey, 1947; Bigeleisen and Mayer, 1947), where isotopic equilibrium fractionation can be treated as an isotope exchange reaction as

(1)\begin{array}{c}{}_{}^{\mathrm{L}}\mathrm{G}{\mathrm{d}}_{\mathrm{a}\mathrm{q}\mathrm{u}\mathrm{e}\mathrm{o}\mathrm{u}\mathrm{s}}^{3+}{+}^{\mathrm{H}}\mathrm{G}{\mathrm{d}}_{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}\mathrm{l}\mathrm{e}\mathrm{x}}^{3+}{=}^{\mathrm{H}}\mathrm{G}{\mathrm{d}}_{\mathrm{a}\mathrm{q}\mathrm{u}\mathrm{e}\mathrm{o}\mathrm{u}\mathrm{s}}^{3+}{+}^{\mathrm{L}}\mathrm{G}{\mathrm{d}}_{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}\mathrm{l}\mathrm{e}\mathrm{x}}^{3+},\end{array}

where H and L denote the heavy and light isotopes, respectively; aqueous represents the aqueous Gd³⁺ phase, and complex refers to the coordinated phase formed by Gd³⁺ and functionalized crown ether resin. The mass-dependent fractionation factor ({\alpha }_{\mathrm{M}\mathrm{D}}) is expressed as the ratio of reduced partition function ratios (RPFR) for the aqueous and complex phases as

(2)\begin{array}{c}{\alpha }_{\mathrm{M}\mathrm{D}}={\displaystyle \frac{\mathrm{R}\mathrm{P}\mathrm{F}{\mathrm{R}}_{\mathrm{a}\mathrm{q}\mathrm{u}\mathrm{e}\mathrm{o}\mathrm{u}\mathrm{s}}}{\mathrm{R}\mathrm{P}\mathrm{F}{\mathrm{R}}_{\mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}\mathrm{l}\mathrm{e}\mathrm{x}}}},\end{array}

(3)\begin{array}{c}\mathrm{R}\mathrm{P}\mathrm{F}\mathrm{R}={\displaystyle \prod _i^{3n-6}{\displaystyle \frac{u_i^{\mathrm{H}}}{u_i^{\mathrm{L}}}}\frac{e^{-u_i^{\mathrm{H}}/2}}{1-e^{-u_i^{\mathrm{H}}}}\frac{{1-e}^{-u_i^{\mathrm{L}}}}{e^{-u_i^{\mathrm{L}}/2}}.}\end{array}

The parameter u is defined as \frac{h\nu }{{\mathrm{k}}_{\mathrm{B}}T}, where h represents the Planck’s constant, ν denotes the harmonic vibrational frequency in Hz, k_B stands for the Boltzmann constant, and T indicates the absolute temperature in Kelvin. For a molecular system comprising n atoms, the total degrees of freedom amount to 3n, which includes three translational along with three/two rotational degrees of freedom for nonlinear/linear molecules. Accordingly, remaining (3n − 6) or (3n − 5) frequencies will be used for calculating the RPFR values.

As we mentioned before, the nuclear volume effect, also known as the nuclear field shift effect, should also be included for theoretical estimations of Gd isotope fractionations due to its prominent relativistic effects. Building upon the work of Almoukhalalati et al. (2016), NVE can be quantitatively characterized through the mass-independent fractionation term lnK_fs (e.g., Bigeleisen, 1996) expressed as

(4)\begin{array}{c}\mathrm{l}\mathrm{n}{K}_{\mathrm{f}\mathrm{s}}=-\frac{1}{{\mathrm{k}}_{\mathrm{B}}T}\frac{Z_{\mathrm{A}}{e}^2}{6{ϵ}_0}{\delta ⟨r_n^2⟩}_{{\mathrm{A}}^{\mathrm{H}}{\mathrm{A}}^{\mathrm{L}}}\left\{\stackrel{-}{\rho }{}_e{}^{\mathrm{A}}\left(\mathrm{A}\mathrm{X}\right)-\stackrel{-}{\rho }{}_e{}^{\mathrm{A}}\left(\mathrm{A}\mathrm{Y}\right)\right\},\end{array}

where Z_A represents the nuclear charge number, e denotes the elementary charge, {ϵ}_0 stands for the vacuum permittivity, and {\delta ⟨r_n^2⟩}_{{\mathrm{A}}^{\mathrm{H}}{\mathrm{A}}^{\mathrm{L}}} corresponds to the mean square difference of nuclear charge radii defined as {\delta ⟨r_n^2⟩}_{{\mathrm{A}}^{\mathrm{H}}{\mathrm{A}}^{\mathrm{L}}}={⟨r_n^2⟩}_{{\mathrm{A}}^{\mathrm{L}}}-{⟨r_n^2⟩}_{{\mathrm{A}}^{\mathrm{H}}} with units of fm², where A^H and A^L refer to the heavy and light isotopes, respectively. The term \stackrel{-}{\rho }{}_e{}^{\mathrm{A}} represents the effective electron density of atom A in molecule AX or AY, abbreviated as effective density with units of bohr^–3. In this context, AX and AY correspond to the complex and aqueous terms in Eq. (1), respectively.

By combining the mass-dependent fractionation factor {\alpha }_{\mathrm{M}\mathrm{D}} with the nuclear volume effect term {\mathit{K}}_{\mathrm{f}\mathrm{s}}, the total isotopic fractionation factors {\alpha }_{\mathrm{t}\mathrm{o}\mathrm{t}} can be obtained:

(5)\begin{array}{c}\mathrm{l}\mathrm{n}{\alpha }_{\mathrm{t}\mathrm{o}\mathrm{t}}=\mathrm{l}\mathrm{n}{\alpha }_{\mathrm{M}\mathrm{D}}+\mathrm{l}\mathrm{n}{K}_{\mathrm{f}\mathrm{s}}.\end{array}

2.2 Computational details

Crown ethers are widely employed for metal ion separation and purification due to their exceptional ability to form highly stable complexes with metal ions. As demonstrated by Boda et al. (2022), functionalized crown ether resins can serve as column chromatographic adsorbents for selectively extracting gadolinium isotopes from aqueous Gd³⁺ solutions while simultaneously determining the corresponding isotopic fractionation factors. During this process, the exchange reaction of Gd³⁺ ions between the complex phase and aqueous solution phase can be represented as

(6)\begin{array}{c}{\left[\mathrm{G}\mathrm{d}{\left({\mathrm{H}}_2\mathrm{O}\right)}_9\right]}^{3+}+\mathrm{C}\mathrm{E}={\left[\mathrm{C}\mathrm{E}-\mathrm{G}\mathrm{d}\right]}^{3+}+9{\mathrm{H}}_2\mathrm{O}.\end{array}

Herein, the crown ether ligands are denoted as CE for simplicity, e.g., PMADB18C6. PMA denotes the polymethyl acrylic acid resin, while DB18C6 represents di-benzo-18-crown-6. Other crown ether analogs, including di-benzo-15-crown-5 (DB15C5), di-benzo-21-crown-7 (DB21C7), and dicyclohexano-18-crown-6 (DCH18C6), will also be employed in such reaction.

Initial structures of water molecules, hydrated gadolinium ions ([Gd(H₂O)₉]³⁺), functionalized crown ether resin ligands, and gadolinium-crown ether complexes were constructed using GaussView 5.0 molecular visualization software (Dennington et al., 2009). All geometry optimizations and frequency analysis were performed using the Gaussian09 quantum chemistry package (Frisch et al. 2009). The initial configurations of crown ether ligands were obtained from the PubChem website (Kim et al. 2016, 2023), with the following entries: DB15C5 was taken from compound CID 373245, DB18C6 from CID 26541, DB21C7 from CID 585932, and DCH18C6 from CID 85955.

Considering the unique characteristics of gadolinium as a lanthanide element, the PBE0 hybrid functional (Perdew et al., 1996; Perdew et al., 1997) was employed in this study. Previous work of Chen et al. (2019) demonstrated the suitability of this functional for lanthanide complexes, who successfully investigated the electronic structures of Pr(III) complexes using PBE0 in combination with the Stuttgart pseudopotential. In the present DFT calculations based on the PBE0 functional, a mixed basis set strategy was adopted. To account for the pronounced relativistic effects of Gd, the MWB53 pseudopotential with relativistic corrections (Dolg et al., 1989a) was applied. For the non-metal atoms (C, H, O, N), the 6-31G(d) basis sets (Ditchfield et al., 1971; Hehre et al., 1972; Hariharan and Pople, 1973) were employed.

To accurately simulate molecular behavior in solution, geometry optimizations of the target system were performed under both gas phase (vacuum) and aqueous phase (solvent model) conditions. For the aqueous phase calculations, the choice of solvent model was carefully considered. Although the SMD (Solvation Model Based on Density) model, which incorporates both polar and nonpolar solvation contributions, is widely recognized as one of the most accurate implicit solvent models, its nonpolar component may occasionally introduce minor numerical perturbations. Such perturbations can result in spurious imaginary frequencies or convergence difficulties in highly symmetric structures. Since the solvation free energy in the SMD framework is partitioned into long-range electrostatic interactions (calculated via the IEF-PCM approach) and short-range nonpolar contributions, and considering that geometry optimization and vibrational frequency analysis are predominantly governed by electrostatic effects (Marenich et al., 2009), this study employed the IEF-PCM implicit solvation model during the geometry optimization and frequency analysis steps to ensure computational stability and consistency. Water was specified as the solvent, and the same level of theory as in the gas phase calculations was maintained. For the subsequent single-point energy calculations, the SMD solvation model was adopted. All optimized structures were verified through vibrational frequency analysis to ensure the absence of imaginary frequencies.

The thermodynamic parameters for Eq. (6) — including internal energy change (ΔU), enthalpy change (ΔH), and Gibbs free energy change (ΔG) — can be determined through theoretical calculations. These parameters are obtained via high-accurate single-point energy calculations coupled with moderate-level frequency analysis, computed using the following expressions:

(7)\begin{array}{c}\mathrm{\Delta }E=9E_{{\mathrm{H}}_2\mathrm{O}}+E_{\mathrm{M}\mathrm{L}}-\left(E_{\mathrm{M}}+E_{\mathrm{L}}\right),\end{array}

where E_ML, E_M, and E represent the energies of the complex [PMA-DB18C6-Gd]³⁺, the hydrated ion [Gd(H₂O)₉]³⁺, and the ligand PMADB18C6, respectively. The numerical values of these thermodynamic parameters serve to evaluate the reaction’s spontaneity and other essential properties.

To obtain more accurate single-point energy data, extra calculations for the gas phase systems were performed at the PBE0/def2-TZVPP (Dolg et al., 1989b; Weigend and Ahlrichs., 2005; Gulde et al., 2012) level. Since single-point energy calculations in the aqueous phase were performed using the SMD model, the M05-2X method (Zhao et al., 2006) was adopted with carefully selected basis sets to accurately evaluate solvation effects (Ho et al., 2010). For systems in the aqueous phase, single-point energy calculations were conducted at M05-2X level with the 6-31G(d) basis set for non-metal atom and the def2-TZVPP basis set for Gd³⁺ ions, ensuring the reliability and precision of the energy results. Furthermore, the Shermo program (Lu and Chen, 2021) was utilized to integrate the zero-point energy (ZPE) obtained from single-point energy calculations with thermal corrections derived from geometric optimization and vibrational frequency analysis. In the frequency analysis of chemical systems, theoretical calculations are typically based on the harmonic approximation, which neglects anharmonic effects as well as inherent errors in the description of the potential energy surface by the chosen method and basis set. These limitations may lead to significant deviations in thermodynamic properties such as fundamental frequencies, zero-point energies, entropies, and enthalpy changes due to temperature variations. Moreover, since the 6-31G(d) basis set was used for the non-metal atoms during geometry optimization and frequency calculations in this study (with only one Gd atom present systems), scaling factors of 0.9726 for zero-point energy and 0.9777 for enthalpy corrections (Merrick et al., 2007) were introduced to reduce these errors. All thermodynamic properties, including Gibbs free energy and enthalpy, were computed at 298.15 K and 1 atm, providing a robust data set for subsequent thermodynamic analysis.

To accurately evaluate the NVE, this study employed the DIRAC19 code (Gomes et al. 2019) to perform all-electron relativistic calculations based on the effective electron density method (Almoukhalalati et al., 2016). However, due to the prohibitively high computational demands for four-component all-electron relativistic calculations on medium-to-large systems (e.g., the [PMADB15C5-Gd]³⁺ complex in this study containing 64 atoms), obtaining direct all-electron relativistic results for functionalized crown ether resin complexes proved unfeasible. To accurately characterize the influence of relativistic effects on the electronic structure, an alternative computational strategy was adopted. First, the contact electron density (referred to as contact density) was calculated for six valence states (0 to + 5) of a single Gd atom using Gaussian09 and Multiwfn 3.8(dev) (Lu and Chen, 2012) at the PBE0/x2c-TZVPPall (Pollak and Weigend, 2017) level, incorporating the fourth-order Douglas–Kroll–Hess scalar relativistic correction with spin–orbit coupling (DKHSO, Nakajima and Hirao, 2012). In parallel, all-electron relativistic calculations were performed on the Gd atom using the DIRAC19 program with the uncontracted Dyall’s triple-zeta all-electron basis set dyall.ae3z (Dyall, 2012; Dyall, 2023) with the four-component Dirac–Coulomb Hamiltonian (Saue, 2011) to obtain the corresponding effective electron density (referred to as effective density). By systematically comparing the results from both methods, a linear correction relationship between the contact density and effective density was established. Second, although DIRAC currently offers the solvation model sole for four-component methods, its application to single atoms or ions lacks clear physicochemical meaning. Therefore, the linear relationship derived under gas phase conditions was extended to solvated environments by assuming the solvation effects on contact/effective densities are negligible. Finally, all-electron relativistic calculations under gas phase conditions were conducted in Gaussian09 for all target systems, including a water molecule, [Gd(H₂O)₉]³⁺, four functionalized crown ether resin ligands, and their Gd³⁺ complexes. The Gd atoms were described with the x2c-TZVPPall all-electron basis set, while other non-metal atoms were treated with def2-TZVPP. All calculations were performed at the PBE0 level with DKHSO relativistic corrections to obtain the contact densities. Using the previously fitted relationship, the scalar-relativistic contact densities from Gaussian09 were corrected into “effective electron densities” equivalent to those from four-component DIRAC calculations.

3 Results

3.1 Structures of functionalized crown ether resins and their Gd³⁺ complexes in gas and aqueous phases

This study obtained the optimized structures of a series of functionalized crown ether resins and their gadolinium complexes through quantum chemical calculations. The optimized structures of the hydrated gadolinium ion [Gd(H₂O)₉]³⁺, the crown ether monomer PMADB15C5 and their corresponding gadolinium ion complexes under both gas phase and aqueous phase are shown in Fig. 1. These structures serve as reference systems and provide a critical benchmark for subsequent stability analysis of the complexes. Additionally, Table 1 summarizes the average Gd-O bond lengths for [Gd(H₂O)₉]³⁺ and the four Gd³⁺–crown ether complexes in both gas and aqueous phases.

3.2 Thermodynamic parameters of Gd³⁺-Crown Ether Complexes

Based on the Eq. (7), Table 2 compiles the thermodynamic parameters for the interaction between Gd³⁺ ions and PMA resins in different phases (gas and aqueous), including: internal energy change (ΔU), enthalpy change (ΔH) and Gibbs free energy change (ΔG).

This study employs solvation free energy \mathrm{\Delta }{G}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{v}} to quantitatively assess the free energy differences between gas phase and solvated conditions. The solvation energy was calculated using the equation: \mathrm{\Delta }{G}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{v}}=E_{\mathrm{a}\mathrm{q}\mathrm{u}}-E_{\mathrm{g}\mathrm{a}\mathrm{s}}. As presented in Table 3, comprehensive calculations were performed to determine the solvation energies for all components, including functionalized crown ether resins, gadolinium complexes, and water molecules.

3.3 Equilibrium Gd isotope fractionations

To achieve theoretical results comparable to the experimental conditions reported by Boda et al. (2022), all subsequent isotope fractionation calculations and analyses in this study were performed using solvation models. The RPFRs and mass-dependent fractionation factor α_MD between Gd³⁺ hydrates and their complexes with functionalized crown ether resins, calculated using Eqs. (2) and (3), are presented in Table 4. The linear fitting results for contact and effective density data of monatomic Gd in the gas phase are shown in Fig. 2.

NVE theoretical calculations were performed using a relativistic scalar approach with all-electron Gaussian basis sets (Table 6). Using the mean-square nuclear charge radius differences {\delta ⟨r_n^2⟩}_{{\mathrm{A}}^{\mathrm{H}}{\mathrm{A}}^{\mathrm{L}}} reported by Angeli and Marinova (2013), we calculated the nuclear volume fractionation factors lnK_fs and total fractionation factors lnα_tot at 298.15 K (Table 7).

The computational results demonstrate that the coordination reaction between PMADB15C5 resin and [Gd(H₂O)₉]³⁺ exhibits the most significant total gadolinium isotope fractionation effect among the studied systems. This theoretical prediction, however, shows a notable discrepancy with the experimental findings reported by Boda et al. (2022), where PMADB18C6 displayed optimal fractionation performance in practical separation applications. To validate their computational results, Boda's research team experimentally determined the fractionation factors of PMADB18C6 for ¹⁵⁵Gd/¹⁶⁰Gd and ¹⁵⁷Gd/¹⁶⁰Gd using column chromatography separation, obtaining values of α(160/155) = 1.00318 and α(160/157) = 1.00187, which are significantly higher than the theoretically predicted values of α_tot(160/155) = 1.000733 and α_tot(160/157) = 1.000423 obtained from our quantum chemical calculations.

4 Discussion

4.1 Structural variations and solvation effects in Gd³⁺-crown ether complexes

To comprehensively evaluate the adsorption performance of various functionalized crown ether resins, it is essential to examine both their thermodynamic properties and structural characteristics. Initially, the strength of adsorption can be preliminarily assessed based on the standard Gibbs free energy change ΔG of the adsorption reaction and the average bond length between the crown ether and Gd³⁺ ions. When the free energy values among different resins are too close (in this study, the free energy change under the solvation model ranges from –64.98 kcal/mol to –61.27 kcal/mol) to allow clear differentiation of their adsorption performance, the resin system with the smallest fractionation effect should be selected for applications in elemental purification, as it imposes lower demands on yields. In scenarios aimed at isotope separation, resins with larger fractionation factors should be prioritized, since a stronger fractionation capacity is more beneficial for enhancing separation efficiency.

The variation in Gd-O bond lengths directly modulates the strength of coordination interactions between Gd³⁺ ions and adsorbents, consequently regulating the thermodynamic properties of the adsorption process. From coordination chemistry perspective, shorter Gd-O bond lengths typically indicate stronger coordination bonds, primarily resulting from enhanced electrostatic interactions between the central metal ion and ligand atoms. Our results show that, with the exception of the [PMADB21C7-Gd]³⁺ complex, all systems exhibit relatively uniform Gd-O bond length distributions with their theoretical coordination numbers (Table 1). Notably, the maximum Gd-O bond length in solvated [PMADB21C7-Gd]³⁺ reaches 3.46 Å—approximately 1 Å longer than other bonds—suggesting this oxygen atom may not form an effective coordination bond with Gd³⁺. This outlier implies potential structural instability in the complex. The selective adsorption capacity of crown ether ligands originates from their unique polar cavity structures. The strong polarity of ether groups and lone electron pairs on oxygen atoms enable intense electrostatic attraction with Gd³⁺ ions, thermodynamically stabilizing the Gd³⁺ at the cavity center. This specific interaction exhibits three key features. 1) Smaller cavities (e.g., DB15C5) impose stronger spatial constraints, yielding shorter average Gd-O bond lengths. 2) Aqueous environments weaken Gd-O interactions, slightly increasing bond lengths. 3) Uniform bond lengths reflect stable coordination environments, while outliers indicate defects.

These structural features provide molecular-level insights into the adsorption performance differences among crown ether resins. Specifically, the anomalously long bond observed in [PMADB21C7-Gd]³⁺ may directly correlate with its relatively low isotopic fractionation efficiency, as the defective coordination geometry likely reduces both thermodynamic stability and isotopic selectivity. The findings highlight how precise control of cavity sizes and coordination environments can optimize Gd³⁺ separation performance.

Gibbs free energy (ΔG), as a thermodynamic function integrating both enthalpy (ΔH) and entropy (ΔS) changes, serves as the crucial indicator for determining the spontaneity of coordination reactions. Our results reveal that except for the gas phase PMADB15C5 system which shows a positive ΔG value, all other systems exhibit negative ΔG values (Table 2), indicating that most Gd³⁺-crown ether coordination reactions are thermodynamically spontaneous at gas phase conditions.

Analysis of solvation effects reveals substantial differences in free energy changes across the studied systems. As shown in Table 3, the solvation free energy ΔG_solv of the [PMADB15C5-Gd]³⁺ complex reaches –451.20 kcal/mol, while the corresponding values for PMADB15C5 and [Gd(H₂O)₉]³⁺ are –26.68 kcal/mol and –406.39 kcal/mol, respectively. This order-of-magnitude difference in solvation energy (|ΔG_solv| > 400 kcal/mol) results in fundamentally opposite trends in reaction spontaneity for the PMADB15C5 system between gas phase (ΔG = 0.32 kcal/mol) and solvated (ΔG = –62.11 kcal/mol) conditions.

Systematic evaluation of computational and physicochemical parameters reveals that the [PMADB15C5-Gd]³⁺ complex exhibits the shortest average Gd–O bond lengths (gas phase: 2.33 Å; aqueous phase: 2.50 Å), indicating the strongest metal-ligand interaction. Furthermore, in the aqueous phase, the free energy differences among the four functionalized crown ether resins are relatively small (ΔG(MAX–MIN) < 3.71 kcal/mol), suggesting that the solvent environment diminishes the thermodynamic distinctions between the systems.

4.2 Gd³⁺-crown ether isotope fractionations

Figure 3 presents the temperature-dependent variations of RPFR and mass-dependent fractionation factors based on data from Table 4 and Table 5. The results reveal significant differences in RPFR values among various systems for both ¹⁶⁰Gd/¹⁵⁵Gd and ¹⁶⁰Gd/¹⁵⁷Gd isotope pairs, following a consistent order of [Gd(H₂O)₉]³⁺ > [PMADB21C7-Gd]³⁺ > [PMADB18C6-Gd]³⁺ > [PMADCH18C6-Gd]³⁺ > [PMADB15C5-Gd]³⁺. This trend demonstrates that hydrated Gd³⁺ ions exhibit higher RPFR values compared to those complexed with functionalized PMA resins, indicating the selective complexation of lighter Gd isotopes by the polar cavities of crown ethers, while the more flexible hydration sphere preferentially accommodates heavier isotopes. A notable observation is the progressive increase in RPFR values with expanding crown ether ring size (from DB15C5 to DB21C7) and increasing number of ether oxygen atoms. This phenomenon likely relates to altered coordination environments resulting from enlarged cavity sizes, where expanded coordination spaces appear more favorable for the enrichment of heavier Gd isotopes. Among all functionalized crown ether resins, the PMADB15C5 complex displays the lowest RPFR values, a characteristic that correlates with its largest mass-dependent fractionation factors. Particularly intriguing is the observation that while the PMADB15C5 complex possesses the lowest complexation free energy (Table 2), it exhibits the largest fractionation factors. This apparent inverse correlation may originate from the smaller polar cavity structure of PMADB15C5, which demonstrates stronger selective complexation capability for lighter Gd isotopes (¹⁵⁵Gd and ¹⁵⁷Gd), consequently leading to significantly reduced RPFR values. The restricted cavity geometry appears to enhance mass-dependent discrimination effects, overriding the thermodynamic stability considerations reflected in the complexation free energies.

Based on all-electron relativistic calculations of Gd atoms in the gas phase using both Gaussian and DIRAC programs, we established a linear correction relationship between the computational results from these two methods. The data reveal that the contact density obtained from Gaussian calculations is approximately 2.8 times greater than that from DIRAC results (Table 5). This systematic discrepancy unequivocally demonstrates the inherent limitations of the scalar relativistic approximation in Gaussian when treating the electronic structure of lanthanide systems, particularly for f-electron containing Gd³⁺ complexes. Notably, while the solvent model implemented in Gaussian simulations introduces additional potential terms to mimic solution environments, the relative contact electron density deviation between solvated and gas phases (RD(Sol/Gas)) consistently remains around 0. This observation highlights a marked deviation from the effective density derived from DIRAC calculations, clearly demonstrating that the magnitude of relativistic corrections far exceeds the electron density variations induced by solvent effects. The systematic differences between Gaussian and DIRAC results primarily stem from inadequate relativistic corrections, whereas solvent-induced electron density variations (\stackrel{-}{\mathrm{\Delta }\rho }{}_e{}^{\mathrm{A}}\left(\mathrm{C}\mathrm{o}\mathrm{m}\mathrm{p}-\mathrm{A}\mathrm{q}\mathrm{u}\right)_max = 6.7) account for merely one-millionth of the relativistic correction magnitude. This orders-of-magnitude disparity unequivocally demonstrates that accurate treatment of relativistic effects constitutes the paramount consideration in lanthanide complex electronic structure calculations, while solvation effects become negligible when proper relativistic corrections are applied. These findings provide a theoretical foundation for simplifying aqueous phase computational protocols: Gas-phase-derived correction relationships remain valid for solvated systems when coupled with rigorous relativistic treatments.

Analysis of Table 6 data reveals that under solvation model conditions, the absolute differences of effective electron density across the systems range from 5.1 to 6.7 Bohr⁻³. This range significantly exceeds twice the standard deviation (SD = 1.8) of the linear fitting curve, confirming the statistical significance of these density differences. However, the variation range of density differences among different complexes is only 1.6 (5.1−6.7). When considering the fitting error (± 1.8), the electron density differences (compared to [Gd(H₂O)₉]³⁺) among [PMADB15C5-Gd]³⁺, [PMADB18C6-Gd]³⁺, [PMADB21C7-Gd]³⁺, [PMADCH18C6-Gd]³⁺ do not show statistical significance. This phenomenon likely stems from the shared Gd³⁺ central ion configuration (constant oxidation state) in all systems, resulting in relatively limited NVE-induced fractionation differences across isotope pairs (e.g., ¹⁶⁰Gd/¹⁵⁵Gd) for different reactions. Notably, nuclear volume effects and mass-dependent fractionation exhibit opposing directions in this case: while mass-dependent fractionation preferentially enriches heavier isotopes, nuclear volume effects promote the enrichment of lighter isotopes. An important observation is that although the [PMADB15C5-Gd]³⁺ complex displays the largest nuclear volume fractionation lnK_fs, it simultaneously possesses the smallest mass-dependent fractionation factors α_MD, ultimately maintaining the highest total fractionation factors α_tot. Further analysis reveals significant variation in the relative contribution of nuclear volume effects (lnK_fs/lnα_MD) for the ¹⁶⁰Gd/¹⁵⁵Gd isotope pair, ranging from 13.82% for [PMADCH18C6-Gd]³⁺ (lowest) to 37.28% for [PMADB21C7-Gd]³⁺ (highest). These results confirm the non-negligible role of nuclear volume effects in Gd isotope separation processes. From the perspective of absolute fractionation magnitudes, the differences in mass-dependent fractionation among functionalized crown ether resin complexes (with a maximum difference of 0.722‰ in 1000lnα_MD) are substantially greater than those of nuclear volume fractionation (showing only a 0.045‰ maximum difference in 1000lnK_fs). This discrepancy primarily originates from coordination environment characteristics: the [Gd(H₂O)₉]³⁺ complex with its larger polar cavity more readily accommodates heavier isotopes, while the smaller cavity of [PMADB15C5-Gd]³⁺ demonstrates stronger selectivity for lighter isotopes. Based on these findings, we recommend prioritizing functionalized crown ether resin complexes with smaller cavity sizes (such as [PMADB15C5-Gd]³⁺) for optimal isotope separation performance. This selection strategy effectively enhances selective enrichment of lighter isotopes through synergistic optimization of both mass-dependent and nuclear volume effects.

As shown in Fig. 3(c) and 3(d), the fractionation factors of the [PMADB18C6-Gd]³⁺ system measured experimentally by the Boda et al. (2022) (red pentagrams; α_tot (160/155) = 1.003180; α_tot (160/157) = 1.001870) differ substantially from the theoretical values obtained in this study (α_tot (160/155) = 1.000733; α_tot (160/157) = 1.000423). The observed differences between theoretical predictions and experimental measurements likely stem from several fundamental factors. First, the theoretical model necessarily incorporates simplifications that may not fully account for dynamic processes in actual separation systems, particularly kinetic effects during the chromatographic process. Second, while the calculations include solvent effects through implicit solvation models, they may not capture all aspects of the complex solvent-resin interactions occurring in real experimental conditions. Third, the interfacial phenomena between the resin and solution phase, which can significantly influence isotope fractionation, are challenging to model accurately.

Despite these quantitative differences, it is important to emphasize that the theoretical calculations successfully reproduced the general fractionation trends observed for Gd isotopes across different crown ether resins. This agreement validates the computational approach as a valuable tool for screening and designing novel separation materials. The systematic comparison between computational predictions and experimental results presented in this study provides crucial insights into the fundamental mechanisms governing lanthanide isotope fractionation at molecular level. These findings establish an important foundation for future studies aiming to optimize separation materials and processes for rare earth elements, while highlighting areas where theoretical models require further refinement to better match experimental observations.

4.3 Rayleigh fractionation model-based control of Gd purification

The measurement of isotope ratios in isotopic analysis is influenced by multiple factors, including the representativeness of sample collection and contamination prevention measures, the accuracy of sample pretreatment (e.g., effectiveness of dissolution and separation, isotopic fractionation), the selection of appropriate analytical methods and standard material calibration, instrument performance (sensitivity, precision, and stability), data processing corrections and error analysis, laboratory environmental conditions, the technical proficiency of operators, as well as sample storage and transportation conditions. Controlling these factors is critical to ensuring the accuracy and reliability of analytical results. Typically, the yield is calculated to assess the loss of the target element during dissolution and separation-purification processes, thereby avoiding significant isotopic fractionation that could affect the analytical results. The yield is defined as

f=\frac{\mathrm{R}\mathrm{e}\mathrm{c}\mathrm{o}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{d}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{m}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{t}\mathrm{y}}{\mathrm{I}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{l}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{m}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{t}\mathrm{y}}\times 100\mathrm{\%}.

In gadolinium isotope separation experiments, a yield ≥ 95% was established as the acceptable standard range. This study defines the isotopic fractionation factor during the adsorption process as K_d (the equilibrium distribution coefficient), with numerical values provided in Table 4.

To investigate isotopic fractionation behavior during the elution process, the following experimental conditions were established. After achieving 100% adsorption efficiency, the complex phase was subjected to elution (desorption). By comparing the isotopic composition differences between the eluate and the initial solution, the elution fractionation factor (denoted as α_elution) was determined. A Rayleigh fractionation model was employed to describe the isotopic fractionation behavior during the elution process:

(8)\begin{array}{c}{R}_{\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{n}}=R_0\cdot {\left(1-f\right)}^{\left(K_{\mathrm{d}}-1\right)}，\end{array}

where R_{\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{n}} represents the isotope ratio (heavy/light) in the complex phase at a given time during the process; R_0 denotes the initial isotope ratio; 1-f indicates the fraction of remaining material (ratio of unseparated mass to initial mass), with \mathrm{f} being the yield; and K_d is the equilibrium coefficient for the Rayleigh fractionation process, defined as K_{\mathrm{d}}={\displaystyle \frac{R_{\mathrm{a}\mathrm{q}\mathrm{u}\mathrm{e}\mathrm{o}\mathrm{u}\mathrm{s}}}{R_{\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{n}}}}, with values provided in Table 4.

According to the principle of mass balance, the fractionation factor for the elution process α_elution is expressed by the following equation:

(9)\begin{array}{c}1000\mathrm{l}\mathrm{n}{\alpha }_{\mathrm{e}\mathrm{l}\mathrm{u}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}=1000\mathrm{l}\mathrm{n}{\displaystyle \frac{1-{\left(1-f\right)}^{K_{\mathrm{d}}}}{f}}‰\end{array}

The elution fractionation factor α_elution was determined by monitoring changes in the Gd isotopic composition of the residual solution. Simulation data demonstrate that at 298.15 K, the fractionation factor α_elution for both ¹⁶⁰Gd/¹⁵⁵Gd and ¹⁶⁰Gd/¹⁵⁷Gd exhibit systematic variations across the 90%–100% yield range, as illustrated in Figs. 4(a) and 4(b), respectively.

The results demonstrate that when employing PMADB15C5-functionalized crown ether resin as the adsorbent under experimental conditions ensuring ≥ 95% yield, the isotopic fractionation effects during the adsorption process can be maintained at minimal levels: < 0.138‰ for the ¹⁶⁰Gd/¹⁵⁵Gd fractionation factor and < 0.080‰ for the ¹⁶⁰Gd/¹⁵⁷Gd fractionation factor.

4.4 Cosmo-chemical and geochemical applications

Di et al. (2024) proposed a theoretical framework that systematically incorporates natural mass-dependent fractionation into radiogenic isotope chronometry. They demonstrated that when parent–daughter elements undergo MDF in natural processes, conventional internal normalization (IN) yields systematically biased radiogenic ratios, thereby affecting isochron slopes and calculated ages. Using the Sm–Nd system in CAIs as an example, they showed how stable isotope fractionation (e.g., δ¹⁴⁶Nd, δ¹⁵²Sm) propagates into derived radiogenic ratios such as ¹⁴³Nd/¹⁴⁴Nd(IN), and emphasized that accurate chronological interpretation requires explicit measurement and parameterization of natural fractionation laws (e.g., the fractionation exponent γ). Although exemplified with Sm–Nd, the framework is broadly applicable to other REE systems.

Besides, Hu et al. (2021, 2025) provided key cosmochemical evidence underscoring the importance of REE, particularly Gd, in early solar system processes. Hu et al. (2021) revealed pronounced light-isotope enrich-ments in Gd, Dy and Er within Allende CAIs, interpreted as signatures of kinetic evaporation–condensation events in the protoplanetary disk. Furthermore, Hu et al. (2025) compared REE isotopic compositions across diverse chondritic components and demonstrated that not all REE display anomalies of equal magnitude, i.e., Gd, Dy, Er and Yb exhibit the strongest deviations. These anomalies were attributed not solely to solar nebular processes but to incomplete mixing of s- and r-process nucleosynthetic components, highlighting that Gd isotopes simultaneously record both astrophysical mixing and disk-level kinetic effects.

Within this context, our theoretical calculations on the adsorption behavior and equilibrium fractionation of ¹⁵⁵Gd/¹⁶⁰Gd and ¹⁵⁷Gd/¹⁶⁰Gd on functionalized crown ether resins establish a molecular-scale link between natural observations and theoretical frameworks. Specifically, 1) the computed differences in bonding energies and vibrational properties among Gd isotopes can be translated into equilibrium fractionation parameters (γ_Gd), providing physically grounded priors for the Di et al. (2024) model; 2) the magnitude and direction of predicted equilibrium shifts furnish baselines against which natural data can be evaluated, aiding in the discrimination of equilibrium versus kinetic fractionations as proposed by Hu et al. (2021); and 3) when viewed alongside Hu et al. (2025), the theoretical results can further constrain the extent to which observed Gd isotopic anomalies reflect nucleosynthetic heterogeneity rather than subsequent nebular processing.

Taken together, the framework of Di et al. (2024), the cosmochemical observations of Hu et al. (2021, 2025), and the present molecular-scale calculations are complementary: the first provides the mathematical and methodological basis, the second supplies direct natural evidence, and the third offers mechanistic insight and quantitative reference values. Integrating these perspectives enables a closed-loop approach—mechanism, parameterization, and correction—that significantly advances the application of REE, particularly Gd isotopes, in geochemistry and planetary science.

5 Conclusions

Based on the revised Bigeleisen–Mayer equation, we explicitly accounted for both mass-dependent fractionation and the nuclear volume effect in calculating the fractionation factors for the adsorption of Gd³⁺ ions onto four functionalized crown ether resins at 298.15 K. The total fractionation factors 1000lnα_tot for the ¹⁶⁰Gd/¹⁵⁵Gd pair were determined to be 0.8786 and 0.7327 for PMADB15C5 and PMADB18C6, respectively, while the corresponding values for the ¹⁶⁰Gd/¹⁵⁷Gd pair were 0.5065 and 0.4224. These results indicate that, compared to the functionalized crown ether resins, [Gd(H₂O)₉]³⁺ tends to enrich heavy isotopes. Furthermore, our findings demonstrate that mass-dependent fractionation and the nuclear volume effect act in opposite directions: in mass-dependent fractionation, [Gd(H₂O)₉]³⁺ enriches heavy isotopes relative to Gd³⁺–crown ether complexes, whereas the nuclear volume effect favors light isotopes. Among the reactions involving the four Gd³⁺–crown ether complexes, fractionation during adsorption is predominantly governed by mass-dependent fractionation, as the oxidation state of gadolinium remains unchanged. The differences in mass-independent fractionation factors 1000lnK_fs due to the nuclear volume effect are relatively small (e.g., for the ¹⁶⁰Gd/¹⁵⁵Gd pair, the maximum is –0.1659 and the minimum is –0.1204). Therefore, despite the presence of nuclear volume effect, [PMADB15C5–Gd]³⁺, which has the largest mass fractionation factor, still shows the highest total fractionation factor. Integrated analysis of thermodynamic properties and bond length data indicates that PMADB15C5 is a more reliable adsorbent for Gd isotope separation due to its higher separation coefficients. In isotopic analysis, when PMADB15C5 is used with a purification yield exceeding 95%, the fractionation induced in the ¹⁶⁰Gd/¹⁵⁵Gd and ¹⁶⁰Gd/¹⁵⁷Gd ratios can be controlled within 0.138‰ and 0.080‰, respectively.

RIGHTS & PERMISSIONS

The Author(s). Published by Higher Education Press. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0)