Proof that a Lorentz-invariant scalar function can only depend on scalar products

Question

If we have a Lorentz-invariant scalar function $f$ of a single four-vector $x^{\mu}$ we can show that $f$ can only depend on $x^2$ (see Argument of a scalar function to be invariant under Lorentz transformations). I am interested in the case where $f$ can depend on many four vectors $f = f(x_1,...,x_n)$. I want to prove that $f$ can only depend on the scalar products $x_i \cdot x_j$ for $i,j \in \{1,...,n \}$. Proceeding as in Argument of a scalar function to be invariant under Lorentz transformations I can show that the Lorentz-invariant condition gives $$ \sum_{j=1}^n \Big (x_j^{\mu} \frac{\partial}{\partial x_j^{\nu}} - x_j^{\nu} \frac{\partial}{\partial x_j^{\mu}} \Big ) f = 0. $$ However, the same argument as in the $n=1$ case does not seem to work here.
Probably it holds the more general statement that for general coordinate transformations $f$ can depend only on combinations like $T_{\mu_1 \mu_2 ...} x_{j_1}^{\mu_1} x_{j_2}^{\mu_2} ...$, where $T$ is an invariant tensor, i.e. $$ T_{\mu_1 \mu_2 ... } \frac{\partial x^{\mu_1}}{ \partial x'^{\nu_1}} \frac{\partial x^{\mu_2}}{ \partial x'^{\nu_2}} ... = T_{\nu_1 \nu_2 ... }? $$ The Lorentz group case follows by noting $g^{\mu \nu}$ is the only tensor that it is invariant under the full Lorentz group.

Answer 1

The proof, even for the Euclidean form, is surprisingly long: See, for example, M. Spivak, A Comprehensive Introduction to Differential Geometry (second edition) Vol. V, pp. 466-481.