We have investigated the field around the radio-quiet $gamma$-ray pulsar, PSR J2021+4026, with a ~140 ks XMM-Newton observation and a ~56 ks archival Chandra data. Through analyzing the pulsed spectrum, we show that the X-ray pulsation is purely thermal in nature which suggests the pulsation is originated from a hot polar cap with $Tsim3times10^{6}$ K on the surface of a rotating neutron star. On the other hand, the power-law component that dominates the pulsar emission in the hard band is originated from off-pulse phases, which possibly comes from a pulsar wind nebula. In re-analyzing the Chandra data, we have confirmed the presence of bow-shock nebula which extends from the pulsar to west by ~10 arcsec. The orientation of this nebular feature suggests that the pulsar is probably moving eastward which is consistent with the speculated proper motion by extrapolating from the nominal geometrical center of the supernova remnant (SNR) G78.2+2.1 to the current pulsar position. For G78.2+2.1, our deep XMM-Newton observation also enables a study of the central region and part of the southeastern region with superior photon statistics. The column absorption derived for the SNR is comparable with that for PSR J2021+4026, which supports their association. The remnant emission in both examined regions are in an non-equilibrium ionization state. Also, the elapsed time of both regions after shock-heating is apparently shorter than the Sedov age of G78.2+2.1. This might suggest the reverse shock has reached the center not long ago. Apart from PSR J2021+4026 and G78.2+2.1, we have also serendipitously detected an X-ray flash-like event XMM J202154.7+402855 from this XMM-Newton observation.