For the experiment using a single GPU with gradient accumulation I divide the loss by the accumulation steps as suggested here 1. (Is this really necessary?)

The linked post focuses on DistributedDataParallel (DDP), which is different from DataParallel (DP). For DDP, every process/rank/device compute its own local loss, while for DP, the loss is global across all devices. I am not sure if this contributes to the discrepancy you saw.

The effective batch size of all 3 experiments is the same.

What do you mean by effective batch size? I assume the input batch size used in every iteration should be the same?