Op-amp integrator response to square wave
I'm not entirely sure how to answer the why question, but simulation reproduces the same result; first the big picture:



But if you zoom in there are tiny little overshoot peaks in there...



which are actually pretty close to those measured. TL081 model is from TI, which is actually a very basic one... yet it predicted the overshoot pretty well.

Alas while I'm not completely sure what actually causes them, thanks to discussion below with the OP, all fingers point to the slew rate of the TL081. If I lower the slopes of the square wave by 3 orders of magnitude (from 1ns to 1us), then the output spikes go away (in simulation), but I also get some DC bias on the output (once it stabilizes):



Yeah, the datasheet slew rate for TL081 is 13 V/us. So in 1ns (first simulation) it would only be able to swing its output by 13mV, not enough to react to 25x greater input swing. Basically the output gets briefly pulled (up or down) by the input via the feedback resistor before the opamp can react to correct that. Actually it almost works like a voltage doubler during those spikes. Another way to check this is to measure the time span of the spikes. In [the first] simulation those spikes last 240ns (=4.15MHz), which is comparable to the 3-4MHz bandwidth of the TL081 (I note that the old datasheet said 4MHz but the new one says 3MHz.)