Vestigial chirality from fluctuating loop currents on the kagome lattice

Introduction

The kagome metals AV$_3$Sb$_5$, where $A$ represents an alkali metal element, have been discovered as a family of topological metals with superconducting ground states [1, 2, 3], combining a commensurate $2 \times2$ charge-density-wave (CDW) transition, strong precursor fluctuations, and reported signatures of broken time-reversal symmetry [4, 5, 6, 7, 8, 9]. Microscopic studies connect the competing multi-$Q$ particle–hole channels to van Hove singularities and sublattice interference [10, 11, 12, 13, 14]. Complementary low-energy classifications and fluctuation-mediated theories show how nonlocal interactions can stabilize charge-bond and loop-current orders [15, 16, 17]. Here loop-current order denotes the imaginary, time-reversal-odd part of a modulated electronic hopping. The candidate loop-current state contains three symmetry-related $M$-point components and therefore breaks both translation and time-reversal symmetry. This multicomponent structure poses a question about the finite-temperature phase structure that cannot be resolved by identifying the leading electronic instability alone: must the two symmetries be restored at the same transition upon heating?

Recent experiments make a separated thermal hierarchy especially timely. Circular-dichroism angle-resolved photoemission on CsV$_3$Sb$_5$ finds a time-reversal-odd response below $T^* \simeq145$–150 K, above the conventional CDW transition near 94 K; it is interpreted as loop-current order that precedes the inverse-Star-of-David real charge-bond reconstruction [18]. A depth-resolved $\mu$SR study likewise reports a near-surface magnetic response in RbV$_3$Sb$_5$ beginning near 175 K, above the bulk charge-order temperature near 110 K [9]. Here “conventional CDW order” denotes the observable real charge/bond and lattice reconstruction; the relative weights of site charge, bond modulation, and atomic displacement depend on the material and probe [19, 20]. It should not be identified automatically with the time-reversal-odd loop-current component. These observations are therefore compatible not only with primary loop-current order preceding the conventional CDW, but also with a vestigial chiral regime in which a uniform time-reversal-odd composite remains ordered after the primary translation-breaking loop-current fields become short ranged. They do not by themselves distinguish these possibilities: the loop-current model used to interpret the photoemission data already doubles the unit cell size, and absence of the real structural reconstruction does not establish restored translations in the loop-current sector.

If the primary translation-breaking field disorders while a uniform composite remains ordered, the intervening state is vestigial [21]. Establishing such a phase requires more than observing a nonzero composite moment on a finite lattice. One must demonstrate long-range order in a symmetry-nontrivial composite and short-ranged primary correlations in the thermodynamic limit. Tsvelik and Sarkar proposed a related phase-only theory that identified a fluctuation-driven composite CDW crossover [22]. Subsequent Monte Carlo simulations of the compact three-phase model found one continuous transition of the primary fields together with that crossover [23]. The phase-locked composite in these works is time-reversal even, so the crossover is not associated with any symmetry breaking. These studies therefore do not establish a thermodynamic phase with long-range time-reversal-odd order and short-ranged primary correlations. In the strongly commensurate limit studied here, the asymptotic phasons are pinned and melting is instead controlled by domain walls [24].

Low-energy and Landau model analyses have classified the relevant three-$Q$ charge, bond, and flux orders [11, 15, 16, 20]; the symmetry-complete homogeneous bond–flux theory is given in Refs. [25, 26]. This classification selects candidate ordered states, but homogeneous theory alone cannot decide how their spatially varying configurations melt. The essential energetic question is whether a domain wall that reverses chirality can remain sufficiently more costly than a wall that changes the translation domain while preserving chirality; see Figure 1.

Image not extracted

We answer this question in two stages. First, we derive an effective fixed-amplitude sign model for the three-$Q$ loop current order, which allows cluster simulations through a linear system size $L = 384$, maps the global phase topology, and converts its multicritical endpoint into a ratio of microscopic wall tensions. Second, a stable amplitude-resolved lattice Ginzburg–Landau (GL) model tests whether the two-stage sequence survives wall-core relaxation. We find that a vestigial chiral phase emerges when a chirality-changing wall is only 12–13% more costly than the same-chirality wall. This phase also survives amplitude fluctuations in a controlled soft-spin parameter family. This is a sufficiency result for the phenomenological theory, not yet a material-specific prediction: a microscopic calculation of the two relaxed wall tensions remains necessary.

Results

Three-$Q$ order and composite chirality

The kagome Bravais lattice is triangular. We place one coarse-grained field $\boldsymbol{\Phi}$ on each kagome unit cell at $\boldsymbol{r}_i = x_i\boldsymbol{a}_1 + y_i\boldsymbol{a}_2$, where $\boldsymbol{a}_1$, $\boldsymbol{a}_2$ are the primitive lattice vectors and $x_i$, $y_i$ are integer coordinates. Its three internal components $\Phi_{ai}$, $a = 1, 2, 3$, are real amplitudes of the symmetry-related $M$-point loop-current patterns; they are not additional spatial sites. Microscopic sublattice and bond structure enters through the transformation laws and current form factors of the components. The triangular lattice is therefore a symmetry-compatible coarse graining of kagome unit cells, not a replacement of the microscopic kagome geometry.

With a convenient ordering at the $M$ points, primitive translations act as

while a threefold rotation cycles the component labels. Time reversal acts as $\mathcal{T}: \boldsymbol{\Phi} \mapsto-\boldsymbol{\Phi}$. Consequently,

is translation invariant and time-reversal odd. Long-range order of $\chi$ with short-ranged correlations of every $\Phi_a$ therefore defines a uniform vestigial chiral phase. It preserves the translations restored by melting the $M$-point order while continuing to break $\mathcal{T}$.

Amplitude-resolved effective Hamiltonian

Real bond and loop-current orders are naturally intertwined because they are the real and imaginary parts of modulated electronic hopping. They generally transform as different irreducible representations, however, and need not be simultaneously critical. Denoting the bond field by $\boldsymbol{B}$, the bond and loop-current combination permits the cubic coupling $-g(B_1\Phi_2\Phi_3 + B_2\Phi_3\Phi_1 + B_3\Phi_1\Phi_2)$ [25, 26]. We assume that the bond field is massive on the scales of interest and integrate it out rather than retain it as an independent Monte Carlo field.

To see the resulting terms, define $X_{1i} = \Phi_{2i}\Phi_{3i}$ and its cyclic permutations. If $A_B$ is the positive quadratic kernel of the massive bond field, its leading contribution to the Hamiltonian is

and integrating out $\boldsymbol{B}$ generates

In the local massive limit, $A_B^{-1} \simeq r_B^{-1}\delta_{ij}$, where $r_B > 0$ is the local bond-field mass, Eq. (4) becomes $-g^2(2r_B)^{-1}\sum_{i,a<b}\Phi_{ai}^2\Phi_{bi}^2$, and is absorbed into the quartic coefficients below. The allowed bond-field cubic invariant generates a term proportional to $\chi_i\chi_j$ at higher order. Thus eliminating a noncritical bond field preserves the loop-current symmetries and motivates, but does not microscopically determine, the effective couplings.

Our production Hamiltonian is consequently the three-field theory

Here $\Phi_i^2 = \sum_a\Phi_{ai}^2$, and the coefficients already include the short-distance renormalizations generated by the eliminated bond field. The interaction proportional to $J_{\chi}$ is sixth order in the primary fields and allowed by every symmetry. It directly changes the relative cost of chirality-changing and translation-related domain walls (Figure 1). We interpret it as an effective interaction generated after shorter-scale electronic, magnetic, or collective modes are integrated out, rather than as a coefficient determined by the homogeneous quartic theory. The on-site sextic term bounds the soft-spin energy. For the triangular bond counting used here, a sufficient strict condition is $c_{\Phi} \ge0$, $|J_{\chi}| < 3c_{\Phi}/2$.

Fixed-amplitude limit and domain wall energetics

Deep in an equal-amplitude three-$Q$ state we write $\Phi_{ai} = A s_{ai}$ with $s_{ai} = \pm1$. Up to constants, Eq. (5) reduces to

with

We denote by $T_s$ the transition at which the primary component signs lose long-range order upon heating and by $T_{\tau}$ the transition at which the composite chirality disorders. The primary phase orders both $s_a$ and $\tau$, and hence breaks translations and time reversal. In the possible vestigial phase, $\tau$ remains ordered while $\langle s_a(r)s_a(0)\rangle$ is short ranged for every component.

The limiting cases provide exact checks. At $q = 0$, the three components are independent ferromagnetic Ising models on a triangular lattice, with a critical temperature $T_c(q = 0) = 4K_1/\ln3$ [27]. For $q \to\infty$, $\tau$ orders at the triangular-Ising scale

Within a fixed-$\tau$ sector the four allowed component states form a ferromagnetic four-state Potts model whose bond contrast is $4K_1$, giving a separate critical temperature [28]

The sign limit turns the phase-topology threshold into a wall-energy target; see Figure 1. A straight same-chirality translation wall flips two component signs. If it cuts $\nu_{\hat n}$ triangular bonds per unit length, its tension is

The cheapest chirality-changing wall flips one component and has tension

For equal orientations,

Thus the numerical endpoint found below, $1.24 \lesssim K_3/K_1 \lesssim1.26$, corresponds to $\sigma_{\mathrm{odd}}/\sigma_{\mathrm{even}} \simeq1.12$–$1.13$. Although the coupling ratio is of order unity, the required wall-tension hierarchy is modest because an odd wall already pays half of the two-component-translation-wall cost at $K_3 = 0$.

We determine the finite-temperature behavior using large-scale cluster and parallel-tempering Monte Carlo simulations on periodic $L \times L$ triangular lattices with $N = L^2$ sites. The sign model is simulated through $L = 384$, while the amplitude-resolved validation reaches $L = 192$. Phase boundaries are located from crossings of dimensionless finite-size ratios; structure-factor scaling and real-space correlations independently test whether the primary fields are short ranged. The update schemes, observable definitions, and uncertainty analysis are given in Appendix Secs. 1–3; run-specific diagnostics and complete transition sequences are collected in the Supplementary Information.

Fixed-amplitude phase diagram

Figure 2(a) summarizes the phase diagram obtained from lattices as large as $L = 384$. At small $q$, the primary and composite orders disappear together. At large $q$, heating first disorders the component signs at $T_s$ while preserving $\langle\tau\rangle\ne0$; time reversal is restored only at a second transition $T_\tau$. The boundaries approach the exact Ising and Potts limits in Eqs. (8) and (9).

Image not extracted

At $q = 2$, the crossings converge to $T_s/K_1 = 5.7131(7)$ and $T_\tau/K_1 = 7.4592(8)$ [Figure 2(b,Figure 2c)]. Their wide separation establishes a robust vestigial interval rather than a finite-size shoulder. Binder ratios provide an independent consistency check, while the predominantly single-peaked energy distributions show no stable coexistence signal at either boundary. The finite-size estimators are defined in Appendix Sec. 2, their joint-bootstrap analysis is described in Appendix Sec. 3, and the full $q = 2$ production record and Binder curves appear in Supplementary Secs. II A and IV.

Near the point where the two boundaries separate, the splitting is unresolved at $q = 1.24$ and first resolved at $q = 1.26$. We therefore obtain the conservative endpoint corridor

The split then widens rapidly with increasing $q$. Resolving the order of the direct transition is a separate question and does not affect this phase topology. The full phase-boundary estimates are tabulated in Supplementary Table S1; the endpoint criterion, tempering diagnostics, and histogram checks are given in Supplementary Sec. II B.

$q$	$T_s$	$T_\tau$	interpretation
0.00	3.634(20)	3.634(20)	direct
1.00	5.080(5)	5.080(5)	direct
1.24	5.382(3)	5.387(3)	unresolved
1.26	5.403(3)	5.414(3)	first resolved
1.28	5.420(3)	5.445(3)	vestigial
2.00	5.7131(7)	7.4592(8)	vestigial
3.00	5.751(20)	10.970(6)	vestigial

Thermal and defect signatures

The specific heat $c_V$ in Figure 3(a) exhibits a dominant sharp peak at $T_s$ and a much weaker anomaly near $T_\tau$. The lower feature accompanies the loss of the three component orders [Figure 3(b)]. The upper feature marks the loss of the composite order and the restoration of time reversal. Its small amplitude explains why a single-size specific-heat curve does not by itself resolve the two transitions reliably; we therefore locate both boundaries from finite-size correlation-ratio crossings. The high-temperature tails in Figure 3(b) decrease with system size as expected for a disordered two-dimensional phase and therefore represent finite-size susceptibility backgrounds, not residual order. The energy-fluctuation definition of $c_V$ is given in Appendix Sec. 2, and the finite-size tail analysis is given in Supplementary Sec. II C.

Image not extracted

The domain wall densities in Figure 3(c) reveal the underlying melting mechanism. At $T_s$, the translation wall density $\rho_s$ rises rapidly while the chirality wall density $\rho_\tau$ remains dilute. The latter grows strongly only near $T_\tau$. The intermediate phase is therefore a plasma of same-chirality translation walls in which time-reversal-changing walls remain suppressed. The wall observables are defined in Appendix Sec. 2.

Real-space characterization

Real-space correlations provide the most direct distinction among the three phases. Figure 4 uses the broad $q=2$ window so that the narrow endpoint corridor does not obscure the limiting behaviors. Both the $s$ and $\tau$ correlation functions $C_s(r)$ and $C_\tau(r)$ approach plateaus in the fully ordered phase, only $C_\tau(r)$ plateaus in the vestigial phase, and both are short ranged above $T_\tau$. The correlation estimator is defined in Appendix Sec. 2, and its sampling details are reported in Supplementary Sec. II A.

Image not extracted

Robustness to amplitude fluctuations

The sign model freezes $|\Phi_a|$. We therefore restore local amplitude fluctuations in Eq. (5) and follow a stable, relatively stiff soft-spin family from the direct-transition regime into the two-transition regime. This family is designed as a controlled bridge from the sign limit, not as a material-specific fit. The parameter family and stability checks are given in Supplementary Sec. III A. The update scheme is described in Appendix Sec. 1, while the production campaign and convergence diagnostics are reported in Supplementary Sec. III B.

The $J_{\chi}/K_{\Phi}=4.75$ case lies on the direct-transition side, while the small separation at $J_{\chi}/K_{\Phi}=5.0$ remains unresolved. At $J_{\chi}/K_{\Phi}=5.5$, the largest size pair gives $T_{\Phi}/K_{\Phi}=2.76918(25)$, $T_{\chi}/K_{\Phi}=2.81310(64)$, and $(T_{\chi}-T_{\Phi})/K_{\Phi}=0.04392(65)$ [Figure 5(a,Figure 5b)]. Agreement among successive large sizes shows that the interval survives amplitude relaxation.

Image not extracted

Separated crossings identify two thermal scales but do not by themselves determine the order between them. Figure 5(c) supplies that test, depicting the normalized zero-momentum structure factors $S_X(0)/L^2$ for the two fields $X\in\{\Phi,\chi\}$. At $T/K_{\Phi}=2.79121$ inside the interval, $S_{\Phi}(0)/L^2$ approaches its short-range $L^{-2}$ scaling with increasing size, whereas $S_{\chi}(0)/L^2$ approaches a nonzero plateau. The amplitude-resolved model therefore has short-ranged primary loop-current correlations but long-range composite chirality. The second-moment estimator and phase-assignment tests are defined in Appendix Sec. 2; the complete crossing sequence and scaling values are given in Supplementary Tables S2–S3 and Sec. III C, with the uncertainty procedure described in Appendix Sec. 3.

Connection to the fixed-amplitude regime

Near the transitions, the measured root-mean-square component amplitude is $A_{\mathrm{rms}}\simeq0.714$. Substitution into Eq. (7) gives $q_{\chi}^{\mathrm{eff}}\simeq(1.23,1.30,1.43)$ for $J_{\chi}/K_{\Phi}=(4.75,5.0,5.5)$. The first and last values fall below the rigid-amplitude endpoint $q_c\simeq1.25$, consistent with their direct and two-transition behaviors. The middle value maps above the endpoint although its thermodynamic splitting remains unresolved, showing that the bulk estimate is not a quantitative threshold.

Wall-core relaxation explains the discrepancy: local suppression of $A$ weakens the chiral contribution more rapidly than the primary stiffness, making chirality reversal easier than the uniform-amplitude estimate suggests. The soft-spin onset is consequently shifted to larger bulk $q_{\chi}^{\mathrm{eff}}$, and fully relaxed even- and odd-wall tensions are the appropriate matching variables. Numerical inputs are given in Supplementary Sec. III D.

Discussion

Symmetry distinction and material criterion

The sign and soft-spin models establish two-stage melting for their specified classical Hamiltonians. Although the numerical location of the endpoint is nonuniversal, the order distinction is exact: the intermediate phase restores the translations broken by the $M$-point fields while retaining a uniform Ising variable that breaks time reversal. The consistent correlation-ratio, structure-factor, and real-space tests distinguish this state from a crossover with a large but finite chiral correlation length.

The attractive $J_{\chi}$ interaction is a symmetry-allowed phenomenological extension. Because the local model with $J_{\chi}/K_{\Phi}=0$ did not resolve the same phase, our result does not imply that it occurs in every kagome loop-current system. The wall mapping in Eq. (12) states the material requirement more directly than the bare coupling: the lowest-energy odd wall must be about 12–13% more costly than the lowest-energy even wall in the rigid limit. Amplitude relaxation changes the numerical threshold but not the matching principle.

A concurrent preprint by Iqbal proposed, on symmetry grounds, that vestigial chiral order is allowed in a doped kagome Dirac spin liquid [29]. This microscopic symmetry analysis is complementary to the present work, which establishes a separate finite-temperature vestigial chiral phase in a commensurate three-$Q$ loop-current model through finite-size scaling and relates its stability to the relative tensions of chirality and translation domain walls.

Microscopic routes to enhanced chirality-wall tension

The required hierarchy is modest enough to have several plausible microscopic origins. First, an odd wall interpolates between domains with opposite orbital magnetization and Berry curvature. Its fermionic reconstruction may suppress the particle–hole gap, redistribute more low-energy spectral weight, or bind wall states, making it more costly than a translation wall connecting domains with the same chirality. The difference can be amplified near the intertwined van Hove singularities that favor loop-current order [13, 14]. Wall-state entropy can partly offset this energetic penalty, however, so the relevant comparison is between finite-temperature free energies rather than zero-temperature band energies alone.

A second route is coupling to a soft time-reversal-odd collective mode, a possibility consistent with fluctuation-mediated loop-current mechanisms [17]. For an orbital magnetization, scalar-spin-chirality, or magnetic field $M_i$ with its Hamiltonian

where $r_M>0$ is the inverse local susceptibility, $J_M$ couples neighboring auxiliary modes, and $\lambda$ is the local $M\chi$ coupling, integrating out $M$ generates, to leading order, an attractive chirality interaction $J_{\chi}^{\mathrm{x}}\sim J_M\lambda^2/r_M^2$. It is enhanced when the auxiliary mode is close to its own instability and directly raises the cost of an odd wall. Microscopic current conservation and the sharing of kagome bonds can instead bind elementary component walls into even-parity pairs: an isolated one-component reversal is then penalized while a two-component translation wall remains comparatively inexpensive. Conversely, sublattice interference, frustrated gradient terms, or proximity to a competing ordering wavevector can soften the primary stiffness and lower the even-wall tension without removing the chirality penalty. Multicomponent ring processes around kagome triangles or hexagons provide another possible source of an effective chirality interaction, especially when flat bands or small strong-coupling denominators enhance higher-order virtual processes.

These mechanisms show why $\sigma_{\mathrm{odd}}>\sigma_{\mathrm{even}}$ is physically plausible, but they do not establish it for a particular compound. A static phonon or ordinary strain is time-reversal even and cannot align $\chi$ linearly, although it can modify amplitudes and primary stiffnesses. Likewise, a homogeneous sixth-order coefficient, or an estimate of $J_{\chi}$ alone, cannot capture the different wall-core relaxations. A material-specific test must compute both fully relaxed wall free energies in the same electronic theory and compare their ratio with the 1.12–1.13 target.

Relation to compact-phase composite order

Refs. [22, 23] study compact three-phase CDW models, including commensuration effects, in which continuous phase fluctuations and vortices govern melting. Their translation-invariant composite is time-reversal even in the proposed intermediate regime, so its development is a crossover rather than a separate symmetry-breaking transition. The Monte Carlo study of Ref. [23] accordingly found one thermodynamic transition and did not establish long-range time-reversal-odd composite order with short-ranged primary correlations.

By contrast, in our strongly commensurate model, the real loop-current amplitudes have discrete translation domains. Their product $\chi=\Phi_1\Phi_2\Phi_3$ is explicitly odd under time reversal, so $\langle\chi\rangle\ne0$ breaks an exact $\mathbb{Z}_2$ symmetry even when the primary $M$-point correlations are short ranged. The distinction is therefore the domain-wall mechanism and a separate thermodynamic vestigial chiral phase, established here by correlation-ratio crossings, structure-factor scaling, and real-space correlations.

Experimental implications

Uniform time-reversal-odd order permits a zero-field anomalous Hall response, Kerr rotation, and weak orbital magnetic fields, with their magnitudes set by the coupling of $\chi$ to the electronic bands [5, 6, 13, 14, 30, 31, 32]. Simultaneously, the primary $2 \times2$ loop-current Bragg correlations are short ranged. The characteristic thermal sequence is therefore a time-reversal-sensitive onset at $T_\chi$, followed only at lower temperature by sharp $M$-point order at $T_\Phi$. Diffuse $M$-point scattering and an enhanced chiral susceptibility should occupy the intermediate window. Opposite-$\chi$ domains can suppress bulk zero-field signals, making field training, hysteresis, and local magnetic probes important diagnostics.

Two recent observations are particularly suggestive of such a hierarchy. The circular-dichroism photoemission response of CsV$_3$Sb$_5$ onsets around 145–150 K, well above the conventional CDW transition, and was attributed to loop-current order [18]. Near-surface $\mu$SR in RbV$_3$Sb$_5$ similarly finds a static magnetic response above the long-range charge-order temperature [9]. Diffuse x-ray scattering also establishes quasi-static, short-range CDW precursors above $T_{\mathrm{CDW}}$ [7, 8]. Taken together, these results are compatible with a time-reversal-odd regime preceding long-range real charge-bond/lattice order. They are not yet a thermodynamic identification of the phase found here: structural CDW scattering does not directly measure the $M$-point loop-current correlator, and the photoemission interpretation permits primary translation-breaking loop-current order rather than only a uniform composite.

The experimental status of time-reversal breaking is also not uniform across probes. Kerr, chiral-transport, and $\mu$SR studies have reported time-reversal-odd responses [6, 9, 31, 32], whereas high-resolution zero-area-loop Sagnac measurements found no spontaneous polar Kerr signal in CsV$_3$Sb$_5$ within a 30 nrad noise floor [33]. A decisive test of vestigial chirality would therefore combine a local or domain-resolved time-reversal-sensitive probe with a wavevector-resolved probe of the primary loop-current correlations in the same sample and depth. The required intermediate regime has long-range $\chi$ order but no resolution-limited primary $M$-point peak.

These are symmetry-allowed signatures, not a prediction of quantized Hall transport. Whether vestigial order produces a Chern metal or insulator depends on the microscopic coupling between $\chi$ and the fermions. In particular, long-range order of the composite does not make the one-electron Hamiltonian equivalent to a static three-$Q$ texture with long-range $\Phi_a$ order.

Relation to vestigial superconductivity and outlook

The symmetry logic parallels charge-$4e$ and charge-$6e$ vestigial superconductors: a composite remains coherent after its root order disorders [34, 35, 36, 37, 38, 39, 40]. Let $\theta$ be the compact phase of an elementary charge-$2e$ pair field, $\Delta_{2e} \propto e^{i\theta}$. Its $n$th-power composite $\Delta_{2ne} \sim(\Delta_{2e})^n \propto e^{in\theta}$ carries charge $2ne$: $n = 2$ and $3$ give charge-$4e$ and charge-$6e$ order, respectively. In a vestigial phase, correlations of $e^{i\theta}$ are short ranged while $e^{in\theta}$ retains phase coherence. Because $\theta$ is compact, the defect spectrum can contain $2\pi/n$ fractional vortices, often tied to discrete defects, whose binding and proliferation help control the transitions.

Our defect mechanism is different. Strong commensuration locks the phasons, leaves real component fields, and replaces the compact-phase defects by even- and odd-parity domain walls. The wall-tension ratio, rather than a fractional-vortex fugacity, therefore controls the intermediate phase. The triangular lattice used here is the Bravais regulator of kagome unit cells and should not be confused with the emergent kagome lattice of phase variables in the charge-$6e$ model of Ref. [38].

We have shown that commensurate three-$Q$ loop-current order can melt in two stages. In the sign model, a direct line gives way near $K_3/K_1 \simeq1.25$ to a phase that restores translations before time reversal. This endpoint requires only a 12–13% hierarchy between the lowest chirality-changing and same-chirality wall tensions. A stable five-size soft-spin calculation retains a nonzero vestigial interval after amplitude relaxation, and structure-factor scaling verifies that the primary loop-current field is short ranged within it. The remaining material question is therefore concrete: determine the two relaxed wall tensions microscopically or experimentally and compare their ratio with the criterion established here.

Note added

After this manuscript was posted on arXiv, we became aware of a concurrent preprint by Iqbal [29], posted shortly before ours, which proposed on symmetry grounds that vestigial chiral order is allowed in a doped kagome Dirac spin liquid. That work addresses a distinct microscopic setting and symmetry classification and is complementary to our finite-temperature demonstration of a thermodynamic vestigial chiral phase and its domain-wall criterion.

Acknowledgments

OpenAI Codex with GPT-5.6 Sol assisted with manuscript organization and language revision. The authors directed its use, checked the cited literature and scientific content, and reviewed and edited all generated text. The authors take full responsibility for the manuscript. This work was supported by the start-up grant from the Institute of Physics, Chinese Academy of Sciences, and CAS Project for Young Scientists in Basic Research (Grant No. YSBR143). S.M. acknowledges financial support from Ministry of Science and Technology (Grant No. 2021YFA1400200), National Natural Science Fund of China (Grants No.12450401 and No. 12025407), and Chinese Academy of Sciences (No. YSBR047).

Data availability

The numerical data supporting the findings of this study are available from the corresponding author upon reasonable request. The simulation inputs required to regenerate the results are included with the archived code described below. The custom code used for the Monte Carlo simulations, data analysis, and figure generation, together with the production input configurations, is openly available under the MIT License on Zenodo [41].

Methods

Monte Carlo sampling

All simulations use periodic $L \times L$ triangular Bravais lattices with $N = L^2$ sites. Angular brackets $\langle\cdots\rangle$ denote thermal Monte Carlo averages. We sample the sign model with embedded Swendsen–Wang cluster updates [42]. Fixed-temperature cluster scans map the broad phase diagram; parallel tempering [43] is used in the narrow direct-transition corridor. For each embedded Swendsen–Wang sub-update, the production GPU implementation batches all coupling points, independent chains, and temperature replicas, and processes their lattice bonds or sites in parallel. We benchmark this cluster updating algorithm against an independent local-update baseline. Writing $\mathbf{r}=x\mathbf{a}_1+y\mathbf{a}_2$, we assign colors to sites by $c(x,y)=(x-y)\bmod3\in\{0,1,2\}$. Every nearest-neighbor move changes $c$ by $\pm1$ modulo three, so no bond connects equal colors. Sites of one color can therefore be updated concurrently with their neighbors fixed; cycling through $c=0,1,2$ gives one conflict-free GPU local sweep.

The soft-spin model combines additive component proposals with two exact embedded-cluster updates. For a one-component reflection of component $a$, the amplitudes and the other two sign fields define a signed conditional Ising model with bond

Satisfied bonds of magnitude $\lvert J_{ij}^{(a)}\rvert$ are activated with probability $1-\exp[-2\lvert J_{ij}^{(a)}\rvert/T]$. This update changes $\chi$ and samples odd walls. A second cluster reflection changes two components together. Its conditional bond is

because every onsite invariant and $\chi_i$ are unchanged by the paired reflection. It therefore moves the primary translation walls without an intermediate chirality wall. Both clusters are needed: at $L=96$, adding the paired update reduces the longest measured autocorrelation time from 1623 to 27 sweeps in an otherwise identical pilot. Four independently seeded chains share a common parallel-tempering ladder.

Order and defect observables

The intensive specific heat plotted in Figure 3(a) is obtained from fluctuations of the extensive energy $E$ (with $k_B=1$),

For the order-parameter amplitudes in Figure 3(b), we first define the three component magnetizations and the scalar chirality magnetization by

The plotted quantities are the thermal root-mean-square amplitudes $\sqrt{\langle m_s^2\rangle}$ and $\sqrt{\langle m_\tau^2\rangle}$. This normalization makes both approach one in a uniform ordered domain, while the squares avoid selecting a particular member of the symmetry-related finite-size ensemble.

For the sign and soft-spin models, we define the common structure-factor notation

where $X\in\{s,\Phi\}$ denotes a three-component field with $n_X=3$, while $X\in\{\tau,\chi\}$ denotes a scalar field with $n_X=1$ and $X_{1i}=X_i$. Let $\mathbf{b}_1,\mathbf{b}_2$ be reciprocal primitive vectors, $\mathbf{a}_\mu\cdot\mathbf{b}_\nu=2\pi\delta_{\mu\nu}$. In the sign-model estimator, $S_X(\mathbf{k}_{\min})$ is shorthand for the average at $\mathbf{b}_1/L$ and $\mathbf{b}_2/L$. In the soft-spin estimator it is averaged over $\mathbf{b}_1/L$, $\mathbf{b}_2/L$, and $(\mathbf{b}_1+\mathbf{b}_2)/L; these are the shortest nonzero reciprocal-lattice directions used in the respective measurements.

For the sign model we use the dimensionless correlation ratio

Crossings of $R_s$ and $R_\tau$ locate the two phase boundaries. To expose their defect content, we also count the number of reversed primary signs on a nearest-neighbor bond,

Thus $h_{ij}=2$ identifies a translation wall that preserves $\tau$, whereas odd $h_{ij}$ reverses $\tau$. Independent random signs give $\rho_s = 3/8$ and $\rho_\tau= 1/2$. These local quantities diagnose wall proliferation but do not locate the transitions. We additionally compute the Binder ratios [44],

The second expression is the scalar convention of Ref. [44]; the factor $5/3$ in the first is the Gaussian fourth-moment ratio for a three-component vector; here $m_s^4 \equiv(m_s^2)^2$. Accordingly, $U_s \to3/5$ and $U_\tau\to1/3$ in their rigid ordered limits, while both approach one in a Gaussian disordered phase. Near the $q = 2$ primary transition, $U_s(T)$ is nonmonotonic and has two successive-size intersections. We therefore use the unique $R_s$ crossing in the main analysis and retain Binder results as Supplementary checks. The monotonic $U_\tau$ curves give a unique crossing consistent with $R_\tau$.

For the real-space diagnostic in Figure 4, we measure

where $r$ is an integer separation and $\hat{e}_d \in\{\mathbf{a}_1, \mathbf{a}_2, \mathbf{a}_1-\mathbf{a}_2\}$ runs over the three positive triangular-lattice directions; the subscripts on $\langle\cdots\rangle$ denote uniform averages over the listed site, component, and direction indices. Both correlators equal one at $r = 0$. For the soft-spin model, the corresponding dimensionless observable is the second-moment ratio

Here $X \in\{\Phi,\chi\}$ and $\xi_X$ is the second-moment correlation length in channel $X$. Successive-size crossings locate the transitions; Binder ratios, susceptibilities, and structure-factor scaling test the resulting phase assignment. No boundary is inferred from a single-size peak.

Convergence and statistical analysis

Time series are blocked at a scale chosen to exceed the integrated autocorrelation time. Bootstrap replicas first resample independent chains and then whole blocks within each selected chain. For parallel-tempering data, each block contains the synchronized measurements from the complete temperature ladder; the same resampled block indices are therefore applied to all temperature slots and primitive observables before nonlinear ratios and crossing temperatures are reconstructed. Fixed-temperature scans are instead resampled independently at each temperature. Different lattice sizes are resampled independently, while the bootstrap curve for a given size is reused in adjacent size-pair crossings. We monitor split-$\hat{R}$ [45], effective sample size, reliable block count, and replica round trips as complementary convergence tests. Near a transport-limited sign-model crossing, the uncertainty is floored at half the native temperature-grid spacing. Ferrenberg–Swendsen/WHAM reweighting of energy histograms, where WHAM denotes the weighted histogram analysis method, tests for phase coexistence [46].

Use of artificial intelligence tools

OpenAI Codex with GPT-5.6 Sol was used interactively to assist with software implementation and debugging, data-analysis workflows, and plotting scripts. The authors specified the physical models, algorithms, and numerical protocols and verified the AI-assisted output through manual code review, automated unit and regression tests, cross-comparisons between independent CPU and Metal/GPU implementations, and equilibrium comparisons between embedded Swendsen–Wang cluster and independent single-site heat-bath samplers. The simulations were additionally benchmarked against the exactly known triangular-lattice Ising and four-state Potts limiting cases. All simulations were executed and all reported numerical results and figures were inspected by the authors.

References

supplementary.pdf

Supplementary Information for ``Vestigial chirality from fluctuating loop currents on the kagome lattice''

Transition summaries

Tables S1 and S2 collect the numerical transition estimates quoted in the main text. They are separated from the main article to keep its Results focused on the phase topology and physical mechanism.

$J_\chi/K_\Phi$	12–18	18–24	24–30	30–36	conclusion
4.75	0.00028(83)	0.00370(130)	−0.00034(62)	0.00096(62)	direct/crossover
5.00	0.00318(112)	0.00370(86)	0.00349(121)	0.00188(72)	thermodynamic limit unresolved
5.50	0.02339(234)	0.02056(303)	0.02451(386)	0.02422(454)	small-size split

Sign-model numerical details

Refined transitions at $q=2$

For the lower boundary, eight independent parallel-tempering ladders were simulated at each of $L = 96, 192, 288, 384$ on a common grid $5.68 \le T/K_1 \le5.78$. At least 4096 measurements were retained per temperature and chain. All sizes satisfy split-$\hat{R} < 1.05$; the $L = 384$ ensemble records 15 aggregate cold–hot–cold round trips, with repeated bidirectional passages in every chain. Joint block bootstrap gives successive $R_s$ crossings 5.71298(21), 5.71226(25), and 5.71306(71). No stable two-peak energy histogram is present.

For the upper boundary, eight independent ladders at each size were run on a 33-slot grid $7.40 \le T/K_1 \le7.52$, retaining 8192 measurements per temperature and chain. Every primitive-moment series has split-$\hat{R} < 1.009$; the $L = 384$ ensemble records 2654 aggregate round trips. The successive $R_\tau$ crossings are 7.4580(7), 7.4588(7), and 7.4592(8). Binder crossings agree within uncertainty, and the energy histograms again show no two-peak candidate.

The real-space correlation data combine four chains and 16 blocks per chain. Four measurements per block are separated by four embedded-cluster sweeps. Uncertainties are the standard errors across chain means.

Endpoint criterion and histogram checks

The endpoint corridor is transport limited: replicas can accept many neighboring exchanges without traversing the full free-energy bottleneck. We therefore classify a split as resolved only when

Here $\sigma_s$ and $\sigma_\tau$ are the one-standard-error uncertainties of $T_s$ and $T_\tau$, respectively. By this criterion, the split is unresolved at $q = 1.24$ and first resolved, at the native-grid precision, at $q = 1.26$. The entries in Table S1 use the independent $L = 192$–288 feedback-tempering campaign, which gives $(T_s,T_\tau)/K_1 = (5.3823, 5.3871)$ at $q = 1.24$ and $(5.4030, 5.4145)$ at $q = 1.26$. No large-size ladder completes a full cold–hot–cold round trip, and the $q = 1.24$ primary estimator has split-$\hat{R} = 1.144$. We therefore retain the half-grid uncertainty floor and do not quote a more precise value of $q_c$.

Energy-histogram scaling does not establish that the direct line is first order. The largest-size distributions are predominantly single peaked, and apparent equal-weight barriers at transport-limited points are not stable against size and chain checks. Determining whether the direct line is continuous or weakly first order is distinct from identifying where the two thermal boundaries separate.

Finite-size tails above the transitions

At $T/K_1 = 7.69 > T_\tau/K_1$, whole-chain bootstrap fits over $L = 48, 96, 192$ give

Both are consistent with the disordered two-dimensional law $L^{-1}$. At fixed $L = 96$, the chiral root-mean-square amplitude decreases from 0.2052(26) at $T/K_1 = 7.69$ to 0.09873(95) at 8.0 and 0.05990(58) at 8.5. The apparent high-temperature tails in the main-text thermal curves are therefore finite-size susceptibility backgrounds rather than residual order.

Amplitude-resolved numerical details

Parameter family and stability

We take $K_\Phi$ as the energy unit and restore amplitude fluctuations using

The common scale of the local coefficients is deliberately large:

These are the cross-component-basis coefficients of the main-text Hamiltonian. Equivalently, the isotropic-plus-cubic-anisotropy convention used by the simulation input has coefficients $\widetilde{u}_{\Phi}=u_{\Phi}-v_{\Phi}=-40K_{\Phi}$ and $\widetilde{v}_{\Phi}=v_{\Phi}=440K_{\Phi}$. The quartic potential is stable both along an equal-component direction, where $(u_{\Phi}-2v_{\Phi}/3)/K_{\Phi}=106.67$, and along a one-component direction, where $u_{\Phi}/K_{\Phi}=400$. For a uniform equal-amplitude state $\Phi_a=A$, the normalized energy per site is

At $J_{\chi}/K_{\Phi}=5.5$, its zero-temperature minimum is $A_0=0.721$, close to the measured transition-region root-mean-square amplitude 0.714. This stiff, bounded family is a controlled bridge from the sign limit rather than a material fit.

Simulation campaign and convergence

The exploratory $L=12$–$36$ runs used a general six-field (three bond fields plus three loop-current fields) implementation whose additional massive block was exactly decoupled from the three loop-current fields. A dedicated three-field implementation reproduces the archived loop-current observables on the same temperature ladder and removes the spectator from updates, reductions, and checkpoints. The promoted $J_{\chi}/K_{\Phi}=5.5$ campaign uses $L=24$, 48, 96, 144, 192, four independent chains, and a common 65-temperature feedback ladder over $2.68\leq T/K_{\Phi}\leq2.84$. Exchanges are attempted every two sweeps, and each chain-temperature pair retains 2048 measurements separated by ten sweeps.

All five sizes pass the preregistered convergence gates. From smallest to largest, the maximum split-$\hat{R}$ values are 1.0015, 1.0011, 1.0022, 1.0057, and 1.0044; the minimum effective sample sizes are 6648, 5230, 3320, 1509, and 697. The $L=192$ run has minimum neighboring acceptance 0.414, a longest autocorrelation time of 71.2 sweeps, and 26 distinct walkers that complete a cold–hot–cold round trip.

Crossings and order-content test

At $J_{\chi}/K_{\Phi}=4.75$, the small-size crossing sequence changes sign and belongs to the direct/crossover side. At $J_{\chi}/K_{\Phi}=5.0$, reasonable correction forms disagree, so that case remains unresolved. For $J_{\chi}/K_{\Phi}=5.5$, the $\xi/L$ sequence in Table S3 converges to the largest-pair values $(T_{\Phi}/K_{\Phi},T_{\chi}/K_{\Phi})=[2.76918(25),2.81310(64)]$, with normalized separation $(T_{\chi}-T_{\Phi})/K_{\Phi}=0.04392(65)$. Correlation-ratio crossings give the consistent separation 0.04358(43); all 4000 joint chain-and-block bootstrap draws have positive separation.

size pair	$T_\Phi/K_\Phi$	$T_\chi/K_\Phi$	$(T_\chi- T_\Phi)/K_\Phi$
24–48	2.7695(4)	2.8000(18)	0.0306(18)
48–96	2.76975(16)	2.81318(71)	0.04343(75)
96–144	2.76973(16)	2.81370(40)	0.04396(43)
144–192	2.76918(25)	2.81310(64)	0.04392(65)

At $T/K_{\Phi}=2.79121$, inside the crossing interval, the normalized primary structure factor is

for $L=24,48,96,144,192$, whereas

Power-law descriptions over the simulated sizes give $S_{\Phi}(0)/L^2\sim L^{-2.11}$ and $S_{\chi}(0)/L^2\sim L^{-0.09}$. The primary structure factor is consistent with saturation at a finite short-range value, while the time-reversal-odd composite remains extensive.

Matching to the rigid limit

The bulk-amplitude estimate

gives $q_{\chi}^{\mathrm{eff}} \simeq1.23, 1.30, 1.43$ for the three focused couplings. Because a wall can suppress the amplitude in its core, this bulk estimate is only a guide: the chiral scale varies as $A^{6}$ while the primary stiffness varies as $A^{2}$. Fully relaxed wall tensions are the appropriate quantitative matching variables.

Binder-ratio consistency

The dedicated $q = 2$ campaigns also measure the Binder ratios

Figure S1 shows that the primary curves develop a pronounced nonmonotonic maximum. Every successive-size difference in panel (b) consequently changes sign twice. The lower roots track the unique $R_s$ crossing used in the main text, whereas the upper roots drift strongly with size. This two-root ambiguity is why the primary Binder ratio is retained as a consistency diagnostic rather than used as the transition locator.

Image not extracted

The chiral Binder curves in Figure S1 are monotonic over the refined upper-transition window. Their successive-size differences in panel (d) have one zero each, providing an unambiguous independent check of the main-text $R_{\tau}$ crossings.