The geomagnetic field is generated by a dynamo process in the Earth’s core and is characterized by a predominant dipole component that has been steadily decreasing in the last few centuries. The physical drivers behind the fluctuations of the geomagnetic dipole field remain poorly understood. One of the possible explanations rely on the interaction between the dipole mode and other multipole terms of the geomagnetic field. To test this hypothesis, we used two millennial scale models based on spherical harmonic fitting of paleomagnetic data, which allowed to reconstruct the geomagnetic field of the past. By performing causality and information statistical analysis, we found significant interactions between the dipole and smaller scale harmonics (quadrupole and octupole) of the geomagnetic field. In particular, both data sets agree that the spherical harmonic $$Y_2^2$$ acts as a source term, whereas the axial dipole term $$Y_1^0$$ consists of the term with least information loss. The results suggest a possible control of core–mantle boundary inhomogeneities on the interaction between the components of the geomagnetic field. Our results also show a net information flux from larger to smaller scales, which is compatible with a direct turbulent cascade view of the geodynamo.