  We determine the magnetic susceptibility of thermal QCD matter by means of
first principles lattice simulations using staggered quarks with physical
masses. A novel method is employed that only requires simulations at zero
background field, thereby circumventing problems related to magnetic flux
quantization. After a careful continuum limit extrapolation, diamagnetic
behavior (negative susceptibility) is found at low temperatures and strong
paramagnetism (positive susceptibility) at high temperatures. We revisit the
decomposition of the magnetic susceptibility into spin- and orbital angular
momentum-related contributions. The spin term -- related to the normalization
of the photon lightcone distribution amplitude at zero temperature -- is
calculated non-perturbatively and extrapolated to the continuum limit. Having
access to both the full magnetic susceptibility and the spin term, we calculate
the orbital angular momentum contribution for the first time. The results
reveal the opposite of what might be expected based on a free fermion picture.
We provide a simple parametrization of the temperature- and magnetic
field-dependence of the QCD equation of state that can be used in
phenomenological studies.
