Notations mostly follow https://en.wikipedia.org/wiki/Mercer%27s_theorem.
Mercer's theorem uses the eigenfunctions $\{e_j\}$ of the integral operator as the expansion function. I wonder if we could do the same with an orthonormal system $\{g_j\}$ at hand, that we do not necessarily know if they are eigenfunctions of an integral operator $K$ associated to a given $\kappa$.
We have $\kappa: [a,b] \times [a,b] \rightarrow \mathbb{R}$, a symmetric and continuous kernel. Further, we assume that $\kappa$ is a positive definite kernel on $[a,b]$, that is, $ \sum_{i=1}^n\sum_{j=1}^n \kappa(x_i, x_j) c_i c_j \geq 0$ holds for for all finite point set $\{x_1, ..., x_n\}\subset[a, b]$ and all choices of $\{c_1,\dotsc,c_n\}\subset\mathbb{R}$.
Consider the integral operator on $(L^2(a,b),\mu)$ associated with $\kappa$, that is, $$ [K \varphi](x) =\int_a^b \kappa(x,s) f(s)\, d\mu(s),\qquad f\in L^2_\mu(a,b),$$ where $L^2_\mu(a,b)$ is the $L^2$ space with respect to a finite measure $\mu$.
Now, we have $L^2_\mu(a,b)$-orthonormal system $\{g_j\}$ at hand.
Question: Does the following hold? $$ \kappa(x,y)=\sum_{j=0}^\infty \gamma_j g_j(x)g_j(y), $$ where $\gamma_j$'s are some constant.
My guess and hope is that it holds with $$ \gamma_j=\int_{[a,b]\times[a,b]} \kappa(y,x) g_j(y) g_j(x) \,d[\mu \otimes \mu](y,x), $$ as it is the case when $\{g_j\}$'s are the eigenfunctions, i.e., $g_j=e_j$.
The thing I cannot be sure is that the expansion is the product form, by which I mean the sum is of the form $\sum_{j=0}^\infty \gamma_j g_j(x)g_j(y)$, as opposed to $\sum_{j=0}^\infty\sum_{k=0}^\infty \gamma_{jk} g_j(x)g_k(y)$, which is what we get if we expand $k(x,y)$ in each term naively. Surely the former expansion gives an positive definite ($K$ is a positive operator and thus $\gamma_i=\langle Kg_j,g_j\rangle\ge0$) kernel $\kappa^g$ but not sure if $\kappa=\kappa^g$.
I tried to use the positive definiteness, or the uniqueness of the reproducing kernel but could not argue properly. If this is a well known result I would like a reference.