"Naturally, this will affect the gradient through the network." this statement is only partially true, let's see why by starting explaining the real aim of batch normalisation.
As the title of the paper suggest, the aim of batch normalisation is to decrease training time by reducing covariance shift. What is covariance shift? We can conceive it as the variation that can happen between between the values of two layers of a network. We are all familiar with the concept that if we have input features with different unite scale, like kilos and euros, most likely lot of values will have different magnitude orders, like thousands for the weight could appear often with hundred of thousands for the money. When applying the activation function to different order of magnitude this discrepancy will remain, causing the values in the first layer to assume a really broad range of values. This is not good, since high fluctuation means more time to converge to stable values. This is why values fed to neural nets are always standardised.
The authors applied the same logic to the hidden layers, arguing that a deep neural net can be conceived as a repetition of itself (every hidden layer is an input layer that send features to another hidden layer) ergo the features should be normalised at every layer. How can we do it? Normalising every batch is the most natural way to procede, but in this way there is the risk to end up transforming the internal representation of a layer, because normalising is not a linear transformation. This is way the authors propose a clever way to perform hidden layers normalisation, which consists in a classic normalisation followed by a linear scaling performed with two trainable parameters $\beta$ and $\gamma$ (which appear in the last step of the batch norm scratch code below).
A really important thing to notice is that the mean and variance used to perform the classic normalisation are mean and variance calculated on the mini batch. I will explain why this is important in a sec, first I want to stress out that the $\beta$ parameter can actually bring to increase overfitting when batch norm is randomly stucked on top of other layer. The reason why is that, as we see in the scale and shift step, $\beta$ is nothing but a biased term, added to perform the shift of the mean of the batch hidden values. So, to avoid having an extra bias that lead to overfit, the bays on the previous layer should be removed, leaving only the classic matrix of weights parameters.
Back to the gradient problem, we can see that in itself doesn't necessarily lead to increased performances, but it does provide an advantage in terms of hidden layer values convergence. The x axis on the two right sub plots of the figure below represent the variation of the hidden values of net trained with and without batch norm. When trained with batch norm the hidden values reach stable ranges after few iteration. This help the net to reach high accuracy in less iterations (first subplot on the left) but we can see that even without batch norm the net reach eventually high accuracy.
The only help provided by batch norm to the gradient is the fact that, as noticed before, the normalisation is firstly performed by calculating the mean and variance on individual batches. This is important because this partial estimation of mean and variance introduces noice. Similarly to drop out, which has a regularisation effect due to the noice generated by randomly deactivating part of the weights, batch norm can introduce regularisation by adding noice due to the larger or smaller mean and variance estimated on individual batches. But still, batch norm was not introduced as a regularisation technique, and the equation you put on the question simply prove that it is possible to calculate the derivatives of the equations applied to perform the batch norm.
