Suppose I have a linear unbounded operator $A:\operatorname{dom}(A) \to \mathscr{H}$, with $\operatorname{dom}(A)$ carefully chosen to make $A$ self-adjoint ($A^\dagger=A$).
Could $A^n$ be automatically self-adjoint $\forall n$ on the same domain of $A$ (in short $\operatorname{dom}(A)=\operatorname{dom}(A^n)$)?
On
A practicing mathematician should answer this question (I'm only a former mathematical physicist), but since no one has answered in a few days, I'll upvote the question and say provisionally that I tend to doubt that ${\rm dom}(A)$ could equal ${\rm dom}(A^2)$ for any unbounded self-adjoint operator $A$ on an infinite-dimensional Hilbert space.
The reason for my doubt is that I'm thinking about the following example. Let ${\bf H}$ be the Hilbert space $L^2({\bf R})$, and let $(D_X, X)$ be the self-adjoint position operator $X$ with domain $D_X$ consisting of all vectors in ${\bf H}$ having, as a representative, a function $f$ with $xf$ square-integrable. (See Geroch, Mathematical Physics, University of Chicago Press, 1985, pp. 330-332.) I don't know for sure what the domain of $X^2$ is, but I think it is clear that in order for $X^2$ to be well defined at all (let alone self-adjoint), its domain will have to be smaller than $D_X$.
So the multiplicative operator $X:f\mapsto xf$ will not serve as an $A$ in your hypothesis. [Strike this: But by the spectral theorem, all self-adjoint operators are multiplicative in some basis, so ] my intuition is that no unbounded self-adjoint $A$ whatsoever can serve in your hypothesis.
(If this answer is nonsense based on faulty memory from 20 years ago, anyone please let me know and I'll quickly delete it.)
No, $\ \text{dom}(A^2)\ $ will always be a proper subset of $\ \text{dom}(A)\ $.
Because $\ A\ $ is unbounded, there exists $ y\not\in\text{dom}(A)\ $, and because $\ A\ $ is self-adjoint, its spectrum lies entirely inside the real line. Therefore $\ A+iI\ $ has a bounded inverse, $\ (A+iI)^{-1}:\mathscr{H}\rightarrow$$\,\text{dom}(A+iI)=$$\,\text{dom}(A)\ $, and $$ z=(A+iI)^{-1}y\in\text{dom}(A)\ . $$ Therefore, \begin{align} Az+iz&=(A+iI)(A+iI)^{-1}y\\ &=y\ . \end{align} Since $\ \text{dom}(A)\ $ is a subspace of $\ \mathscr{H}\ $, then $\ iz\in\text{dom}(A)\ $, and so if $\ Az\ $ were in $\ \text{dom}(A)\ $, then so would be $\ y=$$\,Az+iz\ $, which is not the case. Therefore $\ Az\not\in\text{dom}(A)\ $, so $ z\not\in\text{dom}(A^2)\ $. Thus $\ z\in$$\,\text{dom}(A)\setminus\text{dom}(A^2)\ $, and $\ \text{dom}(A^2)\ $ must be a proper subset of $\ \text{dom}(A)\ $.