As the comments on the original question have stated, this is due to the spherical harmonics of the Earth. This answer is adding more context.
In the two-body problem, one assumes that the central body (in this case Earth) is a perfect sphere with evenly distributed mass. If this is true, its force due to gravity can be modeled as a point source at its geometric center. Of course in real life this isn't true. Since the Earth is rotating, its equatorial radius is larger than its polar radius. This is the J2 perturbation, which is a first order approximation of Earth's spherical harmonics. Also note that the farther you get away from the Earth, the weaker this perturbation gets, since Earth looks more like a point mass from far away.
In the cases of sun synchronous orbits, this is actually extremely useful. One can set up their altitude and inclination such that the J2 perturbation will cause a drift in right ascension that matches the rate of change of Earth's true anomaly around the sun. Here is a plot of the Keplerian orbital elements for a sun synchronous orbit for one year:
Since this sun synchronous orbit is also in LEO, it shows the same oscillations in semi-major axis and eccentricity. The secular decrease in the semi-major axis is numerical error, since we know that gravity is a conservative force, and semi-major axis is inversely proportional to (absolute value) of orbital energy. This plot was made from the results of integration with a RK4/5 solver. Other solvers may have more or less energy drift for these equations of motion
Also note J2 is only a first order estimation. Spherical harmonic models can get much more complex. For example, the GRAIL missions around the Moon produced a 1500x1500 gravity field for the Moon (JGGRX model).
It is possible to represent the exact variation in these orbital elements using the Gaussian Variation of Parameters (VOPs), from Vallado 4th edition, pages 636. Please refer to the ten prior pages for the exact derivation. Note that general astrodynamics tool (like GMAT or Nyx) do not use the Gaussian VOP for their propagation as Keplerian orbital elements are singular (e.g. one cannot use these equations for near-circular or near-equatorial orbits).