Can't understand the proof of the first backpropagation equation in Nielsen's neural network book I am reading Proof of the four fundamental equations chapter of "Neural Networks and Deep Learning" and I think I got the general idea about backpropagation and the math involved, but there is this part which I can't figure out:

Let's begin with Equation (BP1), which gives an expression for the output error, $\delta^L$. To prove this equation, recall that by definition
\begin{eqnarray}
  \delta^L_j = \frac{\partial C}{\partial z^L_j}.
\tag{36}\end{eqnarray}
Applying the chain rule, we can re-express the partial derivative above in terms of partial derivatives with respect to the output activations,
\begin{eqnarray}
  \delta^L_j = \sum_k \frac{\partial C}{\partial a^L_k} \frac{\partial a^L_k}{\partial z^L_j},
\tag{37}\end{eqnarray}
where the sum is over all neurons $k$ in the output layer. Of course, the output activation $a^L_k$ of the $k^{\rm th}$ neuron depends only on the weighted input $z^L_j$ for the $j^{\rm th}$ neuron when $k=j$. And so $\partial a^L_k / \partial z^L_j$ vanishes when $k≠j$. As a result we can simplify the previous equation to
\begin{eqnarray}
  \delta^L_j = \frac{\partial C}{\partial a^L_j} \frac{\partial a^L_j}{\partial z^L_j}.
\tag{38}\end{eqnarray}

I don't understand why, in the second equation, he is summing over all neurons $k$ in the output layer. Intuitively I'd have obtained the third equation directly from the first one by applying the chain rule, considering only the output activation $a^L_j$. I understand that $C$ is a function depending on all the output activations, so maybe the reason is because in its expression $C = \frac{1}{2} \sum_j (y_j-a^L_j)^2$ there is the summation. However, I feel like I am missing something.
 A: That's because the author is considering in that part that the activation function of a neuron $j$ only depends on the value of $z_j^L$. This happens for example with the sigmoid or ReLU activation functions.
However, when using activation functions like Softmax (described in Chapter 3.1.4), we have that:
$$a_j^L = \frac{e^{z_j^L}}{\sum_k e^{z_k^L}}$$
Here we can see that $\partial a_j^L/\partial z_k^L \neq0$ (or equivalently: $\partial a_k^L/\partial z_j^L \neq0$) for $k\neq j$, so the correct computation of $\delta_j^L$ is the equation given by the author:
$$ \delta_j^L=\sum_k\frac{\partial C}{\partial a_k^L}\frac{\partial a_k^L}{\partial z_j^L} $$
So, to sum up, the above equation is the general expresion of $\delta_j^L$, however if we use certain activation functions, it can be simplified to the equation $(38)$ of the question.
A: Yes, in this case, it's obvious. The author is simply writing a general equation for taking derivatives in a multivariate setting. If $f(x,y)$ is a function of $x,y$ and we look for the derivative of $f$ wrt $t$, it's customary to write the following:
$$\frac{\partial f}{\partial t}=\frac{\partial f}{\partial x}\frac{\partial x}{\partial t}+\frac{\partial f}{\partial y}\frac{\partial y}{\partial t}$$
