# American Mineralogist

- © 2002 American Mineralogist

## Abstract

Gd-doped fluorapatite (1.2 ± 0.2 wt% Gd_{2}O_{3}), synthesized from CaF_{2}-rich melts, has been investigated as single crystals and powder samples by using X-band (~9.4 GHz) electron paramagnetic resonance (EPR) spectroscopy at ~295 and 120 K. The well-resolved X-band EPR spectra yielded a previously unreported type of Gd^{3+} center “a” (S = 7/2) and also suggested the possible presence of a second and partly resolved type of Gd^{3+} center “b.” In particular, the single-crystal X-band EPR spectra of center “a” from three orthogonal-rotation planes allowed determination of the spin-Hamiltonian parameters, including the spin terms of type BS (matrix **g**) and S^{2} (matrix **D**) and the parameters associated with the high-spin terms of type S^{4} and S^{6} as well as BS^{3} and BS^{5}. The validity of the parameters has been confirmed by agreement between observed and simulated EPR spectra for both single crystals and powder samples.

The principal values of the matrices **g** and **D** indicate that the local symmetry of center “a” in the X-band EPR spectra is rhombic. The principal axis directions of the **D** suggest that this Gd^{3+} center arises from a substitution of Gd^{3+} ion into the Ca2 type of site. This assignment is supported by the results of a pseudo-symmetry analysis using the S^{4} parameters, e.g., the calculated twofold pseudo-symmetry axis coincides with the twofold rotoinversion axis of the Ca2 site. The local structural environment of this Gd^{3+} ion suggests that the ion is incorporated via the mechanism Gd^{3+} + O^{2−} ↔ Ca^{2+} + F^{−}.

## Introduction

Apatites [Ca_{10}(PO_{4})_{6}(F,OH,Cl)_{2}] are important host minerals for rare-earth elements (REEs, Z = 57 – 71) in igneous, metamorphic, and sedimentary rocks and in biomasses (e.g., Wright et al. 1984; Fleet and Pan 1995a, 1997a and references therein). Numerous studies have applied the REE compositions of apatites as petrogenetic tools in elucidating the sources and evolution of igneous, metamorphic, and sedimentary rocks (e.g., Watson and Green 1981; Fleet and Pan 1997a) and as geochemical tracers in paleoenvironmental reconstruction (e.g., Wright et al. 1984; Grandjean-Lecuyer et al. 1993; Pan and Stauffer 2000). These petrogenetic and paleoenvironmental applications all require understanding of the effects of crystal-chemical and external factors on the uptake of REEs in apatites.

The crystal chemistry of REEs in apatites, as in other Ca-rich minerals, is controlled largely by their accommodation in the Ca sites (Fig. 1a⇓; Fleet and Pan 1995a; Pan and Fleet 1996). Fluorapatite contains two (and only two) types of Ca sites with distinct stereochemistry (Ca1, Fig. 1b⇓; Ca2, Fig. 1c⇓; Hughes et al. 1989), both of which are able to accommodate REEs and a wide variety of other cations (e.g., Hughes et al. 1991; Fleet and Pan 1995b Fleet and Pan 1997a). Previous studies produced contradictory results in the site preference of REEs in apatites (e.g., Ca1, Urusov and Khudolozhkin 1974; Ca2, Borisov and Klevtsova 1963). More recently, high-precision X-ray structure refinement of natural and synthetic apatites revealed that REEs generally prefer the Ca2 site and that the site occupancy ratios (REE-Ca2/REE-Ca1) decrease monotonically with atomic number through the 4f transition series (Hughes et al. 1991; Fleet and Pan 1995b; Fleet et al. 2000a, 2000b). Similarly, Gaft et al. (1997) assigned laser-induced time-delayed luminescence spectra of natural apatites to three distinct types of Eu^{3+} center, two at the Ca2 site and one at the Ca1 site.

The present contribution reports on results of an X-band (~9.4 GHz) electron paramagnetic resonance (EPR) spectroscopic study of single crystals and powder samples of Gd-doped fluorapatite. EPR spectroscopy is a powerful probe for the local structural environments of paramagnetic ions in minerals (Calas 1988; Weil et al. 1994). Previous EPR studies of apatite-group minerals included Cr^{5+} and Mn^{2+} in chlorapatite, fluorapatite, and hydroxyapatite (Warren 1970; Greenblatt 1980; Pifer et al. 1983) and radiation defects in hydroxyapatite (Close et al. 1981; Elliott 1994, and references therein). Synthetic fluorapatite was used in this study because of the generally low abundances of Gd, and the common occurrences of other paramagnetic impurities (e.g., Mn^{2+}) in natural crystals. The well-resolved X-band EPR spectra display a previously unreported type of Gd^{3+} center (electron spin S = 7/2) denoted here by “a”; and also suggest the presence of a second type of Gd^{3+} center, “b,” which is only partly resolved at this frequency but now has been investigated in detail in fluorapatite crystals containing ~100 ppm Gd in a W-band (~94 GHz) EPR study (Chen et al. 2002). Detailed analysis of the single-crystal X-band EPR spectra allows determination of the spin-Hamiltonian parameters of the Gd^{3+} center “a” and its site assignment in the crystal structure of fluorapatite. More importantly, the results of this EPR study provide important insight into the local structural environment of Gd^{3+} ion of the center “a,” which is not available from X-ray structure refinement. In addition, the local structural information yields new insight into the substitution mechanism(s) for the incorporation of Gd^{3+}, as a representative of the trivalent REEs, into fluorapatite.

## Background

### Crystal structure of fluorapatite

Fluorapatite crystallizes in space group *P*6_{3}/*m* and point group 6/*m*, corresponding to proper rotation group *C*_{6} relevant for the EPR spectroscopy. The crystal structure of fluorapatite is made up of Ca1 and Ca2 coordination-polyhedra linked with PO_{4} groups to form a hexagonal network (Fig. 1a⇑). The Ca1 polyhedron is nine-coordinated (i.e., six shorter bonds defining an approximate trigonal prism to oxygen atoms denoted by O1 and O2, and three longer bonds to O3; Fig. 1b⇑). The environment of Ca2 is an irregular CaO_{6}F polyhedron formed by a hemisphere of 6 oxygen atoms capped by one fluoride anion (Fig. 1c⇑; Fleet et al. 2000b). The Ca2 site in end-member fluorapatite is characterized by a horizontal mirror plane normal to the **c** axis (point group *m*), which is magnetically equivalent to a twofold rotation axis parallel to **c** (Weil et al. 1973).

### Spin-Hamiltonian

Gadolinium (Z = 64) includes both even and odd isotopes, but the latter have relatively low natural abundances (i.e., ^{155}Gd = 14.8% and ^{157}Gd = 15.7%, both having nuclear spin of I = 3/2). The spin-Hamiltonian for the even isotopes of Gd (Mombourquette et al. 1986; McGavin 1987; Weil et al. 1994) can be expressed as:

(1)

where the first factors are the coefficients of the double spherical tensor operators *T*(**B**,**S**). The latter can be expanded (Mombourquette et al. 1986) in terms of other (single) spherical tensor operators as follows:

(2)

In Equations 1 and 2, *l* are even integers ranging from |*l*_{1}−*l*_{2}| to *l*_{1} + *l*_{2}; here *l*_{1} ≥ 0 and 0 ≤ *l*_{2} ≤ 2S; *m* is equal to the sum of *m*_{1} and *m*_{2}; here *m*_{k} has values −*l*_{k}, …, + *l*_{k} (k = 1, 2). Terms with *l*_{2} = 0 are of no spectroscopic interest. The term in Equation 1 with *l*_{1} = *l*_{2} = 1 can be expressed as the electronic Zeeman term β_{e}**B**^{T}•**g**•**S** (referred to as the BS term hereafter); here β_{e} is the Bohr magneton, **B**^{T} is the row vector of the external magnetic field (i.e., the Zeeman field), **g** is a 3 × 3 parameter matrix and **S** is the electron-spin column-vector operator. Here the superscript ^{T} indicates transposition of a column vector to the same vector expressed as a row vector, while the symbol • indicates the operation of a scalar product. The terms with *l*_{1} = 0 and *l*_{2} = 2 can be written as **S**^{T}•**D**•**S** (i.e., the S^{2} term); here **S**^{T} is again a row vector and **D** is a traceless 3 × 3 parameter matrix representing the electronic quadrupole effect causing the seven-line fine structure of the Gd^{3+} EPR spectra (Fig. 2a⇓). The major contributions to this quadrupole effect are the electron-electron magnetic dipole interaction, the electronic exchange, and the spin-orbit interaction, which cannot be differentiated experimentally by EPR. Moreover, crystal-field and covalency effects are also involved. Herein, matrices **g** and **D** are taken as being symmetric. The energy levels (and line positions) of the fine structure will be shifted by the magnetic-field-independent high-spin terms of type S^{4} (*l*_{1} = 0, *l*_{2} = 4) and S^{6} (*l*_{1} = 0, *l*_{2} = 6), as well as the magnetic-field-dependent terms of type BS^{3} and BS^{5} with *l*_{2} = 3 and 5 (i.e., *l* = 2, 4 and 4, 6, respectively). A detailed discussion on the form of the “generalized” spin-Hamiltonian can be found in Buckmaster et al. (1972), Al’tshuler and Kozyrev (1974), and Mombourquette et al. (1986).

## Experimental methods

### Synthesis and characterization of fluorapatite crystals

Synthesis of Gd-doped fluorapatite crystals followed the flux method of Prener (1967). Briefly, starting materials of high-purity chemicals (≥99.99%), including CaO obtained from decomposition of CaCO_{3} at 900 °C, and P_{2}O_{5}, CaF_{2}, and Gd_{2}O_{3}, from the Sigma-Aldrich Chemical Company, were weighed and mixed thoroughly to form the composition of fluorapatite having approximately 1 wt% Gd_{2}O_{3}. This mixture (27.5 g) together with flux materials (22.5 g CaF_{2}) were then placed in a tightly covered 50 mL platinum crucible. Synthesis was carried out under atmospheric pressure in a Thermolyne Muffle Furnace 46100 equipped with a programmable controller, at the Department of Geological Sciences, University of Saskatchewan. The mixture was first heated to 1375 °C and held there for 24 hours to ensure complete melting and homogenization, and was then cooled down to 1220 °C at a rate of 2 °C/h and quenched in water. The resulting products typically contained several large crystals (up to 1 cm in length and ≥1 mm in diameter) and numerous smaller grains of fluorapatite, in a quenched melt of CaF_{2}. The CaF_{2} was removed by boiling in an aqueous 20 wt% solution of Al(NO_{3})_{3}·9H_{2}O.

The large crystals of fluorapatite typically are in the form of hexagonal prisms and are terminated by {101̅1}. Most of these contain inclusions of CaF_{2} flux materials, which occur as small rods parallel to the **c** axis. A few large crystals that were clear under microscopic examination (and thus thought to be free of melt inclusions) were selected for further optical characterization and electron microprobe analysis. Back-scattered electron imaging and quantitative analyses on a JEOL JXA-8600 Superprobe equipped with one energy-dispersion and three automated wavelength-dispersion spectrometers allowed the selection of compositionally homogeneous crystals (i.e., 1.2 ± 0.2 wt% Gd_{2}O_{3}) for the EPR measurements.

### EPR experiments

Single-crystal X-band EPR spectra of the Gd-doped fluorapatite were obtained on a Bruker ESP300 spectrometer with a field modulation frequency of 100 KHz, at the Department of Chemistry, University of Saskatchewan. Each selected crystal was mounted on the acrylic three-sided holder (Fig. 3⇓; error in cubicity = ± 0.5°) by using the Laue method (Amoros et al. 1975). The X-ray diffraction analysis also confirmed that the selected crystals of fluorapatite were free of twinning. The holder was then attached to the flat end of a cylindrical Teflon jacket of a goniometer containing an adjustable EPR standard (Chen et al. 1999). This made it possible for simultaneous EPR measurement and standardization by use of the free radical 2, 2-diphenyl-1-picrylhydrazyl (i.e., DPPH). Field-swept EPR experiments were performed at 5° and 10° intervals with the magnetic field **B** in three orthogonal rotation planes by consecutively attaching the three sides of the holder to the goniometer by use of an appropriate vacuum grease (Chen et al. 1999). In these three rotations, the axes **z** and **y** of the ideal experiment system are chosen along the crystallographic axes **c** and **a** (one of three equivalent axes **a**_{1}, **a**_{2}, **a**_{3}) and the direction of the axis **x** is defined to be **y** ⊗ **z** (Figs. 1⇑ and 3⇓). The signs of **c** and **a** are arbitrary. Since the axes of this ideal system are not exactly equal to those of the actual experimental system, the latter are herein referred to as **x**′, **y**′, and **z**′ (i.e., the rotation planes for **B** are x′y′, x′z′, and y′z′; Figs. 4⇓ and 5⇓). All single-crystal X-band spectra were collected at ~295 K. The spectral resolution was 2.0 G (i.e., 1024 field data points over 2000 G in each spectrum), except 0.1 G (1024 points over 100 G) for the three most central (narrowest) lines among the seven. The average microwave frequencies were 9.199(1) GHz for the x′y′ plane, 9.187(2) GHz for the x′z′ plane and 9.205(2) GHz for the y′z′ plane. The numbers of line-position data points collected were 127, 285, and 127 for **B** in planes x′y′, x′z′, and y′z′, respectively.

Powders of the Gd-doped fluorapatite were obtained by grinding single crystals in an agate mortar and were investigated on the Bruker ESP300 X-band spectrometer at both 295 and 120 K. Except for differences in the absolute intensities of signals, the two spectra are indistinguishable in both the peak positions and peak relative intensities.

Single-crystal W-band and powder Q-band (~35 GHz) EPR spectra of synthetic fluorapatite containing approximately 150 ppm Gd were obtained at ~287 K on a Mark II spectrometer and at ~295 K on a Varian E-15 spectrometer, respectively, at the Illinois EPR Research Center (IERC), University of Illinois at Urbana-Champaign. These W-band and Q-band spectra were used in part to test the results from X-band EPR described herein. The results from these experiments are presented in Chen et al. (2002).

## Results and discussion

### EPR spectra

The single-crystal X-band EPR spectra of the Gd-doped fluorapatite mostly consist of seven anisotropic lines from the primary “allowed” transitions (Figs. 2⇑ and 4⇑), consistent with the presence of a Gd^{3+} paramagnetic center, denoted herein by “a,” with electron spin of 7/2. The X-band powder spectra of the Gd-doped fluorapatite at 120 and ~295 K are also consistent with presence of center “a” (Fig. 6⇓). A few single-crystal X-band EPR spectra with **B** close to the **y** axis, however, contain additional lines at the low- and high-magnetic field regions, which have been shown by a separate EPR study at higher frequencies (i.e., Q- and W-bands at IERC) caused by a second Gd^{3+} center, “b” (see Chen et al. 2002). No other Gd^{3+} center was detected. The Q- and W-band spectra showed that center “b” exhibits a relatively large fine structure and hence is generally not detectable in the X-band experiments.

The lines from the “forbidden” transitions of the center “a” (Fig. 2⇑) are not observed in the X-band experiment and have been shown by subsequent computer-based simulation to be very weak (up to only 0.2% of those of the “allowed” transitions). Figure 5⇑ shows the angular dependence of the primary spectra. Crossovers of the EPR transitions occur with **B** in the x′z′ and y′z′ planes but not in the x′y′ plane. A magnetic site splitting into 2 lines was detected with **B** in the plane x′z′ (Figs. 4c⇑ and 5b⇑).

The absence of hyperfine structure in the X-band EPR spectra of the Gd-doped fluorapatite is probably related to the low natural abundances of ^{155}Gd and ^{157}Gd relative to the even isotopes. At natural abundances, the intensity of each individual hyperfine line of ^{155}Gd or of ^{157}Gd is nearly twenty times lower than that of the single EPR line of the even isotopes. Another factor may be peak interference between the ^{155}Gd and ^{157}Gd spectra, since only a small difference is to be expected in the hyperfine parameters of these isotopes (Weil et al. 1994; Table G.4). Furthermore, appreciable nuclear quadrupole effects are known to be present (Chen et al., in preparation). At any rate, simulation of individual peaks at field orientations where no magnetic-site splittings are present yielded well-fitted Lorentzian lines (see below), with no other satellite peaks visible in the difference spectra.

### EPR line shapes and line widths

The basic EPR line shape is of the Lorentzian type, as confirmed by simulation of several of the EPR spectral lines. These were done for the transition between spin-projection eigenvalues *m*_{S}: −1/2 and +1/2 with the Zeeman field **B**||**z**′ and the excitation field **B**_{1}||**x**′**,** at a frequency of 9.200385 GHz and at a temperature of 295.0 K, using a simulation line width (max-min of the 1^{st} derivative) of 6.41 G. The same Lorentzian result was obtained for the transition +3/2 ↔ +1/2, with **B**||**x**′ and **B**_{1}||**z**′ at a frequency of 9.197972 GHz and at a temperature of 293.0 K, as well as for this transition with **B**||**y**′ and **B**_{1}||**z**′ at a frequency of 9.200101 GHz and at a temperature of 294.0 K, using a line width of 11.538 G for both.

Because of the *C*_{6} effective crystal symmetry, Gd^{3+} ions at the six symmetry-related sites (Fig. 1⇑) can contribute as many as 6 “symmetry-related” EPR lines for each transition. For **B**||**z**′**,** approximate superposition of all 6 lines is observed and has also been proved by the single-crystal spectrum simulation with all of the six symmetry-related sites included. Therefore, here the inherent EPR line width (ca., 6.41 G) is detected. The inherent line width may arise from D-strain effects (distributions of the spin-Hamiltonian parameters, especially in **D**, due to crystal imperfections; Wenzel and Kim 1965; Lee and Brodbeck 1977), and from a mosaic-structure effect (small crystallites with their axes **c** slightly away from the average axis **c** of the crystal sample; Shaltiel and Low 1961; Wenzel and Kim 1965; Lee and Brodbeck 1977), and possibly also from superhyperfine interactions between the Gd cation and its nearest- and next-nearest-neighbor nuclides (i.e., ^{19}F and ^{31}P, respectively).

The experimental results indicate that the line widths are a function of the transition (between *m*_{S} ↔ *m*_{S} + 1, −7/2 ≤ *m*_{S} ≤ 5/2) and increase with increasing| *m*_{S}|. If the line-position splitting (due to strain and mosaic effects, etc.) is not sufficiently large to be resolved, then the result is line broadening. This effect can be qualitatively explained by use of second-order perturbation theory (Weil et al. 1994). An expression for the splitting can be derived from the equation of the second-order energies, and shows that the positional splitting (determined by the matrices **g** and **D**) is proportional to (2*m*_{S} + 1)*m*_{S}. Therefore, for a given strength of the strain and mosaic effect, the line widths increase with | *m*_{S} |. This agrees with the X-band observations.

At all field orientations other than **B**||**z**, partial degeneracy (splitting into 2, 3, or 6 lines) is present but not necessarily resolved. Thus, an effective line width greater than the inherent one is generally observed. It has been proved by single-crystal spectrum simulation, in which all symmetry-related sites were included, that there is a weak site-splitting effect even in the plane x′y′. This splitting is caused mainly by the slight deviation of the orientation of the principal direction **D**_{3} from the **c** axis (Table 1⇓).

The effective line width in the plane x′y′ is approximately isotropic. This line width (e.g., for −1/2 ↔ +1/2) in the y′z′ plane has a minimum of 6.4 G when **B**||**z**′ and is 7.4 G when **B**||**y**′. The line shape thus represents an effective line shape exhibiting the contribution from all unresolved site splittings and strains. With **B** in the planes y′z′ and x′z′, the line widths are anisotropic due to the (line-position) contributions arising from the different magnetic sites. The effective line widths in the x′z′ and y′z′ planes were found to be linear with the magnetic-site-splitting line positions predicted by the spin-Hamiltonian obtained in this study. The same effect occurs in the plane x′y′.

### Spin-Hamiltonian optimization and EPR data simulation

The optimization of the spin-Hamiltonian parameters by use of a non-linear least-squares algorithm, determination of the energy eigenvalues by diagonalization, angle corrections for the Zeeman-field directions and the simulation of X-band EPR spectra were all performed by using the software package EPR-NMR (Mombourquette et al. 1996).

Two types of computations were performed: (1) optimization of the spin-Hamiltonian parameters, and (2) calibrations of the crystal alignments in the three rotations. During the process of optimization, both separate fitting for each of the spin-Hamiltonian terms [i.e., the BS term (**g** matrix), the S^{2} term (**D** matrix), and the terms of type BS^{3}, S^{4}, BS^{5}, and S^{6}] and the Zeeman-field directions, as well as simultaneous fitting for all terms were performed iteratively. In each fitting step, iteration was continued until no changes at the last decimal point (the machine accuracy) of the root-mean-sum-of-squares of weighted differences (RMSD) between the calculated and observed line positions occurred. Different weighing factors were adopted for the various transitions, according to the estimated accuracy of each line-position measurement. Thus the data points of the central lines and the lines away from the transition-crossing regions were weighted higher than those of the outer lines and the lines within the transition-crossing regions (Figs. 5b and 5c⇑). Initially an effective *C*_{2} crystal symmetry was assumed during the parameter optimization because of the magnetic site splitting observed in the x′z′ plane (Figs. 4c⇑ and 5b⇑). Afterward, all of the spectrum fittings and simulations were performed using *C*_{6} symmetry to match the actual magnetic crystal symmetry. The final RMSD for line positions is 2.7 G, which is less than half of the average line width for the central line (−1/2 ↔ +1/2) of the observed spectra.

The errors in the alignment of the crystal were estimated after EPR measurement by using the software package EPR-NMR. The normals to the planes x′y′, x′z′, and y′z′ were found to be (𝛉 = 0.362°, φ = 1.050°), (89.973°, 90°), and (90°, 0°), respectively. Note that no error was detected for the y′z′ plane. The final spin-Hamiltonian parameters were expressed in the ideal reference system xyz [Tables 1⇑(polar coordinates) and 2⇓].

As shown below, the six symmetry-related matrices **g** and **D** (effective magnetic point group *C*_{6}) of the center “a” (corresponding to centers denoted by “a^{1}”, “a^{2}”, …, “a^{6}”, respectively) can be assigned as the occupancy of the Gd^{3+} ions at the six symmetry-related Ca2 sites (i.e., the sites Ca2^{1}, Ca2^{2}, …, Ca2^{6}, respectively, Fig. 1a⇑). The matrices **g** and **D** for the center “a^{1}” (corresponding to a Gd^{3+} ion at the site Ca2^{1}; Fig. 1c⇑) are given in Table 1⇑, and the corresponding parameters for the S^{4}, S^{6}, and BS^{3} terms are summarized in Table 2⇑.

The contributions to the parameter optimization from the individual terms of the spin-Hamiltonian were estimated by artificially setting the parameters of each term to be zero, while all the other terms were kept to be those of the optimized results. It was found that the terms BS^{3}, S^{4}, BS^{5}, and S^{6} contribute to improvement of the RMSD by 0.34, 50.31, 0.10, and 0.10 G, respectively. The calculated parameters for the term BS^{3} have a magnitude of 10^{−5} (Table 2⇑). The parameters for the term BS^{5}, of magnitude of 10^{−7}, are not reported herein, and those for the term BS^{7} are expected to be even smaller and therefore were not included in the optimization.

Figures 2a and 2b⇑ show that the calculated energy levels for **B**||**z**′ and **B**||**y**′ are in close agreement with the results in our EPR experiments. The “allowed” transitions (| Δ*m*_{S} | = 1) occur mostly in magnetic-field regions lower than those occupied by the “forbidden” transitions (| Δ*m*_{S} | > 1). The simulated and observed single-crystal spectra, for example with **B**||**z**′ (Fig 7a⇓) and **B**||**y**′ (Fig. 7b⇓), are in excellent agreement. Here the central line positions (i.e., transition +1/2 ↔ −1/2) of the two spectra agree within the resolution limit.

Figure 5⇑ shows that the observed line-position roadmaps of the three orthogonal crystal rotations are in agreement with the calculated ones, and that the magnetic-site splittings in the x′z′ plane are successfully reproduced. Similarly, the validity of the calculated spin-Hamiltonian parameters has been confirmed by the excellent match between the simulated powder spectrum with the observed one (Fig. 6⇑). Furthermore, the spin-Hamiltonian from the present X-band EPR study has been shown to give excellent predictions for the line positions at both Q- and W-bands of the Gd^{3+} center “a” (Chen et al. 2002).

The three principal values of **g** for center “a” are almost isotropic and are slightly smaller than the free-electron value, but in fact are distinct (i.e., rhombic local symmetry). These g values are typical of Gd^{3+} (Al’tshuler and Kozyrev 1974). The principal axis directions of the matrix **g** differ from those of the matrix **D**, also indicative of a relatively low symmetry at the Gd^{3+} center “a” (see below).

In particular, the principal axis **D**_{3} for center “a^{1}” is almost parallel to the **c** axis (Fig. 8⇓), whereas D_{1} and D_{2} (not quite equal) have axes close to the normal of the plane O3^{A}-O2-O3^{J} and the horizontal bisector direction of bond angle O3^{F}-Ca2^{1}-O3^{I} (Fig. 1c⇑), strongly suggesting that the Gd^{3+} ion substitutes into the Ca2 site. The validity of the principal directions of D_{1} and D_{2} has been confirmed by a spin-Hamiltonian optimization artificially assuming D_{1} = D_{2} while varying all other parameters. The RMSD for the line positions here was 3.08 G (i.e., an increase by about 10%). A further optimization, holding D_{1} = D_{2} and forcing the D_{3} axis to be parallel to **z** yielded a relatively large RMSD value of 10.77 G, indicating that the deviation of axis **D**_{3} from **z** of the ideal experimental system is real. This deviation has been confirmed by an optimization assuming **D**_{3}||**z**, with no other constraints imposed. The RMSD value increased considerably, to 10.49 G.

The site assignment is further supported by the values of the reported parameters of the S^{4} terms (Table 2⇑). The relatively large value of B_{4}^{−3} (in the xyz system) indicates the existence of “triad” axes (i.e., 4 threefold pseudo-symmetry axes, see below). In accordance with the common conventions (McGavin 1987), the value of B_{4}^{−3} (or B_{4}^{+3}) will reach a maximum when the primary axis (**z**″) of the reference system is parallel to any “triad.” Therefore, Euler rotations were made (for the parameters of the S^{4} terms) to place the **z**″ axis parallel, in turn, to the 4 threefold pseudo-symmetry axes. The resulting values of B_{4}^{−3} (or B_{4}^{+3}) all increased dramatically, indicating the actual directions of the “triad,” which match the local symmetry of the site Ca2^{1} rather than that of the site Ca1, and thus strengthening the conclusion of the site assignment for center “a.”

Moreover, this site assignment of the center “a” is supported by the fact that the principal directions of its matrix **P** (i.e., the nuclear quadrupole effect), which has been determined from the hyperfine structure spectra of ^{157}Gd-doped fluorapatite, matches the local symmetry at the ideal Ca2 site (Chen et al. in preparation). The local symmetry at the ideal Ca1 site (Fig. 1b⇑), on the other hand, doesn’t match either the principal directions of the matrices **D** and **P** or the pseudo-symmetry of the S^{4} terms (see below) of the center “a.”

### Pseudo-symmetry analysis

The assignment of the Gd^{3+} ion for the centers “a^{i}” (i = 1, 2, …, 6) to the six symmetry-related sites Ca2^{i} has been evaluated by using the pseudo-symmetry method (e.g., Michoulier and Gaite 1972; Gaite 1980, 1987; Gaite et al. 1985a, 1985b; Mombourquette et al. 1986). This method makes use of the parameters of the fourth-degree spin-Hamiltonian terms S^{4}, because of their sensitivity to the immediate environment of paramagnetic ions. Test functions, which were developed to determine the symmetry and/or pseudo-symmetry axes (Michoulier and Gaite 1972), have a value of zero when the axis of the reference system is along a local symmetry axis, or will attain a minimum when the axis is along a local pseudo-symmetry axis. With the estimated local symmetry information, it is possible to assign the paramagnetic ion’s surroundings by comparing the calculated symmetry or pseudo-symmetry elements and their orientations with those of the known sites in the crystal structures. As mentioned above, the Ca1 and Ca2 sites in fluorapatite are characterized by vertical threefold and twofold axes, respectively (Figs. 1b and 1c⇑). Therefore, our pseudo-symmetry analysis involved searching for two- and threefold axes by using the program ROTSTO (Tennant et al. 2000). An initial search with angle interval of 1.0° was attempted to determine the approximate orientations of the symmetry or pseudo-symmetry axes, and was followed by a more precise search with an interval of 0.05°.

The results of the pseudo-symmetry analysis for the Gd^{3+} center “a^{1}” are shown in Figure 9⇓. A twofold pseudo-symmetry axis (𝛉 = 2°, φ = 257°), rather than a threefold axis, is present and matches closely to the vertical twofold rotoinversion axis at 𝛉 = 0° of the Ca2^{1} site (Hughes et al. 1989). Also, other pseudo-symmetry elements have been found but are much weaker than the twofold pseudo-symmetry axis. For example, four weak threefold pseudo-symmetry axes on the upper hemisphere (Fig. 9⇓) correspond approximately to the plane-normals of the O1-O2-O3^{F}, O3^{F}-O2-O3^{A}, O3^{F}-F-O3^{A}, and O3^{F}-O1-F faces of the Ca2^{1} polyhedron (Fig. 1c⇑). Also, a weak fourfold pseudo-symmetry axis is present approximately along the F-Ca2-O2 direction; this apparently corresponds to the arrangement of the four O3 atoms (Fig. 1c⇑), further supporting the assignment of the Gd^{3+} ion to the Ca2 site for center “a.”

The orientations of the calculated pseudo-symmetry axes allow an estimation of the distortion of the Gd-substituted CaO_{6}F polyhedron. For example, the absence of any true symmetry axes in the pseudo-symmetry analysis suggests that the local symmetry of the Gd^{3+} ion for the center “a” is rhombic (i.e., triclinic), consistent with the forms of the matrices **g** and **D**. It is noteworthy that such a situation is well-known for the Ca2 site in chlorapatite and hydroxyapatite, and is attributable to the disordering of the Cl and OH groups respectively along the **c** axis (Elliott et al. 1973; Hughes et al. 1989; Fleet et al. 2000a, 2000b). Fleet et al. (2000a) noted abnormal thermal behavior of the F anions in the X-ray structure refinement of REE-substituted fluorapatites (Fleet and Pan 1995b, 1997b) and attributed it to displacement of F ions along the **c** axis. Another possible cause is distortion of the REE-substituted sites (i.e., partial substitution of F^{−} by O^{2−}; see below).

Admittedly, the basis for the calculated pseudo-symmetry elements of this study may be relatively weak (e.g., compared to those of Mombourquette et al. 1986), because the ideal Ca2 site in fluorapatite deviates significantly from cubic or almost cubic symmetry for which the pseudo-symmetry analysis is best suited (Michoulier and Gaite 1972; Gaite 1980, 1987; Mombourquette et al. 1986). However, our results appear to be rational. Some of the angles between the calculated pseudo-symmetry axes and the bonds and faces of the ideal CaO_{6}F coordination polyhedron as well as the principal directions of D_{1} and D_{2} are probably attributable to minor distortion in the local symmetry of the Gd^{3+}-occupied Ca2 site.

### Substitution mechanism

Fleet and Pan (1995b, 1997a) discussed the following possible mechanisms for the incorporation of REEs into the Ca sites in fluorapatite:

(3)

(4)

(5)

(6)

Equations 3 and 4 can be ignored for the present fluorapatite, because neither Si nor Na is present in any significant quantities. Equation 5 involving a Ca^{2+} vacancy (i.e., □) is also unlikely for “a.” Local charge-balance requirement dictates that the vacancy occurs at a nearby Ca2 site either in the mirror plane or along the c axis, both of which would have caused significant distortion in the local symmetry of the Gd^{3+} ion of the center “a.” According to the superposition model (Brodbeck and Bukrey 1981; Kliava 1982; Gaite et al. 1985a; Mombourquette and Weil 1987), the principal direction associated with D_{3}, which indicates the direction of the largest distortion, is expected to either lie in the mirror plane or be inclined at high angles to the mirror plane, if a vacancy is present in the vicinity of the Gd^{3+} ion. Table 2⇑ and Figure 8⇑ show that the principal direction **D**_{3} is almost parallel to the **c** axis, inconsistent with the presence of such a vacancy.

Therefore, the remaining possible mechanism for the incorporation of Gd^{3+} into the Ca2 site in fluorapatite is Equation 6, which involves a simultaneous substitution of O^{2−} for F^{−} to maintain the local charge balance. A substitution of the larger O^{2−} (*r* = 1.40 Å) for F^{−} (*r* = 1.33 Å) would cause some distortion of the Gd-substituted Ca2 site and, hence, the observed rhombic local symmetry. Such a model of a GdO_{6}O coordination polyhedron is consistent with the absence of any superhyperfine structures in the single-crystal EPR spectra, which might be present for a GdO_{6}F polyhedron, owing to magnetic interaction between Gd^{3+} and ^{19}F^{−} (I = 1/2). The approximately uniaxial feature of matrix **D** (Table 2⇑) also supports this substitution mechanism. A replacement of F^{−} by O^{2−}tends to equalize the bond distance between the substituting Gd and O2 and that between Gd and the O ion at the F site; (Fig. 1c⇑), giving rise to a nearly uniaxial site symmetry (Table 1⇑).

A comparison of the F bond valences (cf., Brown 1981) between the REE-substituted fluorapatites (Hughes et al. 1991; Fleet and Pan 1995b, 1997b) and the end-member fluorapatite (Hughes et al. 1989) further supports this mechanism for the incorporation of trivalent REEs into the Ca2 site. For example, the calculated F bond valences of La-, Nd-, Sm-, and Dy-substituted fluorapatites are 0.950, 0.932, 0.911, and 0.902 v.u., respectively (Fleet and Pan 1995b), whereas that of the end-member fluorapatite is 0.895 (Hughes et al. 1989). The higher values of the F bond valence in the REE-substituted fluorapatites are consistent with partial substitution of F^{−} by O^{2−}. Fleet and Pan (1995b) also noted that the REE site occupancy ratio decreases systematically from La- to Dy-substituted fluorapatite, and shows an inverse correlation with the F bond valence. It is clear now that the partial substitution of O^{2−}for F^{−} contributes to expansion of the Ca2 site and, hence, to increased preference for the larger REEs such as La. Therefore, the present study provides further support for important roles of volatile components in controlling the uptake of REEs in apatites (Fleet and Pan 1997a and references therein). We note that, although the processes of Gd^{3+} uptake in fluorapatite may differ because of changes in the conditions of the fluorapatite crystallization, they do not determine the resulting site symmetry of the Gd-occupied Ca sites for a given substitution mechanism, i.e., our reported result has general significance.

Finally, Mackie and Young (1973) reported that minor amounts of Nd substitute for Ca in only the Ca2 positions in Nd_{2}O_{3}-doped fluorapatite, but in both Ca1 and Ca2 in NdF_{3}- doped crystals, suggesting a possible control of REE site preference by substitution mechanism. The Gd_{2}O _{3}-doped fluorapatite crystals of this study were synthesizd under similar conditions to those for the Nd_{2}O_{3}-doped fluorapatite of Mackie and Young (1973). However, the results of our study show that there are two major Gd^{3+} centers in the Gd_{2}O_{3}-doped fluorapatite and that one of them (the centers “a”) corresponds to the Gd^{3+} ion at the Ca2 site. The second Gd^{3+} center, “b,” has been proved to correspond to occupancy of Gd^{3+} ion at the Ca1 site (Chen et al. 2002).

## Acknowledgments

We thank two anonymous referees and A.M. Hofmeister for constructive reviews and helpful suggestions. We also thank W.C. Tennant for providing the computer program ROTSTO and for helpful comments; M.J. Nilges as well as R.L. Belford, and A. Smirnov for hospitality, collaborative work, and beneficial discussions at IERC; J.-L. Du for assistance with the X-band EPR experiments; and T. Bonli for assistance with the electron microprobe analysis. We gratefully acknowledge financial support from the Natural Sciences and Engineering Research Council (NSERC) of Canada.

## Footnotes

Manuscript handled by Anne M. Hofmeister

- Manuscript Received October 26, 2000.
- Manuscript Accepted August 7, 2001.