[proofplan]
Represent the Wishart matrix as $W=X^\top X$, where the columns of $X$ are independent standard Gaussian vectors in $\mathbb{R}^n$. Apply Gram-Schmidt to these columns and write $X=QR$, with $Q$ having orthonormal columns and $R$ upper triangular with positive diagonal; then $W=R^\top R$, so $A=R^\top$ is the desired Cholesky factor. The distributional content follows from rotational invariance of the standard Gaussian law: at the $j$th step, the Gaussian coordinates in the already constructed orthonormal directions are independent $\mathcal{N}(0,1)$ variables, while the squared norm of the orthogonal residual is chi-squared with $n-j+1$ degrees of freedom and independent of those coordinates.
[/proofplan]
[step:Construct the triangular factor by Gram-Schmidt]
Let
\begin{align*}
X_1,\dots,X_p:\Omega \to \mathbb{R}^n
\end{align*}
denote the column random vectors of $X$. Let $\mathcal{N}_n(0,I_n)$ denote the standard Gaussian probability law on $\mathbb{R}^n$, meaning the law of a random vector whose $n$ coordinates in the standard basis are independent $\mathcal{N}(0,1)$ random variables. Since the entries of $X$ are independent $\mathcal{N}(0,1)$ random variables, the vectors $X_1,\dots,X_p$ are independent and each has distribution $\mathcal{N}_n(0,I_n)$.
We recursively define orthonormal vectors $q_1,\dots,q_p$ and upper triangular entries $r_{ij}$ for $1 \leq i \leq j \leq p$. Set
\begin{align*}
r_{11} := |X_1|, \qquad q_1 := \frac{X_1}{r_{11}}.
\end{align*}
For $j \geq 2$, assuming $q_1,\dots,q_{j-1}$ have been defined, define
\begin{align*}
r_{ij} &:= q_i \cdot X_j, && 1 \leq i < j, \\
Y_j &:= X_j - \sum_{i=1}^{j-1} r_{ij} q_i, \\
r_{jj} &:= |Y_j|, \\
q_j &:= \frac{Y_j}{r_{jj}}.
\end{align*}
Because $X_1,\dots,X_j$ have a joint density on $\mathbb{R}^{nj}$ and $j \leq p \leq n$, the event that $X_1,\dots,X_j$ are linearly dependent has probability zero. Hence $r_{jj}>0$ for every $j$ almost surely, so the construction is valid almost surely.
Let
\begin{align*}
Q:\Omega \to \mathbb{R}^{n \times p}
\end{align*}
be the random matrix with columns $q_1,\dots,q_p$, and let
\begin{align*}
R:\Omega \to \mathbb{R}^{p \times p}
\end{align*}
be the upper triangular random matrix with entries $R_{ij}=r_{ij}$ for $i \leq j$ and $R_{ij}=0$ for $i>j$. The recursion gives
\begin{align*}
X_j = \sum_{i=1}^j r_{ij} q_i
\end{align*}
for each $1 \leq j \leq p$, so
\begin{align*}
X = QR.
\end{align*}
Since the columns of $Q$ are orthonormal, $Q^\top Q=I_p$. Therefore
\begin{align*}
W = X^\top X = R^\top Q^\top Q R = R^\top R.
\end{align*}
Let $E \subset \Omega$ be the full-probability event on which the Gram-Schmidt construction above is valid for every $1 \leq j \leq p$. On $\Omega \setminus E$, define each $q_j$ to be the deterministic standard basis vector $u_j \in \mathbb{R}^n$ and define each $r_{ij}$ with $1 \leq i \leq j \leq p$ to be $0$; this gives everywhere-defined random variables without changing any almost-sure identity or distributional conclusion. Define
\begin{align*}
A:\Omega \to \mathbb{R}^{p \times p}
\end{align*}
by
\begin{align*}
A(\omega) :=
\begin{cases}
R(\omega)^\top, & \omega \in E, \\
I_p, & \omega \in \Omega \setminus E.
\end{cases}
\end{align*}
Then $A$ is a lower triangular random matrix on all of $\Omega$, its diagonal entries are positive on $E$, and
\begin{align*}
W=AA^\top
\end{align*}
on $E$, hence almost surely.
[/step]
[step:Identify the distribution of one Gram-Schmidt step]
Fix $j \in \{1,\dots,p\}$. Let $\mathcal{F}_{j-1}$ denote the $\sigma$-algebra generated by $X_1,\dots,X_{j-1}$; for $j=1$, let $\mathcal{F}_0 := \{\varnothing,\Omega\}$. On the full-probability event where $q_1,\dots,q_{j-1}$ are defined, extend this orthonormal family to an orthonormal basis
\begin{align*}
q_1,\dots,q_{j-1},e_j,\dots,e_n
\end{align*}
of $\mathbb{R}^n$, where each $e_k$ is chosen as an $\mathcal{F}_{j-1}$-measurable vector. To construct the completion measurably, process the fixed deterministic basis $u_1,\dots,u_n$ of $\mathbb{R}^n$ in order, project each $u_m$ onto the orthogonal complement of the span of the vectors already selected, and keep the first vector whose projected norm is nonzero. The ordering gives a fixed tie-breaking rule on the full-probability rank event; on its complement, define the remaining $e_k$ arbitrarily, for instance by the same deterministic basis. Each projection and norm is a Borel function of $q_1,\dots,q_{j-1}$, so the selected vectors are $\mathcal{F}_{j-1}$-measurable.
Define the $\mathcal{F}_{j-1}$-measurable orthogonal matrix
\begin{align*}
O_j:\Omega \to \mathbb{R}^{n \times n}
\end{align*}
whose rows are $q_1,\dots,q_{j-1},e_j,\dots,e_n$ in that order. Since $X_j$ is independent of $\mathcal{F}_{j-1}$ and has distribution $\mathcal{N}_n(0,I_n)$, conditioning on $\mathcal{F}_{j-1}$ fixes the matrix $O_j$. For every deterministic orthogonal matrix $O \in \mathbb{R}^{n \times n}$, the transformed vector $OX_j$ again has law $\mathcal{N}_n(0,I_n)$ by rotational invariance of the standard Gaussian law. Therefore, conditional on $\mathcal{F}_{j-1}$, the coordinates of $X_j$ in this orthonormal basis are independent standard normal random variables. Define these coordinate random variables by
\begin{align*}
Z_i &:= q_i \cdot X_j, && 1 \leq i < j, \\
Z_k &:= e_k \cdot X_j, && j \leq k \leq n.
\end{align*}
Then, conditional on $\mathcal{F}_{j-1}$, the family
\begin{align*}
Z_1,\dots,Z_{j-1},Z_j,\dots,Z_n
\end{align*}
consists of independent $\mathcal{N}(0,1)$ random variables.
By construction, for $1 \leq i<j$,
\begin{align*}
r_{ij}=q_i \cdot X_j=Z_i.
\end{align*}
The residual vector is the [orthogonal projection](/theorems/437) of $X_j$ onto the orthogonal complement of $\operatorname{span}\{q_1,\dots,q_{j-1}\}$:
\begin{align*}
Y_j = \sum_{k=j}^n Z_k e_k.
\end{align*}
Hence
\begin{align*}
r_{jj}^2 = |Y_j|^2 = \sum_{k=j}^n Z_k^2.
\end{align*}
Therefore, conditional on $\mathcal{F}_{j-1}$,
\begin{align*}
r_{1j},\dots,r_{j-1,j}
\end{align*}
are independent $\mathcal{N}(0,1)$ random variables, $r_{jj}^2$ has distribution $\chi^2_{n-j+1}$, and $r_{jj}$ is independent of $r_{1j},\dots,r_{j-1,j}$.
[guided]
Fix the column index $j$. The information from the previous columns is encoded in the $\sigma$-algebra $\mathcal{F}_{j-1}$ generated by $X_1,\dots,X_{j-1}$. Once this information is fixed, the vectors $q_1,\dots,q_{j-1}$ are deterministic orthonormal vectors.
We complete $q_1,\dots,q_{j-1}$ to an orthonormal basis
\begin{align*}
q_1,\dots,q_{j-1},e_j,\dots,e_n
\end{align*}
of $\mathbb{R}^n$. The vectors $e_j,\dots,e_n$ are chosen measurably from the previous data, so conditioning on $\mathcal{F}_{j-1}$ makes the entire basis fixed.
The crucial input is rotational invariance of the standard Gaussian vector. Define the $\mathcal{F}_{j-1}$-measurable orthogonal matrix
\begin{align*}
O_j:\Omega \to \mathbb{R}^{n \times n}
\end{align*}
whose rows are $q_1,\dots,q_{j-1},e_j,\dots,e_n$. Since $X_j$ is independent of $\mathcal{F}_{j-1}$ and $X_j\sim \mathcal{N}_n(0,I_n)$, conditioning on $\mathcal{F}_{j-1}$ fixes $O_j$ while leaving $X_j$ with standard Gaussian law. For each deterministic orthogonal matrix $O \in \mathbb{R}^{n \times n}$, rotational invariance gives $OX_j\sim \mathcal{N}_n(0,I_n)$, so its coordinates are independent $\mathcal{N}(0,1)$ random variables. Thus, conditional on $\mathcal{F}_{j-1}$, the coordinate variables
\begin{align*}
Z_i &:= q_i \cdot X_j, && 1 \leq i < j, \\
Z_k &:= e_k \cdot X_j, && j \leq k \leq n
\end{align*}
are independent standard normal random variables.
The Gram-Schmidt coefficients before the diagonal are exactly the coordinates of $X_j$ in the already constructed directions:
\begin{align*}
r_{ij}=q_i\cdot X_j=Z_i,\qquad 1\leq i<j.
\end{align*}
The residual removes precisely those coordinates, so it lies in the orthogonal complement:
\begin{align*}
Y_j = X_j-\sum_{i=1}^{j-1}(q_i\cdot X_j)q_i
= \sum_{k=j}^n Z_k e_k.
\end{align*}
Taking the Euclidean norm and using orthonormality gives
\begin{align*}
r_{jj}^2=|Y_j|^2=\sum_{k=j}^n Z_k^2.
\end{align*}
A sum of squares of $n-j+1$ independent standard normal variables has distribution $\chi^2_{n-j+1}$. Because the variables $Z_1,\dots,Z_{j-1}$ are independent of $Z_j,\dots,Z_n$, the off-diagonal coefficients $r_{1j},\dots,r_{j-1,j}$ are independent of the diagonal length $r_{jj}$.
[/guided]
[/step]
[step:Promote the conditional distributions to joint independence]
For each $j$, the preceding step gives the conditional distribution of the new block
\begin{align*}
B_j := (r_{1j},\dots,r_{j-1,j},r_{jj})
\end{align*}
given $\mathcal{F}_{j-1}$. This conditional distribution does not depend on the realised values of $X_1,\dots,X_{j-1}$. Moreover, within $B_j$, the variables $r_{1j},\dots,r_{j-1,j}$ and $r_{jj}$ are independent, with
\begin{align*}
r_{ij} &\sim \mathcal{N}(0,1), && 1 \leq i<j, \\
r_{jj}^2 &\sim \chi^2_{n-j+1}.
\end{align*}
Since the conditional law of $B_j$ given $\mathcal{F}_{j-1}$ is fixed and does not depend on $\mathcal{F}_{j-1}$, the block $B_j$ is independent of $\mathcal{F}_{j-1}$. Because $B_1,\dots,B_{j-1}$ are functions of $X_1,\dots,X_{j-1}$, the block $B_j$ is independent of $B_1,\dots,B_{j-1}$. Induction over $j=1,\dots,p$ shows that all entries
\begin{align*}
\{r_{ij}:1\leq i\leq j\leq p\}
\end{align*}
are mutually independent, with the distributions just stated.
[/step]
[step:Translate the upper triangular QR factor into the lower triangular Bartlett factor]
Recall that $A=R^\top$. Hence, for $1 \leq j \leq p$,
\begin{align*}
a_{jj}=r_{jj},
\end{align*}
and for $1 \leq j<i\leq p$,
\begin{align*}
a_{ij}=r_{ji}.
\end{align*}
The independence of the triangular entries of $R$ therefore gives the independence of the lower triangular entries of $A$. The diagonal distributions become
\begin{align*}
a_{jj}^2=r_{jj}^2\sim \chi^2_{n-j+1},
\qquad 1\leq j\leq p,
\end{align*}
and the strictly lower triangular distributions become
\begin{align*}
a_{ij}=r_{ji}\sim \mathcal{N}(0,1),
\qquad 1\leq j<i\leq p.
\end{align*}
Together with $W=AA^\top$ almost surely, this proves the Bartlett decomposition.
[/step]