Your manual approach is not reducing the loss in the first call:
loss = (Y_hat - Y) ** 2

so loss will have the shape [1, 1, 6, 7] and later you will sum this loss.
However, nn.MSELoss uses the reduction='mean' by default, which would yield a lower loss value and thus also gradients with a smaller magnitude.
If you use loss = nn.MSELoss(reduction='sum')(Y_hat, Y), you would get a similar result.