"""Critical temperature (mean-field and Tyablikov RPA) from J(q)."""
from __future__ import annotations
import logging
from typing import Dict
import numpy as np
from mag4.constants import KB_MEV
from mag4.magnon.exchange import build_jq_matrices
log = logging.getLogger(__name__)
[docs]
def compute_tc(evals, S, lam0) -> Dict:
"""Mean-field and Luttinger–Tisza Tyablikov RPA critical temperature.
Both estimates use the Hermitian J(q) eigenvalue spectrum sampled on
the BZ grid.
Returns a dict with keys ``Tc_MF, Tc_RPA, note``.
"""
prefactor = (S + 1) / (3.0 * S)
if lam0 >= 0:
return dict(Tc_MF=0.0, Tc_RPA=0.0, note="lam0 >= 0: no stable order.")
Tc_MF = -prefactor * lam0 / KB_MEV
diffs = lam0 - evals.ravel()
safe = np.abs(diffs) > 1e-10
Tc_RPA = Tc_MF if safe.sum() == 0 else \
-prefactor / (np.mean(1.0/diffs[safe]) * KB_MEV)
return dict(Tc_MF=Tc_MF, Tc_RPA=Tc_RPA, note="")
[docs]
def compute_tc_lswt(bond_data: Dict,
S: float,
ordering_lt: Dict,
q_cart: np.ndarray,
evals_lt: np.ndarray) -> Dict:
"""Method-2 T_c, chosen based on q₀.
q₀ ≈ Γ (FM / type-I AFM)
LSWT is exact in the primitive cell. Build :math:`D^{ord}` with
:math:`s_α = \\mathrm{sign}(v_α)` and apply the Tyablikov RPA to the
acoustic (lowest) magnon branch::
k_B T_c^{LSWT-RPA} = (S+1)/(3S) / ⟨1/ω_0(q)⟩_{ω_0 > 0}
q₀ ≠ Γ (commensurate or incommensurate ordering)
:math:`D^{ord}` in the primitive cell has negative eigenvalues
because :math:`s_α = \\mathrm{sign}(\\cos 2π q₀·r_α)` is not the
exact collinear ground state for q₀ ≠ Γ. Use the LT acoustic-branch
RPA — Tyablikov RPA applied to ``λ_min(q)`` directly::
k_B T_c^{LT-ac-RPA} = (S+1)/(3S) / ⟨1/(λ_min(q₀) − λ_min(q))⟩_{q ≠ q₀}
"""
prefactor = (S + 1) / (3.0 * S)
q0_norm = float(np.linalg.norm(ordering_lt["q0_cart"]))
lam0 = ordering_lt["lam0"]
if q0_norm < 0.01:
spin_dirs = np.sign(ordering_lt["eigvec"])
spin_dirs[spin_dirs == 0.0] = 1.0
log.info(" LSWT T_c (q0≈Γ): spin dirs — %d up, %d down",
int(np.sum(spin_dirs > 0)), int(np.sum(spin_dirs < 0)))
s = spin_dirs
J0_raw = build_jq_matrices(bond_data, np.zeros((1, 3)))[0]
J0 = 0.5 * (J0_raw + J0_raw.T.conj())
onsite = s * (J0.real @ s)
ss = np.outer(s, s)
Jq_raw = build_jq_matrices(bond_data, q_cart)
Jq = 0.5 * (Jq_raw + Jq_raw.swapaxes(1, 2).conj())
N_q = len(q_cart)
eps_tol = 1e-6 # meV — exclude only exact Goldstone zero
acoustic = np.zeros(N_q)
for iq in range(N_q):
Dh = np.diag(onsite) - Jq[iq] * ss
Dh = 0.5 * (Dh + Dh.T.conj())
acoustic[iq] = 2.0 * S * float(np.linalg.eigvalsh(Dh)[0])
acoustic = np.clip(acoustic, 0.0, None)
mask = acoustic > eps_tol
om_acoustic = acoustic[mask]
if len(om_acoustic) == 0:
return dict(Tc_RPA=0.0, method="LSWT-RPA",
note="All acoustic modes ≤ ε_tol — "
"check spin directions.")
om_min = float(om_acoustic.min())
om_mean = float(om_acoustic.mean())
log.info(" LSWT acoustic: ω_min=%.6f meV ω_mean=%.4f meV "
"(%d Goldstone zeros excluded)",
om_min, om_mean, int(np.sum(~mask)))
Tc_RPA = prefactor / (float(np.mean(1.0 / om_acoustic)) * KB_MEV)
return dict(Tc_RPA=Tc_RPA, method="LSWT-RPA", note="",
om_min=om_min, om_mean=om_mean)
else:
eps_ac = evals_lt[:, -1]
diffs = lam0 - eps_ac
safe = diffs < -1e-10
n_excl = int(np.sum(~safe))
if safe.sum() == 0:
return dict(Tc_RPA=0.0, method="LT-ac-RPA",
note="No q-points away from q₀ found.")
Tc_RPA = -prefactor / (float(np.mean(1.0 / diffs[safe])) * KB_MEV)
log.info(" LT acoustic-branch RPA: %d q-points used "
"(%d near-q₀ excluded)",
int(safe.sum()), n_excl)
log.info(" Acoustic branch: λ_min range [%.4f, %.4f] meV.S2",
float(eps_ac.min()), float(eps_ac.max()))
return dict(Tc_RPA=Tc_RPA, method="LT-ac-RPA", note="",
lam_min=float(eps_ac.min()),
lam_mean=float(eps_ac.mean()))