A Statistical Thermodynamic Lattice Model for Cu–In Coupled Substitution in Sphalerite

Abstract

Version 2.0. This version corrects the equilibration argument of Section 4.7. Replica-exchange sampling reaches a state 3.1 eV over the box below anything the version 1 annealing protocol attained, showing that the two-trajectory test reported there was not a sufficient test of equilibrium. Re-measuring both Hamiltonians shows the effect is confined to the full model, which shifts the Cu-In enrichment factor at 300 degrees C from 24.48 to 25.00, about 2 per cent; the electrostatics-only calculation reproduces its published values exactly. That contrast identifies the hypothetical covalent term, not the charge-compensation term, as the origin of the sampling difficulty, and the mechanism in Section 5.2 is revised accordingly. No conclusion of version 1 is altered. Preprint of a study deriving, rather than merely documenting, the copper-indium coupling in sphalerite. Source code, numerical results and figures are openly available at https://github.com/Ruqing1963/sphalerite-lattice-mc Indium in sphalerite is almost invariably accompanied by copper in a near-1:1 molar ratio, a relation reproduced across volcanogenic massive sulphide, Mississippi Valley-type, vein and skarn deposits and universally attributed to the coupled substitution of two Zn(II) by one Cu(I) plus one In(III). The relation has been documented extensively but never derived: no existing framework predicts its magnitude, its temperature dependence, or the conditions under which it should fail. This matters beyond mineral chemistry because indium is essentially never mined for its own sake, the overwhelming majority of primary production being recovered as a by-product of zinc refining, so that the global indium endowment is controlled not by a distinct ore-forming process but by the trace-element chemistry of sphalerite itself. The work constructs a four-term configurational Hamiltonian for the sphalerite cation sublattice, comprising a local charge-compensation penalty defined on the sulphur coordination tetrahedra, an elastic misfit energy in the Eshelby continuum approximation, a substitution energy referenced to the ore fluid, and a cation-vacancy term. The Hamiltonian is sampled by Metropolis Monte Carlo with Kawasaki exchange dynamics in the canonical ensemble, so that composition is conserved exactly. The central technical device is an exact topological property of the zinc-blende structure: two nearest-neighbour cations share precisely one bridging anion, and higher neighbours share none. This allows the tetrahedral electroneutrality functional to be reduced, without approximation, to a composition constant plus a nearest-neighbour pair interaction, whose coefficient is then fixed by the measured dielectric constant of ZnS rather than fitted, giving lambda between 0.227 and 0.368 eV. The resulting Cu-In pair binding energy is minus two lambda, approximately minus 0.45 to minus 0.73 eV, an order of magnitude larger than the thermal energy at hydrothermal temperatures. The ground state of the charge-compensation term at a copper to indium ratio of one to one is shown analytically to be exactly the roquesite cation ordering of CuInS2, recovered without any input describing that structure, and is degenerate across all orderings satisfying the two-copper-two-indium tetrahedral rule. Comparing the two documented charge-balancing routes at fixed indium content gives an electrostatic preference of eight lambda, about 2.4 eV, for the copper-compensated route over the vacancy-compensated one. Monte Carlo simulation at two atomic percent substitution confirms these predictions. Electrostatics alone, with all chemical interaction terms switched off, enriches Cu-In nearest-neighbour bonds sixteenfold over the random solid solution at 300 degrees Celsius and removes 96 percent of the local charge imbalance; adding a hypothetical covalent term raises the enrichment to 24.5-fold. Heterovalent bonds are enriched 2.7 times more strongly than homovalent ones, the direct signature of the charge-product coupling. Elastic terms, retained at full strength in a control calculation, produce no detectable ordering whatever, the largest strain-mediated pair interaction being some forty times smaller than the thermal energy at ore-forming temperatures. The ordering crossover lies at a reduced temperature three to seven times above the entire hydrothermal window, so short-range order is saturated wherever sphalerite forms and carries no thermometric information. Below the crossover the solute population condenses into a single connected nanodomain, and an ergodicity test with quenched and pre-annealed trajectories confirms that this state is thermodynamic rather than an artefact of kinetic arrest. Three consequences follow for ore geology. Indium tenor in sphalerite should be limited by the availability of a monovalent charge compensator rather than by indium supply, so that prospectivity should be assessed against the coupled availability of copper and indium. Vacancy compensation should be diagnostic of copper-poor systems and detectable as cation deficiency, which constitutes the sharpest near-term falsification test of the model. And the ideal dilute-mixing entropy assumed in conventional partition-coefficient models is substantially overestimated at ore-forming temperatures. The limitations are stated explicitly rather than obscured. At four thousand cation sites and four atomic percent total solute, the solute population is exhausted by a single condensed domain, so these calculations cannot resolve the solid-solution versus nanoinclusion boundary; the saturation of the lambda scan above 0.10 eV and the seed-to-seed dispersion above 0.15 eV are finite-size artefacts. The acceptance ratio collapses to about ten to the minus five below 1300 kelvin, placing simple Kawasaki dynamics at its limit. Iron is absent from the model, the treatment is rigid-lattice, Monte Carlo sweeps are not physical time, and all calculations are at a single composition at the upper end of natural indium contents. Settling the solid-solution question would require boxes of order fifty cells cubed, enhanced sampling such as parallel tempering or cluster moves, a composition scan, and first-principles parameterisation of the chemical term. All calculations use fixed seeds and are reproducible. The accompanying code includes a verification suite that checks lattice topology, energy bookkeeping, every analytic limit to machine precision, and the random-solution limit of the order parameters. This study was funded by the Hezhou Municipal Scientific Research and Development Program, Project Numbers 2024143, 2024141 and 2024104.