In a continuum mechanical context, let $\mathbf{A}$ be a second order tensor function.
We consider our operations in an orthonormal context
$$\mathbf{A} = A_{ij}\, \mathbf{e}_i\otimes\mathbf{e}_j\tag{1}$$
How do I calculate the derivatives
$$\cfrac{\partial\mathbf{A}}{\partial\mathbf{A}} = \cfrac{\partial A_{ij}}{\partial A_{kl}} \tag{2}$$
and
$$\cfrac{\partial\mathbf{AA}}{\partial\mathbf{A}} = \cfrac{\partial A_{ij}\,A_{jk}}{\partial A_{lm}}\tag{3}$$
or namely, is there an algorithmic basis on which I can take the derivatives of arbitrary indicial expressions such as above?
The thing is, I already know the answers from literature. For (1), the answer is simply given as $\delta_{ik}\delta_{jl}$, but I don't know exactly how that is derived.
For (2), according to the literature, the product rule should be utilized such as
\begin{align*} \cfrac{\partial A_{ij}\,A_{jk}}{\partial A_{lm}}&=\cfrac{\partial A_{ij}}{\partial A_{lm}}A_{jk}+A_{ij}\cfrac{\partial A_{jk}}{\partial A_{lm}}\\ &=\delta_{il}\delta_{jm}A_{jk}+A_{ij}\delta_{jl}\delta_{km}\\ &=(\delta_{il}A_{mk}+A_{il}\delta_{km})\,\,\, \text{with the basis}\,\,\, \mathbf{e}_i\otimes\mathbf{e}_k\otimes\mathbf{e}_l\otimes\mathbf{e}_m\tag{4} \end{align*}
But how do we prove that the product rule yields a valid result in the context of these multi-dimensional objects? Can I define an algorithm for the differential operator that doesn't expand using the product rule but somehow using the indices?
Just to give a different perspective, let's assume temporarily that $\mathbf{A}$ is a simple tensor consisting of the dyadic product of two vector functions
$$\mathbf{A} = \mathbf{a} \otimes \mathbf{b} = a_i b_j \,\,\mathbf{e}_i\otimes\mathbf{e}_j\tag{5}$$
Then (3) would be of the form
$$\cfrac{\partial a_i b_j a_j b_k}{\partial a_l b_m}\tag{6}$$
What types of strategies can be used when approaching this? Note that I am assuming that these objects are completely separate from geometry, i.e. just multi-dimensional arrays of scalar functions.