What's the commutator of $|\mathbf{\hat{x}}|$ and $|\mathbf{\hat{p}}|$? Straight to the point: what's the result of the commutator of the magnitude of the position and  the momentum operators and how can I approach it, i.e., $[|\mathbf{\hat{x}}|,|\mathbf{\hat{p}}|]=$ ?
My efforts: (1) trying to use $|\mathbf{\hat{x}}|=\sqrt{\sum \hat{x}_i\hat{x}_i}$ doesn't seem to help because of the square root; (2) identifying $|\mathbf{\hat{x}}|$ as the radial position operator $\hat{r}$ and using the "radial momentum" operator $\hat{p}_r$ was of no help because $\hat{p}_r \neq |\mathbf{\hat{p}}|$ and $\hat{p}_r$ is the square root of a sum relating $|\mathbf{\hat{p}}|^2$ and the square of the angular momentum operator.
ADDENDUM (Oct 19th, 2020)
I see there may be a number of technical difficulties defining $|\mathbf{\hat{x}}|$ and $|\mathbf{\hat{p}}|$ because these are square roots of operators $\hat{x}_i$ and $\hat{p}_i$. To be clearer, I'm looking for a formal expression for $[|\mathbf{\hat{x}}|,|\mathbf{\hat{p}}|]$ that might (but not necessarily should) overlook technical issues regarding square root of operators. For instance, I wonder if that's the case of the well-known result $[x_i,F(p_x,p_y,p_z)] = i\hbar \frac{\partial F}{\partial p_i}$ if we take $F=|\mathbf{\hat{p}}|$ so that $[x_i,|\mathbf{\hat{p}}|] =i\hbar\frac{p_i}{|\mathbf{\hat{p}}|}$.
 A: Here are some incomplete untested ideas.

*

*Define operators
$$\begin{align}J_-~:=&~\frac{1}{2}{\bf x}^2,\qquad
J_+~:=~\frac{1}{2}{\bf p}^2,\cr
J_z~:=&~\frac{1}{4}\sum_{j=1}^3 \{x^j, p_j\}_+.\end{align}\tag{1}$$
One may prove that they form an $sl(2,\mathbb{R})$ Lie algebra
$$ [J_-,J_+]~=~2i\hbar J_z, \qquad [J_z,J_{\pm}]~=~i\hbar J_{\pm}.\tag{2}$$
In this language OP wants to calculate the commutator $2[\sqrt{J_-},\sqrt{J_+}]$.


*Define normalization
$$\sigma_{\pm}~:=~\frac{J_{\pm}}{\sqrt{2}\hbar}, \qquad \sigma_z~:=~\frac{J_z}{i\hbar}.\tag{3}$$
Then
$$ [\sigma_+,\sigma_-]~=~\sigma_z, \qquad [\sigma_z,\sigma_{\pm}]~=~\sigma_{\pm}.\tag{4}$$
The $2\times 2$ Pauli matrices satisfy the same Lie algebra (although the underlying associative algebra is different).


*It might be possible to define a  notion of square root operators similar to this Phys.SE post.


*It might be possible to adapt a representation  a la Holstein-Primakoff or Dyson-Maleev, cf. this Phys.SE post, or perhaps some squeezed state methods.
A: In quantum mechanics, there are different possibilities to express an operator. Let's consider a matrix of the operator $\hat{A} = [|\hat{\bf{x}}|, |\hat{\bf{p}}|]$ in coordinate reperesentation:
$$
A(\bf{x},\bf{x'}) = \langle \bf{x} | \hat{A}|\bf{x'}\rangle = (|\bf{x}| - |\bf{x'}|) \langle \bf{x} | |\hat{\bf{p}}||\bf{x'}\rangle
$$
According to the spectral representation of $\bf{p}$, we have
$$
\langle \bf{x} | |\hat{\bf{p}}||\bf{x'}\rangle = {\rm \frac1{(2\pi\hbar)^3}} \int |\bf{p}| e^{\frac{i}\hbar\bf{p}(\bf{x}-\bf{x'})} \bf{d^3 p} {\rm\ \equiv \frac1{(2\pi\hbar)^3} f(}{\bf x - x'}{\rm )}
$$
Further, I shall treat the function
$$
f(\bf{R}) = \int |\bf{p}| e^{\frac{i}\hbar\bf{p}\bf{R}} \bf{d^3 p}
$$
as a generalized function. Then we have
$$
f(\bf{R}) = -\hbar^2\Delta_{\bf{R}} \int \frac1{|\bf{p}|} e^{\frac{i}\hbar\bf{p}\bf{R}} \bf{d^3 p} = {\rm -4\pi\hbar^4} \Delta_{\bf{R}}{\rm \frac1{{\bf R}^2}}\quad {\rm (1)}
$$
It is well known that
$$
-\Delta \frac1{|\bf{R}|} = {\rm 4\pi\delta(}\bf{R}{\rm )}
$$
If analogous local representation exists for the function (1), then it might be possible to express $\hat{A}$ as a simple combination of $\hat{\bf x}$ and $\hat{\bf p}$ operators. Otherwise, I suppose there is no simple expression for $\hat{A}$.
A: I'll provide a tentative solution based on a result derived by M. K. Transtrum and J.-F. S. Van Huele, J. Math. Phys. 46, 063510 (2005). They derived a general expression for the commutator of functions $f(A,B)$ and $g(A,B)$ of noncommuting operators $A$ and $B$:
\begin{equation}
\left[f(A,B),g(A,B)\right] = \sum_{k=1}^\infty \frac{(-c)^k}{k!} \left( \frac{\partial^k g}{\partial A^k} \frac{\partial^k f}{\partial B^k} - \frac{\partial^k f}{\partial A^k} \frac{\partial^k g}{\partial B^k}\right),
\quad \text{where} \quad
c=[A,B].
\end{equation}
Tentative solution
I consider the particular case $f=f(A)$ and $g=g(B)$:
\begin{equation}
\left[f(A),g(B)\right] = - \sum_{k=1}^\infty \frac{(-c)^k}{k!} \frac{\partial^k f}{\partial A^k} \frac{\partial^k g}{\partial B^k}
= \left[ - \sum_{k=1}^\infty \frac{(-c)^k}{k!} \frac{\partial^k }{\partial A^k} \frac{\partial^k }{\partial B^k} \right] f(A)g(B),
\end{equation}
where I believe the last step is not problematic as long we understand its meaning: derivative $\partial_A\equiv\frac{\partial}{\partial_A}$ acts on $f(A)$ and derivative $\partial_B\equiv\frac{\partial}{\partial_B}$ acts on $g(B)$. Finally we simplify the result to
\begin{equation}
[f(A),g(B)] = \left( 1-e^{-c\partial_A \partial_B} \right)f(A)g(B) \quad \text{or} \quad g(B)f(A)=e^{-c\partial_A \partial_B}f(A)g(B).
\end{equation}
Heading to $[|\mathbf{\hat{x}}|,|\mathbf{\hat{p}}|]$, I'll now ommit the "hats" for simplicity and use the notation $\mathbf{x}=(x,y,z)$ and $\mathbf{p}=(p_x,p_y,p_z)$ for the position and momentum operators respectively. My approach is writing |\mathbf{x}| and |\mathbf{p}| as power series:
\begin{equation}
|\mathbf{x}| = \sum_{abc}A_{abc} x^a y^b z^c \quad \text{and} \quad |\mathbf{p}| = \sum_{uvw}B_{uvw} p_x^u p_y^v p_z^w.
\end{equation}
Then
\begin{equation}
[|\mathbf{x}|,|\mathbf{p}|] = \sum_{abc}\sum_{uvw} A_{abc}B_{uvw} (x^a y^b z^c p_x^u p_y^v p_z^w - p_x^u p_y^v p_z^w x^a y^b z^c).
\end{equation}
The last term can be recast as
\begin{equation}
p_x^u p_y^v p_z^w x^a y^b z^c = (p_x^u x^a) (p_y^v y^b) (p_z^w z^c) = (e^{-i\hbar\partial_x \partial_{p_x}} x^a p_x^u) (e^{-i\hbar\partial_y \partial_{p_y}} y^b p_y^v) (e^{-i\hbar\partial_z \partial_{p_z}} z^c p_z^w) = e^{-i\hbar(\partial_x \partial_{p_x} + \partial_y \partial_{p_y} + \partial_z \partial_{p_z})} x^a p_x^u y^b p_y^v z^c p_z^w = e^{-i\hbar\partial_\mathbf{x}\cdot\partial_\mathbf{p}} x^a p_x^u y^b p_y^v z^c p_z^w,
\end{equation}
where
\begin{equation}
\partial_\mathbf{x}\cdot\partial_\mathbf{p} \equiv \sum_i \frac{\partial}{\partial x_i} \frac{\partial}{\partial p_i}.
\end{equation}
Finally,
\begin{equation}
[|\mathbf{x}|,|\mathbf{p}|] = \sum_{abc}\sum_{uvw} A_{abc}B_{uvw} (1-e^{-i\hbar\partial_\mathbf{x}\cdot\partial_\mathbf{p}})x^a y^b z^c p_x^u p_y^v p_z^w = (1-e^{-i\hbar\partial_\mathbf{x}\cdot\partial_\mathbf{p}})|\mathbf{x}| |\mathbf{p}|
\end{equation}
or, explicitly,
\begin{equation}
[|\mathbf{x}|,|\mathbf{p}|] = - \sum_{n=1}^\infty \frac{(-i\hbar)^n}{n!}(\partial_\mathbf{x} \cdot \partial_\mathbf{p})^n |\mathbf{x}| |\mathbf{p}|.
\end{equation}
A few remarks:

*

*The one-dimensional version of this equation recovers the expected expression given in the beginning of this answer.


*The r.h.s. is rotationally invariant as expected considering the l.h.s.
Prospects:

*

*Maybe a nice, closed form for the above result could be achieved writing $\partial_\mathbf{x} \cdot \partial_\mathbf{p}$ in spherical polar coordinates?


*This derivation seems to work without modifications for any $[|\mathbf{x}|^n,|\mathbf{p}|^m]$. If the last remark is met with success, we can check if the proposed answer recovers some commutators that can be computed easily -- e.g., one with $n=m=2$.
I'm investigating these prospects and I will modify this answer accordingly.
A: Here I'll overlook technical issues that may be relevant when defining operators $|\mathbf{x}|$ and $|\mathbf{p}|$ and head toward a formal expression for $[|\mathbf{x}|,|\mathbf{p}|]$. I'll omit the "hat" above operators for simplicity.
I start defining the symmetrized radial momentum operator,
$$p_r \equiv \frac{1}{2}\left( \frac{\mathbf{x}}{|\mathbf{x}|} \cdot \mathbf{p} + \mathbf{p} \cdot \frac{\mathbf{x}}{|\mathbf{x}|}
\right).$$
This operator came into my attention on Alvarez & González, Am. J. Phys. 57, 923 (1989) and Liboff, Nebenzahl & Fleischmann, Am. J. Phys. 41, 976 (1973) but I learned it is standard content in QM textbooks. In position space represention $\mathbf{p}=-i\hbar\mathbf{\nabla}$ and spherical coordinates, straightforward calculation reveals
$$p_r = -i\hbar \left(\frac{\partial}{\partial r} + \frac{1}{r} \right),$$
where $r\equiv|\mathbf{x}|$. In this last form, the following result is easily verified:
$$[|\mathbf{x}|,p_r]=i\hbar.$$
From the relation $[A,f(B)]=[A,B]\frac{\partial f}{\partial B}$ valid if $[A,[A,B]]=0$ we conclude
$$[|\mathbf{x}|,f(p_r)]=i\hbar \frac{\partial f}{\partial p_r}.$$
Now I set $f(p_r) = |\mathbf{p}| = \sqrt{p_r^2 + \frac{L^2}{r^2}}$, with the last equality coming from the splitting of $|\mathbf{p}|$ into radial and angular parts, and with $L$ the angular momentum. From the previous commutator, we get
$$[|\mathbf{x}|,|\mathbf{p}|]=i\hbar \frac{1}{|\mathbf{p}|}\left[p_r + \frac{1}{2}\frac{\partial}{\partial p_r}\left(\frac{L^2}{r^2}\right) \right] = i\hbar \frac{p_r}{|\mathbf{p}|}.$$
The last derivative vanished, as I see it, because $L^2$ comprises only of angular pieces of $|\mathbf{p}|$ but no radial piece and because $\partial/\partial p_r$ doesn't act on $r$ as we can check applying it to $[|\mathbf{x}|,p_r]=i\hbar$. Finally, with the definition of $p_r$ above, we get
$$\boxed{
[|\mathbf{x}|,|\mathbf{p}|] = \frac{i\hbar}{2}\frac{1}{|\mathbf{p}|} \left( \frac{\mathbf{x}}{|\mathbf{x}|} \cdot \mathbf{p} + \mathbf{p} \cdot \frac{\mathbf{x}}{|\mathbf{x}|}
\right) = i\hbar \left( \frac{\mathbf{p}}{|\mathbf{p}|} \cdot \frac{\mathbf{x}}{|\mathbf{x}|} + \frac{i\hbar}{|\mathbf{p}| |\mathbf{x}|}
\right).
}$$
A few remarks:

*

*This result is rotationally invariant as expected.


*For the cases I checked, the derivation presented here predicts results derived by other means. For instance, $[|\mathbf{x}|,|\mathbf{p}|^2]=[|\mathbf{x}|,\sum p_i p_i] = \sum\left( [|\mathbf{x}|,p_i] p_i + p_i [|\mathbf{x}|,p_i]\right) = i\hbar\left( \frac{\mathbf{x}}{|\mathbf{x}|}\cdot\mathbf{p} + \mathbf{p}\cdot\frac{\mathbf{x}}{|\mathbf{x}|}\right)$ and the derivation presented here gets the same, $[|\mathbf{x}|,|\mathbf{p}|^2]=i\hbar\frac{\partial |\mathbf{p}|^2}{\partial p_r} = 2i\hbar p_r = i\hbar \left( \frac{\mathbf{x}}{|\mathbf{x}|} \cdot \mathbf{p} + \mathbf{p} \cdot \frac{\mathbf{x}}{|\mathbf{x}|}
\right).$
