In a paper I have been reading they are talking about a function of two variables $\xi(r,t)$ ($r$ is the radial coordinate in a cylindrical coordinate system and $t$ is actually time, with $0 < r_1 < r < r_2$ and $0 \leq t$).
The authors state that as we are working in cylindrical coordinates, the eigenfunctions of the regular Sturm-Liouville operator defined by the Bessel equation and some appropriate boundary conditions (which they provide) in $r_1$ and $r_2$ can be used to write $\xi$ as \begin{equation} \xi(r,t) = \displaystyle\sum_{k=1}^{\infty} a_n(t) \psi_n(r) \text{,} \end{equation} with $\psi_n$ the eigenfunctions.
I am not an expert in pure mathematics and I have never seen a linear combination like this where the coefficients themselves are variables. I understand that if $t$ isn't there the $\psi_n$ form a complete set and a function can be written as a linear combination of the $\psi_n$, but when there is a second coordinate ($t$ in this case) can you just make the coefficients depend on that coordinate and that's it? Or is there some theory that states this?
It's easier than you think. For fixed $t$, $\xi(t,r)$ is a function of one variable (namely, $r$). Hence for each such $t$, it can be expanded in terms of eigenfunctions:
$\xi(t,r)=\sum_1^\infty a_n(t)\psi_n(r)$.
So yes your formula follows in a straightforward manner from the completeness of the eigenfunctions (for the purposes of this thread -- and since it seems related more to physical problems -- im assuming pointwise convergence of the series)
At this point the coefficients $a_n(t)$ seem entirely arbitrary however we can improve this situation: We know exactly how to get those coefficients $a_n(t)$. If you multiply both sides of the above equation by $r\psi_m(r)$ (as I recall, the bessel eigenfunctions are orthogonal w.r.t to this weighted inner product) and integrate w.r.t r you get:
$\int r\xi(t,r) \psi_m(r)=\int\sum_1^\infty a_n(t)\psi_n(r)\psi_m(r)$
We can interchange the sum and the integral under fairly general conditions (for instance, if the integrand on the right hand side is absolutely integrable), and thus we get:
$\int r\xi(t,r) \psi_m(r)=\sum_1^\infty \int a_n(t)\psi_n(r)\psi_m(r) = \sum_1^\infty a_n(t)\delta_{nm} = a_m(t)$
Hence $a_m(t)=\int r\xi(t,r) \psi_m(r)$.
So not only can you do such an expansion but the coefficients can be explicitly calculated and the coefficients will be well behaved if $\xi(t,r)$ is. For instance the coefficients will be continuous or differentiable if $\xi$ is $C^0$ or $C^1$ respectively.
Edit: I can see where your confusion lies. So let me perhaps go into a bit more detail why eigenfunctions are independent of time:
Let us forget about the time variable for now and suppose you have two single-variable functions $g(r)$ and $h(r)$ defined on $0<r_1<r<r_2$ (the domain on which any function can be written in terms of your eigenfunctions). Then as we know we can expand each of them in terms of eigenfunctions:
$g(r)= \sum_1^\infty a_n(g)\psi_n(r)$
$h(r)= \sum_1^\infty a_n(h)\psi_n(r)$
Where the coefficients $a_n(g)$ and $a_n(h)$ are meant to indicate that they are the Bessel series coefficients for $g$ and $h$ respectively. Are you ok with the above steps? If so then observe that the eigenfunctions $\psi_n$'s are the same for both $g$ and $h$, what differs are the coefficients (this should be familiar to you, it's the same thing you do in Fourier series, you always expand in terms of sines and cosines, what you must find are the coefficients). The $\psi's$ are the same for both functions. This follows directly from the fact that the eigenfunctions "span" the space.
Now introduce the time variable. Let $t_1$ and $t_2$ be two different points in time, then similar to the above $\xi(t_1,r)$ ($t_1$ fixed) and $\xi(t_2,r)$ ($t_2$ fixed) are two different (in general) functions on the domain $0<r_1<r<r_2$. So similar to the above we expand them:
$\xi(t_1,r)=\sum_1^\infty a_n(t_1)\psi_n(r)$.
$\xi(t_2,r)=\sum_1^\infty a_n(t_2)\psi_n(r)$.
But notice that the above can be written for an arbitrary fixed $t$:
$\xi(t,r)=\sum_1^\infty a_n(t)\psi_n(r)$.