  Supernovae (SN), the most energetic stellar feedback mechanism, are crucial
for regulating the interstellar medium (ISM) and launching galactic winds. We
explore how supernova remnants (SNRs) create a multiphase medium by performing
3D hydrodynamical simulations at various SN rates, $S$, and ISM average
densities, $\bar{n}$. The evolution of a SNR in a self-consistently generated
three-phase ISM is qualitatively different from that in a uniform or a
two-phase warm/cold medium. By travelling faster and further in the low-density
hot phase, the domain of a SNR increases by $>10^{2.5}$. Varying $\bar{n}$ and
$S$, we find that a steady state can only be achieved when the hot gas volume
fraction $f_{\rm{V,hot}}\lesssim 0.6 \pm 0.1 $. Above that level, overlapping
SNRs render connecting topology of the hot gas, and the ISM is subjected to
thermal runaway. Photoelectric heating (PEH) has a surprisingly strong impact
on $f_{\rm{V,hot}}$. For $\bar{n}\gtrsim 3 \cm-3 $, a reasonable PEH rate is
able to suppress the thermal runaway. Overall, we determine the critical SN
rate for the onset of thermal runaway to be $S_{\rm{crit}} = 200
(\bar{n}/1\cm-3)^k (E_{\rm{SN}}/10^{51}\erg)^{-1} \kpc^{-3} \myr-1$, where $k =
(1.2,2.7)$ for $\bar{n} \leq 1$ and $> 1\cm-3 $, respectively. We present a
fitting formula of the ISM pressure $P(\bar{n}$, $S$), which can be used as an
effective equation of state in cosmological simulations. Despite the 5 orders
of magnitude span of $(\bar{n},S)$, the average Mach number varies little:
$\mathcal{M} \approx \ 0.5\pm 0.2, \ 1.2\pm 0.3,\ 2.3\pm 0.9$ for the hot, warm
and cold phases, respectively.
