Now I have a function $f=f(x,y)$, smooth and symmetric(i.e. $f(x,y)=f(y,x)$ everywhere), with arguments defined on a compact set: $(x,y)\in[0,1]\times[0,1]$.
I'd wish to know if $f$ can be expanded as $f(x, y) = \sum_{l=1}^\infty \lambda_l \xi_l(x)\xi_l(y)$, where $\xi_l(\cdot)$ are functions that are norm 1 and orthogonal to each other -- here I define the needed inner product between two functions $u, v$ on $[0,1]$ to be $<u, v> := \int_0^1 u(t)v(t)dt$.
A reasonably easy to follow reference will be good enough for my question. Just I didn't know what the correct key words to use -- I tried "functional spectral theory" etc and got something too abstract for me to understand(and they were unnecessarily general for my purpose of application). Thank you so much in advance!
For simplicity, assume $f(x,y)=f(y,x)$ is a real function. Because $f$ is in $L^{2}([0,1]\times[0,1])$, then the integral operator $K$ given by $$ Kg = \int_{0}^{1}f(x,y)g(y)\,dy $$ is a selfadjoint Hilbert-Schmidt integral operator on $L^{2}[0,1]$. So there is an orthonormal basis $\{\varphi_{n}\}_{n=1}^{\infty}$ consisting of real eigenfunctions of $K$ with corresponding real eigenvalues of $K$. The eigenfunctions with non-zero eigenvalues are smooth because they inherit smoothness from $k$: $$ \varphi_{n}(x)=\frac{1}{\lambda_{n}}\int_{0}^{1}f(x,y)\varphi_{n}(y)\,dy. $$ To show that $\{ \varphi_{n}(x)\varphi_{m}(y)\}_{n,m=1}^{\infty,\infty}$ is an orthonormal basis of $L^{2}([0,1]\times[0,1])$, notice that this set is an orthonormal subset of $L^{2}([0,1]\times[0,1])$, and Parseval's equality for $L^{2}[0,1]$ gives Parseval's equality for every $g \in L^{2}([0,1]\times[0,1])$: $$ \begin{align} &\sum_{n,m}\left|\int_{0}^{1}\int_{0}^{1}g(x,y)\varphi_{n}(x)\varphi_{n}(y)\,dx\,dy\right|^{2} \\ & = \sum_{n}\sum_{m}\left|\int_{0}^{1}\left(\int_{0}^{1}g(x,y)\varphi_{n}(x)\,dx\right)\varphi_{m}(y)\,dy\right|^{2} \\ & = \int_{0}^{1}\sum_{n}\left|\int_{0}^{1}g(x,y)\varphi_{n}(x)\,dx\right|^{2}\,dy \\ & = \int_{0}^{1}\int_{0}^{1}|g(x,y)|^{2}\,dx\,dy. \end{align} $$ So, the following sum converges in $L^{2}([0,1]\times[0,1])$: $$ \begin{align} K(x,y) = \sum_{n,m}\left[\int_{0}^{1}\int_{0}^{1}f(x,y)\varphi_{n}(y)\varphi_{m}(x)\,dydx\right]\varphi_{n}(x)\varphi_{m}(y) \end{align} $$ Because $\varphi_{n}$ is an eigenfunction of $K$ with eigenvalue $\lambda_{n}$, then $$ \begin{align} K(x,y) & = \sum_{n,m}\left[\int_{0}^{1}\left[\int_{0}^{1}k(x,y)\varphi_{n}(y)\,dy\right]\varphi_{m}(x)\,dx\right]\varphi_{n}(x)\varphi_{m}(y) \\ & = \sum_{n,m}\left[\lambda_{n}\int_{0}^{1}\varphi_{n}(x)\varphi_{m}(x)\,dx\right]\varphi_{n}(x)\varphi_{m}(y) \\ & = \sum_{n}\lambda_{n}\varphi_{n}(x)\varphi_{n}(y). \end{align} $$ All functions $\varphi_{n}$ in the final sum are smooth because $0$ eigenvalue terms drop out. The convergence of the sum is in $L^{2}([0,1]\times[0,1])$.