Mean-Field Self-Consistency

A self-consistent mean-field calculation solves a nonlinear equation in which an effective Hamiltonian depends on the state it produces. Iteration is only one possible solver for that equation. A physically meaningful stationary solution can repel a naive iteration, while a converged numerical branch need not minimize the relevant energy or free energy.

Helpful background. Hartree–Fock computation and BdG numerics supply concrete matrix examples.

Fixed points and linearized iteration

Collect densities, fields, and any constrained variables into $x$. A direct self-consistent cycle defines $x_{\rm out}=F(x)$, and the desired equation is

$$R(x)=F(x)-x=0.$$

Near a fixed point $x_*$, an error evolves under direct iteration as $\delta x_{n+1}\approx J\delta x_n$, where $J=F'(x_*)$. For a continuously differentiable map, spectral radius $\rho(J)<1$ guarantees local attraction, although a nonnormal $J$ can produce transient growth. A value $\rho(J)>1$ gives an unstable direction. At $\rho(J)=1$, the linear test is inconclusive: higher-order terms can produce slow convergence or divergence. For example, $F(x)=x-x^3$ attracts sufficiently small nonzero $x$ toward zero despite $F'(0)=1$.

Linear mixing replaces direct iteration by

$$x_{n+1}=x_n+\alpha R(x_n),$$

whose linearized eigenvalues are $1-\alpha+\alpha\lambda_j(J)$. Reducing $\alpha$ can stabilize negative or oscillatory modes, but cannot fix every unstable eigenvalue. In particular, a real $\lambda_j>1$ remains above one for every positive $\alpha$.

Preconditioning acts on residual directions with different physical scales. Long-wavelength charge fluctuations in a large Coulomb system can converge much more slowly than local fluctuations. A scalar damping factor treats both equally and can be inefficient.

A scalar example that isolates mixing

For $F(x)=b-2x$, the fixed point is $x_*=b/3$. Direct iteration multiplies the error by $-2$ and diverges. With mixing,

$$\delta x_{n+1}=(1-3\alpha)\delta x_n.$$

Convergence requires $0<\alpha<2/3$; at $\alpha=1/3$, this linear problem reaches the fixed point in one step. The example demonstrates an algorithmic instability without any ambiguity about the unique solution of the equation.

A Newton step instead solves $[F'(x)-I]\delta x=-R(x)$. Approximate Jacobians, quasi-Newton methods, and DIIS can accelerate more realistic systems. DIIS chooses a linear combination of prior iterates or effective Hamiltonians by minimizing their residual combination subject to coefficients summing to one. A nearly dependent residual history needs regularization or pruning.

Constraints and physical stability

Enforce particle number, Hermiticity, normalization, and any imposed symmetries in the representation. If a chemical potential is adjusted to fix number, include the number residual in the convergence test and define its relative weight. A small field residual with a wrong particle number is not a solved constrained problem.

For the illustrative mean-field equation

$$m=\tanh[\beta(Jm+h)],$$

with $J>0$ and $\beta>0$, the zero-field solution $m=0$ has fixed-point slope $\beta J$. For $\beta J<1$, it attracts direct iteration; above one, two symmetry-broken solutions appear. In this special equilibrium model, the bifurcation also has thermodynamic meaning. That coincidence does not make iteration stability a general substitute for an energy Hessian or free-energy comparison.

Run different initial states, scan parameters in both directions, and compare stationary branches using the correct thermodynamic potential. A hysteresis loop may reflect metastability, a first-order transition, or failure to converge; distinguish these with residuals and independent starts.

Stopping and error allocation

Monitor a norm of $R$, change in the objective, constraint errors, and target observables. Scale residual components to avoid letting large units dominate. Tighten the inner eigenproblem as outer convergence improves; oversolving early steps wastes work, while loose final inner solves set a residual floor.

After convergence, remove artificial level shifts or continuation fields unless they define the intended model, and solve again. Refine the basis and integration grid separately. Physical mean-field error persists even for an exact root of the represented equation.

Exercises

For a fixed-point Jacobian eigenvalue $\lambda=-5$, determine the allowed positive mixing interval.

Solution

The mixed eigenvalue is $1-\alpha+\alpha(-5)=1-6\alpha$. Requiring its magnitude below one gives $0<\alpha<1/3$. This condition applies to that mode; every other mode must also satisfy its own stability condition.

References

• Nocedal, J., and S. J. Wright, Numerical Optimization, 2nd ed., Springer, 2006, doi:10.1007/978-0-387-40065-5 — nonlinear equation and optimization methods.
• Pulay, P., “Convergence acceleration of iterative sequences. The case of SCF iteration,” Chemical Physics Letters 73, 393–398 (1980), doi:10.1016/0009-2614(80)80396-4 — DIIS acceleration.