The following question is not as well formulated as I would like, but here goes.
How should one think about unbounded linear operators on a Hilbert space $H$? Even though I have read a little on the subject (Conway, Reed and Simon, etc.), I feel like I'm still not able to manipulate them as I would with bounded linear operators. Multiplication $AB$ of unbounded operators is weird and the topology of on the space of unbounded operators is tricky. Maybe I'm not getting the big picture, so some general tips or maybe a nice reference to review the subject would be appreciated.
To put my general question into context, I am mostly interested in dealing with creation/annihilation operators $b^*, b$ on the bosonic Fock space $\mathscr{F}$ and quadratic Hamiltonians of the form $H = \sum_{k=1}^\infty \epsilon_k b_k^* b_k$ where $\epsilon_k \in \mathbb{R}$ are "well-behaved". A particular example would be that I would like to calculate $$ \langle b_k^* b_k \rangle = \operatorname{tr} (b^*_k b_k \rho) $$ where $\rho = e^{-H}/Z$ is a nonnegative linear operator on $\mathscr{F}$ with $\operatorname{tr} (\rho)=1$ ($Z$ is a normalization factor so that the trace is 1). In this particular example, even though $\rho$ is of trace-class, I am not sure why $\operatorname{tr} (b^*_k b_k \rho)$ is well-defined.
TL;DR: I would like some references that deal with self-adjoint unbounded operators and how they "interact" with the trace function, bounded linear operators, etc. Examples of application in physics (especially second quantization) would be appreciated.