Is there a geometric reason for why every map has a divergence-free Jacobian cofactor matrix?

$\newcommand{\Cof}{\text{Cof}}$

Consider the following claim:

Let $f:\mathbb{R}^n \to \mathbb{R}^n$ be a $C^2$ map. Then $$ \text{div} (\Cof df)=0, $$ where $\Cof df$ is the cofactor matrix of $df$, and the divergence is taken row-by-row.

In other words $$ \sum_{j=1}^n \frac{\partial(Cof(Du))_{kj}}{\partial x_j} = 0,$$ for every $1 \le k \le n$.

This is proved in Evan's PDE book, in section 8.1 of the Calculus of Variation.

This identity is "universal", that is, it is satisfied by any smooth (or $C^2$) map $\mathbb{R}^n \to \mathbb{R}^n$.

I follow his proof, but is there an intuitive/geometric way to see why this should be true? I seems like an arbitrary identity - how would someone know this?

$\newcommand{\Cof}{\operatorname{Cof}}$ $\newcommand{\R}{\mathbb{R}}$

My previous answer was insufficient, and did not really prove the required identity. The following is an updated, completely new answer.

Here is a geometric way to think about this result. A convenient way for describing it is to use differential forms:

First, a lemma:

Lemma 1: Let $f_0, f_1:M\rightarrow N$ smooth homotopic maps between smooth manifolds $M$ and $N$. Suppose $M$ is compact with no boundary. Then, for every closed form $\omega\in \Omega^m(N)$ (where $m=\dim M$), $$\int_{M}f_0^*\omega=\int_Mf_1^*\omega.$$

The proof is essentially Stokes theorem.

This lemma implies that the functional $E:C^{\infty}(M,N) \to \mathbb{R}$ defined by $$ E(f)=\int_{M}f^*\omega$$ is a null-Lagrangian, i.e. every smooth map is a critical point of it.

In, particular we have the following corollary:

Let $\Omega \subseteq \R^n$. Then $E:C^{\infty}(\Omega,\R^n) \to \mathbb{R}$ defined by $$ E(f)=\int_{\Omega} f^*\text{Vol}_{\R^n}=\int_{\Omega} \det df \text{Vol}_{\R^n}$$ is a null-Lagrangian.

Now, it is easy to prove that the Euler-Lagrange equation of $E$ is exactly $\text{div} \Cof (df)=0,$ when the divergence is taken row-wise. Indeed,

Let $f_t=f+tV$ for some smooth variation field $V:\Omega \to \R^n$. Then $$ \left. \frac{d}{dt} \right|_{t=0} E(f_t)=\int_{\Omega} \left. \frac{d}{dt} \right|_{t=0} \det df_t=\int_{\Omega} \left. \frac{d}{dt} \right|_{t=0} \det (df+tdV)=\int_{\Omega} \langle \Cof df,dV \rangle=\int_{\Omega} \langle \text{div} \Cof df,V \rangle,$$

where we used two facts:

1. The cofactor is the gradient of the determinant.
2. The divergence is the adjoint of the gradient $d$. (This is essentially Stokes theorem or integration by parts).

Thus, a map is critical if and only if $$\text{div} \Cof df=0. \tag{1}$$

Now, since any map $f:\Omega \to \R^n$ is critical (by lemma 1) this proves that identity $(1)$ holds universally.

Further discussion:

It turns out that there is nothing special in the Euclidean case; everything here can be generalized to arbitrary Riemannian manifolds:

Let $M,N$ be oriented $d$-dimensional Riemannian manifolds. Then $$E(f)=\int_{M}f^*\text{Vol}_N$$ is a null-Lagrangian. Its Euler-Lagranges equation is $\delta (\Cof df)=0$ which is satisfied by any smooth map $f:M \to N$.

Here $\delta$ and $\Cof df$ are natural generalizations of the cofactor matrix and the divergence operator to Riemannian settings:

1. $\Cof df:TM \to f^*(TN)$ is the cofactor map of $df$ defined by $$ \Cof df= (-1)^{d-1} \star_{f^*TN}^{d-1} (\wedge^{d-1} df) \star_{TM}^1. $$

$\Cof df$ represents the action of $df$ on $d-1$ dimensional cubes. ($\star$ denotes the Hodge-dual operator, and the subscripts indicate w.r.t which metric it is taken).

2. Note that $\Cof df \in \Omega^1(M,f^*TN)$. $$\delta:\Omega^1(M,f^*TN) \to \Omega^0(M,f^*TN)=\Gamma(f^*TN)$$ is the adjoint of the pullback connection $f^*\nabla^{TN}$ of the Levi-Civita connection of $N$.

The full details of the derivation of the E-L equation in the Riemannian setting can be found in my paper, in lemma 2.9.

A non-variational approach for deriving this identity (as well as more information on the variational approach) can be found in this paper of mine.