We introduce an efficient Langevin method to study bilinear fermionic Hamiltonians interacting with classical fields. Our approach is orders of magnitude faster than previous methods when applied to very large systems with high accuracy requirements. 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] 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,Evaluat...