  Evaluating the linear response of a driven system to a change in environment
temperature(s) is essential for understanding thermal properties of
nonequilibrium systems. The system is kept in weak contact with possibly
different fast relaxing mechanical, chemical or thermal equilibrium reservoirs.
Modifying one of the temperatures creates both entropy fluxes and changes in
dynamical activity. That is not unlike mechanical response of nonequilibrium
systems but the extra difficulty for perturbation theory via path-integration
is that for a Langevin dynamics temperature also affects the noise amplitude
and not only the drift part. Using a discrete-time mesh adapted to the
numerical integration one avoids that ultraviolet problem and we arrive at a
fluctuation expression for its thermal susceptibility. The algorithm appears
stable under taking even finer resolution.
