Both nonzero temperature and chemical potentials break the Lorentz symmetry present in vacuum quantum field theory by singling out the rest frame of the heat bath. This leads to complications in the application of thermal perturbation theory, including the appearance of novel infrared divergences in loop integrals and an apparent absence of four-dimensional integration-by-parts (IBP) identities, vital for high-order computations. Here, we propose a new strategy that enables the use of IBP techniques in the evaluation of Feynman integrals, in particular vacuum or bubble diagrams, in the limit of vanishing temperature T but nonzero chemical potentials μ. The central elements of the new setup include a contour representation for the temporal momentum integral, the use of a small but nonzero T as an IR regulator, and the systematic application of both temporal and spatial differential operators in the generation of linear relations among the loop integrals of interest. The relations we derive contain novel inhomogeneous terms featuring differentiated Fermi-Dirac distribution functions, which severely complicate calculations at nonzero temperature, but are shown to reduce to solvable lower-dimensional objects as T tends to zero. Pedagogical example computations are kept at the one- and two-loop levels, but the application of the new method to higher-order calculations is discussed in some detail.