How does batch normalization work for convolutional neural networks

I am trying to understand how batch normalization (BN) works in CNNs. Suppose I have a feature map tensor $T$ of shape $(N, C, H, W)$

where $N$ is the mini-batch size,

$C$ is the number of channels, and

$H,W$ is the spatial dimension of the tensor.

Then it seems there could a few ways of going about this (where $gamma, beta$ are learned parameters for each BN layer)

Method 1: $T_{n,c,x,y} := gamma*frac {T_{c,x,y} – mu_{x,y}} {sqrt{sigma^2_{x,y} + epsilon}} + beta$ where $mu_{x,y} = frac{1}{NC}sum_{n, c} T_{n,c,x,y}$ is the mean for all channels $c$ for each batch element $n$ at spatial location $x,y$ over the minibatch, and

$sigma^2_{x,y} = frac{1}{NC} sum_{n, c} (T_{n, c,x,y}-mu_{c})^2$ is the variance of the minibatch for all channels $c$ at spatial location $x,y$.

Method 2: $T_{n,c,x,y} := gamma*frac {T_{c,x,y} – mu_{c,x,y}} {sqrt{sigma^2_{c,x,y} + epsilon}} + beta$ where $mu_{c,x,y} = frac{1}{N}sum_{n} T_{n,c,x,y}$ is the mean for a specific channel $c$ for each batch element $n$ at spatial location $x,y$ over the minibatch, and

$sigma^2_{c,x,y} = frac{1}{N} sum_{n} (T_{n, c,x,y}-mu_{c})^2$ is the variance of the minibatch for a channel $c$ at spatial location $x,y$.

Method 3: For each channel $c$ we compute the mean/variance over the entire spatial values for $x,y$ and apply the formula as

$T_{n, c,x,y} := gamma*frac {T_{n, c,x,y} – mu_{c}} {sqrt{sigma^2_{c} + epsilon}} + beta$, where now $mu_c = frac{1}{NHW} sum_{n,x,y} T_{n,c,x,y}$ and $sigma^2{_c} = frac{1}{NHW} sum_{n,x,y} (T_{n,c,x,y}-mu_c)^2 $

In practice which of these methods is used (if any) are correct for?

The original paper on batch normalization , https://arxiv.org/pdf/1502.03167.pdf , states on page 5 section 3.2, last paragraph, left side of the page:

For convolutional layers, we additionally want the normalization to
obey the convolutional property – so that different elements of the
same feature map, at different locations, are normalized in the same
way. To achieve this, we jointly normalize all the activations in a
minibatch, over all locations. In Alg. 1, we let $mathcal{B}$ be the set of all
values in a feature map across both the elements of a mini-batch and
spatial locations – so for a mini-batch of size $m$ and feature maps of
size $p times q$, we use the effective mini-batch of size $m^prime = vert mathcal{B} vert = m cdot pq$. We learn a pair of parameters $gamma^{(k)}$ and $beta^{(k)}$ per feature map,
rather than per activation. Alg. 2 is modified similarly, so that
during inference the BN transform applies the same linear
transformation to each activation in a given feature map.

I’m not sure what the authors mean by "per feature map", does this mean per channel?