Model of a crystal and deriving the expected diffraction image
In single-crystal X-ray crystallography, the diffraction image is due to scattering of electrons in a crystalline sample. A convenient way of describing the electron density is to first specify crystal symmetry and mean atomic positions (coordinates), and then to describe in more or less detail how atoms deviate from these mean positions (modeled by atomic displacement factors). Combining this with a spherical description of the electrons "belonging" to each atom (modeled by atomic form factors), you arrive at an electron density.
This description is useful because in the model, the electron density is the convolution of the crystal lattice points with the coordinates in the unit cell with the distribution of atoms around mean positions with the distribution of electrons around the atoms. The convolution theory states that the Fourier transform of a convolution is the product of the Fourier transforms. So in this model, you can separately Fourier-transform the different levels of the model and then multiply them to arrive at the expected diffraction image (via the structure factors).
Atomic displacement and Debye-Waller factor
The atomic displacement around a mean position may be described through a probability density function $p_k(\mathbf{u})$, where $\mathbf{u}$ is the position (3D-vector) with the mean position as origin. The Fourier transform of this function defines the Debye-Waller factor (source: equation 1.4.8):
$$T_k(\mathbf{h}) = \int p_k(\mathbf{u}) e^{2 \pi i \mathbf{h} \cdot \mathbf{u}}d^3 \mathbf{u}\tag{3}$$
Here, $\mathbf{h}$ represents coordinates in reciprocal space (just as Fourier transform of a time series gives you information in frequency space, the Fourier transform connects real space - coordinate $\mathbf{u}$ - with reciprocal space - coordinate $\mathbf{h}$).
If the displacement is modeled by a Gaussian, the integral can be determined as (source: equation 1.4.10)
$$T_k(\mathbf{h}) = e^{-2 \pi^2 \langle(\mathbf{h} \cdot \mathbf{u})^2\rangle}\tag{4}$$
This reflects that the Fourier transform of a Gaussian is a Gaussian again, and the width of one Gaussian is the inverse of the other (i.e. if the atoms are widely distributed, the signal will diminish a lot with increasing resolution). Again using the time series as an example, if a sound is played for a very short time, the pitch will be ill-defined (or if an NMR FID signal decays rapidly, the peaks will be broad).
Displacement parameters $\mathbf{U}$
The displacement parameters $\mathbf{U}$ or $\mathbf{U^{jl}}$ in the OP's equation (1) and (2) are derived from the $\langle(\mathbf{h} \cdot \mathbf{u})^2\rangle$ expression in equation (3), either direction-independent in the case of an isotropic approximation (spherical distribution around the mean) or separately along the three lattice axes (leading to six anisotropic displacement factors). I am not explicitly showing this - it gets a bit complicated because the derivation requires transformation into the sometimes non-orthogonal coordinate system based on the crystal lattice (see here). In any case, the expression is a mean square displacement (analogous to a variance of a one-dimensional function). For that reason, the dimensions are length squared and the typical units are $Å^2$. As stated in the IUCR nomenclature document:
the elements of the tensor U have dimension (length)2 and can be directly associated with the mean-square displacements of the atom considered in the corresponding directions
and
If the atomic pdf is assumed to be a trivariate Gaussian, the characteristic function corresponding to this pdf - by definition, its Fourier transform - can be described by the second moments of the pdf, which in the present context are called anisotropic mean-square displacements.
Second moments (variance) have dimensions that are the square of the dimension of the quantity described by the probability distribution function (pdf). In this case, we are looking at the distribution of the position (dimensions length), and the second moment will have dimensions of length squared.