How to see from definition that determinant is just volume? The determinant can be defined recursively by Laplace expansion (LE), or equivalently, directly defined by
Lebniz formula (LF)  using parity of permutation.
My question is:

How do we see that these calculations indeed carry out the volumn of the parallelopiped expanded by row vectors of the determinant?

I found some explanations, all of which are not satisfactory for me.

*

*This answer utilizes diagnolization of matrix and then we see that the product of the eigenvalues is indeed the volume.

*This interactive textbook first defines the determinant to be an abstract multilinear function, and then work out the calculation formula, and finally proves by: The only function satisfying the defining properties of determinant is the determinant itself, and the volume of vectors satisfies these properties.

These explanations, without exception, all circumvents the calculation formulas (LE, LF). But I'm wondering this. Determinant is first discovered when solving system of linear equations and defined by LF. How do the earliest researchers, without knowledges of diagnalization, abstract multilinear functions, etc., sees that this definition (LE or LF) is just the volume of vectors?
One more thing. Please do not use the definition of cross product of vectors. Because this operation is also best calculated by formal determinant (as I see it). In fact, this question first came to me when I was thinking: Why the cross product, or the vector orthogonal to two given vectors with length equal to the area of the parallelogram, can be calculated by
$$ \mathbf{v}_1 \times \mathbf{v}_2
= \begin{vmatrix} 
\mathbf{i} & \mathbf{j} & \mathbf{k} \\
x_1 & y_1 & z_1 \\
x_2 & y_2 & z_2
\end{vmatrix}$$
 A: 
How do the earliest researchers, without knowledges of diagnalization,
abstract multilinear functions, etc., sees that this definition (LE or
LF) is just the volume of vectors?

Here's a slide I made that gives a simple explanation for the determinant of a 2 by 2 matrix.

A similar but more laborious argument could be made for the determinant of a 3 by 3 matrix. At that point, we can guess that this interpretation holds for the determinant of an $n \times n$ matrix
A: 
TLDR; You can use row operations on a matrix to convert it to a diagonal matrix, without changing the determinant. Geometrically, this corresponds to shear operations on a parallelepiped, transforming it to cuboid with sides parallel to axes. In the end, the volume of a cuboid is the product of side lengths and the determinant of a diagonal matrix is the product of diagonal values.

Let's assume that we agree with the following statement about $n$-dimensional volume of parallelepiped on vectors $x_1,\ldots,x_n$. If $x_1 = h+e$, where for $h\perp x_i$, $i=2\ldots n$ and $h\perp e$, then
$$V(x_1,\ldots,x_n) = |h|\cdot V(x_2,\ldots,x_n).$$
In other words, the volume is the product of height and base area.

On the other hand, that means that if you add to $x_1$ any linear combination of $x_2,\ldots,x_n$, the volume won't change, since the volume doesn't depend on the vector $e$.
In matrix form it corresponds to the row elimination process:
$$\det M =
\det\begin{pmatrix}x_1\\x_2\\x_3\\\vdots\\x_n\end{pmatrix} = 
\det\begin{pmatrix}x_1\\x_2+\alpha_2x_1\\x_3\\\vdots\\x_n\end{pmatrix}
=\det\begin{pmatrix}x_1\\x_2+\alpha_2x_1\\x_3+\alpha_3x_1\\\vdots\\x_n\end{pmatrix}
=\ldots=\det\begin{pmatrix}x_1\\x_2+\alpha_2x_1\\x_3+\alpha_3x_1\\\vdots\\x_n+\alpha_nx_1\end{pmatrix} = \det M'. \!\,\!\,
$$
Assume, $x_{11}\neq 0$, then you can choose $\alpha_2,\ldots,\alpha_n$, so $x'_{21}=x'_{31}=\ldots=x'_{n1}=0$:
$$       
M' = \begin{pmatrix}
x_{11}&x_{12}&\ldots&x_{1n}\\
0   &x'_{22}&\ldots&x'_{2n}\\
0   &x'_{32}&\ldots&x'_{3n}\\
\ldots&\ldots&\ldots&\ldots\\
0   &x'_{n2}&\ldots&x'_{nn}
\end{pmatrix}
$$
Now it's clear that on the one hand:
$$
\det M' = x_{11}\det \begin{pmatrix}
x'_{22}&\ldots&x'_{2n}\\
x'_{32}&\ldots&x'_{3n}\\
\ldots&\ldots&\ldots\\
x'_{n2}&\ldots&x'_{nn}
\end{pmatrix}
$$
On the other hand, in parallelepiped $M'$, $|h|=|x_{11}|$, so the volume is too $V(x_1,x'_2,\ldots,x'_n)=|x_{11}|V(x'_2,\ldots,x'_n)$. The formulas are the same up to sign. So both for determinant and for volume, the step of induction is the same. And since the base of the induction is also the same (|\det(a)| = |a|), we can conclude that determinant is a signed volume. “Signed” part needs more rigour from our side, but I hope you can follow how volume appears.
