# Hermiticity of Momentum Operator (matrix) Represented in Position Basis

I read that the momentum operator, $\hat P$ should be Hermitian (some would say by QM postulate). That makes perfect sense when it is represented in the momentum (p) basis: $\hat P =\int_{-\infty}^{\infty} \! p \,\, |p \rangle \langle p | \, \mathrm{d}p$, which certainly is Hermitian. ("$p$" is momentum eigenvalue)

It states that the momentum operator, in the position (x) basis, is $\hat P =-i\hbar\int_{-\infty}^{\infty} \! \partial/ \partial x \,\, |x \rangle \langle x | \, \mathrm{d}x$. That does not seem Hermitian, as the diagonal contains complex values.

Further, sometimes the momentum operator is used in the interaction Hamiltonian ($\hat H_{int}$), such as $\hat H_{int} = kA\hat P$ (where A is a variable of one system, $\hat P$ is the momentum operator of the other, and $k$ is a coupling constant. That could make $\hat H_{int}$ non-Hermitian, if the momentum operator is represented in the position basis (I also read that the Hamiltonian should be Hermitian, except in special cases, such as radioactive decay.)

Could someone please clarify how $\hat P =-i\hbar\int_{-\infty}^{\infty} \! \partial/ \partial x \,\, |x \rangle \langle x | \, \mathrm{d}x$ makes sense along with the statement that $\hat P$ and $\hat H$ should be Hermitian?

The following is a copy-paste from Dirac(1)

This gives $$\bigg< \phi\dfrac{\mathrm d \hphantom{q}}{\mathrm d q}\bigg| q'\bigg>=\boldsymbol{-}\dfrac{\mathrm d \phi\left(q'\right)}{\mathrm d q'}, \nonumber$$ showing that $$\bigg< \phi\dfrac{\mathrm d \hphantom{q}}{\mathrm d q}=\boldsymbol{-}\bigg<\dfrac{\mathrm d\phi}{\mathrm dq}. \tag{016}\label{eq016}$$ Thus $\:\mathrm d /\mathrm d q\:$ operating to the left on the conjugate complex of a wave function has the meaning of $\mathrm{minus}$ differentiation with respect to q.

And two paragraphs later

The conjugate complex of the linear operator $\:\mathrm d /\mathrm d q\:$ can be evaluated by noting that the conjugate imaginary of $\:\mathrm d /\mathrm d q.\psi\rangle\:$ or $\:\mathrm d\psi /\mathrm d q\rangle\:$ is $\:\langle\mathrm d \overset{\boldsymbol{\_\!\_}}{\psi}/\mathrm d q$, or $\:\boldsymbol{-}\langle\overset{\boldsymbol{\_\!\_}}{\psi}\,\mathrm d /\mathrm d q\;$ from \eqref{eq016}. Thus the conjugate complex of $\:\mathrm d /\mathrm d q\:$ is $\:\boldsymbol{-}\mathrm d /\mathrm d q\:$, so $\:\mathrm d /\mathrm d q\:$ is a pure imaginary linear operator.

(1) The Principles of Quantum Mechanics, P.A.M.Dirac, 4th edition 1967- $\S\;$22. Schröndiger's representation.

Let $\;p\;$ an operator in a finite $n-$dimensional Hilbert space represented by the finite $n\times n$ matrix $\mathrm{p}_{\mu\nu}$. The hermitian conjugate operator $\;p^{\boldsymbol{\dagger}}\;$ is represented by the matrix $\mathrm{p}^{\boldsymbol{\dagger}}_{\mu\nu}$, the hermitian conjugate of $\mathrm{p}_{\mu\nu}$ which is produced by transposing and complex conjugating $$\mathrm{p}^{\boldsymbol{\dagger}}_{\mu\nu}=\overset{\boldsymbol{\_\!\_}}{\mathrm{p}}_{\nu\mu} \tag{01}\label{eq01}$$ The operator $\;p\;$ is hermitian if its matrix $\mathrm{p}_{\mu\nu}$ is hermitian $$p^{\boldsymbol{\dagger}}=p \quad \Longleftarrow\!\Longrightarrow \quad \mathrm{p}^{\boldsymbol{\dagger}}_{\mu\nu}=\mathrm{p}_{\mu\nu} \tag{02}\label{eq02}$$

Now, if $\;q\;$ is the position operator for the $x-$axis and $\;p\;$ the operator of the $x-$component of the linear momentum in the $\;q-$basis, that is with respect to the complete set of the eigenstates of $\;q$, then
$$p=-i\hbar\dfrac{\mathrm d\hphantom{q}}{\mathrm dq} \tag{03}\label{eq03}$$ If we want to check the hermiticity of this momentum operator by transposing and complex conjugating the elements of a matrix then equation \eqref{eq03} is not convenient, as it is, since it's far away from being an equation between matrices.

In the following we'll see a matrix representation of \eqref{eq03} by which we'll prove the hermiticity of $\;p\;$ exactly by transposing and complex conjugating of the elements of a matrix. This is achieved by the basic quantum-mechanical equation $$qp-pq =i\hbar\, I \tag{04}\label{eq04}$$ and the Dirac $\delta-$function. It is known that there don't exist matrices $\;\mathrm{q}_{\mu\nu},\mathrm{p}_{\mu\nu},\mathrm{I}_{\mu\nu}\;$ representations of $\;q,p,I\;$ satisfying equation \eqref{eq04} with $\mu,\nu$ taking values from a finite or even countable set of values (any effort to succeed this would convince us). So these matrices must be of infinite uncountable dimensions.

Now, our first step would be to find a representation for the operator $\;q\;$ since we must express the momentum operator $\;p\;$ in the $\;q-$basis. This problem is equivalent to the determination of the eigenstates of $\;q$.

Before examine the infinite continuous case we'll try to build an example for a finite discrete case and see the correspondence. So, suppose that a particle could be only on a finite set of points on the $\;x-$axes, let $x_1=1,x_2=2, \cdots,x_n=n$. The space of states is the $\;n-$dimensional Hilbert space $\;\mathbb{C}^{n\;}$ so it's reasonable to take as eigenstates of the position operator the $n$ state vectors
$$\left(\mathbf{e}_\mu\right)_\nu=\delta_{\mu\nu}\quad \text{that is : } \mathbf{e}_1= \begin{bmatrix} 1\\ 0\\ \vdots\\ 0\\ 0 \end{bmatrix} \quad \mathbf{e}_2= \begin{bmatrix} 0\\ 1\\ \vdots\\ 0\\ 0 \end{bmatrix} \quad \cdots\cdots\cdots \quad \mathbf{e}_n= \begin{bmatrix} 0\\ 0\\ \vdots\\ 0\\ n \end{bmatrix} \tag{05}\label{eq05}$$ and since the operator $\;q\;$ has eigenvalues $\;1,2,\cdots,n\;$ it is represented by the diagonal matrix $$\mathrm q= \begin{bmatrix} 1 & \hphantom{n} &\hphantom{n} &\hphantom{n}&\hphantom{n} \\ \hphantom{n} & 2 &\hphantom{n} &\hphantom{n}&\hphantom{n} \\ \hphantom{n} & \hphantom{n} & \ddots &\hphantom{n}&\hphantom{n}\\ \hphantom{n} & \hphantom{n} &\hphantom{n} &\ddots&\hphantom{n}\\ \hphantom{n} & \hphantom{n} &\hphantom{n} &\hphantom{n} & n \end{bmatrix} \tag{06}\label{eq06}$$ expressed via the Kronecker $\delta$ $$\mathrm q_{\mu\nu}=\mu\delta_{\mu\nu} \tag{07}\label{eq07}$$ In \eqref{eq07} the Einstein's convention for repeated indices is not valid.

A general state $\;f\;$ is a $\;n-$dimensional complex vector represented with respect to the basis \eqref{eq05} by the one column matrix $$\mathrm f= \begin{bmatrix} f(1)\\ f(2)\\ \vdots\\ \vdots\\ f(n) \end{bmatrix} \tag{08}\label{eq08}$$

Application of the $\;q\;$ operator on $\;f\;$ produces a new state $\;g\;$ $$g=qf \tag{09}\label{eq09}$$ represented by the matrix $$\mathrm g= \begin{bmatrix} 1\cdot f(1)\\ 2\cdot f(2)\\ \vdots\\ \vdots\\ n\cdot f(n) \end{bmatrix} \tag{10}\label{eq10}$$ that is $$g(\mu)=\sum\limits_{\nu}\mathrm q_{\mu\nu}f(\nu)=\sum\limits_{\nu}\mu\delta_{\mu\nu}f(\nu)=\mu f(\mu) \tag{11}\label{eq11}$$

Now, returning to the infinite continuous case, a particle could be at any point on the real $\;x-$axis so we expect the $\;q\;$ to be represented in the basis of its own eigenstates by a diagonal $^\prime$improper$^\prime$ matrix with all real numbers as eigenvalues on the main diagonal as below
$$\mathrm q= \begin{bmatrix} \boldsymbol{-\infty} & \hphantom{n} &\hphantom{n} &\hphantom{n}&\hphantom{n}\vphantom{\dfrac{a}{b}} \\ \hphantom{n} & \nwarrow &\hphantom{n} &\hphantom{n}&\hphantom{n}\vphantom{\dfrac{a}{b}} \\ \hphantom{n} & \hphantom{n} & \mathbb{R} &\hphantom{n}&\hphantom{n}\vphantom{\dfrac{a}{b}}\\ \hphantom{n} & \hphantom{n} &\hphantom{n} &\searrow&\hphantom{n}\vphantom{\dfrac{a}{b}}\\ \hphantom{n} & \hphantom{n} &\hphantom{n} &\hphantom{n} & \boldsymbol{+\infty}\vphantom{\dfrac{a}{b}} \end{bmatrix} \tag{12}\label{eq12}$$ But the row and column indices $\;\mu,\nu\;$ which take discrete values in \eqref{eq07} must be replaced by the indices $\;x_{\boldsymbol{r}},x_{\boldsymbol{c}}\;$ which take continuous values in $\;\mathbb{R}$ as shown below \begin{align} &\begin{matrix} \hphantom{\cdots\cdots} &\hphantom{\!\!\! \cdots\cdots} &\hphantom{\!\!\! \cdots\cdots} & \:\:\:x_{\boldsymbol{c}} &\hphantom{\!\!\!\! \cdots\cdots} \end{matrix} \nonumber\\ &\begin{matrix} \hphantom{\cdots\cdots} &\hphantom{\!\!\! \cdots\cdots} &\hphantom{\!\!\! \cdots\cdots} & \:\:\:\:\downarrow &\hphantom{\!\!\!\! \cdots\cdots} \end{matrix} \nonumber\\ \mathrm q \left( x_{\boldsymbol{r}}, x_{\boldsymbol{c}}\right)\boldsymbol{=}\quad \begin{matrix} \vphantom{\dfrac{a}{b}} \\ x_{\boldsymbol{r}} \longrightarrow\vphantom{\dfrac{a}{b}} \\ \vphantom{\dfrac{a}{b}} \\ \vphantom{\dfrac{a}{b}} \\ \vphantom{\dfrac{a}{b}} \end{matrix}\!\! &\begin{bmatrix} \hphantom{n} & \hphantom{n} &\hphantom{n} &\vdots&\hphantom{n} \\ \cdots\cdots\!\!\!\! & \cdots\cdots\!\!\!\! &\cdots\cdots\!\!\!\! &\cdots\vdots\cdots &\!\!\!\! \cdots\cdots \vphantom{\dfrac{a}{b}}\\ \hphantom{n} & \hphantom{n} &\hphantom{n} &\vdots&\hphantom{n}\vphantom{\dfrac{a}{b}} \\ \hphantom{n} & \hphantom{n} &\hphantom{n} &\vdots&\hphantom{n} \\ \hphantom{n} & \hphantom{n} &\hphantom{n} &\vdots&\hphantom{n} \end{bmatrix} \tag{13}\label{eq13} \end{align} For notation reasons in \eqref{eq13} the row and column indices $\;x_{\boldsymbol{r}},x_{\boldsymbol{c}}\;$ are not placed as subscripts as $\;\mu,\nu\;$ do in $\;\mathrm q_{\mu\nu}\;$ but are placed in parentheses so that $\;\mathrm q \left( x_{\boldsymbol{r}}, x_{\boldsymbol{c}}\right)\;$ looks like a function of two real variables. Here the eigenstates of the position operator $\;q\;$ that correspond to those of the finite case, see equation \eqref{eq05}, are $$\left(\mathbf{e}_{ x_{\boldsymbol{r}}}\right)_ {x_{\boldsymbol{c}}}\boldsymbol{=}\delta\left( x_{\boldsymbol{r}}\!\boldsymbol{-}\!x_{\boldsymbol{c}}\right)\qquad x_{\boldsymbol{r}},x_{\boldsymbol{c}} \in \mathbb{R} \tag{14}\label{eq14}$$ where $\;\delta\left(x\right)\;$ the Dirac $\;\delta-$function.

The position operator $\;q\;$ is represented by the $^\prime$improper$^\prime$ matrix $$\mathrm q \left( x_{\boldsymbol{r}}, x_{\boldsymbol{c}}\right)\boldsymbol{=} x_{\boldsymbol{r}}\delta\left( x_{\boldsymbol{r}}\!\boldsymbol{-}\!x_{\boldsymbol{c}}\right) \tag{15}\label{eq15}$$ an expression that corresponds to equation \eqref{eq07} of the finite case.

The matrix \eqref{eq15} is hermitian due to the fact that the Dirac $\;\delta-$function is real and even, that is transposing and complex conjugating $$\mathrm{q}^{\boldsymbol{\dagger}}\left( x_{\boldsymbol{r}}, x_{\boldsymbol{c}}\right)\boldsymbol{=}\overset{\boldsymbol{\_\!\_}}{\mathrm{q}}\left( x_{\boldsymbol{c}}, x_{\boldsymbol{r}}\right)\boldsymbol{=} x_{\boldsymbol{r}}\delta\left( x_{\boldsymbol{c}}\!\boldsymbol{-}\!x_{\boldsymbol{r}}\right)\boldsymbol{=} x_{\boldsymbol{r}}\delta\left( x_{\boldsymbol{r}}\!\boldsymbol{-}\!x_{\boldsymbol{c}}\right)\boldsymbol{=}\mathrm q \left( x_{\boldsymbol{r}}, x_{\boldsymbol{c}}\right) \tag{16}\label{eq16}$$ Note that the identity operator $\;I\;$ is represented by $$\mathrm I \left( x_{\boldsymbol{r}}, x_{\boldsymbol{c}}\right)\boldsymbol{=} \delta\left( x_{\boldsymbol{r}}\!\boldsymbol{-}\!x_{\boldsymbol{c}}\right) \tag{17}\label{eq17}$$ A general state $\;f\;$ is a function $\;f\left(x\right)\;$ of the $\;x-$coordinate and application of the position operator $\;q\;$ produces an other state $\;g\;$ expressed in matrix form as $$g\boldsymbol{=}qf\boldsymbol{=}\!\!\!\!\int\limits_{\boldsymbol{-\infty}}^{\boldsymbol{+\infty}}\!\!\!\!\mathrm q \left( x_{\boldsymbol{r}}, x_{\boldsymbol{c}}\right)f\left(x_{\boldsymbol{c}}\right)\mathrm dx_{\boldsymbol{c}}\boldsymbol{=}\!\!\!\!\int\limits_{\boldsymbol{-\infty}}^{\boldsymbol{+\infty}}\!\!\!\!x_{\boldsymbol{r}}\delta\left( x_{\boldsymbol{r}}\!\boldsymbol{-}\!x_{\boldsymbol{c}}\right)f\left(x_{\boldsymbol{c}}\right)\mathrm dx_{\boldsymbol{c}}\boldsymbol{=}x_{\boldsymbol{r}}f\left(x_{\boldsymbol{r}}\right) \tag{18}\label{eq18}$$ that is $$g\left(x\right)\boldsymbol{=}xf\left(x\right) \tag{19}\label{eq19}$$ Equations \eqref{eq18}-\eqref{eq19} correspond to \eqref{eq11} of the finite case where the summation $\;\sum_\nu\;$ with respect to the index $\;\nu\;$ is replaced by the integral $\;\int\mathrm dx_{\boldsymbol{c}}\;$ with respect to the real index $\;x_{\boldsymbol{c}}$.

The final step will be to find the matrix representation of the linear momentum operator $\;p$, let it be $\;\mathrm{p}\left( x_{\boldsymbol{r}}, x_{\boldsymbol{c}}\right)$. So, let a general state $\;f$, that is a function $\;f\left(x\right)$. Application of the basic quantum-mechanical equation \eqref{eq04} on $\;f\;$ yields \begin{align} &\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\left(qp-pq\right)f =i\hbar\, I f \quad \Longrightarrow \nonumber\\ &\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\int\limits_{\boldsymbol{-\infty}}^{\boldsymbol{+\infty}}\!\!\!\!\mathrm q \left( x_{\boldsymbol{r}},u\right)\left(\:\int\limits_{\boldsymbol{-\infty}}^{\boldsymbol{+\infty}}\!\!\!\!\mathrm p \left( u,x_{\boldsymbol{c}}\right)f\left(x_{\boldsymbol{c}}\right)\mathrm dx_{\boldsymbol{c}}\right)\mathrm du\boldsymbol{-}\!\!\!\!\int\limits_{\boldsymbol{-\infty}}^{\boldsymbol{+\infty}}\!\!\!\!\mathrm p \left( x_{\boldsymbol{r}},v\right)\left(\:\int\limits_{\boldsymbol{-\infty}}^{\boldsymbol{+\infty}}\!\!\!\!\mathrm q \left( v,x_{\boldsymbol{c}}\right)f\left(x_{\boldsymbol{c}}\right)\mathrm dx_{\boldsymbol{c}}\right)\mathrm dv \nonumber\\ &\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\boldsymbol{=}i\hbar\!\!\!\int\limits_{\boldsymbol{-\infty}}^{\boldsymbol{+\infty}}\!\!\!\!\mathrm I \left( x_{\boldsymbol{r}}, x_{\boldsymbol{c}}\right)f\left(x_{\boldsymbol{c}}\right)\mathrm dx_{\boldsymbol{c}} \quad \stackrel{\eqref{eq15},\eqref{eq17}}{=\!=\!=\!\Longrightarrow} \nonumber\\ &\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\int\limits_{\boldsymbol{-\infty}}^{\boldsymbol{+\infty}}\!\!\!\!x_{\boldsymbol{r}}\delta\left( x_{\boldsymbol{r}}\!\boldsymbol{-}\!u\right)\left(\:\int\limits_{\boldsymbol{-\infty}}^{\boldsymbol{+\infty}}\!\!\!\!\mathrm p\left( u,x_{\boldsymbol{c}}\right)f\left(x_{\boldsymbol{c}}\right)\mathrm dx_{\boldsymbol{c}}\right)\mathrm du\boldsymbol{-}\!\!\!\!\int\limits_{\boldsymbol{-\infty}}^{\boldsymbol{+\infty}}\!\!\!\!\mathrm p \left( x_{\boldsymbol{r}},v\right)\left(\:\int\limits_{\boldsymbol{-\infty}}^{\boldsymbol{+\infty}}\!\!\!\!v\delta\left( v\!\boldsymbol{-}\!x_{\boldsymbol{c}}\right)f\left(x_{\boldsymbol{c}}\right)\mathrm dx_{\boldsymbol{c}}\right)\mathrm dv \nonumber\\ &\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\boldsymbol{=}i\hbar\!\!\!\int\limits_{\boldsymbol{-\infty}}^{\boldsymbol{+\infty}}\!\!\!\!\delta\left( x_{\boldsymbol{r}}\!\boldsymbol{-}\!x_{\boldsymbol{c}}\right)f\left(x_{\boldsymbol{c}}\right)\mathrm dx_{\boldsymbol{c}} \quad \Longrightarrow \nonumber\\ &\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! x_{\boldsymbol{r}}\!\!\!\int\limits_{\boldsymbol{-\infty}}^{\boldsymbol{+\infty}}\!\!\!\!\mathrm p \left( x_{\boldsymbol{r}}, x_{\boldsymbol{c}}\right)f\left(x_{\boldsymbol{c}}\right)\mathrm dx_{\boldsymbol{c}}\boldsymbol{-}\!\!\!\!\int\limits_{\boldsymbol{-\infty}}^{\boldsymbol{+\infty}}\!\!\!\!\mathrm p \left( x_{\boldsymbol{r}}, x_{\boldsymbol{c}}\right)x_{\boldsymbol{c}}f\left(x_{\boldsymbol{c}}\right)\mathrm dx_{\boldsymbol{c}}\boldsymbol{=}i\hbar\!\!\!\int\limits_{\boldsymbol{-\infty}}^{\boldsymbol{+\infty}}\!\!\!\!\delta\left( x_{\boldsymbol{r}}\!\boldsymbol{-}\!x_{\boldsymbol{c}}\right)f\left(x_{\boldsymbol{c}}\right)\mathrm dx_{\boldsymbol{c}} \quad \Longrightarrow \nonumber\\ &\int\limits_{\boldsymbol{-\infty}}^{\boldsymbol{+\infty}}\!\!\!\!\Bigl[\left( x_{\boldsymbol{r}}\!\boldsymbol{-}\!x_{\boldsymbol{c}}\right)\mathrm p \left( x_{\boldsymbol{r}}, x_{\boldsymbol{c}}\right)\!\boldsymbol{-}\!i\hbar\delta\left( x_{\boldsymbol{r}}\!\boldsymbol{-}\!x_{\boldsymbol{c}}\right)\Bigr]f\left(x_{\boldsymbol{c}}\right)\mathrm dx_{\boldsymbol{c}}=0 \tag{20}\label{eq20} \end{align} Since \eqref{eq20} must be valid for any state-function $\;f\left(x\right)\;$ the expression in brackets under the integral must be identically zero, so $$\bbox[yellow] {\:\:\:\mathrm p \left( x_{\boldsymbol{r}}, x_{\boldsymbol{c}}\right)\boldsymbol{=} i\hbar \dfrac{\delta\left( x_{\boldsymbol{r}}\!\boldsymbol{-}\!x_{\boldsymbol{c}}\right)}{\left( x_{\boldsymbol{r}}\!\boldsymbol{-}\!x_{\boldsymbol{c}}\right)}\boldsymbol{=}\boldsymbol{-}i\hbar\delta^{\,\prime}\!\left( x_{\boldsymbol{r}}\!\boldsymbol{-}\!x_{\boldsymbol{c}}\right)\boldsymbol{=} \boldsymbol{-}i\hbar \dfrac{\mathrm d\delta\left( x_{\boldsymbol{r}}\!\boldsymbol{-}\!x_{\boldsymbol{c}}\right)}{\mathrm d\left( x_{\boldsymbol{r}}\!\boldsymbol{-}\!x_{\boldsymbol{c}}\right)}\:\:\:\vphantom{\dfrac{\dfrac{a}{b}}{\dfrac{a}{b}}} } \tag{21}\label{eq21}$$ where $$\delta^{\,\prime}\!\left( x\right)\boldsymbol{=} \dfrac{\mathrm d\delta \left(x\right)}{\mathrm dx}\boldsymbol{=} \boldsymbol{-}\dfrac{\delta \left(x\right)}{x} \tag{22}\label{eq22}$$ the first derivative of the $\;\delta-$function.

Now, the hermiticity of the momentum operator is proved by transposing and complex conjugating the matrix in \eqref{eq21}

$$\mathrm{p}^{\boldsymbol{\dagger}}\left( x_{\boldsymbol{r}}, x_{\boldsymbol{c}}\right)\boldsymbol{=}\overset{\boldsymbol{\_\!\_}}{\mathrm{p}}\left( x_{\boldsymbol{c}}, x_{\boldsymbol{r}}\right)\boldsymbol{=} \boldsymbol{-}i\hbar \dfrac{\delta\left( x_{\boldsymbol{c}}\!\boldsymbol{-}\!x_{\boldsymbol{r}}\right)}{\left( x_{\boldsymbol{c}}\!\boldsymbol{-}\!x_{\boldsymbol{r}}\right)}\boldsymbol{=}i\hbar \dfrac{\delta\left( x_{\boldsymbol{r}}\!\boldsymbol{-}\!x_{\boldsymbol{c}}\right)}{\left( x_{\boldsymbol{r}}\!\boldsymbol{-}\!x_{\boldsymbol{c}}\right)}\boldsymbol{=}\mathrm p \left( x_{\boldsymbol{r}}, x_{\boldsymbol{c}}\right) \tag{23}\label{eq23}$$ The matrix \eqref{eq21} is hermitian due to the fact that the first derivative $\;\delta^{\,\prime}\!\left(x\right)\;$ of the Dirac $\;\delta-$function is real and odd and to the antihermiticity of the imaginary unit $\;i$.

Note also that equation \eqref{eq21} is a matrix representation of \eqref{eq03} since application of the momentum operator $\;p\;$ on a state $\;f\;$ yields $$pf=\!\!\!\!\int\limits_{\boldsymbol{-\infty}}^{\boldsymbol{+\infty}}\!\!\!\!\mathrm p \left( x_{\boldsymbol{r}}, x_{\boldsymbol{c}}\right)f\left(x_{\boldsymbol{c}}\right)\mathrm dx_{\boldsymbol{c}}\boldsymbol{=}\boldsymbol{-}i\hbar\!\!\!\int\limits_{\boldsymbol{-\infty}}^{\boldsymbol{+\infty}}\!\!\!\!\delta^{\,\prime}\!\left( x_{\boldsymbol{r}}\!\boldsymbol{-}\!x_{\boldsymbol{c}}\right)f\left(x_{\boldsymbol{c}}\right)\mathrm dx_{\boldsymbol{c}}\boldsymbol{=}\boldsymbol{-}i\hbar\dfrac{\mathrm df}{\mathrm dx} \tag{24}\label{eq24}$$

(Sorry I can't write comments!) it is simply because the $\partial/\partial x$ is not hermitian itself: Consider (in dirac notation) where f,g are square integrable functions (represented by bras and kets below)

$$\left<f \left| \frac{\partial}{\partial x}\right|g\right>:=\int_{-\infty}^{\infty}f \frac{\partial}{\partial x}[g] \,dx \; \overset{IBP}{=} [fg]_{-\infty}^{\infty} - \int_{-\infty}^{\infty}g \frac{\partial}{\partial x}[f] \,dx =: -\left<g \left| \frac{\partial}{\partial x}\right|f\right>$$ So it picks up a minus sign. But the extra minus sign(s) from the momentum operator get absorbed by the integral, giving you the result you need. I'll try answer the next bit below