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

Abstract

This record contains the complete reference implementation, the numerical results underlying every figure and table, and the manuscript for a statistical thermodynamic lattice model of copper-indium coupled substitution in sphalerite. 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. The repository is organised as follows. The src directory holds the core module, which contains the lattice construction, the Hamiltonian, the Monte Carlo engine, the order parameters and the verification suite; an experiment driver that reproduces all four figures and the results table; and an analytic script that locates the ground states of the charge-compensation term by annealed minimisation. The data directory holds the numerical results underlying all figures together with the full console log of the production run. The figures directory holds Figures 1 to 4 in both vector PDF and raster PNG form. The paper directory holds the manuscript source and compiled PDF in English and in Chinese translation. All runs use fixed seeds and are reproducible. The verification suite checks lattice topology including the sharing lemma, energy bookkeeping to a deviation of 1.8 times ten to the minus fourteen electron volts, every analytic limit to machine precision, the random-solution limit of the order parameters over twenty-four independent configurations, and conservation of composition. The software requires Python 3.10 or later with NumPy, matplotlib for figures, and optionally Numba, which gives roughly a fiftyfold speed-up; the full experiment suite completes in about eight minutes on a single core. Licensing is mixed and is stated here for clarity. The source code in the src directory is released under the MIT License. The data and figures are released under the Creative Commons Attribution 4.0 International License. The manuscript remains the copyright of the authors. This study was funded by the Hezhou Municipal Scientific Research and Development Program, Project Numbers 2024143, 2024141 and 2024104.