Expanding two-variable function $f(x,y)$ over the complete sets $\{ g_{i}(x) \}$ and $\{ h_{j}(y) \}$ Quite often (see, for example, this PDF, 50 KB) when discussing the Born-Oppenheimer approximation the following assertion is made: any well-behaved function of two independent variables $f(x,y)$ can always be expanded over the complete set of functions $\{ g_{i}(x) \}$ and $\{ h_{j}(y) \}$ in the following way
$$
f(x,y) = \sum\limits_{i} \sum\limits_{j} c_{ij} g_{i}(x) h_{j}(y) \, , \quad (1)
$$
or by defining
$$
c_{j}(x) = \sum\limits_{i} c_{ij} g_{i}(x) \, ,
$$
as
$$
f(x,y) = \sum\limits_{j} c_{j}(x) h_{j}(y) \, .
$$
Sometimes an argument in favor of the statement above is that it is equivalent to expanding a function of one variable $f(y)$ over the complete set of functions $\{ h_{j}(y) \}$
$$
f(y) = \sum\limits_{j} c_{j} h_{j}(y) \, ,
$$
with the difference that coefficients $c_{j}$ in the former case carry the $x$-dependence.
For me it seems like this argument is a bit far-fetched. So my question is: how do we know that (1) is true? Is it a theorem, or an axiom, or something?
 A: Start with $f(x,y)$. Fix a certain $x$, let's call it $x_0$. Obviously then $f(x_0, y)$ is a function of just one variable, so it has an expansion in the complete set,
$$f(x_0, y) = \sum_j c_j(x_0) h_j(y)$$
While the expansion coefficients now depend on $x_0$, since we get a different $y$-dependent function for each value of $x_0$ that we fix, the ability to expand doesn't depend on that. So we might as well call it just $x$ again,
$$f(x,y) = \sum_j c_j(x) h_j(y)$$
Next step: Since $c_j(x)$ is a single-variable function of $x$, we can expand it:
$$c_j(x) = \sum_i c_{ij} g_i(x)$$
Putting it all together,
$$f(x,y) = \sum_{i}\sum_j c_{ij} g_i(x) h_j(y)$$
EDIT: In the true fashion of a physicist, I did assume that the functions are sufficiently well behaved to do all these steps. Technically we have to worry about various notions of convergence of functions, but especially for square-integrable Hilbert spaces we should be fine :)
A: OP is asking a very important mathematical question (v1), which is used in many parts of mathematical physics, e.g. in solving PDEs by separation of variables.
Example: Consider functions on the positive real half line $\mathbb{R}_{>0}$ with a basis $(g_n)_{n\in\mathbb{Z}}$ of integer power functions $g_n(x):=x^n$.
Similarly $h_n(y):=y^n$, $n\in\mathbb{Z}$.
Let there be given a function $f: \mathbb{R}_{>0}^2\to \mathbb{C}$ on the product space as
$$ f(x,y)~:=~\frac{1}{x+iy}~=~\left\{\begin{array}{rcl}\frac{1}{x}\sum_{n=0}^{\infty}\left(\frac{y}{ix}\right)^n &\text{for} & 0<y<x, \\
\frac{1}{iy}\sum_{n=0}^{\infty}\left(\frac{ix}{y}\right)^n &\text{for} & 0<x<y.
\end{array}\right.
$$
Note that one cannot use the same series expansion in the whole domain  $\mathbb{R}_{>0}^2$. End of example.
To answered OP's question properly, one should know the precise definition of the pertinent function spaces and their topology. E.g. are the infinite sums supposed to be pointwise convergent, uniform convergent, or  convergent in ${\cal L}^{p}$-norm, or something different?
A surprisingly efficient and fairly broad setting is if the pertinent function space can be viewed as a Hilbert space with appropriate inner product. Often the basis of functions are orthonormal eigenfunctions for some Hermitian operator.
Say e.g., that we have two Hilbert spaces $H$ and $K$ with orthonormal base $(e_i)_{i\in I}$ and $(f_j)_{j\in J}$. Then it's a theorem that $(e_i\otimes f_j)_{(i,j)\in I\times J}$ is an orthonormal basis for the product Hilbert spaces $H\otimes K$.
Finally, if OP is more interested in continuous functions and uniform convergence, he might be interested in product space applications of Stone–Weierstrass theorem.
