Square root of number operator for quantum harmonic oscillator Let $a$, $a^{\dagger}$ denote the standard annihilation and creation operators for the quantum harmonic oscillator, with $[a, a^{\dagger}] = \mathbb{I}$. The number operator is then defined as $a^{\dagger}a$. As far as I am aware, $a^{\dagger}a$ is an unbounded self-adjoint operator in the complex Hilbert space for the oscillator.
My first question: does the square root operator $\sqrt{a^{\dagger} a}$ exist in a rigorous mathematical sense? By the square root $A = \sqrt{B}$ of some operator $B$, I mean an operator that satisfies $A^2 = B$. I guess this boils down to the existence of square roots for general unbounded self-adjoint operators in complex Hilbert spaces? I would be thankful if someone could provide a reference containing a proof.
My second question: what are the commutation relations for $\sqrt{a^{\dagger} a}$? In particular, I am interested in the commutators $[\sqrt{a^{\dagger} a}, a^{\dagger}a]$, $[\sqrt{a^{\dagger} a},  a]$ and $[\sqrt{a^{\dagger} a},  a^{\dagger}]$. I think the first one should be zero, but I am not quite sure what a rigorous proof would look like. For functions $f(A)$ of some operator such that $f(A)$ contains only integer non-negative powers of $A$, I believe that $[A, f(A)] = 0$ since any operator commutes with itself and thus also with products of itself, but I am not immediately convinced that this translates to roots $A^{1/k}$ (if they exist) since those cannot just be written as products of operators.
Since $a a^{\dagger} = a^{\dagger} a + \mathbb{I}$, I assume that the answers to the above questions could also be applied to $\sqrt{a a^{\dagger}}$ with only slight modifications.
 A: The square root exists and it is  defined by standard functional calculus for every selfadjoint operator $A: D(A) \to H$ with non negative spectrum (which is equivalent to $\langle x, Ax\rangle \geq 0$ for all $x\in D(A)$):
$$\sqrt{A}:= \int_{\sigma(A)}\sqrt{\lambda} dP^{(A)}(\lambda)\:.$$
Above $P^{(A)}$ is the spectral measure of $A$.
It is possible to prove that $\sqrt{A}$ above is the unique positive selfadjoint operator $B$ such that $BBx= Ax$ for all $x\in D(A)$.
Concerning commutativity of $A$ with measurable functions of $A$, what matters, and give rise to subtleties,  are the domains of the compositions $Af(A)$ and $f(A)A$ when $A$ is not bounded and/or $f$ is not bounded as a function.
However the following fact holds
$$Af(A)x = f(A)Ax$$
if $x\in H$ is such that both sides are defined. There are a number of subcases where well definedness of one side implies the same for the other side. All that is standard spectral theory. As a reference I can quote my couple of books 1, 2, but there are very well written modern textbooks on the subject. For instance the recent book by Schmudgen on spectral theory and the Reed-Simon series (however they cover a so wide range of subjects that it is sometime difficult to point out a spot where to find the precise result one is interested in).
So $$\sqrt{a^\dagger a} a^\dagger a x = a^\dagger a \sqrt{a^\dagger a} x$$ for the vectors $x\in H$ such that both sides are defined.
(Notice that the standard creation operator is not the adjoint of $a$ but just a restriction of it and thus $a^\dagger a$ has to be interpreted as a positive selfadjont extension of that product. A theorem by von Neumann proves that, using the closures of the indicated operators,  it exists and it is unique so that we can safely apply the procedures of spectral calculus.)
I stress that the use of Taylor series is just euristic and it does not guarantee the validity of any found identity. It has some technical utility in finite dimensional spaces or bounded operators, but it is not the way to rigorously tackle the operator theory in a general setting, except for very special cases (domains of analytic vectors).
Regarding a close form for $[\sqrt{a^\dagger a}, a]$ and $[\sqrt{a^\dagger a}, a^\dagger]$, I do not know, because $a$ and $a^\dagger$ are not selfadjoint, nor normal and they are not a function of
$a^\dagger a$ in spectral sense. So the formula you are looking for, if any,  has not a spectral origin, but it depends on other properties of the written operators.
A: Asking for such references may belong to the MSE, instead of this site. Most sensible, well-meaning, physicists work out the requisite relations in Fock space themselves, extending the humdrum Sylvester formula for finite matrices, and obtain relations which mostly work.
So it is generally first assumed that functions of $a^\dagger a$ are expandable in powers of $a^\dagger a$, for which commutation rules are straightforward, so, obviously,
$$
[f(a^\dagger a), a^\dagger a ]=0, \\
[f(a^\dagger a), a]= a(f(a^\dagger a+1)- f(a^\dagger a)),\\
 [f(a^\dagger a), a^\dagger ]= (f(a^\dagger a)-f(a^\dagger a+1)   ) a^\dagger,
$$
etc, and then one looks of singularities or breakdowns in more recondite cases.
Quantum optics delight in investigating marginal twisted cases, represented in several questions or else on this site, but the rule is you first work out your relations yourself before hitting sourcebooks and compilations that are likely to confuse you.
