  We analyze new XMM-Newton and archival Chandra observations of the
middle-aged $\gamma$-ray radio-quiet pulsar J1957+5033. We detect, for the
first time, X-ray pulsations with the pulsar spin period of the point-like
source coinciding by position with the pulsar. This confirms the pulsar nature
of the source. In the 0.15--0.5 keV band, there is a single pulse per period
and the pulsed fraction is $\approx18\pm6$ per cent. In this band, the pulsar
spectrum is dominated by a thermal emission component that likely comes from
the entire surface of the neutron star, while at higher energies ($\gtrsim0.7$
keV) it is described by a power law with the photon index $\Gamma \approx 1.6$.
We construct new hydrogen atmosphere models for neutron stars with dipole
magnetic fields and non-uniform surface temperature distributions with
relatively low effective temperatures. We use them in the spectral analysis and
derive the pulsar average effective temperature of $\approx(2-3)\times10^5$ K.
This makes J1957+5033 the coldest among all known thermally emitting neutron
stars with ages below 1 Myr. Using the interstellar extinction--distance
relation, we constrain the distance to the pulsar in the range of 0.1--1 kpc.
We compare the obtained X-ray thermal luminosity with those for other neutron
stars and various neutron star cooling models and set some constraints on
latter. We observe a faint trail-like feature, elongated $\sim 8$ arcmin from
J1957+5033. Its spectrum can be described by a power law with a photon index
$\Gamma=1.9\pm0.5$ suggesting that it is likely a pulsar wind nebula powered by
J1957+5033.
