can someone point me in the right direction at least to diagonalize the modified Hamiltonian?
This answer explains how to diagonalize the Hamiltonian and how to derive the $n$-dependent term in $M_n$. That's the interesting part, because the $n$-independent part is just an overall constant term in the Hamiltonian, which has no observable consequences in a non-gravitational theory.
Eliminating the linear terms
Sometimes a little extra generality can help. Start with
$$
H = \frac{1}{2}\int dx\ H(x)
\tag{1}
$$
where
$$
\newcommand{\pl}{\partial}
H(x) = \pi^2+(\phi')^2+V(\phi)
\tag{2}
$$
with $\phi'\equiv\pl\phi/\pl x$. Suppose $\phi_c$ solves
$$
2\phi_c'' - V'(\phi_c)=0
\tag{3}
$$
with $V'\equiv \pl V/\pl\phi$. Replace $\phi\to\phi_c+\phi$ in (2) to get
$$
H(x) = \pi^2+(\phi')^2+2\phi'\phi_c'+V'(\phi_c)\phi
+\frac{V''(\phi_c)}{2}\phi^2 + \cdots
\tag{4}
$$
where the omitted terms are either independent of $\phi$ or are of order $\phi^3$ and higher. After subsituting this back into (1), we can integrate-by-parts to convert the term $2\phi'\phi_c'$ to $-2\phi''\phi_c$. Then equation (3) says that the terms linear in $\phi$ add up to zero, which leaves
$$
H = \frac{1}{2}\int dx\
\left(\pi^2+(\phi')^2 +\frac{V''(\phi_c)}{2}\phi^2 + \cdots\right).
\tag{5}
$$
As before, the terms represented by "$\cdots$" are either independent of $\phi$ or are of order $\phi^3$ and higher.
The $x$-dependent "mass" term
The Hamiltonian (5) looks like a harmonic oscillator but with an $x$-dependent coefficient
$$
\frac{V''(\phi_c)}{2}
\tag{6}
$$
for the $\phi^2$ term (the "mass" term). The function $V(\phi)$ in the question is
$$
V(\phi)=-\mu^2\phi^2+\frac{\lambda}{2}\phi^4,
\tag{7}
$$
which gives
$$
\frac{V''(\phi_c)}{2}=-\mu^2+3\lambda\phi_c^2.
\tag{8}
$$
The question also provides this expression for $\phi_c$:
$$
\phi_c=\frac{\mu}{\sqrt{\lambda}}
\tanh\left(\frac{\mu x}{\sqrt{2}}\right).
\tag{9}
$$
Use (9) in (8) to get
$$
\frac{V''(\phi_c)}{2}
=\mu^2\left(3\tanh^2\left(\frac{\mu x}{\sqrt{2}}\right)-1\right).
\tag{10}
$$
The goal is to diagonalize the Hamiltonian (5) with $V''(\phi_c)$ given by (10).
Diagonalization
Here's a summary of what we've derived so far: after omitting higher-than-quadratic terms, the Hamiltonian is
$$
H = \frac{1}{2}\int dx\
\left(\pi^2+(\phi')^2 +f(x)\phi^2\right)
\tag{5b}
$$
with
$$
f(x)
=\mu^2\left(3\tanh^2\left(\frac{\mu x}{\sqrt{2}}\right)-1\right).
\tag{10b}
$$
Integrate-by-parts to put (5b) in the form
$$
H = \frac{1}{2}\int dx\
\Big(\pi^2+\big(-\phi'' +f(x)\phi\big)\phi\Big).
\tag{11}
$$
This shows that we can finish diagonalizing the Hamiltonian by solving the eigenvalue equation
$$
-\phi''+f(x)\phi=E\phi,
\tag{12}
$$
because solutions of (12) with different eigenvalues $E$ are mutually orthogonal.
Some intuition about the spectrum
The function $f(x)$ approaches a positive constant $2\mu^2$ as $x\to\pm\infty$ but dips down to $-\mu^2$ at $x=0$. Since (12) has the familiar form of a time-independent Schrödinger equation for a particle in the potential well $f(x)$, we expect that equation (12) will have a finite number of solutions with discrete eigenvalues. These eigenvalues are determined below.
Calculation of the spectrum
To save some writing, define $y\equiv \mu x / \sqrt{2}$ so that equation (12) becomes
$$
\newcommand{\pl}{\partial}
(-\pl_y^2+6\tanh^2 y-2)\phi=2(E/\mu^2)\phi.
\tag{13}
$$
To save even more writing, define $t\equiv \tanh y$. The identity
$$
-\pl_y^2+6t^2-2 = (-\pl_y+2t)(\pl_y+2t)
\tag{14}
$$
immediately gives us one eigenfunction, namely the function
$$
\phi_2(y) \propto \exp\left(-\int_0^y dz\ \tanh(z)\right),
\tag{15}
$$
which is annihilated by $\pl_y+2t$ and which approaches $0$ for $y\to\pm\infty$. The corresponding eigenvalue of (14) is $0$. Use this in (13) to get $E=0$. We can generate more eigenvalues using this identity:
$$
(-\pl_y+nt)(\pl_y+nt)+n
=
\big(\pl_y+(n+1)t\big)\big(-\pl_y+(n+1)t\big)-(n+1),
\tag{16}
$$
which holds for every integer $n$. To use this, let $\phi_n$ be the function annihilated by $\pl_y+nt$. Using (16), we can verify that
$$
(-\pl_y+2t)\phi_1
\tag{17}
$$
is an eigenfunction of (14) with eigenvalue $3$. Use this in (13) to get $E=(3/2)\mu^2$. Again using (16), we can verify that
$$
(-\pl_y+2t)(-\pl_y+t)\phi_0
\tag{18}
$$
is an eigenfunction of (14) with eigenvalue $4$. Use this in (13) to get $E=2\mu^2$. This pattern cannot be continued, because the only function annihilated by $\pl_y+0t=\pl_y$ is a constant function, which is not normalizable unless it is zero. Thus we have found exactly three discrete values for $E$, namely
$$
E = 0,\ \frac{3}{2}\mu^2,\ 2\mu^2,
\tag{19}
$$
which can be summarized as
$$
E = \frac{n(4-n)}{2}\mu^2
\tag{20}
$$
with $n=0,1,2$. For each of these values of $E$, equation (11) gives a corresponding harmonic-oscillator mode with energy $\sqrt{E}$. This reproduces the $n$-dependent term in $M_n$ shown in the question, which is the interesting part.
References
This general approach to diagonalizing the Hamiltonian is consistent with section 2 of ref 1. The approach I used to solve equation (12) with the potential (10b) was inspired by sections 2 and 3 in ref 2. The same equation is also treated in problem 39 in ref 3, starting on page 94. For more information about solvable special cases of the one-dimensional Schrödinger equation, see refs 4 and 5.
Polyakov (1976), "Isomeric states of quantum fields," Sov. Phys. JETP, 41:988-995 (http://jetp.ac.ru/cgi-bin/dn/e_041_06_0988.pdf)
Bougie et al (2012), "Supersymmetric Quantum Mechanics and Solvable Models," Symmetry 4:452-473 (https://www.mdpi.com/2073-8994/4/3/452/htm)
Flügge (1999), Practical Quantum Mechanics, Springer (https://books.google.com/books?id=VpggN9qIFUcC&pg=PA94)
Alvarez-Castillo (2007), "Exactly Solvable Potentials and Romanovski Polynomials in Quantum Mechanics" (https://arxiv.org/abs/0808.1642)
Dereziński and Wrochna (2010), "Exactly solvable Schrödinger operators" (https://arxiv.org/abs/1009.0541)
