Consider the position and momentum operator on $L^2(\mathbb{R})$ defined on a dense domain, say $D(X) = D(P)=\mathscr{S}(\mathbb{R})$ with the position operator being defined as $$ X\phi=x\phi(x) \qquad \forall \phi \in \mathscr{S}(\mathbb{R}) $$ and the momentum operator is defined by $$ P\phi=-i \frac{d}{dx}\phi(x) \qquad \forall \phi \in \mathscr{S}(\mathbb{R}) $$ Then we have the canonical commutation relations $$ [X,P]\phi=i\phi \qquad \forall \phi \in \mathscr{S}(\mathbb{R}) \tag{1}$$ For any unitary operator $U$ (for example the Fourier transform), it is easy to see that $U^{\dagger}PU$ and $U^{\dagger}XU$ satisfy the same commutation relations.
Question: If $2$ self-adjoint, densely defined linear operators satisfy (1), are they unitarily equivalent to the position and momentum operator?
In other words: Do the canonical commutation relations already determine position and momentum operator up to unitary equvialence?
I am asking this question out of curiousity, since it stuck in my head for quite some time. I have thought about (formally) taking operator exponentials, since $\exp(-iX)\phi=e^{ix}\phi(x)$ and $\exp(iaP)\phi=\phi(x+a)$ and then use the Stone-von Neumann theorem for the Heisenberg group, however, the issue with the domain remains and also taking operators exponential is non-trivial for this reason. Feel free to modify my domain to e.g. $D(\cdot)=W^{1,2}(\mathbb{R}),C_c^{\infty}(\mathbb{R}),..$ or whatever you see fit.
I am somewhat familiar with the basic spectral theory for unbounded operators and I have already encountered the concept of a rigged Hilbert space, so feel free to use them without elaborating every detail.
Here's a starting point that may help: $$ \left[X,P\right]=iI \\ (XP-PX)=iI \\ X(P-\lambda I)-(P-\lambda I)X = iI \\ (P-\lambda I)^{-1}X-X(P-\lambda I)^{-1}=i(P-\lambda I)^{-2}=\frac{d}{d\lambda}i(P-\lambda I)^{-1} \\ [(P-\lambda I)^{-1},X]=i\frac{d}{d\lambda}(P-\lambda I)^{-1}. $$ Then you can use Stone's Formula for the Spectral Measure $E$ associated with $P$: $$ s-\lim_{v\downarrow 0}\frac{1}{2\pi i}\int_{a}^{b}(P-(u+iv)I)^{-1}-(P-(u-iv)I)^{-1}du = \frac{1}{2}(E[a,b]+E(a,b)) $$ (The $s-\lim$ refers to convergence in the strong operator topology.)