Contexts and definitions of the $LU$ method:
Consider a matrix $A \in \mathbb{R}^{n \times n}$:
$$ A = \begin{bmatrix}a_{11} & \dots & a_{1n} \\ \vdots & \ddots & \vdots \\ a_{n1} & \dots & a_{nn} \end{bmatrix} $$
We define: $$ \ell_{i,j} := \frac{a_{ij}}{a_{jj}}, \,\,\, m_j = \begin{bmatrix} 0 \\ \vdots \\ 0 \\ \ell_{j+1, j} \\ \ell_{j+2,j} \\ \vdots \\ \ell_{n,j} \end{bmatrix}, e_j = \begin{bmatrix} 0 \\ \vdots \\ 0 \\ 1 (j_{th} \text{ line}) \\ 0 \\ \vdots \\ 0 \end{bmatrix}, M_j = I - m_je^T_j$$
$M = M_{n-1}P_{n-1} \dots M_1P_1 $, where each $ P_k$ is a permutation matrix and $P = P_{n-1} \dots P_1$.
The $LU$ method with pivoting says: $$ PA = LU \\ U = MA \text{ is upper triangular}\\ L = PM^{-1} \text{ is unit lower triangular}$$
I want to proof the following: $$ L = PM^{-1} \text{ is unit lower triangular}$$
Reference: this is an exercice left to the reader of the book "Numerical Mathematics" by Alfio Quarteroni, second edition (page 89).
Given that we are talking about LU factorisation, I assume that the permutation matrix $P_i$ exchanges line $j$ and $k$ with $j, k \ge i$ (Typically $j=i$ and $k>i$).
With this assumption, we can rewrite the matrix $M$ in the following form:
$$ M = M_{n-1}P_{n-1} M_{n-2}P_{n-2} \dots M_2P_2M_1P_1 = L_{n-1}L_{n-2}\dots L_2L_1 \cdot P_{n-1}P_{n-2}\dots P_2P_1$$
with $$L_{n-1} = M_{n-1}, \quad L_{n-2} = P_{n-1}M_{n-2}P_{n-1}^{-1},\quad L_{n-3} = P_{n-1}P_{n-2}M_{n-3}P_{n-2}^{-1}P_{n-1}^{-1}, \quad \text{etc}$$
thus having $M_i$ equal to $L_i$, but with the sub-diagonal nonzero entries permuted by $\dots P_{i+2}P_{i+1}M_i$. This is a result of the assumption that $P_i$ does not change the order of the rows before row $i$. If this hypothesis is not respected, then $L$ is no longer unit triangular. It follows that the matrices $L_i$ are lower triangular and unit diagonal, and thus so is their product $L^{-1}$.
$$ M = \underbrace{L_{n-1}L_{n-2}\dots L_2L_1}_ {L^{-1}} \cdot \underbrace{P_{n-1}P_{n-2}\dots P_2P_1}_{P}$$
The inverse of a unit diagonal triangular matrix is also a unit triangular matrix which implies that $L$ is unit triangular.
Complementary note N°1: Here is an example (with n = 4) showing that the $M$ is actually equal to $L_{n-1}L_{n-2}\dots L_2L_1 \cdot P_{n-1}P_{n-2}\dots P_2P_1$, as it might be useful to convince yourself:
Using the definition for $L_i$:
$$L_3 L_2 L_1 \cdot P_3 P_2 P_1 = M_3 (P_3 M_2 P_3^{-1})(P_3 P_2 M_1 P_2^{-1} P_3^{-1}) \cdot P_3 P_2 P_1$$
cancelling the permutation matrices gives the expected matrix $M$:
$$L_3 L_2 L_1 \cdot P_3 P_2 P_1 = M_3 P_3 M_2 P_2 M_1 P_1 = M$$
Complementary note N°2: Example showing why $L_i$ and $M_i$ have the same structure. Let's have a look at $L_2$ in the same example: $$ L_2 = P_3 M_2 P_3^{-1}$$ with $M_2$ being of the form \begin{pmatrix} 1 & & & \newline & 1 & & \newline & x & 1& \newline & x && 1 \end{pmatrix}
($x$ means some nonzero entry). $P_3$ will shuffle rows 3 and 4, thus $P_3 L_2$ will be of the form: \begin{pmatrix} 1 & & & \newline & 1 & & \newline & x & 1 \text{ or } 0 & 1 \text{ or } 0\newline & x & 1 \text{ or } 0 & 1 \text{ or } 0 \end{pmatrix}
$P_3^{-1}$ will then shuffle columns 3 and 4 in an opposite way, resetting the 2x2 lower right bloc back to an identity sub-matrix.
Edit: if you want more details, I would suggest reading lectures 20 and 21 of the book "Numerical linear algebra" by Lloyd N. Trefethen & David Bau III, the whole algorithm for LU is explained and it's surprisingly pleasant to read