# Canonical transformation for the "nearly diagonal" Hamiltonians

The Hamiltonians in question are the Hamiltonians that have all the non-diagonal equal to zero except for the first row and the first column (assuming for simplicity that $$U_i\in\mathbb{R}$$ are real): $$H_2 = \begin{bmatrix} \epsilon_0 & U_1 \\ U_1 & \epsilon_1 \end{bmatrix}$$ $$H_3 = \begin{bmatrix} \epsilon_0 & U_1 & U_2\\ U_1 & \epsilon_1 & 0 \\ U_2 & 0 & \epsilon_2 \end{bmatrix},$$ $$H_4 = \begin{bmatrix} \epsilon_0 & U_1 & U_2 & U_3\\ U_1 & \epsilon_1 & 0 & 0 \\ U_2 & 0 & \epsilon_2 & 0 \\ U_3 & 0 & 0 & \epsilon_3 \end{bmatrix},$$ and so on.

Relevance to physics
Such Hamiltonians are quite ubiquitous, and many people could come with the examples typical for their field. I will mention just a few:

• Coupled sublattice sites (order $$n$$ of $$H_n$$ is typically 3 or 4, see also this question)
• A level coupled to a band (order $$n$$ is very large)
• Anderson impurity (with some complications added)

Solvability
Many particular cases of this Hamiltonian are solvable: e.g., when all the couplings are the same ($$U_1=U_2=...=U_{n_1}$$) and the number of levels is infinite (broad-band limit for a resonant level) or if $$\epsilon_1=\epsilon_2=...=\epsilon_{n-1}\neq \epsilon_0$$, when the characteristic equation can be reduced to a quadratic equation.

Moreover, in cases $$n=3$$, $$n=4$$ the characteristic equation is exactly solvable in principle, although the results may not look pretty.

Question
What is the general form of a canonical transformation for diagonalizing such Hamiltonians?

For $$n=3$$ and $$n=4$$ this would amount to a simply convenient parametrization, whereas for arbitrary $$n$$ it is a transformation that takes into account the greater simplicity of this Hamiltonians in comparison to a general $$n\times n$$ Hamiltonian.

To suggest a specific line of reasoning: the canonical transformation matrix for a real Hamiltonian is an orthogonal matrix with $$n(n-1)/2$$ independent parameters (see here, e.g.). However, the zeros in the Hamiltonian impose additional $$(n-1)(n-2)/2$$ constraints, thus leaving only $$n-1$$ independent parameters.

Update
I am looking for a closed analytical expression parametrized in terms of $$n-1$$ parameters, which could be unknown. For example, this could be $$n-1$$ eigenvalues - since we are not sure that they can be found analytically for $$n>4$$ (a definitive statement to this end would be also handy). I might be already providing an answer here...

• Just for clarification: are you asking whether there is a closed analytic expression for the eigenvectors of this class of symmetric matrices? Commented Feb 19, 2021 at 16:41
• @TomášBrauner I hoped for some nice, intuitive symmetric expression, like the rotation matrix for 2-by-2 case... but yes, technically you have pointed out a solution, although not a very pretty one. Commented Feb 19, 2021 at 17:31
• Note that a solution to the problem you describe can be used recursively diagonalise any symmetric matrix. Based on our understanding of the complexity of the latter problem this constrains the complexity of finding a solution to your specific problem. Commented Feb 19, 2021 at 18:58
• @Vadim unfortunately ytlu's method cannot work. The characteristic polynomials of your matrices $H_n$ do not factorise into quadratic equations, and so cannot be transformed to $2 \times 2$ problems by a matrix whose elements are linear combinations the matrix elements of $H_n$. Commented Feb 20, 2021 at 0:15
• This should have a nice interpretation in terms of perturbation theory, since the matrix is only every coupling the 1st state to one of the others, and back. And perturbation theory can be interpreted as decoupling the two diagonal blocks to any given order. Commented Feb 22, 2021 at 16:58

The working matrix $$(n+1)\times (n+1)$$: $$\mathbf{H} = \begin{bmatrix} \epsilon_0 & U_1 & U_2 & U_3 & ... &U_{n-1}&U_n\\ U_1 & \epsilon_1 &0&0&...&0&0\\ U_2& 0 & \epsilon_2 &0&...&0&0\\ U_3& 0 & 0 & \epsilon_3 &...&0&0\\ ...&...&...&...&...&....&...\\ U_{n-1} &0&0&0&...&\epsilon_{n-1}&0\\ U_n &0&0&0&...&0&\epsilon_n \end{bmatrix}$$ Assum that all $$U_i \ne = 0$$, and are real numbers. The eigen equation for this matrix:

$$\tag{1} det\left( \mathbf{H} - \lambda \mathbf{I} \right) = 0$$ Since $$\mathbf{H}$$ is a symmetric real matrix, it guarantees that there exists $$n+1$$ real eigen values and the corresponding $$n+1$$ orthogonal eigen vectors.

A closed form of the eigen equation can be obtained by Gauss-Jordan elimination. First, add to the first row, the $$n$$th row mutiplied by $$-\frac{U_n}{\epsilon_n -\lambda}$$. This will leave $$H_{0n}=0$$, and adds a term to $$H_{00}$$, $$H_{00} = (\epsilon_0 - \lambda) -\frac{U_n^2}{\epsilon_n -\lambda}$$. Repeat this process for $$i = n-1, n-2, ... 1$$. After $$n$$ eliminations, the upper triangular elements of $$\mathbf{H} - \lambda \mathbf{I}$$ becomes all zeros, and the element $$H_{00}$$ collects $$n$$ more terms from the elimination process:

$$\tag{2} H_{00} = (\epsilon_0 - \lambda) - \sum_{i=1}^n \frac{U^2_i}{\epsilon_i - \lambda}.$$

and the determinant is the product of all diagonal elements:

$$det\left(\mathbf{H} - \lambda \mathbf{I} \right) = \left\{ (\epsilon_0 - \lambda) - \sum_{i=1}^n \frac{U^2_i}{\epsilon_i - \lambda} \right\} \prod_{i=1}^n (\epsilon_i - \lambda).$$ $$\tag{3} = \prod_{i=0}^n (\epsilon_i - \lambda) - \sum_{i=1}^n U^2_i \prod_{j=1, j\ne i}^n (\epsilon_j - \lambda).$$

Case 1. Degenerate diaginal elements

In case that there are some energies are equal, saying 3-fold degeneracy: $$\epsilon_a = \epsilon_b = \epsilon_c$$ for $$a$$, $$b$$, $$c$$ are arbitrary three integers in $$\{1,2,3, ..., n\}$$.

Under this case, a term $$(\epsilon_a - \lambda )^2$$ can be factored out of the Eq. (3). The energy $$\epsilon_a$$ is therefore an eigen value of degeneracy $$3-1=2$$. The corresponding eigen vectors can be calculated from the $$4\times 4$$ submatrix:

$$\tag{4} \begin{bmatrix} \epsilon_0-\epsilon_a & U_a & U_b & U_c\\ U_a & 0 & 0 & 0\\ U_b & 0 & 0 & 0\\ U_c & 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} x_0\\ x_a\\ x_b\\ x_c \end{bmatrix} = \begin{bmatrix} 0\\ 0\\ 0\\ 0 \end{bmatrix}$$

The equation for eigen vectors: $$x_0 = 0;$$ $$x_a U_a + x_b U_b + x_c U_c = 0.$$

Two orthonomal eigen vectors can be constructed from these two conditions. It would be very interesting to observe how these two domension block-diaginalized and decoupled from the rest $$(n-1)\times (n-1)$$ dimensions, and how the coupling elements is changed for the remaining state of the 3 degenerate levels.

We conclude that for $$p$$-fold degeneracy energy $$\epsilon_a$$, ther existed $$(p-1)$$ degenerate eigen vectors with eigen value $$\lambda = \epsilon_a$$. This $$(p-1)$$ basis vectors span an $$(p-1)$$ subspace. The motion of the subspace is decoupled from the rest of the system. One example realized this matrix is a coupled oscillator system, which all the oscillators are coupled only to the $$0$$th oscillator (see Figure). The dynamical matrix for the normal mode will resemble this matrix form. Then, we can select a subgroup of oscillation modes, that they are coupled only within the subgroup.

Another worthy mentioned situation is that all $$\epsilon_i = \epsilon_1$$ for all $$i=1, 2, 3,..,n$$. Then we immediately has $$n-1$$-fold eigen values $$\lambda = \epsilon_1$$. Their eigen vectors span the $$(n-1)\times (n-1)$$ subspace.

$$\vec{v}_\lambda = [x_0, x_1, x_2, ..., x_{n-1}, x_n]$$

Conditions for eigen vector:

$$x_0 = 0;$$ $$\sum_{i=1}^n x_i U_i = 0.$$

Since we have $$(n+1)\times (n+1)$$ symmetric matrix, there are two more eigen values. From the determinant, after factoring out $$(\epsilon_1-\lambda)^{n-1}$$, we find that the last two eigen values are the root of the equation: $$(\epsilon_0-\lambda) (\epsilon_1-\lambda) - \sum_{i=1}^n U_n^2 = 0.$$ Which is two-level matrix with coupling strength $$U =\sqrt{ \sum_{i=1}^n U_n^2}$$. The eigen values are readily given: $$\lambda_\pm = \frac{\epsilon_0+\epsilon_1}{2} \pm \sqrt{ \frac{(\epsilon_0-\epsilon_1)^2}{2} + U^2}$$ The total $$n+1$$ eigen values consist of $$(n-1)$$-fold of $$\epsilon_1$$ in the middle and two $$\lambda_\pm$$ in higher and lower values.

Non-degenrate case

we may now assume that all $$\epsilon_i$$ are different. All eigen values will then solely determined by zeros of the element $$H_{00}$$ in Eq.(2)': {2} $$H_{00} = (\epsilon_0 - \lambda) - \sum_{i=1}^n \frac{U^2_i}{\epsilon_i - \lambda} = 0.$$ $$= (\epsilon_0 - \lambda) - \sum_{i=1}^n \frac{U^2_i}{(\epsilon_i-\epsilon_0 ) -( \lambda -\epsilon_0 ) } = 0.$$

Setting $$\epsilon_0 = 0$$ without lossing generality: $$= -\lambda - \sum_{i=1}^n \frac{U^2_i}{\epsilon_i - \lambda } = 0.$$

For better understanding of this equation, lets assum $$\epsilon_i = i$$. We can now examine the value of $$H_{00}(\lambda)$$ by scanning $$\lambda$$ from $$-\infty$$ to $$\infty$$.

The value of $$H_{00}(\lambda)$$ is positive when $$\lambda \ll 0$$, but the term of summation is negative. Just before $$\lambda \to 0^-$$, the smaller negative summation catches up with $$-\lambda$$. It is the smallest eigen value (the $$\lambda_-$$ in the total degenerate case.)

Then the $$H_{00}$$ approach $$-\infty$$ at $$\epsilon_1^-$$ due to divergence in $$-\frac{U_1^2}{\epsilon_1 - \lambda}$$. After $$\lambda = \epsilon_1^+$$, the function emerges with large positive value. And approach $$-\infty$$ again at $$\epsilon_2^-$$, behaves like the $$\cot$$ function. Therefore, it locates first eigen value at $$\lambda_0 <0$$, and $$n-1$$ eigen values at $$\epsilon_i < \lambda_i < \epsilon_{i+1}$$, a final eigen value locates at $$\lambda_n > \epsilon_n$$, because when $$\lambda \to \infty$$, $$H_{00} \to -\infty$$, as shown in Fig.(2).

We can see this figure topologically, as $$\epsilon_1$$,..$$\epsilon_n$$ are drawn closer, and finally into one energy $$\epsilon_1$$. All the levels in between, are forced collapse into a single degenerate level, but the first and the last, stay away from this squeezing process.

pertubation for small $$U_i$$'s

In case that $$U_i \ll |\epsilon_a - \epsilon_b|$$ for any two diagonal element, $$\epsilon_a$$, and $$\epsilon_b$$. The eigen value and eigen vector can be obtained by a approximation method similar to the second order perturbation in quantum mechanics:

$$\lambda_0 = \sum_{i=0}^n \frac{U_i^2}{\epsilon_i}.$$ Where is replace the $$\lambda$$ in the summation by $$\epsilon_0 = 0$$.

For $$i \ge 1$$ $$\lambda_i = \epsilon_i + \frac{U_i^2}{\epsilon_i - \epsilon_0} = \epsilon_i + \frac{U_i^2}{\epsilon_i}.$$ Where we replace the $$\lambda$$ in the summation by $$\epsilon_i$$, and made expansion for small $$U_i$$s.

Eigen vectors and tranformation

Each eigen value renders an eigen vector $$\vec{v}_\lambda = \left[ x_0, x_1, x_2, ...,x_n \right]$$: $$\vec{v}_\lambda = A_1 \begin{bmatrix} 1, \frac{U_1}{\lambda - \epsilon_1}, \frac{U_2}{\lambda - \epsilon_2}, ..., \frac{U_n}{\lambda - \epsilon_n} \end{bmatrix}$$ Leave $$A_1$$ as the normalization constant. There are $$n$$ ortho-normal vectors, providing the bases for the orthogonal transformation $$\mathbf{R}$$ to the diagonalisation of the Hamiltonian.

$$\mathbf{R} = \begin{bmatrix} A_0 & A_1 & A_2 & A_3& ... & A_n\\ \frac{A_0 U_1}{\lambda_0 - \epsilon_1} &\frac{A_1 U_1}{\lambda_1 - \epsilon_1} & \frac{A_2 U_2}{\lambda_2 - \epsilon_2} &\frac{A_3 U_1}{\lambda_3 - \epsilon_1} & ... & \frac{A_{n} U_1}{\lambda_{n} - \epsilon_1}\\ \frac{A_0 U_2}{\lambda_0 - \epsilon_1} &\frac{A_1 U_2}{\lambda_1 - \epsilon_2} & \frac{A_2 U_2}{\lambda_2 - \epsilon_1} &\frac{A_3 U_2}{\lambda_3 - \epsilon_2} & ... & \frac{A_{n} U_2}{\lambda_{n} - \epsilon_2}\\ ....&... & ... & ... & ... \\ \frac{A_0 U_{n}}{\lambda_0 - \epsilon_1} & \frac{A_1 U_{n}}{\lambda_1 - \epsilon_{n}} & \frac{A_2 U_{n}}{\lambda_2 - \epsilon_{n}} & \frac{A_3 U_{n}}{\lambda_3 - \epsilon_{n}} &... & \frac{A_{n} U_{n}}{\lambda_{n} - \epsilon_{n}} \end{bmatrix}.$$

Since the $$\mathbf{H}$$ is Hermitian, the $$n+1$$ eigen vectors are mutually orthogonal:

$$\mathbf{R} \mathbf{R}^T =\mathbf{R}^T \mathbf{R} = \mathbf{I}$$ $$\mathbf{R}^T \mathbf{H}\mathbf{R} =\mathbf{D}$$ Where $$\mathbf{D}$$ is the diagonal matrix with all $$\lambda_i$$ as the elements.

Lets consider a matrix equation: $$\mathbf{H} \vec{v} = \vec{w}$$

Transform it to the eigen coordinate: $$\mathbf{R}^T \mathbf{H} \left(\mathbf{R} \mathbf{R}^T \right) \vec{v} = \mathbf{R}^T\vec{w}$$ $$\mathbf{D} \left( \mathbf{R}^T \vec{v} \right) = \left(\mathbf{R}^T\vec{w}\right)$$ $$\mathbf{D} \vec{v}' = \vec{w}'$$

In $$\vec{v}'$$ and $$\vec{w}'$$ primed coordinate, the matrix $$\mathbf{D}$$ is diagonalized, thus each components of $$\vec{v}'$$ and $$\vec{w}'$$ is decoupled.

The transformation $$\mathbf{R}^T \vec{w} = \vec{w}' ;$$ and $$\mathbf{R} \vec{w}' = \vec{w} ;$$

• Case 1 is just the general case - in terms of canonical transformation one still has $n$ unknown eigenvalues and $n-1$ normalization constants the canonical transformation. It is not clear what the case 2 is differznt ftom case 3. In case 3 - what are $H_1$ and $H_2$? It is not clear to me, if you are really onto something new... Commented Feb 21, 2021 at 7:14
• There are $n$ roots from the Eq. (1), and reders $n$ eigen vectors.
– ytlu
Commented Feb 21, 2021 at 11:27
• @Wolpertinger equation 1 is the standard eigenvalue equation $|H-\lamba I|=0$. It is easily written for $n=2,3,4$ and generalized to higher order by induction. Commented Feb 21, 2021 at 13:15
• @Wolpertinger what you are saying is correct for an arbitrary n-by-n Hamiltonian, it doesn't bring us closer to understanding the Hamiltonians in question. General polynomial equation of order $n$ is not solvable, but there are many particular cases which are solvable. E.g., the degenerate case $\epsilon_1=\epsilon_2=...\epsilon_{n-1}$ is exactly solvable for arbitrary $n$, $\epsilon_0$, and $U_1, U_2,...,U_{n-1}$ - as pointed out in the question. Commented Feb 21, 2021 at 17:53
• @ytlu Getting Eq. (1) by Laplace expansion is trivial. Also, if you insist on LU decomposition - once your matrix has upper or lower triangular elements zero, the determinant is simply the product of the diagonal elements. Commented Feb 22, 2021 at 11:12

this is not a complete solution just some idea

$$\mathbf A=\left[ \begin {array}{ccc} \epsilon _{{1}}&U_{{1}}&U_{{2}} \\ U_{{1}}&\epsilon _{{2}}&0\\ U_{ {2}}&0&\epsilon _{{3}}\end {array} \right]$$

Step I

Transform $$~\mathbf A~$$ to

$$\mathbf B=\mathbf J^T\,\mathbf A\,\mathbf J= \left[ \begin {array}{ccc} 1&{\frac {U_{{1}}}{\sqrt {\epsilon _{{1}}} \sqrt {\epsilon _{{2}}}}}&{\frac {U_{{2}}}{\sqrt {\epsilon _{{1}}} \sqrt {\epsilon _{{3}}}}}\\ {\frac {U_{{1}}}{\sqrt { \epsilon _{{1}}}\sqrt {\epsilon _{{2}}}}}&1&0\\ { \frac {U_{{2}}}{\sqrt {\epsilon _{{1}}}\sqrt {\epsilon _{{3}}}}}&0&1 \end {array} \right] \mapsto \left[ \begin {array}{ccc} 1&U_{{1}}&U_{{2}}\\ U_{{ 1}}&1&0\\ U_{{2}}&0&1\end {array} \right]$$

with: $$\mathbf J=\left[ \begin {array}{ccc} {\frac {1}{\sqrt {\epsilon _{{1}}}}}&0&0 \\ 0&{\frac {1}{\sqrt {\epsilon _{{2}}}}}&0 \\ 0&0&{\frac {1}{\sqrt {\epsilon _{{3}}}}} \end {array} \right]$$

Step II

transform $$\mathbf B$$ to diagonal shape. Ansatz: $$\mathbf B_d=~\mathbf T^T\,\mathbf B\,\mathbf T-I_3=\mathbf 0$$

with :

$$\mathbf T=\left[ \begin {array}{ccc} 1&1&1\\ {\it T21}&1&1 \\ {\it T31}&{\it T32}&1\end {array} \right]$$

you get three equations

$$eq_1=1+U_{{1}}+{\it T32}\,U_{{2}}+{\it T21}\,U_{{1}}+{\it T21}+{\it T31}\,U _{{2}}+{\it T31}\,{\it T32} =0$$ $$eq_2=1+{\it T21}\,U_{{1}}+{\it T31}\,U_{{2}}+U_{{1}}+{\it T21}+U_{{2}}+{ \it T31} =0$$ $$eq_3=2+2\,U_{{1}}+{\it T32}\,U_{{2}}+U_{{2}}+{\it T32}=0$$

for the three unknows $$~T21~,T31~,T32$$

Results:

$$\mathbf T=\left[ \begin {array}{ccc} 1&1&1\\ -{\frac {U_{{1}} +1-{U_{{2}}}^{2}}{U_{{1}}+1}}&1&1\\ -U_{{2}}&-{ \frac {2+2\,U_{{1}}+U_{{2}}}{U_{{2}}+1}}&1\end {array} \right]$$

$$B_d(1,1)=-{\frac {2\,{U_{{1}}}^{3}+ \left( 2-{U_{{2}}}^{2} \right) {U_{{1}}}^{2 }+ \left( -2+2\,{U_{{2}}}^{2} \right) U_{{1}}-2-{U_{{2}}}^{4}+3\,{U_{{ 2}}}^{2}}{ \left( U_{{1}}+1 \right) ^{2}}}$$ $$B_d(2,2)=4\,{\frac {{U_{{1}}}^{2}}{ \left( U_{{2}}+1 \right) ^{2}}}+{\frac { \left( -2\,{U_{{2}}}^{2}+4\,U_{{2}}+10 \right) U_{{1}}}{ \left( U_{{2 }}+1 \right) ^{2}}}+{\frac {-3\,{U_{{2}}}^{2}+4\,U_{{2}}+6-2\,{U_{{2}} }^{3}}{ \left( U_{{2}}+1 \right) ^{2}}}$$ $$B_d(3,3)=3+2\,U_{{1}}+2\,U_{{2}}$$

and $$~B_d(i,j)=0~,i\ne j$$

for

$$\mathbf B=\left[ \begin {array}{cccc} 1&U_{{1}}&U_{{2}}&U_{{3}} \\U_{{1}}&1&0&0\\ U_{{2}}&0&1&0 \\ U_{{3}}&0&0&1\end {array} \right]$$

$$\mathbf T=\left[ \begin {array}{cccc} 1&1&1&1\\ {\frac {U_{{ 1}}+1-{U_{{3}}}^{2}-{U_{{2}}}^{2}}{U_{{1}}+1}}&1&1&1 \\ -U_{{2}}&-{\frac {2+2\,U_{{1}}+U_{{2}}-{U_{{3}}}^ {2}}{U_{{2}}+1}}&1&1\\ -U_{{3}}&-U_{{3}}&-{\frac {3+ 2\,U_{{1}}+2\,U_{{2}}+U_{{3}}}{U_{{3}}+1}}&1\end {array} \right]$$

you can try now to find a pattern between the transformation matrix $$~\mathbf T_{n=3}~$$ and $$~\mathbf T_{n=4}~$$ if you find a pattern you can obtain the transformation matrix $$~\mathbf T_{n}~$$

• I've been thinking about this approach as well. In principle, once you have a matrix with equal diagonal elements, you csn find all the eigenvalues (see my own answer). The problem is that your transformation $J$ is not a unitary/orthogonal: $J^T\neq J^{-1}$... Commented Feb 22, 2021 at 17:25
• I've been thinking about this approach as well. In principle, once you have a matrix with equal diagonal elements, you csn find all the eigenvalues (see my own answer). The problem is that your transformation $J$ is not a unitary/orthogonal: $J^T\neq J^{-1}$... Commented Feb 22, 2021 at 17:25

Disclaimer: Due to the disagreement (in the comments) about what should be considered an answer to my question I post this here: it refines the existing answer and raises more precise questions about what needs to be understood for a complete treatment of this problem.

Standard approach to eigenvalue problems includes

1. Solving the characteristic equation for the eigenvalues, $$|H-\lambda I|=0,$$ where $$I$$ is the identity matrix.
2. Solving the eigenvectors equations, given known eigenvalues, $$H\mathbf{x}=\lambda \mathbf{x}.$$

The characteristic equation
Evaluating determinants is is trivial for this type of the Hamiltonians. In particular we obtain:
$$n=2$$
$$|H-\lambda I| = \left|\begin{matrix} \epsilon_0-\lambda & U_1 \\ U_1 & \epsilon_1 - \lambda \end{matrix}\right| =(\epsilon_0 - \lambda)(\epsilon_1-\lambda) - U_1^2 = 0$$ $$n=3$$ $$|H-\lambda I| = (\epsilon_0 - \lambda)(\epsilon_1-\lambda)(\epsilon_2-\lambda) - U_1^2(\epsilon_2-\lambda) - U_2^2(\epsilon_1-\lambda) = 0$$ $$n=3$$ $$|H-\lambda I| = (\epsilon_0 - \lambda)(\epsilon_1-\lambda)(\epsilon_2-\lambda)(\epsilon_3-\lambda) - U_1^2(\epsilon_2-\lambda)(\epsilon_3-\lambda) - U_2^2(\epsilon_1-\lambda)(\epsilon_3-\lambda) -U_3^2(\epsilon_1-\lambda)(\epsilon_2-\lambda) = 0$$ One can then show by induction that for an arbitrary order we get: $$|H-\lambda I| = (\epsilon_0 - \lambda)\prod_{i=1}^{n-1} (\epsilon_i - \lambda) - \sum_{i=1}^{n-1}U_i^2 \prod_{j=1, j\neq i}^{n-1} (\epsilon_i - \lambda) = \\ \left(\epsilon_0 - \lambda - \sum_{i=1}^{n-1}\frac{U_i^2}{\epsilon_i-\lambda}\right)\prod_{i=1}^{n-1} (\epsilon_i - \lambda)=0$$ The last form of the determinant obviously makes sense only if the roots are not equal to any of $$\epsilon_i, i>0$$. It is also easily seen that, if any two energies are equal $$\epsilon_i=\epsilon_j$$, one can factor out $$\epsilon_i-\lambda$$, obtaining one root of the equation.

Q1: The non-trivial question on the level of the characteristic equation is: can it be solved or can it be shown that it cannot be solved for $$n>0$$. Note that while the latter statement is known to be true for a general polynomial equation, it is not clear a priori that this is also the case for this particular form.

Solving for eigenvectors
Assuming that the eigenvalues $$\lambda_i, i=1,..,n$$ are known, we can solve for the eigenvectors (Note that this is again a standard procedure for an arbitrary matrix Hamiltonian). For simplicity we assume that none of the roots is equal to any of the energies $$\epsilon_j$$, although we will show that the solution can be represented in a general form, suitable for a degenerate case.

The solution for this types of the Hamiltonians is trivial (formally one could call it a trivial case of Gaussian elimination). The eigenvector corresponding to eigenvalue $$\lambda_i$$ takes form: $$\mathbf{x}_i = A_i\begin{bmatrix} 1\\\frac{U_1}{\lambda_i - \epsilon_1}\\\frac{U_2}{\lambda_i - \epsilon_1}\\...\\\frac{U_{n-1}}{\lambda_i - \epsilon_{n-1}}\end{bmatrix},$$ where $$A_i = \left[1 + \sum_{j=1}^{n-1}\frac{U_j^2}{(\lambda_i - \epsilon_j)^2}\right]^{-\frac{1}{2}}.$$ Note that by replacing $$1$$ in the first element in the eigenvector with $$\prod_{j=1}^{n-1}(\lambda_i-\epsilon_j)$$ we could remove singularities from thz denominators, obtaining the answers that are also valid for a non-degenerate case.

The canonical transformation diagonalizing the Hamiltonian ($$S^\dagger H S = \Lambda$$, where $$\Lambda_{ij} = \delta_{i,j}\lambda_i$$) can be written as $$\mathbf{S} = \begin{bmatrix} A_1 & A_2 & ... & A_3\\ \frac{A_1 U_1}{\lambda_1 - \epsilon_1} & \frac{A_2 U_1}{\lambda_2 - \epsilon_1} & ... & \frac{A_{n} U_1}{\lambda_{n} - \epsilon_1}\\ \frac{A_1 U_2}{\lambda_1 - \epsilon_2} & \frac{A_2 U_2}{\lambda_2 - \epsilon_2} & ... & \frac{A_{n} U_2}{\lambda_{n} - \epsilon_2}\\ ... & ... & ... & ... \\ \frac{A_1 U_{n-1}}{\lambda_1 - \epsilon_{n-1}} & \frac{A_2 U_{n-1}}{\lambda_2 - \epsilon_{n-1}} & ... & \frac{A_{n} U_{n-1}}{\lambda_{n} - \epsilon_{n-1}} \end{bmatrix}.$$ This canonical transformation is parametrized by $$n$$ eigenvalues $$\lambda_i$$, which can be considered as parameters: even if we fail to solve the characteristic equations the canonical transformation will have this form.

Q2: Contrary to what was suggested in the OP, the transformation diagonalizing an arbitrary n-by-n Hamiltonian is parametrizable in terms of its $$n$$ eigenvalues. So the question is about the special constrains and the number of independent parameters that we have in case of the simpler Hamiltonians considered in this question. Given the presence of the additional constrains - shouldn't there exist a relation between the eigenvalues, meaning that some cases with $$n>4$$ are actually exactly solvable?

Simpler form of the canonical transformation
The canonical transformation above can be written as $$\mathbf{S} = \begin{bmatrix} 1 & 0 & ... & 0 \\ 0 & \frac{U_1}{U} & ... & 0 \\ ... & ... & ... & ... \\ 0 & 0 & ... & \frac{U_{n-1}}{U} \end{bmatrix} \begin{bmatrix} A & A & ... & A\\ \frac{A U}{\lambda_1 - \epsilon_1} & \frac{A U}{\lambda_2 - \epsilon_1} & ... & \frac{A U}{\lambda_{n} - \epsilon_1}\\ \frac{A U}{\lambda_1 - \epsilon_2} & \frac{A U}{\lambda_2 - \epsilon_2} & ... & \frac{A U}{\lambda_{n} - \epsilon_2}\\ ... & ... & ... & ... \\ \frac{A U}{\lambda_1 - \epsilon_{n-1}} & \frac{A U}{\lambda_2 - \epsilon_{n-1}} & ... & \frac{A U}{\lambda_{n} - \epsilon_{n-1}} \end{bmatrix} \begin{bmatrix} \frac{A_1}{A} & 0 & ... & 0 \\ 0 & \frac{A_2}{A} & ... & 0 \\ ... & ... & ... & ... \\ 0 & 0 & ... & \frac{A_n}{A} \end{bmatrix}$$ This decomposition is more obvious, if we take $$U=A=1$$, but defining these factors as $$U=\left(\prod_{i=1}^{n-1}U_i\right)^{\frac{1}{n-1}},\\ A=\left(\prod_{i=1}^{n-1}A_i\right)^{\frac{1}{n}}$$ preserves the dimensionality and assures that we are dealing with a product of three canonical transformations (three unitary/orthogonal matrices).

Q3: The ultimate holy grail here would be to have a representation of the canonical transformation that is manifestly unitary and "symmetric" in the way the well-known rotation matrices are for $$n=2$$ and $$n=3$$. At present these properties are hidden in the eigenvalues $$\lambda_i$$.

Other possible approaches

• Canonical transformation with an ad-hoc parametrization of the transformation matrix (e.g., using Cayley transform)
• Formal perturabtion expansion to all orders of magnitude, e.g., using the resolvent operator
• Using an auxhiliary transformation, as suggested by @Eli
• From eigen vectors to the answer of Q3 is a process of elementary physical math. It is orthogonal (unitary for complex matrix), but not symmetric. If your symmetry means dual transformation, then again, it is a buit-in property of an orthogonal matrix. All these features have been granted in the form of the Hamiltonian - a hermitian matrix.
– ytlu
Commented Feb 22, 2021 at 11:14
• @Wolpertinger I appreciate your efforts, but think you at_ytlu simply do not understand the problem - probably because you lack experience with relevant physical problems. Please let others contribute. Commented Feb 22, 2021 at 15:19
• @Wolpertinger Should I consider it as a threat or revenge? Moderators would normally point out that you do not close the question that you answer and vice versa. You in fact not only tried to answer it, but also tried to promote another answer. Commented Feb 22, 2021 at 18:01