  Starting with the solution to the Bethe-Salpeter equation for the pion, in a
beyond rainbow-ladder truncation to QCD's Dyson-Schwinger equations (DSEs), we
determine the pion's $l_z=0$ and $|l_z|=1$ leading Fock-state light-front wave
functions (LFWFs) [labeled by $\psi_{l_z}(x,k_T^2)$]. The leading-twist
time-reversal even transverse momentum dependent parton distribution function
(TMD) of the pion is then directly obtained from these LFWFs. A key
characteristic of the LFWFs, which is driven by dynamical chiral symmetry
breaking, is that at typical hadronic scales they are broad functions in the
light-cone momentum fraction $x$. The LFWFs have a non-trivial $(x,k_T^2)$
dependence and in general do not factorize into separate functions of each
variable. The $l_z=0$ LFWF is concave with a maximum at $x=1/2$, whereas
orbital angular momentum effects causes the $|l_z|=1$ LFWF to have a slight
{\it double-humped} structure for quark transverse momentum in the range $0.5
\lesssim k_T^2 \lesssim 5\,$GeV$^2$. For $k_T^2 \lesssim 1\,$GeV$^2$ the
$k_T^2$ dependence of the LFWFs is well described by a Gaussian, however for
$k_T^2 \gtrsim 10\,$GeV$^2$ these LFWFs behave as $\psi_0 \propto x(1-x)/k_T^2$
and $\psi_1 \propto x(1-x)/k_T^4$, and therefore exhibit the power-law behavior
predicted by perturbative QCD. The pion's TMD inherits many features from the
LFWFs, where for $k_T^2 \lesssim 1\,$GeV$^2$ the $k_T^2$ dependence is well
described by a Gaussian, and for large $k_T^2$ the TMD behaves as $f_\pi^q
\propto x^2(1-x)^2/k_T^4$. At the model scale we find the average transverse
momentum, defined by a Bessel-weighted moment with $b_T=0.3\,$fm, to equal
$\langle k_T^2\rangle = 0.19\,$GeV$^2$. The TMD evolution of our result is
studied using both the $b^*$ and $\zeta$ prescriptions which allows a
qualitative comparison with existing Drell-Yan data.
