Solid solutions

A pure solid has  : it is its own standard state, its chemical potential does not depend on how much of it there is, and it either precipitates or it does not. A solid solution is a single crystalline phase whose composition varies continuously between end-members, and that changes the thermodynamics qualitatively: the activity of an end-member depends on the composition of the phase, so the phase can absorb a component gradually instead of appearing all at once.

This matters in cement more than anywhere else. C-S-H is not a compound with a formula but a phase of variable Ca/Si; the AFm phases exchange sulfate, carbonate and hydroxide on the same interlayer site; alkalis partition into hydrates rather than staying in solution. A model that only knows pure phases predicts sharp appearances and disappearances that are not observed.

1. Ideal mixing: the activity is the mole fraction

Put moles of each end-member into one phase, with mole fractions  . Even with no interaction at all, mixing lowers the Gibbs energy, because there are more ways to arrange a mixture than a pure solid. Per mole of phase,

which is entropy alone — the of mixing. Differentiating for one end-member gives its chemical potential and hence its activity:

That is IdealSolidSolutionModel, and it is the default. Note what it already buys: an end-member at   has   , so its saturation index is shifted by decades relative to the pure phase. A trace component is stabilized by being diluted in a host, which is why a solid solution can take up an ion that would never precipitate as its own phase.

2. Non-ideal mixing: the excess Gibbs energy

Real end-members interact. Everything beyond the ideal term is collected into the excess Gibbs energy , and the activity coefficient is its partial molar derivative:

and then   . This is exactly what the code does: _solid_solution_lna! writes log(x[k] + ϵ) + _excess_ln_gamma(model, k, x, T) into the log-activity vector, so the two terms are visibly separate and a model supplies only the second.

3. RegularSolutionModel — one interaction energy per pair

The simplest non-ideal form keeps only pairwise contacts, with one energy per pair of end-members:

which is the symmetric multicomponent Margules expression, and it is _excess_ln_gamma(::RegularSolutionModel, …) line for line. For a binary it collapses to the form usually quoted:

(J/mol) is the energy of an – contact relative to the average of – and –, so its sign is physical:

•   — unlike neighbors are unfavorable.  , the phase resists mixing, and above a threshold it unmixes;

•   — unlike neighbors are favorable.  , mixing is stabilized beyond ideal, and the phase may order.

W must be square and symmetric, which a validating inner constructor enforces; the reason it is inner rather than outer is recorded in the source, and it is a mistake worth not repeating.

The threshold is  , and it is checkable

For a symmetric binary, the molar Gibbs energy of mixing is       . A phase is unstable to unmixing where that function is concave, and

which at   is  : negative as soon as  . So   — about 5 kJ/mol at 25 °C — is the critical point of a symmetric regular solution: below it one homogeneous phase, above it a miscibility gap that widens as grows. The threshold is evaluated, and the sign of the second derivative tabulated on either side of it, in Solid solution models, in numbers.

This matters in practice because nothing detects it. A solid solution is entered as one phase, the activity expression goes on returning values past the threshold, and those values describe a metastable single phase. Checking against 2 is the caller's job.

4. RedlichKisterModel — a binary that need not be symmetric

A regular solution is symmetric by construction: swapping the end-members leaves unchanged. Real binaries often are not, so the excess is expanded in a polynomial of the composition difference. As implemented, for two end-members,

with in J/mol. Read the roles off the formulas: is the symmetric term and is exactly a regular solution's ; and are the asymmetric corrections, so setting them to zero must reproduce the regular model. That is an identity between two independently written methods, and it is checked rather than trusted in Solid solution models, in numbers, where the two agree to the last digit.

RedlichKisterModel requires exactly two end-members, and SolidSolutionPhase refuses the combination at construction rather than letting a ternary reach an expression written for a binary. For three or more end-members the choices are the ideal model or RegularSolutionModel with a full matrix.

5. What the three models do to an activity

Evaluated over the whole composition range in Solid solution models, in numbers, with    kJ/mol — inside the stable range of §3, so that the numbers describe a phase that stays homogeneous. Two things to read off it.

A positive pushes the activity of a dilute end-member above its mole fraction: the host is rejecting it. A negative pulls it below: the host is stabilizing it. And at   all three models converge, because   as the phase becomes pure — the standard state of an end-member is the pure end-member, so they agree there by construction.

The ideal case is the one to keep in mind for cement. An end-member at   has  , so its saturation index sits three decades below the pure phase: a trace component is stabilized simply by being diluted in a host. That is how a solid solution takes up an ion which would never precipitate as a phase of its own, and it is the whole reason C-S-H and the AFm phases can absorb alkalis, sulfate and carbonate continuously instead of in jumps.

6. The spinodal, the common tangent, and why one amount is not enough

Sections 3 and 4 gave models whose excess term can be strong enough to make the molar Gibbs energy of mixing concave over an interval. This section is about what that means for the equilibrium, because it is not a detail of the model: it changes what the answer is.

Convexity is the hypothesis under everything else

For a binary write the molar Gibbs energy of mixing, in units of , as

The ideal part is convex everywhere — its second derivative is     — so unmixing is always the excess term's doing. Where

the phase is inside its spinodal, and there the straight line joining two compositions lies below the curve between them. A system at an intermediate overall composition therefore lowers its energy by separating into those two, and the minimum of is not a point but a pair.

For the symmetric regular solution of section 3 the criterion is exactly    at  , so unmixing begins at   — the threshold section 3 already quoted, here derived from the same inequality. spinodal_interval evaluates the second derivative on a grid rather than in closed form, which is why the same routine covers the asymmetric Redlich-Kister case without a separate derivation.

The answer inside a gap is the common tangent

Which pair? The two compositions   at which the chemical potentials of both end-members agree:

which says geometrically that one straight line is tangent to the curve at both points — the common tangent. Between and the equilibrium state is a mixture of the two, in the proportions the lever rule gives, and the energy follows the tangent line rather than the curve.

The common tangent construction is wider than the spinodal: the binodal contains the spinodal, and between the two the phase is metastable rather than unstable. A minimization sees only the tangent.

Computing the pair: what PHREEQC does, and what it costs

The pair is not found by minimizing. It is computed from the model alone, by common_tangent: two equations in two unknowns, solved by Newton with ForwardDiff supplying the derivatives. Nothing about the rest of the chemical system enters — not the aqueous solution, not the other phases, not the element budget — which is why it costs microseconds and can be used as the starting point of a full equilibrium rather than as its result.

This is the construction PHREEQC uses for binary solid solutions, after (Glynn and Reardon, 1990). It was validated here against a case with a closed form: for a symmetric model the pair must be symmetric about  , and it comes out at with a residual of   and the symmetry exact to the last bit.

miscibility_split then applies the lever rule. Inside the gap the two compositions are fixed and only their proportions move with the overall composition — which is what makes a miscibility gap flat in a phase diagram, and what a single-composition answer cannot reproduce at any resolution.

Where the three codes agree, and where this one adds a step

The criterion is the same object in GEM-Selektor and here, arrived at independently from the same KKT conditions — its phase stability index   is term for term what phase_split_measure computes (Kulik et al., 2013). The construction of the pair is PHREEQC's, after Glynn & Reardon. Neither this package nor the others invented either.

What differs is only where the duplication comes from. GEM-Selektor's users get the right answer inside a gap because CEMDATA18 ships the AFm and AFt binaries under two names each, so the declaration is already doubled in the database; here it is asked for by a keyword. Both are sound, and the database route has the advantage of being the published one.

The step this package adds is the refusal: SolidSolutionPhase evaluates the second derivative at construction and declines a model that unmixes, naming the interval — and OptimaSolver's certificate applies the tangent-plane test to phases that are present, not only to absent ones. Assuming convexity and leaving the duplication to the caller is a reasonable design for a general-purpose code; checking it is worth the few lines here, because a phase sitting inside its own spinodal otherwise certifies on the stationarity of its members alone and says nothing.

Why a formulation with one amount per species cannot hold it

The composition vector of a ChemicalSystem has one entry per species, so a species belongs to exactly one phase and a phase has exactly one composition. A common-tangent pair is two compositions of one substance, and there is no way to write that down with one amount each.

So the formulation has to be given the substance twice. That is what SolidSolutionPhase(...; instances = 2) does: ChemicalSystem builds a second copy of each end-member under a derived symbol (monosulphate12#2) sharing the same thermodynamic record, and each copy belongs to its own phase. The groups stay disjoint, the duplicated columns of the conservation matrix are copies of columns already there — so the row rank, and with it conservation, is untouched — and the minimization is free to populate either lobe or both.

It is refused for a convex model, and the reason is worth stating because it is not conservatism. Two instances of a convex phase are degenerate: the energy is the same however the amount is split between them, so the minimum becomes a flat manifold and the optimizer is asked to pick a point on it arbitrarily. Inside a spinodal the common-tangent pair is unique and no such direction exists.

Michelsen's tangent-plane distance, which is what detects it

Representation is one problem; knowing that it is needed is another. A phase sitting inside its own spinodal satisfies every first-order condition — its end-members are stationary, the element balance closes, nothing absent is supersaturated — so stationarity alone certifies a non-minimum.

The test that separates them is Michelsen's: for a trial composition of the phase, the tangent-plane distance

measures how far the Gibbs surface at lies above the tangent plane drawn at the reported composition . If is negative anywhere, some other composition lies below that plane and the reported state is not a minimum. At   the distance is exactly zero and is a stationary point of the search, which is why the measure has to be started from the end-member corners — the other lobe is where the negative value lives.

OptimaSolver's phase_split_measure computes this for every present mole-fraction phase and folds the result into the certificate. Concretely: check_convexity = false stops being a silent loss of the proof — a concave declaration that does unmix now fails to certify, with the phase named.

The same criterion, in another code

GEM-Selektor's PhaseSelection computes, for every phase, a stability index

the activity from the dual solution divided by the activity coefficient from the primal one — an estimate of the mole fraction. That is term for term what this package computes, in _repair_start's   (the ideal case  ) and in phase_split_measure's log-sum-exp of    . The GEMS3K paper calls "a generalization of the saturation index" and notes that it derives from the KKT conditions (Kulik et al., 2013).

So the criterion is one object, arrived at independently. What differs is only the declaration: CEMDATA18 ships the AFm and AFt binaries under two names each, so a GEMS user represents a gap by declaring the binary twice in the database; instances = 2 asks for the same thing with a keyword. Neither code splits a phase by itself.

The executed counterpart of this section is the miscibility-gap page, which runs one cement three ways and reports the certificate each time.

7. How a solid solution is declared

A SolidSolutionPhase names its end-members and carries a model:

phase = SolidSolutionPhase("CSHQ", [sp["CSHQ-JenD"], sp["CSHQ-JenH"],
                                    sp["CSHQ-TobD"], sp["CSHQ-TobH"]];
                           model = IdealSolidSolutionModel())

Every end-member must be AS_CRYSTAL, and the phase is checked against its model's arity and against the convexity of its mixing energy at construction; instances is how a declaration asks for more than one coexisting composition, which section 6 is about. Sets of end-members with their interaction parameters are read from a TOML through build_solid_solutions, which is also the only interaction-parameter file format in the package — data/solid_solutions.toml is the shipped example.

Two consequences to keep in mind when reading results:

• a ChemicalSystem groups end-members into ss_groups, and every aqueous activity model must call the solid-solution branch or the end-members come back with  , i.e. treated as pure phases;

• an end-member's saturation_indices entry is relative to its current mole fraction, since its activity is . For an end-member sitting at the solver's lower bound that is a statement about a vanishing phase, not about whether the solid solution would form.

See also

• Thermochemistry — where enters

• Activity models — the aqueous half

• Chemical Equilibrium — declaring and solving with solid solutions