Assume $\mathbf{x}$ to be a lightlike 4-vector in reference Frame $F$ and $\mathbf{x'}$ the corresponding vector in reference frame $F'$ where both frames share the same origin. Then due to the principle of relativity:
$$(x_0)^2-(x_1)^2-(x_2)^2-(x_3)^2 = (x_0')^2-(x_1')^2-(x_2')^2-(x_3')^2 = 0$$
According to Dodson/Poston (page 16), from this "it follows fairly easily that there is a positive number $a$ such that for any [vector] labelled $(x_0, x_1, x_2, x_3)$ by one system and $(x_0', x_1', x_2', x_3')$ by the other, not just [lightlike vectors], we have
$$\tag{*}(x_0)^2-(x_1)^2-(x_2)^2-(x_3)^2 = a((x_0')^2-(x_1')^2-(x_2')^2-(x_3')^2)$$"
I just changed the nomenclature a bit and set $c=1$ but otherwise cited the original text. The idea is to derive some basic property of Minkowski space / Lorentz transformations directly from the constancy of the speed of light.
I tried to restate the problem like this:
- Let $f: \mathbb{R}^4 \rightarrow \mathbb{R}: \mathbf{x} \mapsto (x_0)^2 - (x_1)^2 - (x_2)^2 - (x_3)^2$ (the Lorentz form)
- Let S be the set of all mappings $\mathbb{R}^4 \rightarrow \mathbb{R}^4$ such that for every $L \in S$, $L(0) = 0$ and $f(\mathbf{x}) = 0 \Rightarrow f(L(\mathbf{x})) = 0$ (constancy of the speed of light in different reference frames)
- It should then be possible to derive from (1) and (2) that for every $\mathbf{x} \in \mathbb{R}^4$ and a given $L \in S$: $$f(\mathbf{x}) = a f(L(\mathbf{x}))$$ with $a$ some constant positive real number.
EDIT: As has been pointed out in the comments, the restrictions I place on $L$ are not sufficient to derive $(*)$. I cannot think of more conditions that I could directly extract from the text, though. Anyone can help to come up with a minimum of conditions that would allow the deduction proposed by Dodson/Poston?
Instead of the quadratic form $f$, use its related bilinear form $B({\mathbf x},{\mathbf y}) = x_0y_0 - x_1y_1 - x_2y_2 - x_3y_3$, which has more flexibility since it has two vector inputs instead of one.
Theorem. Let $B$ be a nondegenerate bilinear form on a finite-dimensional real vector space $V$. The following conditions on a linear transformation $L \colon V \rightarrow V$ are equivalent:
a) $B({\mathbf x},{\mathbf y}) = 0 \Longrightarrow B(L{\mathbf x},L{\mathbf y}) = 0$. That is, $L$ preserves orthogonality relative to $B$.
b) There is a real number $c$ such that $B(L{\mathbf x},L{\mathbf y}) = cB({\mathbf x},{\mathbf y})$ for all ${\mathbf x}$ and ${\mathbf y}$ in $V$.
Proof: See Theorem 3.20 at http://www.math.uconn.edu/~kconrad/blurbs/linmultialg/bilinearform.pdf, which is carried out with vector spaces over a general scalar field, and you can take the scalars to be the real numbers.
In this theorem $L$ is assumed to be linear, and that is missing from the setup in your problem. The geometric content of this theorem is that a linear transformation preserves orthogonality relative to $B$ exactly when it changes $B$-values by a uniform scaling factor. (Your $a$ is my $1/c$.)
To connect this theorem to your question, let $Q({\mathbf x}) = B({\mathbf x},{\mathbf x})$ be the quadratic form associated to $B$. (So $Q = f$ in your notation.) Then condition a implies that $Q({\mathbf x}) = 0 \Longrightarrow Q(L{\mathbf x}) = 0$ by using ${\mathbf y} = {\mathbf x}$ in a, but that implication about where $Q$ takes the value $0$ doesn't imply condition a in general.
Condition b in the theorem, however, is equivalent to $Q(L{\mathbf x}) = cQ({\mathbf x})$ for all ${\mathbf x}$ in $V$. In one direction, setting ${\mathbf y} = {\mathbf x}$ in condition b gives us $Q(L{\mathbf x}) = cQ({\mathbf x})$ for all ${\mathbf x}$. In the other direction, by the formula $$ B({\mathbf x},{\mathbf y}) = \frac{1}{2}(Q({\mathbf x} + {\mathbf y}) - Q({\mathbf x}) - Q({\mathbf y})), $$ if $Q(L{\mathbf x}) = cQ({\mathbf x})$ for all ${\mathbf x}$ then condition b is true by using $L{\mathbf x}$ in place of $\mathbf x$ and $L{\mathbf y}$ in place of $\mathbf y$ in the formula for $B({\mathbf x},{\mathbf y})$ in terms of values of $Q$; we need the additivity of $L$ for that, which is a special instance of being linear.
On a related note, a function $V \rightarrow V$ that preserves all $B$-values (not just $B$-value $0$) must be linear. Here is that result.
Theorem. Let $B$ be a nondegenerate bilinear form on a finite-dimensional vector space $V$. If $\sigma \colon V \rightarrow V$ satisfies $B(\sigma({\mathbf x}),\sigma({\mathbf y})) = B({\mathbf x},{\mathbf y})$ for all ${\mathbf x}$ and ${\mathbf y}$ in $V$, then $\sigma$ is linear and must be injective and surjective.
Proof: See my answer at Isometry without injection and surjection.
There is definitely no theorem about the condition $Q(\sigma{\mathbf x}) = Q({\mathbf x})$ for all ${\mathbf x}$ implying $\sigma$ is linear, since such a condition is satisfied if $\sigma$ is a random permutation in each level set of $Q$ (the sets where $Q$ takes a common value).