{'fraction_non_alphanumeric': 0.03945013245507267, 'fraction_numerical': 0.017827736808190737, 'mean_word_length': 4.420256111757858, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Efficient Langevin Simulation of Coupled Classical Fields and Fermions

We introduce an efficient Langevin method to study bilinear Fermionic Hamiltonians interacting with classical fields. Our method is suitable for very large systems and offers high accuracy. To demonstrate the method, we study complex non-coplanar chiral spin textures on the triangular Kondo lattice model. We also explore non-equilibrium mesoscale physics such as chiral domain coarsening and Z2 vortex annihilation.

Lattice models of fermions interacting with classical fields encompass a wide range of physics. Popular examples in condensed matter include Kondo lattice (KL) models of itinerant electrons interacting with localized magnetic moments [1], Falicov-Kimball models of metal-insulator transitions in rare-earth materials [2] and Bogoliubov-De Gennes equations for superconductivity [3]. The Hubbard-Stratonovich transformation is another path to obtaining bilinear fermionic systems coupled to an auxiliary classical field [4,5]. This broad class of models poses a notoriously difficult numerical challenge: Monte Carlo (MC) sampling of the classical field requires repeated diagonalization of the single-particle fermion matrix.
Several MC methods have been developed to more efficiently sample the classical field [6][7][8][9]. Spurred by colossal magnetoresistance (CMR) [10][11][12], these methods have largely been applied to the ferromagnetic transition in KL models at large coupling. This transition is relatively easy to study using moderate temperatures and small system sizes.
Recent interest has shifted to exotic spin-textures, which would occur in KL models at small to moderate couplings. Skyrmion lattices have recently been observed with spatial modulations up to 0.1µm [13,14]. Chiral textures lead to an anomalous Hall effect associated with huge (∼ 10 5 T ) effective magnetic fields, as predicted in the KL model with triangular lattice [15,16] and experimentally observed in Pr 2 Ir 2 O 7 [17] and UCu 5 [18]. Compared to ferromagnetism, these spin-textures can be very challenging to study. High precision and large system sizes may be needed to capture the physics of low temperatures and effective long-range interactions. State of the art numerical methods are often impractical.
In this Letter we introduce a suitable Langevin sampling method that is very efficient-the cost scales linearly with system size-at high accuracy. Our method is based on a non-trivial gradient transformation [19] of the kernel polynomial method (KPM) [20]. Below we outline our method and then demonstrate it with a study of the triangular KL model. Our lattices are large enough to uncover interesting non-equilibrium effects such as chiral domain coarsening and Z 2 vortex dynamics. In this way, we bridge the gap between quantum atomic scale and mesoscale physics.
Our method applies to a general bilinear fermionic Hamiltonian coupled to continuous, classical degrees of freedom φ, with sparse matrix A. We work at fixed temperature β −1 and chemical potential µ. The partition function is a trace over classical and fermionic degrees of freedom, The energy contains effective long-range many-body interactions of the classical field.
A key difficulty in MC sampling the classical field is the calculation of ∆F in response to changes in φ. KPM estimates the density of states ρ( ) using a series of Chebyshev polynomials T m ( ) truncated at order M [20,21]. The cost of KPM is linear in system size N for sparse matrix A. We state directly the recursive KPM procedure to construct an unbiased stochastic estimate of F (φ), Here, r is a random column vector whose components satisfy r * i r j = δ ij . We draw complex r i from a uniform distribution |r i | 2 = 1. The coefficients is chosen to damp Gibbs oscillations yet retain high accuracy. KPM requires that the eigenvalues of A have magnitude less than 1, which can usually be achieved without loss of generality as a rescaling of energy. There are two independent sources of error in KPM: (1) truncation at finite order M and (2) stochastic estimation by averaging over finitely many random vectors r. Both are well controlled and will be discussed below.
To sample fields {φ} from the Boltzmann distribution, , we apply the overdamped Langevin equation. In discretized form, where η i (t) are uncorrelated Gaussian random variables with unit variance and t is a fictitious time. The Langevin approach simultaneously updates all components φ i . Efficient and accurate estimation of the gradient ∂F/∂φ i is crucial. We exclude inertial terms from the Langevin equation because they would amplify errors in the gradient estimate. The technique of automatic differentiation with "reverse accumulation" [19] ensures that, by careful application of the chain rule, we can transform the KPM procedure to estimate F , Eqs. 3-5, into one that estimates ∂F/∂φ i at the same cost. We perform this transformation analytically and state only the final result, The row vectors β m are given by reverse recursion, from m = M − 2 down to m = 0, The sequence of vectors α m are the same as in the original KPM, but are here required in reverse order. We recalculate them as needed using α m = 2Aα m+1 − α m+2 . The recursion begins with α M −1 and α M −2 , which are available at the end of the original KPM procedure.
Note that the procedure to estimate all components of ∂F/∂φ i , Eqs. 7-8, has computational cost equivalent to the original KPM procedure to calculate F .
The gradient calculation also inherits the approximation errors of KPM, controlled by two parameters: the truncation order M of the Chebyshev series, and the dynamical stochastic error z ≡ ∆t/Q, where Q is the number of of KPM random vectors used per time step. The KPM estimated density of states ρ( ) is resolved to order ∆ /M , where ∆ = max − min is the span of extremal eigenvalues. The parameter M should be chosen large enough to resolve the physically relevant features of the density of states. Finite stochastic error z > 0 acts much like an additional noise term in the Langevin dynamics, effectively rescaling its magnitude by an amount T → T + ∆T eff . For matrices of the form A(Jφ), where J 1 is a small coupling constant, the estimate of the Langevin force term ∆t∂F/∂φ i includes a stochastic error that scales like ∆tJ/ √ Q = √ J 2 z∆t. Comparison with Eq. 6 suggests modeling this stochastic error as an additional Langevin noise term with a temperature that scales as ∆T eff ∼ J 2 z. The parameter z should be chosen small enough that ∆T eff T for the smallest relevant temperature scale T .
Our Langevin sampling remains efficient at high accuracy: the cost to integrate the Langevin equation one  [7]. In simulations presented below with d = 2, N = 100 2 , M = 1000, and z = 0.02, our Langevin approach outperforms existing methods by 2 orders of magnitude or more. Another linear cost algorithm is based on hybrid Monte-Carlo [6], but a direct comparison is difficult because the cost of precision at low temperatures is unclear.
We apply our method to the triangular KL model defined by the Hamiltonian, where c † jσ (c jσ ) is the creation (annihilation) operator of an electron with spin σ on site j, S j is a classical Heisenberg spin with |S j | = 1, and σ µν = (σ x µν , σ y µν , σ z µν ) is a vector of Pauli matrices. The hopping coefficients are t ij = t when i and j are nearest neighbor sites on the triangular lattice, and t ij = 0 otherwise. In the following, we fix the energy scale by taking t → 1, and the spatial scale by taking the lattice spacing to 1.
Martin and Batista argue, by perfect nesting of the Fermi surface, that the chiral 3q configuration (Fig. 1c) is the ground state at 3/4 electron filling fraction with small coupling [15]. This state is of special interest, as it exhibits a spontaneous quantum Hall effect. Variational calculation on the 2 × 2 plaquette also predicts stability of the 3q state [23]. However, unconstrained Monte-Carlo study of this phase at 3/4 filling has not yet been achieved due to severe numerical difficulties. Very low temperatures, T 0.001 and small couplings J 0.3 are required to stabilize 3q. The numerical method must be very accurate to resolve the small gap in the density of states (of width ∼ J 2 ). Furthermore, the 3q state is stabilized by a susceptibility that diverges like log 2 N , so very large lattices sizes are required (N ≈ 100 2 ). Due to these challenges, this system offers a rigorous test of our Langevin method.
In Fig. 2(a) we observe melting of the chiral 3q ground state. The mean chirality χ = S i × S j · S k of triangular plaquettes [ijk] abruptly disappears at a first order phase transition at T ≈ 0.0010. This transition, however, is not to a paramagnetic phase. To distinguish the phases 3q, 2q, and 1q we consider the (unordered) set of nearest-neighbor spin-spin correlations C = { S x · S x+ 1,0 , S x ·S x+ 1, metry of the triangular lattice). Fig. 2(b) plots the three elements of C as a function of temperature. At T = 0 we find C = C 3q as expected. We now observe three first order transitions at temperatures T = 0.0010, 0.0017, and T = 0.0029 to the 2q, 1q, and paramagnetic phases, respectively. The 2q phase is identified by its correlation set C, which has two zero elements and one negative element. In the 1q phase, C has one positive element and two negative (symmetric) elements. To avoid equilibration issues in the above data, we used initial conditions with explicitly broken chiral symmetry, χ > 0, We now investigate the dynamical, non-equilibrium process by which 3q chiral symmetry breaking occurs at low temperatures. We use our Langevin dynamics to study the phase ordering kinetic of chiral domains following a quench from infinite to zero temperature.
Langevin "time", properly speaking, is fictitious. However, in the spirit of time-dependent Ginzburg-Landau (TDGL) models, we expect the energetic relaxation of overdamped Langevin dynamics to qualitatively capture the large-scale aspects of phase ordering [24]. A point of comparison is Model A dynamics in the Hohenberg and Halperin classification [25], the prototypical TDGL model for phase-ordering of a non-conserved scalar order parameter. Interestingly, we find that the KL model has effective effective long-range many-body interactions that introduce dynamical features not present Model A.
First we consider the case of ∼ 1/4 filling with J = 3 and µ = −3.2. Previous work found a robust 3q phase at system sizes up to N = 16 2 [16]. We use our Langevin dynamics to study the ordering dynamics at N = 200 2 with accuracy parameters M = 500, z = 0.005. Figure 3(a) shows the evolution of chirality ranging from red (positive) to blue (negative). The coarsening of chiral domains is analogous to Model A, but we observe strong anisotropy of domain walls. The dynamics slows at large times, consistent with a characteristic length scale that grows as ∼ t 1/2 [26].
Next we consider ∼ 3/4 filling, with J = 0.2 and µ = 1.947. We use accuracy parameters M = 500 and z = 0.02. The ordering dynamics is displayed in Fig. 3(b). A new dynamical feature appears: the chiral domains evolve into a remarkable pattern which is a very long lived metastable state. This is a reproducible phenomenon. At higher temperatures this metastable pattern could be annealed to the pure 3q phase, but in experimental practice it is easier to explicitly break chiral symmetry with an applied external magnetic field B z [17]. We introduce orbital coupling into our model, Eq. 9, by applying a non-uniform phase θ ij = B zẑ · x i × x j /2 to the hopping coefficients t ij = t exp(−iθ ij ), with x i the position of lattice site i. The smallest magnetic field consistent with periodic boundaries, B z = 8π/ √ 3N , causes the system to rapidly reach the uniform chiral 3q phase, shown in Fig. 3(c).
Topological defects, visible as small white dots, are apparent at both 1/4 and 3/4 filling. These defects are Z 2 vortices associated with winding of the SO(3) topological manifold, and predicted to have fractional charge [27]. A Z 2 vortex is enlarged in Fig. 3(c), second and third panels. One of the four 3q spin sub-lattices are shown. We can understand this defect by constructing a closed Burger's circuit (draw in green) that encircles it. Every point on the green circuit is identified with an element in SO (3), plotted as a trajectory in the fourth panel. The SO(3) manifold, in the axis-angle representation, is a filled-sphere with antipodal points identified. This Burger's circuit has winding number 1 because it wraps SO(3). The vortex is Z 2 because the only homotopically distinct winding numbers are 0 and 1. Consequently, any pair of Z 2 vortices may annihilate. In the ordering dynamics of Fig. 3(a), we observe many vortices annihilating with each other and with domain walls.
In conclusion, we have introduced a numerical method to study the broad class of Hamiltonians that couple fermions to classical degrees of freedom. Our method is highly accurate and efficient, enabling the study of complex systems at unprecedented size. Large system sizes may be necessary to resolve logarithmic divergences associated with nesting of the fermi surface. In the triangular Kondo lattice model at 3/4 filling, we found that lattices of size N = 100 2 with six digits of precision are required. Large system sizes also allow us to bridge the gap between quantum and mesoscopic physics. With systems of size N = 200 2 we are able to probe chiral domain dynamics, metastable trapping, and Z 2 vortex dynamics, effects typically inaccessible in standard approaches.
We thank Ivar Martin and Cristian Batista for useful discussions. This work was carried out under the auspices of the NNSA of the U.S. DOE at LANL under Contract No. DE-AC52-06NA25396 and supported by the LANL/LDRD Program. The calculations presented were performed using the CCS-7 Darwin cluster.