  We compute rigorously the scaling limit of multi-point energy correlations in
the critical Ising model on a torus. For the one-point function, averaged
between horizontal and vertical edges of the square lattice, this result has
been known since the 1969 work of Ferdinand and Fischer. We propose an
alternative proof, in a slightly greater generality, via a new exact formula in
terms of determinants of discrete Laplacians. We also compute the main term of
the asymptotics of the difference $\mathbb{E}(\epsilon_{V}-\epsilon_{H})$ of
the energy density on a vertical and a horizontal edge, which is of order of
$\delta^{2}$, where $\delta$ is the mesh size. The observable
$\epsilon_{V}-\epsilon_{H}$ has been identified by Kadanoff and Ceva as (a
component of) the stress-energy tensor.
  We then apply the discrete complex analysis methods of Smirnov and Hongler to
compute the multi-point correlations. The fermionic observables are only
periodic with doubled periods; by anti-symmetrization, this leads to
contributions from four "sectors". The main new challenge arises in the doubly
periodic sector, due to the existence of non-zero constant (discrete) analytic
functions. We show that some additional input, namely the scaling limit of the
one-point function and of relative contribution of sectors to the partition
function, is sufficient to overcome this difficulty and successfully compute
all correlations.
