We present a novel simulation prescription for thermal quantum fields on a lattice that operates directly in imaginary frequency space. By distinguishing initial conditions from quantum dynamics it provides access to correlation functions also outside of the conventional Matsubara frequencies ωn = 2πnT . In particular it resolves their frequency dependence between ω = 0 and ω1 = 2πT , where the thermal physics ω ∼ T of e.g. transport phenomena is dominantly encoded. Real-time spectral functions are related to these correlators via an integral transform with rational kernel, so their unfolding is exponentially improved compared to Euclidean simulations.We demonstrate this improvement within a 0 + 1-dimensional scalar field theory and show that spectral features inaccessible in standard Euclidean simulations are quantitatively captured.Modern experiments ranging from heavy-ion collisions at RHIC and LHC [1, 2] to ultracold quantum gases [3] explore the physics of strongly-correlated quantum systems across vastly separated temperature scales (T ∼ 10 12 ..10 −9 K) [4]. A unifying aspect of these studies are thermal phenomena, such as the transport of conserved charges [5][6][7][8][9][10][11][12][13][14][15][16] or the in-medium modification of heavy bound states [17][18][19][20][21][22][23]. Connecting these to fundamental interactions, such as QED and QCD, requires a thorough theoretical understanding of the equilibrium properties of quantum fields. Unfortunately in realistic experimental environments many approximate techniques, such as perturbation theory, are inapplicable and a non-perturbative numerical computation is required.Conventional simulations of quantum fields, e.g. lattice QCD, are based on non-perturbative Monte-Carlo techniques relying on an analytic continuation of the realtime axis into the complex plane [24,25], for real-time approaches to spectral functions see e.g. [8][9][10][26][27][28][29][30][31][32][33][34].Two conceptual challenges plague the necessary reconstruction of real-time correlation functions from those simulated in imaginary time: the crucial one is the finite extent of the Euclidean time axis with τ ∈ [0 , 1/T ]. It leads to discrete Matsubara frequencies ω n = 2πnT , which limits the resolution of imaginary frequencies. In turn, increasing the number of temporal points N τ simply increases the maximum frequency available while keeping the spacing fixed. Thus, a large N τ does not improve the access to the relevant regime ω ∼ T , where thermal physics is dominantly encoded in G(ω). The signal to noise ratio for thermal contributions actually decays rapidly above ω ≈ T and indeed, already G(ω 1 = 2πT ) may only be marginally relevant for thermal physics.In order to resolve this conceptual issue, we present a simulation algorithm, which operates in imaginary frequency space with arbitrary resolution, given by the number of frequency points N ω . In particular this approach resolves frequencies between ω 0 and ω 1 .The second issue concerns the unfolding of spectral functions ρ(µ) f...