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 \gamma -ray binary LS I + 61 ^ { \circ } 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 \sim 26.4960 -day orbital period .
Here we show that , during the time of our observations , the \gamma -ray emission of LS I + 61 ^ { \circ } 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 \gamma -ray observations to study the outflows of massive stars in eccentric binary systems .