The equation $(3)$ in the OP is formally1 correct, and it is in fact one of the main ingredients in the functional integral formulation of quantum field theory. See Ref.1 for the explicit construction.
For completeness, we sketch the derivation here. We use a notation much closer to the one in the OP than to that of Ref.1, but with some minors modifications (to allow for a more general result).
The setup.
Let $\{\phi_a,\pi^a\}_a$ be a set of phase-space operators, where $a\in\mathbb R^{d-1}\times \mathbb N^n$ is a DeWitt index (i.e., it contains a continuous part, corresponding to the spatial part of spacetime $\mathbb R^d$ and a discrete part, corresponding to a certain vector space $\mathbb N^n$ whose base is spacetime). Note that we are taking these operators to be in the Schrödinger picture. They are assumed conjugate:
$$\tag1
[\phi_a,\pi^b]=i\delta_a^b
$$
where $\delta$ is a Dirac-Kronecker delta. Here, $[\cdot,\cdot]$ denotes a commutator (we assume $\phi,\pi$ to be Grassmann even; we could consider the general case here by keeping tracks of the signs, but we won't for simplicity). The rest of commutators are assumed to vanish. We take the phase-space operators to be hermitian (or otherwise, we double the dimension and split them into their real and imaginary parts).
As $[\phi_a,\phi_b]=[\pi^a,\pi^b]=0$, we may simultaneously diagonalise them:
\begin{equation}
\begin{aligned}
\phi_a|\varphi\rangle&=\varphi_a|\varphi\rangle\\
\pi^a|\varpi\rangle&=\varpi^a|\varpi\rangle
\end{aligned}\tag2
\end{equation}
where $\varphi_a,\varpi^b$ are $c$-numbers. After normalising them, if necessary, these eigenvectors are orthonormal:
\begin{equation}
\begin{aligned}
\langle\varphi|\varphi'\rangle&=\prod_a\delta(\varphi_a-\varphi'_a)\equiv\delta(\varphi-\varphi')\\
\langle\varpi|\varpi'\rangle&=\prod_a\delta(\varpi^a-\varpi'^a)\equiv\delta(\varpi-\varpi')
\end{aligned}\tag3
\end{equation}
and, as per $[\phi_a,\pi^b]=\delta_a^b$, we also have
$$\tag4
\langle \varphi|\varpi\rangle=\prod_a\frac{1}{\sqrt{2\pi}}\mathrm e^{i\varphi_a\varpi^a}
$$
As the sets $\{\phi_a\}_a$ and $\{\pi^a\}_a$ are both assumed complete, we also have
\begin{equation}
\begin{aligned}
1&\equiv \int\prod_a|\varphi\rangle\langle\varphi|\,\mathrm d\varphi_a\\
1&\equiv \int\prod_a|\varpi\rangle\langle\varpi|\,\mathrm d\varpi^a
\end{aligned}\tag5
\end{equation}
as OP anticipated.
As mentioned in the comments, and for future reference, we note that there is another identity operator that is fundamental in a quantum theory, to wit, the resolution in terms of energy eigenstates:
\begin{equation}\tag6
1\equiv \int|E\rangle\langle E|\,\mathrm dE
\end{equation}
where $\mathrm dE$ is the counting measure in the case of discrete eigenvalues. Typically, we assume that the Hamiltonian is non-negative and that its ground state energy is zero; furthermore, we assume that there is a non-zero mass gap so that $E=0$ is a regular eigenvalue of $H$ (as opposed to a singular value).
It is at this point convenient to switch into the Heisenberg picture, where $\{\phi_a,\pi^a\}_a$ become time-dependent:
\begin{equation}
\begin{aligned}
\phi_a(t)\equiv\mathrm e^{iHt}\phi_a\mathrm e^{-iHt}\\
\pi^a(t)\equiv\mathrm e^{iHt}\pi^a\mathrm e^{-iHt}
\end{aligned}\tag7
\end{equation}
with eigenstates
\begin{equation}
\begin{aligned}
|\varphi;t\rangle&=\mathrm e^{iHt}|\varphi\rangle\\
|\varpi;t\rangle&=\mathrm e^{iHt}|\varpi\rangle
\end{aligned}\tag8
\end{equation}
such that
\begin{equation}
\begin{aligned}
\phi_a(t)|\varphi;t\rangle&=\varphi_a|\varphi;t\rangle\\
\pi^a(t)|\varpi;t\rangle&=\varpi^a|\varpi;t\rangle
\end{aligned}\tag9
\end{equation}
The Heisenberg picture eigenstates satisfy the same completeness and orthonormality relations as the Schrödinger picture eigenstates (inasmuch as the transformation is unitary).
The functional integral. Time slicing.
From this, and by the usual arguments (time slicing), Ref.1 derives the phase-space functional integral representation of the transition amplitude, to wit,
\begin{equation}
\begin{aligned}
&\langle\varphi_\mathrm{in};t_\mathrm{in}|\mathrm T\left(O_1(\pi(t_1),\phi(t_1)),\dots,O_n(\pi(t_n),\phi(t_n))\right)|\varphi_\mathrm{out};t_\mathrm{out}\rangle\equiv\\
&\hspace{10pt}\int_{\varphi(t_\mathrm{in})=\varphi_\mathrm{in}}^{\varphi(t_\mathrm{out})=\varphi_\mathrm{out}}\left(O_1(\varpi(t_1),\varphi(t_1)),\dots,O_n(\varpi(t_n),\varphi(t_n))\right)\cdot\\
&\hspace{25pt}\cdot\exp\left[i\int_{t_\mathrm{in}}^{t_\mathrm{out}}\left(\sum_a\dot \varphi_a(\tau)\varpi^a(\tau)-H(\varphi(\tau),\varpi(\tau))\right)\mathrm d\tau\right]\mathrm d\varphi\,\mathrm d\varpi
\end{aligned}\tag{10}
\end{equation}
where $O_1,\dots,O_n$ is any set of operators; $\mathrm T$ denotes the time-ordering symbol; and $\mathrm d\varphi,\mathrm d\varpi$ denote the measures
$$\tag{11}
\mathrm d\varphi\equiv\prod_{\tau,a}\mathrm d\varphi_a(\tau),\qquad \mathrm d\varpi\equiv\prod_{\tau,a}\frac{1}{2\pi}\mathrm d\varpi^a(\tau)
$$
The time-slicing procedure is standard. We only consider the case $O_i=1$ here. We begin by considering the case where $t_\mathrm{in}$ and $t_\mathrm{out}$ are infinitesimally close:
$$\tag{12}
\langle\varphi';\tau+\mathrm d\tau|\varphi;\tau\rangle=\langle\varphi';\tau|\exp\left[-iH\mathrm d\tau\right]|\varphi;\tau\rangle
$$
where $H=H(\phi(\tau),\pi(\tau))$ is the Hamiltonian (the generator of time translations, essentially defined by this equation). We take the convention that all the $\phi$ must always be moved to the left of $\pi$. In this case, and up to factors of order $\mathrm O(\mathrm d\tau)^2$, we may replace $\phi$ by its eigenvalue, to wit, $\varphi$. To deal with $\pi$, we insert the identity $1$ in the form of the completeness relation:
$$\tag{13}
\langle\varphi';\tau+\mathrm d\tau|\varphi;\tau\rangle\overset{(5)}=\int\exp\left[-iH(\varphi',\varpi)\mathrm d\tau+i\sum_a(\varphi'_a-\varphi_a)\varpi^a\right]\ \mathrm d\varpi
$$
where each $\varpi^a$ is integrated over $\mathbb R$ unrestrictedly.
To find the transition amplitude over a finite interval, we just compose an infinite number of infinitesimal transition amplitudes: we break up $t'-t$ into steps $t,\tau_1,\tau_2,\dots,\tau_N,t'$, with $\tau_{k+1}-\tau_k=\mathrm d\tau=(t'-t)/(N+1)$. With this, and inserting the identity $1$ in the form of the completeness relation at each $\tau_k$, we get
\begin{equation}
\begin{aligned}
\langle\varphi';t'|\varphi;t\rangle&\overset{(5)}=\int\langle \varphi';t'|\varphi_N;\tau_N\rangle\langle \varphi_N;\tau_N|\varphi_{N-1};\tau_{N-1}\rangle\cdots\langle \varphi_1;\tau_1|\varphi;t\rangle\ \prod_{k=1}^N\mathrm d\varphi_k\\
&\overset{(13)}=\int\left[\prod_{k=1}^N\prod_a\mathrm d\varphi_{k,a}\right]\left[\prod_{k=0}^N\prod_a\frac{\mathrm d\varpi^a_k}{2\pi}\right]\cdot\\
&\hspace{20pt}\cdot\exp\left[i\sum_{k=1}^{N+1}\left(\sum_a(\varphi_{k,a}-\varphi_{k-1,a})\varpi_{k-1}^a-H(\varphi_k,\varpi_{k-1})\mathrm d\tau\right)\right]
\end{aligned}\tag{14}
\end{equation}
where $\varphi_0\equiv\varphi$ and $\varphi_{N+1}\equiv\varphi'$. By taking the formal limit $N\to\infty$, we indeed obtain the claimed formula. The generalisation of the proof to include insertions is straightforward.
Vacuum-to-vacuum transition amplitude. Feynman's $\boldsymbol{+i\epsilon}$ prescription.
In non-relativistic quantum mechanics, the functional integral as written above is the natural object to work with. On the other hand, when dealing with particle physics in the relativistic regime, one usually works with $S$ matrix elements, that is, one considers the transition amplitude, not from eigenstates of $\phi_a$, but of the creation and annihilation operators. As is well-known from the LSZ theorem, it suffices to consider the vacuum-to-vacuum transition amplitude,
\begin{equation}\tag{15}
\langle 0|\mathrm T\mathrm e^{iJ^a\phi_a}|0\rangle
\end{equation}
and from which all $S$-matrix elements can be computed. We thus want to obtain the functional integral representation of this transition amplitude.
In Ref.1 there is a rather clean derivation of such object which, unfortunately, is only worked out for a scalar boson field. The outcome is that the vacuum-to-vacuum transition amplitude is given by the functional integral over all field configurations, and the correct boundary conditions are enforced by Feynman's $+i\epsilon$ prescription. The higher spin case is non-trivial because the propagators (and ground-state wave-functionals) are gauge-dependent. Indeed, where to put the $+i\epsilon$ in an arbitrary gauge theory is a very complicated matter (e.g., in the axial gauge it is far from clear what to do with, say, $k^4$: should it be $(k^2+i\epsilon)^2$? Should it be $k^4+i\epsilon$?), and the general case has not been worked out to the best of my knowledge. It is nevertheless possible to argue that, at least in the 't Hooft-Feynman gauge $\xi=1$, the propagators and ground-state wave-functional are identical to the scalar case (up to a unit matrix in colour space), so that the derivation in the reference holds for a bosonic field of arbitrary spin. The fermionic case requires Grassmann integration, but a similar analysis (in the $\xi=1$ case) is possible. Gauge invariance then implies the general case.
For completeness, we will prove the claim by an alternative method which, although admittedly not nearly as clean, works for a field of arbitrary spin. The trick is to use the completeness relation in terms of energy eigenstates instead of $\phi_a$ eigenstates.
We proceed as follows. We want to calculate
$$\tag{16}
\lim_{t\to-\infty}\langle \varphi_\mathrm{in};t|\overset{(6)}=\lim_{t\to-\infty}\int\mathrm e^{-iEt}\langle \varphi_\mathrm{in}|E\rangle\langle E|\,\mathrm dE
$$
If we send $t$ to $-\infty$ in a slightly imaginary direction, then all excited states acquire a real and negative part in the exponential factor, which vanishes in the large $t$ limit. We are thus left with the ground state only:
$$\tag{17}
\lim_{t\to-\infty+i\epsilon}\langle \varphi_\mathrm{in};t|=\langle \varphi_\mathrm{in}|0\rangle\langle 0|
$$
Similarly,
$$\tag{18}
\lim_{t\to+\infty+i\epsilon}|\varphi_\mathrm{out};t\rangle=|0\rangle\langle 0|\varphi_\mathrm{out}\rangle
$$
With this, the matrix element $\langle\varphi_\mathrm{in};-\infty|O|\varphi_\mathrm{out};+\infty\rangle$ can be written as
\begin{equation}
\begin{aligned}
&\langle \varphi_\mathrm{in}|0\rangle\langle 0|O|0\rangle\langle 0|\varphi_\mathrm{out}\rangle
=
\int_{\varphi(-\infty)=\varphi_\mathrm{in}}^{\varphi(+\infty)=\varphi_\mathrm{out}}O\cdot\\
&\hspace{20pt}\cdot\exp\left[i\int_{(1+i\epsilon)\mathbb R}\left(\sum_a\dot \varphi_a(\tau)\varpi^a(\tau)-H(\varphi(\tau),\varpi(\tau))\right)\mathrm d\tau\right]\mathrm d\varphi\,\mathrm d\varpi
\end{aligned}\tag{19}
\end{equation}
Integrating both sides with respect to $\mathrm d\varphi_\mathrm{in}\,\mathrm d\varphi_\mathrm{out}$, we get the vacuum-to-vacuum transition amplitude in its standard form, where the integral over $\mathrm d\varphi$ is unrestricted:
\begin{equation}\tag{20}
\langle 0|O|0\rangle
=N^{-1}
\int O\exp\left[iS(\varphi,\varpi)\right]\mathrm d\varphi\,\mathrm d\varpi
\end{equation}
where
\begin{equation}\tag{21}
S(\varphi,\varpi)\equiv\int_{(1-i\epsilon)\mathbb R}\left(\sum_a\dot \varphi_a(\tau)\varpi^a(\tau)-H(\varphi(\tau),\varpi(\tau))\right)\mathrm d\tau
\end{equation}
is the classical action, and where
\begin{equation}\tag{22}
N\equiv \int\langle \varphi_\mathrm{in}|0\rangle\langle 0|\varphi_\mathrm{out}\rangle\,\mathrm d\varphi_\mathrm{in}\,\mathrm d\varphi_\mathrm{out}
\end{equation}
is an inconsequential normalisation constant (it is the norm of the ground-state wave-functional). This proves the claim: the vacuum-to-vacuum transition amplitude is given by the standard functional integral, but over all field configurations (unrestricted integral over $\mathrm d\varphi$); and the correct boundary conditions are essentially those that result from a Wick rotation $\tau\to-i\tau_\mathrm{E}$.
The configuration space functional integral. The Lagrangian.
Finally, it bears mentioning how the configuration space functional integral is obtained. Ref.1 considers the case where the Hamiltonian is a quadratic polynomial in $\varpi$:
\begin{equation}\tag{23}
H(\varphi,\varpi)=\frac12\varpi^a A_{ab}(\varphi)\varpi^b+B_a(\varphi)\varpi^a+C(\varphi)
\end{equation}
In this case, the integral over $\mathrm d\varpi$ is gaussian and so its stationary phase approximation is in fact exact. The stationary point $\varpi^\star$ is easily computed to be
\begin{equation}\tag{24}
\dot\varphi_a=\frac{\partial H}{\partial\varpi^a}\bigg|_{\varpi\to\varpi^\star}
\end{equation}
which agrees with the classical canonical relation. Therefore, the Hamiltonian at $\varpi^\star=\varpi^\star(\dot\varphi)$ is nothing but the Lagrangian $L=L(\varphi,\dot \varphi)$, and therefore
\begin{equation}\tag{25}
\langle 0|O|0\rangle
\propto
\int O\exp\left[iS(\varphi)\right]\mathrm d\varphi
\end{equation}
where
\begin{equation}\tag{26}
S(\varphi)\equiv \int_{(1-i\epsilon)\mathbb R}L(\varphi,\dot\varphi)
\end{equation}
and where $\mathrm d\varphi$ implicitly includes the determinant of the Vilkovisky metric,
\begin{equation}\tag{27}
\mathrm d\varphi\to\sqrt{\det(A)}\,\mathrm d\varphi
\end{equation}
which is required for covariance in configuration space (or, equivalently, unitarity). If the metric is flat the determinant can be reabsorbed into $N$.
If $H$ is not a quadratic polynomial, we may nevertheless use the stationary phase approximation, but the measure will acquire higher order corrections:
\begin{equation}\tag{28}
\mathrm d\varphi\to\left(\sqrt{\det(A)}+\mathcal O(\hbar)\right)\,\mathrm d\varphi
\end{equation}
which can be computed order by order in perturbation theory. This proves that the configuration space functional integral can always be made both covariant and unitary by carefully taking care of the integration measure.
In any case, there is a subtlety that must be mentioned: the integral over $\mathrm d\varpi$ as written above is only valid if $O$ depends on $\varphi$ only; for otherwise the integral is not gaussian. In other words, we may not use the configuration space functional integral to compute matrix elements of derivatives of $\phi$. In pragmatical terms, this is easy to understand: the time ordering symbol $\mathrm T$ does not commute with time-derivatives. The resolution involves the introduction of the so-called covariant time-ordering symbol which is defined such that it commutes with time derivatives (cf. this PSE post).
References.
- Weinberg, QFT, Vol.I, chapter 9.
1: Formal in the sense that this is not a rigorous statement, inasmuch as the whole functional-integral formalism is not rigorous itself. It seems hard to formalise the sum over all fields, but one may argue that the picture is at least consistent.
Best Answer
You have gone awry at the end of your question, utilizing ambiguous/meaningless notation.
Here, it is best to use the caret $\hat {\mathbf v}$ to signify a unit vector, as in Sakurai-Napolitano (3.6.22,23), $$ \langle {\mathbf x}|nlm\rangle = R_{nl}(r) Y^m_l(\theta,\phi)= R_{nl}(r) \langle \hat {\mathbf x}|l,m\rangle , $$ and not an operator, as in the punchline of your question: use another symbol for operators, instead.
Consequently, by inspection, $$ \langle r\rangle =\int r^2 dr ~~ r (R_{nl}(r))^2 , $$ since the $d\Omega$ integral factored out to 1.
The linear power of r in this expression is just just a c-number, the radial value (length) of the c-number vector ${\mathbf x}= r \hat {\mathbf x}$, since you are working in the coordinate representation.
If you must think of an operator with eigenvalue r, it is $\large \mathbb r$ s.t. $$ {\large \mathbb r} ~ |{\mathbf x}\rangle= r |{\mathbf x}\rangle , $$ and you could apply the machinery of Dirac's ket notation in the above, which, however, is superfluous, and, evidently, confusing here.