Probing the Local Atomic Structure of In and Cu in Sphalerite by XAS Spectroscopy Enhanced by Reverse Monte Carlo Algorithm

Abstract

The distortion of atomic structure around In and Cu dopants in sphalerite ZnS was explored by extended X-ray absorption fine structure (EXAFS) spectroscopy enhanced by reverse Monte Carlo (RMC) simulation (RMC-EXAFS method). These data were complemented with quantum chemical Density Functional Theory (DFT) calculations and theoretical modeling of X-ray absorption near edge spectroscopy (XANES) spectra. The RMC-EXAFS method showed that in the absence of Cu, the In-bearing solid solution is formed via the charge compensation scheme 3Zn²⁺↔2In³⁺ + □, where □ is a Zn vacancy. The coordination spheres of In remain undistorted. Formation of the solid solution in the case of (In, Cu)-bearing sphalerites follows the charge compensation scheme 2Zn²⁺↔Cu⁺ + In³⁺. In the solid solution, splitting of the interatomic distances in the 2nd and 3rd coordination spheres of In and Cu is observed. The dopants’ local atomic structure is slightly distorted around In but highly distorted around Cu. The DFT calculations showed that the geometries with close arrangement (clustering) of the impurities—In and Cu atoms, or the In atom and a vacancy—are energetically more favorable than the random distribution of the defects. However, as no heavy In atoms were detected in the 2nd shell of Cu by means of EXAFS, and the 2nd shell of In was only slightly distorted, we conclude that the defects are distributed randomly (or at least, not close to each other). The disagreement of the RMC-EXAFS fittings with the results of the DFT calculations, according to which the closest arrangement of dopants is the most stable configuration, can be explained by the presence of other defects of the sphalerite crystal lattice, which were not considered in the DFT calculations.

1. Introduction

In our previous works [1,2], we studied the local atomic environment and electronic state of dopants (In, Cu, Ag, Au, Cd) in natural and synthetic ZnS crystals by means of X-ray absorption spectroscopy (XAS). Analysis of the X-ray absorption near edge structure (XANES)/Extended X-ray absorption fine structure (EXAFS) spectra revealed that in substitutes for Zn in the structure of sphalerite. The “formal” oxidation state of in equals to 3⁺ and does not depend on the presence of other impurities. The difference in the “formal” oxidation states of In³⁺ and Zn²⁺ implies the substitution scheme 3Zn²⁺↔2In³⁺ + □, where □ is a Zn vacancy. In most natural sphalerites, enrichment in In is associated with elevated concentrations of Cu. According to our previous data, in the presence of In, the oxidation state of Cu is 1+, which means that the formation of the In-Cu-bearing sphalerite solid solution follows the scheme 2Zn²⁺↔Me⁺ + In³⁺. In the solid solution state, the local atomic environment of Cu derived by EXAFS spectra fitting is close to the one of pure ZnS. However, we found that the 2nd and 3rd coordination shells of Cu are considerably distorted—the atoms in each shell form two subshells, which implies significant deformation of ZnS crystal structure near Cu impurity. According to quantum chemical Density Functional Theory (DFT) calculations, such deformation can be observed in the case of the close arrangement of impurities (In and Cu), or the atom of In and the vacancy. Nonetheless, the conventional method of EXAFS spectra approximation, realized in the IFEFFIT software package [3], usually provides less accurate results for the disordered or amorphous materials than for the structures of high symmetry [4,5]. To improve the quality of the EXAFS spectra fittings, we employed the reverse Monte Carlo (RMC) method, which is highly informative in exploring the distortion of crystal structure around impurities [6,7].

In the present work, we studied the local atomic geometry of dopants in (i) sphalerite with one dopant—In, and (ii) In-Cu-bearing sphalerites. For this purpose, we used the RMC algorithm as realized in EvAX code [8].In contrast to the conventional method of EXAFS spectra approximation, the algorithm allows more efficient exploration of disorder effects beyond the first coordination shell when multiple-scattering (MS) effects are not negligible. As a result of the RMC fittings of the experimental EXAFS spectra, the atomic coordinates are derived, which can be used for further processing. The RMC fits of EXAFS spectra were supplemented by the DFT calculations and theoretical XANES spectra simulations.

2. Materials and Methods

The methods of synthesis of the doped sphalerites and reference materials, analytical methods, and the methods of XAS measurements and preliminary spectra treatment (normalization, background removing) are described in detail in our previous publication (Section 2 in Trofimov et al. [2], this issue). The sample IDs are identical throughout both the contributions.

2.1. Sample Preparation and Characterization

In the present work, we investigated three samples of synthetic sphalerite. The compositions of the samples are given in

Table 1. Formation energy of sphalerite solid solutions (E_EF) derived by DFT calculations.

System	EEF, eV
ZnS + Cu + In (In and Cu atoms are in neighboring sites)	0.18
ZnS + Cu + In (In and Cu atoms are placed far from each other)	0.28
ZnS + 2In, (2 In atoms are placed far from each other	2.99
ZnS + 2In, (2 In atoms are in the nearest positions)	3.06
ZnS + 2In + □, (In atoms and vacancy are placed far from each other)	0.89
ZnS + 2In + □, (In atoms and Zn vacancy are placed in the nearest positions)	0.35

in Trofimov et al. [2]. One sample was doped with one element—In (Sample 3757—0.7 wt.% In), and the other two samples were doped with two elements—In + Cu (Sample 4108—0.1 wt. Cu and 0.09 wt.% In; 4186—2.7 wt.% Cu and 4.9 wt.% In). The reference materials include pure sphalerite ZnS and roquesite CuInS₂.

2.2. X-ray Absorption Spectroscopy (XAS) Measurements

The Zn, In, and Cu K-edge spectra were recorded at the Rossendorf Beamline BM20 of the European Synchrotron Radiation Facility (ESRF) in Grenoble, France. Energy calibration was performed using the K-edge excitation energy of Zn (9659 eV), In (27,940 eV), and Cu (8979 eV) metal foil. The spectra of reference substances were collected in transmission mode, while the spectra of sphalerite samples were recorded in total fluorescence yield mode using a 13-element high-throughput Ge-detector. The detected intensity was normalized to the incident photon flux. Normalization and extraction of the χ(k) signal were performed using the ATHENA program (IFEFFIT software package [3], see details in Trofimov et al. [2]).

2.3. Density Functional Theory (DFT) Calculations

In our previous studies [1,2], it was demonstrated by the analysis of the X-ray absorption spectra of In-Cu-bearing sphalerites that the formation of the sphalerite solid solution takes place according to the charge compensation schemes: 3Zn²⁺↔2In³⁺ + □ and 2Zn²⁺↔In³⁺ + Cu⁺. In the present study, the formation energy of In- and In-Cu-bearing sphalerite solid solutions was calculated by means of Density Functional Theory (DFT). Besides, inthis way, the distortion of the local atomic geometry around dopants was explored.

The software package QUANTUM ESPRESSO (version 6.1) [9] was used for the calculations. We employed a projector-augmented wave description of the electron–ion interactions [10,11] with the Perdew-Burke-Ernzerhof (PBE) exchange–correlation functional. The self-consistent field (SCF) method with a 70 Ry kinetic energy cut-off for the plane waves, a 1000 Ry charge density cut-off, and a SCF tolerance >10⁻⁹ Ry was applied in the electronic structure calculations. The optimizations of the crystal structure and supercell parameters were performed using the Broyden-Fletcher-Goldfarb-Shanno algorithm for the atomic coordinates with convergence threshold 10⁻³ Ry/a.u. for the forces and 10⁻⁴ Ry for the energy. The relaxation of the atomic positions and cell parameters were applied for a 3 × 3 × 3 supercell, which contained 108 Zn and 108 S atoms, with periodic boundary conditions. In all cases, the large unit cell allowed gamma point approximation to be employed. Our previous investigations [1,12,13] demonstrated the high accuracy of this approach to DFT calculations in the case of sphalerite-type structures. The discrepancy between the calculated interatomic distances and EXAFS results was 0.02, 0.01, and 0.04 Å for the 1st, 2nd, and 3rd coordination shells, respectively.

The formation energies of sphalerite solid solutions (E_FE, [12,14]) were obtained using the calculated total energies. We considered the following models: two In atoms substitute for two Zn atoms in ZnS crystal structure; two In atoms substitute for two Zn atoms in ZnS crystal structure with the formation of a vacancy at Zn position; one In and one Cu atom substitute for two Zn atoms in ZnS crystal structure. Accordingly, the E_FE values per impurity atom after geometry relaxation were calculated as defined by the following equations:

2E_FE = E(Zn₁₀₆S₁₀₈In₂) + 2E(ZnS) + E(S) − E(Zn₁₀₈S₁₀₈) − E(In₂S₃)

(1)

2E_FE = E(Zn₁₀₅S₁₀₈In₂) + 3E(ZnS) − E(Zn₁₀₈S₁₀₈) − E(In₂S₃)

(2)

2E_FE = E(Zn₁₀₆S₁₀₈InCu) + 2E(ZnS) − E(Zn₁₀₈S₁₀₈) − E(CuInS₂)

(3)

where E is the total energy of the ZnS, ZnS-In, or ZnS-Cu-In supercells; E(ZnS), E(CuInS₂), and E(In₂S₃) are the energies calculated by DFT for pure sphalerite, roquesite, and In₂S₃ crystal structures; E(S) is calculated for the isolated S atom placed in a vacuum box. The structures with lower E_EF are more stable.

2.4. Reverse Monte Carlo (RMC) EXAFS Spectra Fitting

To explore distortion of the local atomic structure caused by the presence of dopants, we employed a reverse Monte Carlo (RMC) method [15,16]. The evolutionary algorithm (EA) of the RMC method was used for the EXAFS spectra fitting as implemented in EvAX code [8]. An important advantage of the RMC method is the automatic handling of multiple scattering paths (MS). The inclusion of MS paths in the conventional EXAFS fitting procedure is usually associated with a drastic increase in the number of fitted parameters, which often makes the obtained results unreliable.

The RMC fitting aimed to obtain the best agreement between the experimental and the calculated spectra by small variations of the atomic coordinates. The maximum displacement of atoms from their initial positions was fixed at 0.4 Å; this value is enough to describe the thermal movement of atoms and distortion of the sphalerite crystal structure. A comparison of experimental and theoretical spectra was made in the Wavelet transform (WT) space [17,18] (see Trofimov et al. [2] for details of the WT calculation). The experimental EXAFS spectra were fitted in the R-range between 1.3 and 4.5 Å; the k-range was chosen for every spectrum depending on the signal/noise ratio. The E₀ shift was determined for each EXAFS spectrum by fitting the first peak of Fourier transform magnitude.

Calculations of the scattering amplitude and phase functions were performed by the FEFF8 [19] code. We used atomic structures derived from ZnS crystal by substituting the central Zn atom by In or Cu as a starting geometry for EXAFS refinement. Such clusters contain a central absorbing atom (In or Cu) and three coordination shells with 4 S, 12 Zn, and 12 S atoms (see Figure S1). The radii of coordination shells around central Cu or In atoms were corrected to match the results of the conventional EXAFS spectra fitting [2]. To speed up the calculations, in the calculations of doped sphalerites, we used isolated clusters without periodic boundary conditions. EXAFS spectra at the Cu K-edge were fitted using a cluster of ZnS structures, with and without In in the second coordination shell, to examine the possibility of the close arrangement of Cu and In dopants. The In K-edge EXAFS spectra were fitted without Cu in the second coordination shell since Cu has the same scattering amplitudes and phase functions as Zn. As will be shown below, the 5 Å cluster radius is enough to fit the EXAFS spectra in the R range from 1.3 to 4.5 Å. It is important to have at least a few thousand absorbing atoms involved in the fitting procedure to ensure good quality of the calculated spectra, Radial Distribution Function (RDF) curves, and reliable results of the spectra approximation [20]. For the RMC procedure, we used 128 ensembles of 432 clusters placed at significant distance from each other to mimic uniform distribution of dopants in sphalerite.

To obtain the reference RDFs, we performed a fit of the EXAFS spectrum of pure ZnS measured at the Zn K-edge. The RMC-EXAFS procedure for pure sphalerite was applied for a 4 × 4 × 4 supercell with the periodic boundary condition. In all the calculations, the amplitude reduction factor S₀² was fixed at values estimated for the references in [1]: S₀²(Zn) = 0.85, S₀²(Cu) = 0.75, S₀²(In) = 0.95.

As a result of the RMC optimization, the atomic coordinates were derived. Based on these data, the RDFs were built for each pair of atoms, and the structural parameters (interatomic distances) were obtained. Gaussian functions were used to approximate the calculated RDF curves. Coordination numbers associated with Gaussian decomposition were calculated by integration of the RDFs. Analysis of the RDF shape provided information which cannot be derived by means of the conventional EXAFS fitting procedure, including detection of small distortion of the local atomic structure.

2.5. XANES Spectra Calculation

Two series of XANES spectra calculations were performed. The first series included In and Cu K-edge XANES spectra for which the models of the local atomic geometry around the absorbing atom were obtained by means of DFT geometry optimization. The second series (only Cu K-edge XANES) was based on the atomic cluster geometries derived by means of RMC-EXAFS fitting. To model the RMC-EXAFS spectra, about 100 random geometries obtained by RMC-EXAFS fitting were averaged.

Theoretical calculations of Cu and In K-edge XANES spectra were performed using the finite difference method (FDM) employing the FDMNES program [21,22]. Relativistic self-consistent field FDMNES calculations were carried out with the exchange–correlation part of potential in a local density approximation [23]. The final electronic states were calculated in a full core-hole screening. Atomic clusters inside the spheres with radii of 5 Å were used in the calculations. Many-body effects and core-hole lifetime broadening were accounted for in the frame of the arctangent convolution model [24]. Parameters of the XANES spectra convolution were chosen to obtain the best co-incidence of the experimental and calculated data.

3. Results

3.1. DFT Calculations

The optimization of crystallographic parameters of the pure ZnS obtained for a 3 × 3 × 3 supercell yields the lattice parameter a_DFT = 5.439 Å, which is slightly bigger than the experimental XRD value of a_XRD = 5.410 Å, with relative deviation of 0.5%. This accuracy level is high enough to model the local atomic geometry and compare it to the structural parameters derived from EXAFS spectra.

In case of substitution of Zn atoms by dopants, the atomic structure of sphalerite is preserved, but the interatomic distances are relaxed to match the In and Cu ionic radii. The close arrangement of dopants and/or Zn vacancies results in the significant distortion of the coordination shells, which is translated to the splitting of the interatomic distances. Table S1 (Supplementary Materials) shows the interatomic distances of the optimized structure of ZnS around the In and Cu dopants. Localization of dopants (Me = In, Cu) and vacancy in the nearest cationic positions causes splitting of the Me-Zn distances up to 0.12 and 0.16 Å in the 2nd and 3rd coordination shells, respectively. This splitting is sufficiently high to be detected by EXAFS spectroscopy. The range of variation of the interatomic distances in the 1st coordination shell, which consists of four S atoms, is less pronounced (<0.05 Å). It is difficult to determine the contributions of individual atoms, which are characterized by small variations of interatomic distances within the 1st coordination shell, by means of EXAFS spectra analysis. Accordingly, the effect of the close arrangement of defects in sphalerite structure—In and a Zn vacancy, and In and Cu atoms—results in important splitting of interatomic distances in the 2nd and 3rd coordination shells, which can be detected and characterized by means of EXAFS.

Calculated values of the formation energy of sphalerite solid solutions are presented in Table 1. The positive sign of E_EF obtained in all calculations implies endothermic character of the substitution reactions. Embedding two In atoms in ZnS without a vacancy creation does not obey the charge compensation scheme and is characterized by high value of E_EF. Therefore, this configuration is not stable. Comparison of E_EF calculated for different arrangements of In atom and a vacancy in the cationic sublattice shows that the location of In and a vacancy in the nearest cationic positions is energetically favorable. The same conclusion can be reached for Cu and In atoms: the close position of these dopants is more stable because of low values of E_EF.

To summarize the results of the DFT calculations, the close arrangement of the impurities—In and Cu atoms, or In atom and a vacancy—is energetically more favorable than the random distribution of the defects. The close arrangement of the defects causes distortion of the coordination polyhedra of the 2nd and 3rd coordination shells which is high enough to be detected by means of EXAFS. These results will be compared with the data derived from the experimental EXAFS spectra treatment.

3.2. EXAFS Spectra Analysis (RMC-EXAFS)

The main structural information which is obtained from the RMC-EXAFS spectra fitting is the shape and position of RDF peaks, which are calculated using atomic coordinates produced by the fit. In the primary stage of the calculations, the Zn K-edge EXAFS spectrum of pure ZnS was fitted using the RMC method. Figure 1 compares the experimental and calculated spectra of ZnS. The RDF consists of three distinct peaks which correspond to the 1st, 2nd, and 3rd coordination shells. These peaks are accurately fitted by Gaussian functions (optimized Gaussian functions are shown by solid curves in Figure 1). The Gaussian curves that correspond to the 2nd coordination shell (Zn-Zn RDF, red peak in Figure 1) will be further compared with RDFs derived by the analysis of In and Cu K-edge EXAFS spectra fitting. The centroid positions and variances of Gaussian functions, which fit coordination shells around Zn in the ZnS structure, are presented in

Table 2. Structural parameters derived from RMC-EXAFS fitting (all RDF peaks are fitted with one Gaussian function). Coordination numbers of sphalerite structure were fixed for all samples. Results of the present study (RMC) are compared with values derived by conventional EXAFS analysis performed using IFFEFIT software [2].

Experiment	Fitting Procedure	Coordination Spheres Around Absorbing Atoms
		S		Zn		In		S
		R,Å	σ2*103, Å2	R,Å	σ2*103, Å2	R,Å	σ2*103,Å2	R,Å	σ2*103,Å2
ZnS, Zn K-edge	RMC	2.34	5.3	3.83	14.4	-	-	4.50	16.1
	IFEFFIT, [2]	2.34	5.0	3.85	17.0	-	-	4.46	15.0
Sample 3757, In K-edge	RMC	2.45	3.3	3.91	13.0	-	-	4.48	13.7
	IFEFFIT, [2]	2.45	4.0	3.91	14.0	-	-	4.48	14.0
Sample 4108, In K-edge	RMC	2.46	2.0	3.91	13.7	-	-	4.49	12.3
	IFEFFIT, [2]	2.46	3.0	3.91	15.0	-	-	4.49	11.0
Sample 4186, In K-edge	RMC	2.45	3.4	3.91	13.8	-	-	4.49	14.2
	IFEFFIT, [2]	2.46	4.0	3.91	16.0	-	-	4.47	12.0
Sample 4108, Cu K-edge	RMC, with In in 2nd shell	2.28	4.3	3.84	21.3	3.87	8.8	4.47	14.6
Sample 4108, Cu K-edge	RMC, no In in 2nd shell	2.28	4.3	3.84	21.0	-	-	4.47	14.7
	IFEFFIT, [2]	2.30	5.0	3.81	7.0	-	-	4.44	11.0
				4.01
Sample 4186, Cu K-edge	RMC, with In in 2nd shell	2.29	7.5	3.84	22.5	3.88	24.3	4.46	13.4
Sample 4186, Cu K-edge	RMC, no In in 2nd shell	2.30	7.9	3.84	24.0	-	-	4.47	22.5
Sample 4186, Cu K-edge	IFEFFIT, [2]	2.31	6.0	3.76	11.0	-	-	4.31	9.0
				3.92		-	-	4.52

.

The next step consisted of fitting EXAFS spectra measured at K-edges of dopants in Samples 3757 (In K-edge), and 4186 and 4108 (In K-edge, Cu K-edge). The results of the fits are shown in Figure 2 (In K-edge spectra) and Figure 3 (Cu K-edge spectra). The positions of the Gaussian curves, which were used to approximate the RDFs, are listed in Table 2. Visual examination of the images in Figure 2 and Figure 3 shows good agreement between the experimental data and model spectra. The RDFs are described by three distinct peaks which correspond to coordination shells of the probed cation with 4 S, 12 Zn, and 12 S atoms. Images of Cu K-edge EXAFS fitting with In in the 2nd coordination shell and corresponding RDFs are shown in Supplementary Materials (Figures S2 and S3). The coordination shell radii derived from EXAFS fitting for different samples measured at the same edge are very close to each other, which demonstrate high stability of the RMC fitting procedure. Adding In to the 2nd coordination shell doesn’t improve the description of the experimental EXAFS data (compare the fit results in Figure 3 and Figure S2), which is in agreement with the performed wavelet analysis applied to the collected EXAFS data [2].

The first peak in the RDF plots, which represents the mean bond length between the central atom (Cu or In) and its nearest neighbors (4 S atoms), is of asymmetrical and sharp shape. The peak asymmetry can be explained by strong distortion of the nearest atomic coordination or limited k-range of measured EXAFS spectra, especially in the case of the Cu-K edge. We believe that the limited k-range of the spectra is the more probable reason for the observed asymmetrical peak shape. The second peak of RDF functions at ~3.9 Å (In K-edge) and 3.84 Å (Cu K-edge) corresponds to the 2nd coordination shell, which was approximated by two models—12 Zn or (only for Cu K-edge) 11 Zn + 1 In atoms. The last peak in RDFs corresponds to 12 S atoms located in the 3rd coordination shell. As follows from data presented in Table 2, the mean radii of the three coordination shells determined by RDF peak fittings with Gaussian functions are in good agreement with the results of conventional EXAFS spectra fits (Trofimov et al. [2]).

The most important information concerning the atomic arrangement around a dopant can be obtained by exploring the 2nd coordination shell. The peaks, which correspond to the metal atoms located in the 2nd coordination shells of Cu and In, have asymmetrical shape for all fitted spectra and cannot be accurately approximated by single Gaussian function (Figure 2 and Figure 3). The asymmetry can be accounted for by decomposition of the peaks into two Gaussian components for the In-Zn RDFs of Samples 3757 and 4108, and three Gaussian components in the case of the In-Zn RDF of Sample 4186 and Cu-Zn RDFs. The results of the decomposition are given in

Table 3. Structural parameters derived from Cu-Zn and In-Zn RDF decomposition. Results of the present study (RMC) are compared with values derived by conventional EXAFS analysis performed using IFFEFIT software [2].

Parameter	3757 In K-Edge		4108 In K-Edge		4186 In K-Edge		4108 Cu K-Edge		4186 Cu K-Edge
	RMC	[2]	RMC	[2]	RMC	[2]	RMC	[2]	RMC	[2]
Centroid position, Å	3.90 4.05	3.91	3.88 4.04	3.91	3.73	3.91	3.66	3.81 4.01	3.61	3.76 3.92
					3.89		3.82		3.79
					4.03		4.00		3.97
Gaussian variance × 103, Å2	11.2 1.7	14.0	9.9 5.3	15.0	3.7	16.0	8.8	11.0	5.4	7.0
					4.4		7.2		7.6
					2.8		6.5		8.9
Coordination number	11.8 0.2	12	11.0 1.0	12	1.4	12	2.8	6.4 5.6	1.0	5.4 6.6
					7.8		6.2		5.6
					2.8		3.0		5.3

and are shown in Figure 4.

The In-Zn RDF of Sample 3757 (In is the only dopant) derived from RMC fitting of the In K-edge EXAFS is very close to Zn-Zn RDF in pure sphalerite, but bimodal convolution shows appearance of a small shoulder at 4.05 Å. This small component of bimodal Gaussian function can be assigned to the interaction between In impurity and Zn vacancy. The In-Zn RDFs of Samples 4108 and 4186 (In-Cu-bearing sphalerites) have much more intense features at nearly 4.04–4.06 Å. The increase in the intensity of the RDF shoulder (or the second component of the Gaussian function) can be caused by the presence of Cu impurity atoms in ZnS crystal structure. Splitting of the 12 Zn atoms between subshells demonstrates an increase in ZnS crystal structure distortion around In atom in a row 3757→4108→4186. We note that in the first two samples (3757, 4108), the coordination numbers associated with the distorted Zn atoms in the 2nd coordination shell of In do not exceed N = 1 (Table 3). This means that the distortion of the 2nd coordination shell of In is almost negligible. Sample 4186 contains a high amount of Cu and In impurities (nearly 2 mol.% of Cu and In [2]), and in this case, the 2nd coordination shell around In is much more distorted.

The positions of the three Gaussian functions, which were used to fit the Cu-Zn RDFs in Samples 4108 and 4186 (Cu K-edge EXAFS), are very stable and only slightly change from one sample to another. The significant splitting of the 2nd coordination shell of Cu results in the formation of three groups of Zn atoms, two of which are shifted by 0.16–0.2 Å from the average position of 3.84 Å (the average Cu-Zn distance is close to the Zn-Zn distance of 3.83 Å in pure sphalerite, Table 2).

All Cu K-edge EXAFS spectra were accurately fitted without In in the 2nd coordination shell. Addition of an In atom does not improve the agreement between the experimental and calculated spectra. The shapes of the RDFs derived from the EXAFS spectra fitting with and without In are close to each other (Figure S4). Since the WTs (Figure S2) do not reveal any heavy atoms in the 2nd coordination shell of Cu, we can conclude that the 2nd coordination shell consists mostly of Zn atoms.

The results of RMC-EXAFS spectra fittings demonstrate that the structures of the local atomic environments of In and Cu dopants are not identical. The local atomic symmetry of In depends on the concentration of Cu. The local atomic structure of In in sphalerite with only In impurities is close to one of Zn in pure sphalerite. An increase in Cu and In concentration causes distortion of the local atomic structure of In. In contrast, the interatomic distances in the 2nd coordination sphere of Cu exhibit considerable splitting, which doesn’t depend on the amount of impurities in sphalerite. This splitting cannot be explained by the close arrangement of the dopants (In and Cu) as (i) no heavy In atoms were detected in the second shell of Cu (Trofimov et al. [2]), (ii) inclusion of In atoms in the model does not improve the quality of the RMC-EXAFS fitting, and (iii) the 2nd shell of In in Sample 4108, with strong distortion of atomic structure around Cu, is almost undistorted. Therefore, despite the splitting of interatomic distances in the 2nd shell of Cu, RMC-EXAFS results argue for random (or at least, not close) arrangement of the dopants—In and Cu. The disagreement with the results of DFT calculations, according to which the closest arrangement of dopants is the most stable configuration, can be explained by the presence of other defects of the sphalerite crystal lattice, which were omitted during the DFT calculations.

3.3. XANES Spectra Simulation

3.3.1. In K-Edge

Since XANES spectra of Samples 3757, 4108, and 4186 are identical [2], we used the experimental XANES spectrum of Sample 3757 (In-bearing sphalerite) as a reference to compare with calculated In K-edge XANES spectra. The simulations of In K-edge XANES spectra were performed using atomic configurations optimized by means of DFT. The In-S distance calculated by DFT differs slightly from the one derived by EXAFS analysis. Therefore, XANES calculations were performed with atomic coordinates corrected to match the In-S distance obtained by EXAFS. The calculation of XANES spectra were performed for (i) sphalerite crystal structure with two Zn atoms substituted by In and Cu atoms (Figure 5a); (ii) sphalerite crystal structure with two Zn atoms substituted by two In atoms and one Zn vacancy (Figure 5b); (iii) sphalerite crystal structure with two Zn atoms substituted by two In atoms (this structural model is beyond the charge compensation scheme) (Figure 5c). Calculations for the three models listed above were performed with dopants and vacancy placed in the nearest positions and far from each other to check the influence of the crystal structure defects arrangement on XANES spectra. As can be seen in Figure 5a–c, all calculated spectra reproduce correctly all features of the experimental spectrum. The difference between the calculations performed for the different structural models is negligible, which can be explained by the significant natural broadening of In K-edge spectra.

3.3.2. Cu K-Edge

The Cu K-edge XANES spectrum of Sample 4186 (In-Cu-bearing sphalerite) was used as a reference to compare with the calculated spectra. The results of the Cu K-edge XANES calculations are shown in Figure 6. The figure compares the results of simulations which are based on atomic geometries obtained by DFT optimization and RMC-EXAFS fitting. In both cases, calculations were performed on structures with and without In in the 2nd coordination shell. The experimental XANES spectrum contains spectral features A, B, C, and D. All the structures used and both the calculation methods (DFT and RMC-EXAFS) reproduce the overall shape of the experimental spectrum. However, the geometries obtained by RMC-EXAFS fitting provide better agreement with the experimental data regarding the energy position and intensity of spectral features B and C. Another useful conclusion follows from the comparison of calculated spectra of structures with and without In in the 2nd coordination shell: influence of one In atom on the shape of the XANES spectrum is negligible.

4. Discussion

In the present study, the local atomic geometry of In and Cu in sphalerite ZnS was determined by the reverse Monte Carlo (RMC) method applied to EXAFS spectra approximation. The EXAFS spectra of synthetic In- and In-Cu-bearing sphalerites were fitted by means of RMC with determination of the interatomic distances (for the 2nd shell atoms) in the first three coordination shells of dopants. The advantage of the RMC-EXAFS method in comparison with conventional EXAFS spectra fitting (realized, for example, in the IFEFFIT software) is the possibility to obtain reliable information about small distortion of coordination shells. In general, the results of the present study are in agreement with our data obtained previously by conventional EXAFS spectra fittings. Incorporation of In into the cationic position of ZnS crystal structure causes expansion of the nearest-to-In coordination polyhedra. The bond lengths increase (relatively to pure ZnS) from 2.34 to 2.45 Å in the 1st coordination shell (N_S = 4), from 3.85 to 3.91 Å in the 2nd coordination shell (N_Zn = 12), and are close to the Me-S distance in the pure sphalerite in the 3rd coordination shell (N_S = 12, R_In-S = 4.48 Å). The In Radial Distribution Function (RDF) consists of three distinct peaks, which correspond to the first three coordination shells. We observe no splitting of the coordination shells of In into subshells when In is the only dopant. The presence of Cu in sphalerite distorts the local atomic surrounding of In, and the distortion extent of the In-Zn coordination shell depends on the concentration of impurity atoms. In Sample 4108 (concentration of In and Cu is less than 0.1 mol%), the shape of In-Zn RDF is close to the one of In-Zn for the sphalerite with only In dopants. An increase in concentrations of In and Cu causes the splitting of In-Zn RDF into three contributions.

The Cu impurity in the In-Cu-bearing sphalerite also presents in the solid solution state substituting for Zn. Since the ionic radii of Zn and Cu are close to each other, only small changes of the average coordination shell sizes take place around Cu compared to the pure sphalerite. However, analysis of Cu RDF shows significant deformation of the 2nd coordination shell. The Cu-Zn RDF peak, which corresponds to the 2nd coordination shell of the cation, can be decomposed into three Gaussian contributions with centroids at ~3.64, ~3.80, and ~3.99Å, which means considerable deformation of ZnS crystal structure near Cu impurity.

The In and Cu K-edge XANES spectra were simulated using structural models obtained by the DFT and RMC-EXAFS methods. It was found that the deformed coordination polyhedra around Cu are better reproduced by the RMC-EXAFS method, which provides better agreement between the simulated and the experimental spectra. The XANES spectra calculations confirm conserving of the of ZnS structure near In impurity and its deformation near Cu impurity.

Substitution In³⁺→Zn²⁺ results in the formation of a positively charged site located on the dopant. On the contrary, the Cu⁺→Zn²⁺ substitution results in the formation of a negatively charged site. The DFT calculations show that the close arrangement of In and Cu is more energetically favorable than random distribution of dopants, which can be explained by the local charge balance. The DFT shows significant distortion in ZnS crystal structure around the dopants (In and Cu). The splitting of interatomic distances, calculated by means of DFT, reaches 0.10 Å for the 2nd shell around Cu and 0.12 Å for the In 2nd shell. This coordination shell splitting should be observed by means of EXAFS, even in the samples with low concentrations of In and Cu. Analysis of the experimental EXAFS spectra shows considerable deformation of ZnS geometry around Cu, which doesnot depend on the concentration of Cu or In. In the case of In, the local atomic distortion depends on the concentration of impurities, and its extent correlates with the amount of dopants in sphalerite. This implies that Cu and In impurities in sphalerite with low concentration of dopants are placed at a significant distance from each other. The apparent disagreement between the results of the DFT calculations and experimental EXAFS spectra fittings can stem from the oversimplified model of doped sphalerite structure used in the DFT calculations. In the calculations, we employed the simplest models with defects localized at the atoms of In, Cu, or vacancy. The real structure of doped sphalerite can be more complicated and contain defects, which were not considered in our calculations.

5. Conclusions

The reverse Monte Carlo (RMC) simulation was applied to interpret the In and Cu K-edge EXAFS spectra of In- and In-Cu-bearing sphalerites. In the studied samples, both In and Cu occur in the solid solution state, substituting for Zn. The local atomic structure of In in ZnS without Cu is undistorted. The presence of both dopants in ZnS causes variations of the interatomic distances, which is clearly seen in the splitting of the 2nd coordination spheresof In and Cu. The degree of distortion of the 2nd coordination sphere of In correlates with the concentration of In + Cu, but is weak in general. The degree of distortion of the 2nd coordination sphere of Cu shows no correlation with the concentration of the dopants, but is much more pronounced than that of In. The results of RMC-EXAFS fittings argue for a random (or at least not close) arrangement of the dopants—In and Cu. In contrast to the RMC-EXAFS method, the Density Functional Theory (DFT) calculations showed that the close arrangement of In and Cu is energetically favorable. We attribute the disagreement between the RMC-EXAFS and DFT methods to the presence of defects or vacancies in ZnS structure in close vicinity of the dopant sites, which were not taken into account in the DFT calculations. The nature and types of such defects are the questions of further studies.