In gradient descent's update rule, why do we use $\sigma(z^{l-1})\frac{\delta C_0}{ \delta w^{l}}$ instead of $\frac{\delta C_0}{\delta w^{l}}$?
Okay - the answer is here  https://explained.ai/matrix-calculus/#sec6.2 and it is pretty involved. Basically, there is a difference when you derive the equation for one neuron and when you have to do practically for a set of neurons. The answer is matrix calculus. Here goes from what I could make out. Feel free to correct if I am wrong
Gradient Vector/Matrix/2D tensor of Loss function wrto Weight
$$
C = \frac{1}{2} \sum_j (y_j-a^L_j)^2
$$
Assuming a neural net with 2 layers, we have the final Loss as
$$
C = \frac{1}{2} \sum_j (y_j-a^2_j)^2
$$
Where
$$
a^2 = \sigma(w^2.a^1)
$$
We can then write
$$
C = \frac{1}{2} \sum_j v^2 \quad \rightarrow (Eq \;A)
$$
Where
$$
v= y-a^2 
$$
Partial Derivative of Loss function wrto Weight
For the last layer, lets use Chain Rule to split like below
$$
\frac {\partial C}{\partial w^2} = \frac{\partial v^2}{\partial v} * \frac{\partial v}{\partial w^2} \quad \rightarrow (Eq \;B)
$$
$$
 \frac{\partial v^2}{\partial v} =2v \quad \rightarrow (Eq \;B.1)
$$
$$
\frac{\partial v}{\partial w^2}=  \frac{\partial (y-a^2)}{\partial w^2} = 
0-\frac{\partial a^2}{\partial w^2} \quad \rightarrow (Eq \;B.2)
$$
$$
\frac {\partial C}{\partial w^2} = \frac{1}{2} *2v(0-\frac{\partial a^2}{\partial w^2}) \quad \rightarrow (Eq \;B)
$$
 
Now we need to find $\frac{\partial a^2}{\partial w^2}$
Let
$$
a^2= \sigma(sum(w^2 \otimes a^1 )) = \sigma(z^2) 
$$
$$
z^2 = sum(w^2 \otimes a^1 )
$$
$$
z^2 = sum(k^2) \; \text {where} \; k^2=w^2 \otimes a^1 
$$
We now need to derive an intermediate term which we will use later
$$
\frac{\partial z^2}{\partial w^2} =\frac{\partial z^2}{\partial k^2}*\frac{\partial k^2}{\partial w^2}
$$
$$
=\frac {\partial sum(k^2)}{\partial k^2}* \frac {\partial (w^2 \otimes a^1 )} {\partial w^2}
$$
$$
\frac{\partial z^2}{\partial w^2} = (1^{\rightarrow})^T* diag(a^1) =(a^{1})^T \quad \rightarrow (Eq \;B.3)
$$
How the above is, you need to check this  in https://explained.ai/matrix-calculus/#sec6.2
Basically though these are written like scalar here; actually all these are partial differention of vector by vector, or vector by scalar; and a set of vectors can be represented as the matrix here.
Note that the vector dot product $w.a$ when applied on matrices becomes the elementwise multiplication $w^2 \otimes a^1$ (also called Hadamard product)
Going back to  $Eq \;(B.2)$
$$
\frac {\partial a^2}{\partial w^2} = \frac{\partial a^2}{\partial z^2} * \frac{\partial z^2}{\partial w^2}
$$
Using $Eq \;(B.3)$ for the term in left
$$
=  \frac{\partial a^2}{\partial z^2} * (a^{1})^T
$$
$$
=  \frac{\partial \sigma(z^2)}{\partial z^2} * (a^{1})^T
$$
$$
\frac {\partial a^2}{\partial w^2} =   \sigma^{'}(z^2) * (a^{1})^T \quad \rightarrow (Eq \;B.4)
$$
Now lets got back to partial derivative of Loss function wrto to weight
$$
\frac {\partial C}{\partial w^2} = \frac {1}{2}*2v(0-\frac{\partial a^2}{\partial w^2}) \quad \rightarrow (Eq \;B)
$$
Using $Eq \;(B.4)$ to substitute in the last term
$$
= v(0- \sigma^{'}(z^2) * (a^{1})^T)
$$
$$
= v*-1*\sigma^{'}(z^2) * (a^{1})^T
$$
$$
= (y-a^2)*-1*\sigma^{'}(z^2) * (a^{1})^T
$$
$$
\frac {\partial C}{\partial w^2}= (a^2-y)*\sigma^{'}(z^2) * (a^{1})^T
$$