is Batch Norm a little bit stochastic by default?
While I am no expert in batch normalization, I have noticed that in the paper they include a constant, epsilon (ε) to the minibatch variance for numerical stability (in the normalize step from this paper, page 3). So while your math might be right and verified by others, the actual computation could still be wrong.

This issue appears here. The solution on the bottom is to add variance clipping.

And with a stab at why you might be seeing it more with smaller learning rates. Well that might be just because you're taking a lot more steps to get to the same accuracy so that's just many more steps of numerical instability.
EDIT: I now realize when you say full-batch gradient descent you mean you are doing the entire dataset before updating the weights (just standard gradient descent). This means that with a smaller learning rate you will have to be doing more steps... so I'm sure that explains the learning rate difference.