  We theoretically investigate the origin of the square-type skyrmion crystal
in centrosymmetric itinerant magnets, motivated from the recent experimental
finding in GdRu$_2$Si$_2$ [N. D. Khanh $\mathit{et}$ $\mathit{al.}$, Nat.
Nanotech. $\mathbf{15}$, 444 (2020)]. By simulated annealing for an effective
spin model derived from the Kondo lattice model on a square lattice, we find
that a square skyrmion crystal composed of a superposition of two spin helices
is stabilized in a magnetic field by synergy between the positive biquadratic,
bond-dependent anisotropic, and easy-axis anisotropic interactions. This is in
stark contrast to triangular skyrmion crystals which are stabilized by only one
of the three, suggesting that the square skyrmion crystal is characteristic of
itinerant magnets with magnetic anisotropy. We also show that a variety of
noncollinear and noncoplanar spin textures appear depending on the model
parameters as well as the applied magnetic field. The present systematic study
will be useful not only for identifying the key ingredients in GdRu$_2$Si$_2$
but also for exploring further skyrmion-hosting materials in centrosymmetric
itinerant magnets.
