Maxwell's multiple pole theory provides the basis for a convenient means of determining the spherical harmonics of geomagnetic fields from directional paleomagnetic data. The relationship between the Maxwell poles and axes and the corresponding spherical harmonics was given by Maxwell. We show that the distribution of virtual geomagnetic poles (VGPs) is symmetrical about the Maxwell axes and converge to their poles. Utilizing this property of VGPs and a complimentary property of the distribution of equatorial virtual poles (EVPs), which are defined as points 90¿ from the VGP on the site meridian, leads to a means of obtaining the Gaussian spherical harmonic coefficients from paleomagnetic data. This VGP method involves the minimization of the horizontal components of the magnetic fields at VGPs and the radial components of the magnetic field at EVPs to yield the best fitting spherical harmonic coefficients. A hybrid variant of the VGP method involving this minimization and fitting the available mean vectors measured at sites has also been developed. Tests of the VGP and hybrid methods on model fields derived from the International Geomagnetic Reference Field (IGRF) 1995 demonstrate that they are effective means of determining the Gauss coefficients for any field requiring no a priori assumptions about the field. The hybrid method has been applied in a preliminary analysis of the time-averaged paleomagnetic data for the past 5 Myr and for the Bruhnes chron. The resulting mean fields were dominantly dipolar, but there were also persistent second order features, suggesting some long-term smaller-scale control over the dynamo process. The hybrid method was also applied to data from the Brunhes-Matuyama reversal. Four models were obtained: (1) a mean field model for reversed polarity, (2) a mean field model for normal polarity, (3) a mean field model for the entire reversal, and (4) a time sequence field model. The results were consistent with the Americas being a persistent site of inward field lines in the nondipole field and suggested that the reversal is initiated by decreases in strength of dipole, quadrupole, and octupole terms. The true dipole path of the transitional field tracks across eastern Asia, but a strong radially inward flux bundle moves over Africa. The paucity of demonstrably reliable data in these reversal records, however, requires that these results be interpreted cautiously. ¿ 1999 American Geophysical Union