What motivates that strange definition of the determinant?

So historically, Leibniz first came up with the following definition of determinant while trying to derive a solvability condition for system of linear equations.
Afterwards, Laplace came up with the Laplace expansion formula.
Here, I am going to derive the Laplace expansion formula by solving a system of 3x3 equations "inducting" on the determinant definition for a 2x2 system.

We all know that determinant for a 2x2 matrix \begin{bmatrix} a & b \\ c & d \end{bmatrix} the determinant is $ad-bc$. This definition naturally arises when you try to row-reduce the system.
The inverse matrix (and hence solution to the associated system) is given by: $$ \frac{1}{ad-bc} \begin{bmatrix} d & -c \\ -b & a \end{bmatrix} $$ Let's try to find the definition of determinant for 3x3 matrix by row-reduction and re-using this result here.
So we start with: $$ M = \begin{bmatrix} a & b & c\\ d & e & f\\ g & h & i\\ \end{bmatrix} $$ We can apply row reduction. At each step, i am going to give the resulting matrix and also the "elementary matrix" which encodes the row reduction operation. Namely, I give you a matrix \(E\) such that when you perform \(EM\) you get the result of the row-reduction.

Multiply row 2 with \(a\) and write to row 2. First I give you the resulting matrix followed by the elementary matrix: $$ \begin{bmatrix} a & b & c\\ ad & ae & af\\ g & h & i\\ \end{bmatrix} \qquad \begin{bmatrix} 1 & 0 & 0\\ 0 & a & 0\\ 0 & 0 & 1\\ \end{bmatrix} $$

new row 3 = old row 3 * a. Remember that I am giving you the "cumulative" matrix and the "individual" elementary matrices. $$ \begin{bmatrix} a & b & c\\ ad & ae & af\\ ag & ah & ai\\ \end{bmatrix} \qquad \begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & a\\ \end{bmatrix} $$

new row 3 = old row 3 - ((row1) * g) $$ \begin{bmatrix} a & b & c\\ ad & ae & af\\ 0 & ah-bg & ai-cg\\ \end{bmatrix} \qquad \begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ -g & 0 & 1\\ \end{bmatrix} $$

new row 2 = old row 2 - ((row1) * d) $$ \begin{bmatrix} a & b & c\\ 0 & ae-bd & af-cd\\ 0 & ah-bg & ai-cg\\ \end{bmatrix} \qquad \begin{bmatrix} 1 & 0 & 0\\ -d & 1 & 0\\ 0 & 0 & 1\\ \end{bmatrix} $$ Now you are seeing some "determinant-like" expressions. These are exactly the determinant of adjugate matrices. For example, I can define \(M_{33}\) as the determinant of the matrix obtained by removing row 3 and column 3 from the original matrix \(M\): \(M_{33} = ae - bd\).

Using this new "determinant of adjugate sub-matrices" definition, we have row-reduced our original matrix \(M\) to \begin{bmatrix} a & b & c\\ 0 & M_{33} & M_{32}\\ 0 & M_{23} & M_{22}\\ \end{bmatrix}

Our cumulative elementary operation matrix is \begin{bmatrix} 1 & 0 & 0\\ -d & a & 0\\ -g & 0 & a\\ \end{bmatrix}

So now, we can think of the the right bottom sub-matrix (call it \(S\)) and use the 2x2 inverse formula to get the following elementary operation matrix (bad name since it is not really elementary): \begin{bmatrix} 1 & 0 & 0\\ 0 & \frac{1}{\det(S)}M_{22} & \frac{-1}{\det(S)}M_{32}\\ 0 & \frac{-1}{\det(S)}M_{23} & \frac{1}{\det(S)}M_{33}\\ \end{bmatrix}

This will take the original matrix further to a simple \begin{bmatrix} a & b & c\\ 0 & 1 & 0\\ 0 & 0 & 1\\ \end{bmatrix}

at this point, if we "left-apply" the elementary matrix \begin{bmatrix} 1/a & -b/a & -c/a\\ 0 & 1 & 0\\ 0 & 0 & 1\\ \end{bmatrix} you will row-reduce \(M\) to identity.

but now we multiply all the elementary matrix to get the whole complete inverse. $$ \begin{bmatrix} 1/a & -b/a & -c/a\\ 0 & 1 & 0\\ 0 & 0 & 1\\ \end{bmatrix} \begin{bmatrix} 1 & 0 & 0\\ 0 & \frac{1}{\det(S)}M_{22} & \frac{-1}{\det(S)}M_{32}\\ 0 & \frac{-1}{\det(S)}M_{23} & \frac{1}{\det(S)}M_{33}\\ \end{bmatrix} \begin{bmatrix} 1 & 0 & 0\\ -d & a & 0\\ -g & 0 & a\\ \end{bmatrix} $$

Lots of really nice cancellations happen, for example \(gM_{32}-dM_{22} = -aM{12}\). Very soon, you should see something like this $$ \frac{a}{\det(S)}\begin{bmatrix} M_{11} & -M_{21} & M_{31}\\ -M_{12} & M_{22} & -M_{32}\\ M_{13} & -M_{23} & M_{33}\\ \end{bmatrix} $$ Simplify the scalar coefficient you get the 3x3 inverse and the 3x3 determinant formula in Laplace expansion form. This obviously can be generalized inductively to n-by-n.