@ptrblck thanks for pointing out the channel dimension, I revisited the docstring. I only had vectorial data (N x L input size) in mind, and did not think about the use case of convolution with (N x C x L input size).
The current behaviour of batch norm on tensors (N x C x L) is designed to conform with what should be used with 1D convolutional outputs. That is, each channel will get its own mean and var, right? And all elements in the channel will be standardised with the same mean, var pair.
I can then reframe my issue:
The default behaviour of BatchNorm on tensors (instead of matrices), does not broadcast the matrix operation (element wise batch norm), but completely switches to a different kind operation (channel wise batch norm).
I appreciate the fact, that this may be a sound default for the many users working with CNNs (including me  ). I however see many cases, where an actual broadcasting of the BatchNorm layer to higher order tensors is desired. The broadcasting would then be identical to how linear layers cope with higher dimensional tensors (only apply to last dimension, broadcast else).
Maybe this is too didactic, but it could be helpful, if the library explicates the difference of applying an elementwise batch normalization to a channel wise batch normalization by setting a default flag.