Consider a separable, complex Hilbert space $H$ and a self-adjoint operator $A: D(A) \longrightarrow H$, where $D(A)$ is a dense subspace. Assume that $A$ admits an orthonormal basis $\{\varphi_n\}_{n\in \mathbb N}$ as eigenvectors: $A \varphi_n =a_n\varphi_n$. I want to prove the following:
$$x\in \mathcal D(A) \iff (a_n \langle \varphi_n,x\rangle)_{n\in \mathbb N} \in \ell^2\quad . $$
Proof attempt. If $x\in D(A)$, then we can expand $Ax \in H$ in the eigenbasis of $A$, which yields that $$Ax=\sum\limits_{n\in \mathbb N} \langle \varphi_n,Ax\rangle \varphi_n = \sum\limits_{n\in \mathbb N} a_n \langle \varphi_n,x\rangle \varphi_n \quad .$$ But as $Ax\in H$, we have that $(a_n \langle \varphi_n,x\rangle)_{n\in \mathbb N} \in \ell^2$.
Conversely, if $(a_n \langle \varphi_n,x\rangle)_{n\in \mathbb N} \in \ell^2$ holds for some $x\in H$, then there exists a unique $y\in H$ such that
$$ y= \sum\limits_{n\in \mathbb N} a_n \langle \varphi_n,x\rangle \varphi_n \quad .$$
Expanding $y$ in the eigenbasis of $A$, on the other hand, gives $$y= \sum\limits_{n\in \mathbb N} \langle \varphi_n,y\rangle \varphi_n \quad , $$
which due to the uniqueness of the expansion coefficients yields
$$\forall n\in \mathbb N: a_n \langle x,\varphi_n\rangle = \langle x,A\varphi_n\rangle= \langle y,\varphi_n\rangle \quad . $$
Recalling that $$D(A^*):=\left\{ x\in H | \exists y\in H: \forall \varphi \in D(A): \langle x,A\varphi\rangle = \langle y,\varphi\rangle\right\} $$ and, as $A$ is self-adjoint, $A=A^*$ and $D(A)=D(A^*)$, it thus suffices to show that $x \in D(A^*)$ in the sense above. Then also $A^*x =Ax=y$.
To do so, expand $\varphi\in D(A)$ in the eigenbasis of $A$, which shows that $$\langle y,\varphi\rangle = \sum\limits_{n\in \mathbb N} \langle y,\varphi_n\rangle \langle \varphi_n,\varphi \rangle =\sum\limits_{n\in \mathbb N} a_n \langle x,\varphi_n\rangle\langle \varphi_n,\varphi\rangle =\sum\limits_{n\in \mathbb N} \langle x,A\varphi_n\rangle\langle \varphi_n,\varphi\rangle = \langle x,A\varphi\rangle \quad , $$
for all $\varphi \in D(A)$. Thus $x\in D(A^*)=D(A)$, concluding the proof.
We've used the continuity of the inner product, the self-adjointness (and thus symmetry) of $A$ as well as (variants of) the Riesz-Fischer theorem multiple times.
Is this proof correct? If so, can it be improved?
I think the last part should be
$$\sum\limits_{n\in \mathbb N} a_n \langle x,\varphi_n\rangle\langle \varphi_n,\varphi\rangle =\sum\limits_{n\in \mathbb N} \langle x,\varphi_n\rangle\langle \bar{a_n}\varphi_n,\varphi\rangle = \sum\limits_{n\in \mathbb N} \langle x,\varphi_n\rangle\langle a_n\varphi_n,\varphi\rangle = \sum\limits_{n\in \mathbb N} \langle x,\varphi_n\rangle\langle A\varphi_n,\varphi\rangle = \sum\limits_{n\in \mathbb N} \langle x,\varphi_n\rangle\langle \varphi_n,A\varphi\rangle = \langle x,A\varphi\rangle$$
Using that the $a_n$ must be real and $A$ is self-adjoint. Other than that looks ok to me.
[I'm assuming you're using the Physics convention of conjugate-linearity on the left which is consistent with your formulas.]