  The goal of our work is to study the mode changes of the radio-quiet
gamma-ray pulsar PSR J2021+4026 with improved detail. By accurately
characterizing variations in the gamma-ray spectrum and pulse profile, we aim
to relate the Fermi-LAT observations to theoretical models and interpret the
mode changes in terms of variations in the structure of a multipolar
dissipative magnetosphere. We continually monitored the rotational evolution
and the gamma-ray flux of PSR J2021+4026 using more than 13 years of Fermi-LAT
data with a binned likelihood approach. We clearly detect the previous mode
changes and confirm a more recent mode change that occurred around June 2020.
We investigated the features of the phase-resolved spectrum and pulse profile,
and we inferred the macroscopic conductivity, the electric field parallel to
the magnetic field, and the curvature radiation cutoff energy. These physical
quantities are related to the spin-down rate and the gamma-ray flux and
therefore are relevant to the theoretical interpretation of the mode changes.
We computed the relative variations in the best-fit parameters, finding typical
flux changes between 13% and 20%. Correlations appear between the gamma-ray
flux and the spectral parameters, as the peak of the spectrum shifts by about
10% toward lower energies when the flux decreases. The analysis of the pulse
profile reveals that the pulsed fraction of the light curve is larger when the
flux is low. We introduced a simple magnetosphere model that combines a dipole
field with a strong quadrupole component. We simulated magnetic field
configurations to determine the positions of the polar caps for different sets
of parameters, and we conclude that some configurations could explain the
observed multiwavelength variability.
