Gamma-ray binaries are stellar systems for which the spectral energy distribution (discounting the thermal stellar emission) peaks at high energies. Detected from radio to TeV gamma rays, the γ-ray binary LS I +61°303 is highly variable across all frequencies. One aspect of this system's variability is the modulation of its emission with the timescale set by the ~26.4960 day orbital period. Here we show that, during the time of our observations, the γ-ray emission of LS I +61°303 also presents a sinusoidal variability consistent with the previously known superorbital period of 1667 days. This modulation is more prominently seen at orbital phases around apastron, whereas it does not introduce a visible change close to periastron. It is also found in the appearance and disappearance of variability at the orbital period in the power spectrum of the data. This behavior could be explained by a quasi-cyclical evolution of the equatorial outflow of the Be companion star, whose features influence the conditions for generating gamma rays. These findings open the possibility to use γ-ray observations to study the outflows of massive stars in eccentric binary systems.
