A two-dimensional, global-scale hybrid simulation is carried out to study the kinetic structure of the quasi-parallel bow shock. In the simulation the bow shock forms by the interaction of the supersonic solar wind and the geomagnetic field. Strong temporal electromagnetic waves occur in the shock transition and foreshock regions due to backstreaming and reflected ions at the bow shock. This self-consistent, global-scale simulation shows the formation of diamagnetic cavities in the foreshock regions of the parallel and quasi-parallel shocks, where ion beams interact with the incoming solar wind plasma. The cavities are crater-like and saturate at a width ~1--2 RE, with a low-density and low-magnetic field center bounded by a rim of high density and high magnetic field. The bulk flow speed decreases slightly in the center of the cavity, while the ion temperature often, but not always, increases. The craters convect downstream with the solar wind flow. When the interplanetary magnetic field (IMF) lies nearly parallel to the solar wind flow, the craters develop into elongated phase-standing spatial structures along the field lines, both upstream and downstream, with alternate increases and decreases in the density and corresponding in-phase variations in the magnetic field strength. In the general cases with an oblique IMF the foreshock cavities can convect toward the Earth as transient compressional structures and impinge on the magnetopause. The generation and structure of the diamagnetic cavities are compared with those of hot flow anomalies. The simulation results are also compared with satellite observations of the foreshock of quasi-parallel shocks.