  We numerically investigate the stabilization mechanisms of skyrmion crystals
under thermal fluctuations and external magnetic field in itinerant
centrosymmetric tetragonal magnets. By adopting an efficient steepest descent
method with a small computational cost, we systematically construct the
magnetic field-temperature phase diagrams of the effective spin model derived
from the itinerant electron model on a two-dimensional square lattice. As a
result, we find that a square-type skyrmion crystal is stabilized by either or
both of the high-harmonic wave-vector interaction and the biquadratic
interaction under an external magnetic field. Especially, we discover that the
former high-harmonic wave-vector interaction can stabilize the skyrmion crystal
only at finite temperatures when its magnitude is small. In addition to the
skyrmion crystal, we also find other stable multiple-$Q$ states in the phase
diagram. Lastly, we discuss the correspondence of the phase diagrams between
the effective spin model and the skyrmion-hosting material GdRu$_2$Si$_2$. The
present results suggest a variety of multiple-$Q$ states could be driven by
thermal fluctuations and external magnetic fields in centrosymmetric itinerant
magnets.
