Waveguides and Maxwell's Equations Background
In many textbooks on electromagnetism, the sequence in which topics are introduced are generally as follows -- electrostatics, magnetostatics, electrodynamics, wave propagation in free space, and wave propagation in confined space, in that order. As an example, consider Griffiths' 3rd edition of Introduction to Electrodynamics. Similarly, consider Orfanidis' Electromagnetic Waves and Antennas.
By using the discussion on wave propagation in free space as a launchpad, the discussion on wave propagation in a confined space usually introduces some additional boundary conditions that the solutions to Maxwell's equations for wave propagation in free space did not have.
If a proposed solution can be plugged into Maxwell's equations and not produce any contradictions, then they are valid. The goal is then to maintain the monochromatic wave solutions that worked for free space wave propagation with some "patches" to account for the fact that the propagating EM energy will be confined by the hollow waveguide. The reason monochromatic (one-wavelength) wave solutions are desirable is that, by superposition, any signal can be expressed as a sum of one or more monochromatic waves, and so any signal's behavior in a hollow waveguide can be understood.
One of these "patches" is the boundary condition that $E_{\parallel}=0$. This is because the net electric field in an ideal conductor, such as the material that makes up the confining hollow waveguide, must be $0$ because any externally applied electric field will influence electrons to rush to counter the externally applied electric field. If there is no net electric field (meaning $\vec{E}=\vec{0}$), then there is no net electric field in the tangential direction ($E_{\parallel}=0$) along and inside the metal.
If there is no net electric field along the metal, then, by Faraday's law, there must be no net changing magnetic field $\frac{\partial \vec{H}}{\partial t}$. This means that each of the components $\frac{\partial H_x}{\partial t}$, $\frac{\partial H_y}{\partial t}$, and $\frac{\partial H_z}{\partial t}$ must all be $0$ too. If there is no initial magnetic field (meaning $\vec{H}=\vec{0}$, then there must be no magnetic field at any later point. Therefore $\vec{H}=\vec{0}$ and $H_{\perp}=0$.
So these are the two additional boundary conditions that must be considered. An implicit assumption is that the material on the interior of the waveguide is isotropic, homogeneous, and linear. This means that $\vec{B}=\mu\vec{H}$ and therefore $\vec{B}=\vec{0}$ and $B_{\perp}=0$ at the boundary of the metal + hollow interior interface as well.
Problem
Frequently, the monochromatic wave solution to Maxwell's equations given the boundary conditions of a hollow waveguide are presented as
$$
\vec{F}\left(x,y,z,t\right)=\left[\text{spatial component}\right]\left[\text{temporal component}\right]
$$
For example from Orfanidis,
$$
\vec{E}\left(x,y,z,t\right)=\vec{E}\left(x,y\right)e^{j\omega t - jkz}
$$
and from Griffiths,
$$
\vec{E}\left(x,y,z,t\right)=\widetilde{\vec{E}}_0\left(x,y\right)e^{i\left(kz - \omega t\right)}
$$
Further, $\vec{E}$ is a 3-dimensional vector with components $E_x\left(x,y\right)$, $E_y\left(x,y\right)$, and $E_z\left(x,y\right)$.
Bottom line: Why is the spatial component only a function of $x$ and $y$? Further, why does spatial component have $x$, $y$, $z$ components [e.g. $E_x$, $E_y$, and $E_z$] (where $E_x$, $E_y$, and $E_z$ are each functions of only $x$ and $y$)?
 A: To begin with, Maxwell's equations are vector partial differential equations, which makes their solution intrinsically more difficult than that of scalar partial differential equations.
To help obviate this difficulty, an ansatz was long ago made (by Hertz himself, perhaps?) who proposed two useful categories of wave solution: transverse electric (TE), and transverse magnetic (TM).  The term 'Transverse' here means 'transverse to the waveguide axis', which is conventionally the 'z' direction, and the presumed direction of propagation. 
These categories are somewhat self-explanatory.  TM waves have no z directed magnetic field component, and TE waves have no z directed electric field component.  Each of these is derivable as the curl of an appropriate vector potential, magnetic or electric.  If the idea of an electric vector potential is troublesome to you, try the book by Harrington 'Time Harmonic Electromagnetic Fields.' (With no denigration of Griffiths, I advise that you widen the scope of your reading.) In either case the vector potential is assumed to have a single, z-directed, component, comprising a scalar function of the form f(x,y)exp(ikz). This entire scalar function is presumed to satisfy the scalar Helmholtz equation (with propagation constant k), which insures that taking its curl (and the curl of the resultant) will lead to Maxwell's equations.
In essence, TE and TM waves are assumed; these ideas came out of someone's brain.  Their justification is that they lead us (as noted) to Maxwell's equations, and from there to innumerable useful solutions and applications. Inasmuch as a waveguide is designed to propagate a monochromatic wave of some chosen frequency, the simplest form of traveling wave solution is exp(ikz), as noted.  As far as why one wants a monochromatic wave, one answer would be that (at a guess) >90 % of communication electronics involves the modulation of a carrier wave, with the frequency bandwidth of modulation typically small in comparison to the carrier frequency.
Sorry to go on at such length. I hope this helps.
A: About your "bottom line": there is translational invariance along $z$ in the problem. so the electric field can be expanded into a Fourier integral with respect to $z$. So there is dependence on $z$, but it can be expressed as a superposition of exponents.
Some waveguide modes do have nonzero $z$-component of electric field (although this component is zero on the boundaries of the ideally conducting waveguide).
A: Akhmeteli is absolutely spot on. Since he/she has thrust to the problem's very heart with his / her keen dagger, just a few more words to clean up the lacerations on the wound!
It IS all about z-axis shift invariance: so this symmetry requires the existence of solutions which are eigenfunctions of $\partial_z$, namely functions of the form $f(x,y) \exp(i\,k_z\,z)$. Note that this form gainsays your assertion that the spatial component is a function of $x$ and $y$ only: there is the $\exp(i\,k_z\,z)$ as well. As you see, the $z$ variation is so simple that you can often overlook it (I have made this very mistake many times myself).
To understand the solution to this problem at a practical level, imagine we are solving this problem numerically and represent the transverse field as some numerical grid: we can rewrite a 2D grid into a 1D column vector $\psi$ with $N = n\times m$ complex components where the transverse 2D grid resolution is $n$ points by $m$ points. Maxwell’s equations are linear so that the action of a length $z$ of waveguide can be represented by a big square $N\times N$ matrix $U(z)$, i.e. $\psi \mapsto U\,\psi$. Now we think about concatenating lengths of waveguide. Owing to our translational invariance, we have:
$$U(z_1)\,U(z_2) = U(z_2)\,U(z_1) = U(z_1 + z_2)$$
for any lengths $z_1,\,z_2\in\mathbb{R}$, and the only continuous solution to this functional equation is $U(z) = \exp(H\,z)$, for some constant square $N\times N$ matrix $H$. Now if we set up our column vector right, its $\ell^2$-length (i.e. its norm) $\sqrt{\psi^\dagger\psi}$ is proportional to the total propagating power through any transverse cross section of the waveguide. If the waveguide is lossless, this norm cannot change, so that $U(z)$ is also unitary (preserves lengths of and angles between vectors) and this can be shown to be equivalent to the condition $H^\dagger = -H$, i.e. $H$ is skew Hermitian. The spectral theory for finite normal operators (i.e. ones that commute with their adjoint, so square skew Hermitian matrices are such operators) is very well studied and easy: such a matrix can always be diagonalised by a unitary matrix of $N$ orthonormal column vectors, and so the matrix exponential  $U(z) = \exp(H\,z)$ can be expressed in the form:
$$U(z) = \Omega \, \operatorname{diag}[e^{i\,k_1\,z}, \, e^{i\,k_2\,z}, \, e^{i\,k_3\,z}, \,\cdots] \Omega^\dagger$$
where now the $ e^{i\,k_1\,z}$ are all scalar exponentials and the $k_j$ are all real (this together with the unitary decomposition is equivalent to skew-Hermitianhood of $H$). So now think what the eigenvectors – the columns of $\Omega$ mean, given that this is a numerical problem. They stand for certain special field configurations in the transverse plane, i.e. they are discretised versions of the eigenfunction variation $f_j(x, y)$. So, through this spectral factorization, you can see all numerical solutions to the problem are linear superpositions of the kind:
$$\sum_j\alpha_j\,f_j(x,\,y)\, \exp(i\,k_j\,z)$$
(for the numerical problem, you can think of $x$ and $y$ as discrete grid indices passing to continuous transverse co-ordinates in the limit and the $f_j$ are our "grid" eigenfields) and hopefully you can accept that this numerical representation will approach the real one as the grid spacing gets smaller and smaller, so the above is an explanation for the kinds of solutions you are asking about. The numerical solution of this problem seldom works like this - $N = n\times m$ is huge (think of a problem with a $1024\times1024$ transverse discretization) and there are much more efficient ways of representing the action of propagation through an axial step $\Delta z$ than to work out a matrix $U(\Delta z)$ or $H$ (which would hold $10^{12}$ complex numbers in the $1024\times1024$ grid example!), but you can see the above is an in-principle description: there is a square unitary grid operator and it can always be factorized in the way I describe. So the numerical procedure is always equivalent to the above: one just finds clever algorithms that are equivalent in their action but do not need to build all of the huge $U(\Delta z)$ in memory at once. Finite difference methods, for example, assume, accurately, that over small steps, the field at a point in the updated grid is only influenced by the fields at neighbouring points (for first order PDE problems) or at neighbouring and next-to-neighbouring points (for second order PDE problems, as Maxwell's equations are). So we only need to represent in memory the leading diagonal together with two neighbouring stripes of the matrices $H$ and $U(z)$: this is roughly $5\times 1024\times 1024$ points in my example, or about $16\times 5 = 80$ megabytes on a 64-bit machine representing complex numbers as two one-word (8 bytes) reals, i.e. readily workable in today's technology. The co-efficients in the $H$ matrix come straight from discretized versions of Maxwell's equations and tridiagonal (three stripe) pentadiagonal (five stripe) matrix factorization numerical procedures are very well worked out, tested and numerically sound.
Incidentally, the above discussion is wholly analogous to the theory behind the general Schrödinger picture in quantum mechanics: see my answers to the Physics SE questions A Simple Explanation for the Schrödinger Equation and Model of Atom?, Reference frame involved in the Schrödinger's equation and Why can the Schroedinger equation be used with a time-dependent Hamiltonian?.
In non-numerical simulation words: in a lossless waveguide, the $z$-shift operation is indeed unitary (no energy loss) so the "infinitesimal" shift $\partial_z$ is skew-Hermitian i.e. the eigenvalues $i\,k_z$ are purely imaginary. You need a bit of theory - a vector analogue of Sturm-Liouville theory - to see why such eigenfunctions can be added together to match arbitrary boundary conditions on a transverse cross-section. Actually, given the vector nature and further conditions of solenoidal fields ($\nabla\cdot\vec{E}=\nabla\cdot\vec{H}=0$), there's a bit more to a proof than standard Sturm-Liouville theory. Standard SL theory will work of course for two out of three Cartesian components of the $\vec{E}$ or $\vec{H}$ field (or Cartesian component of the Lorenz-gauged potential four-vector). For a metal waveguide with a compact cross-section, the spectrum of eigenfunctions is wholly discrete: for an optical fibre, there is a continuous component to the spectrum - the radiation field - as well as the discrete, losslessly guided bound eigenfields.
The full story is to be found in Snyder and Love, "Optical Waveguide Theory", particularly chapters 30 and 31 for groundings and earlier chapters for the radiation field of optical waveguides (I understand from John Love that the ANU will make this book freely downloadable sometime) and R. E. Collins, "Field Theory of Guided Waves". However, neither of these books is particularly up front about the mathematical foundations of the "completeness" of the vector eigenfunctions and I don't know of any reference that discusses this idea in simple, modern terms. 
