  We propose a variational modelling method with differentiable temperature for
canonical ensembles. Using a deep generative model, the free energy is
estimated and minimized simultaneously in a continuous temperature range. At
optimal, this generative model is a Boltzmann distribution with temperature
dependence. The training process requires no dataset, and works with arbitrary
explicit density generative models. We applied our method to study the phase
transitions (PT) in the Ising and XY models, and showed that the
direct-sampling simulation of our model is as accurate as the Markov Chain
Monte Carlo (MCMC) simulation, but more efficient. Moreover, our method can
give thermodynamic quantities as differentiable functions of temperature akin
to an analytical solution. The free energy aligns closely with the exact one to
the second-order derivative, so this inclusion of temperature dependence
enables the otherwise biased variational model to capture the subtle thermal
effects at the PTs. These findings shed light on the direct simulation of
physical systems using deep generative models
