Is there a similarity transformation rendering all diagonal elements of a matrix equal? I'm especially interested in SL$(2,\mathbb C)$, i.e. $2\times2$ matrices with determinant one, in which case I'm looking for a transformation from $\begin{pmatrix}a&b\\c&d\end{pmatrix}$ to $\begin{pmatrix}\frac{a+d}2&x\\ y&\frac{a+d}2\end{pmatrix}$ (the trace is conserved). Does such a similarity transformation exist? What about general $n\times n$ matrices?
Bonus points for an analytical formula (even if only for the 2x2 case).
 A: I will present the final algorithm first, with follow up descriptions. Use the matrix
$$\pmatrix{d_0 & x \\ y & d_1}$$
To obtain the complex $c$ and $s$ values, follow these steps (the text hopefully becomes clear with the later descriptions):
\begin{align}
  y_s &= y - \bar{x} \quad&\text{ from the skew part}\\
  \delta &= d_1 - d_0 \\
  \delta_s &= \operatorname{Imaginary}\{\delta \} \\
  \Delta &= \delta_s^2 + y_s\bar{y}_s &  \text{the discriminant in the formulation for the sine} \\
  s_0 &= j(\delta_s + \sqrt{\Delta}) & \text{the scaled sine value to diagonalize the skew part} \\
  x_0 &= \delta y_s s_0 + y_s^2 x - s_0^2 y & \text{the first rotation's (scaled) result}\\
  k_2c_2  &= j|x_0|y_s + x_0 \bar{s}_0 & \text{ the final (scaled) rotation values} \\
  k_2s_2 &= j|x_0|s_0 - x_0\bar{y}_s \\
\end{align}
$y_s$ and $\delta_s$ are all the information needed from the skew Hermitian portion of the matrix, needed to diagonalize the
skew Hermitian component to give a result of zero (for the new $x$ and $y$) in the skew Hermitian part of the result. This means that $x=\bar{y}$ in
the result (first step).
$y_s$ is not divided by $2$ as it would normally be in the calculation of the skew portion, since it is an unnecessary scale factor that
is removed in the final scaling for $c$ and $s$.
Solving a particular quadratic (see this previous question) gives $k_0c_0=y_s$ and $k_0s_0 = j(\delta_s + \sqrt{\Delta})$, as described. The $k_0$ need never
be found; it is a consistent scale factor that is removed in the final result (thus not included as a variable).

After $c_0$ and $s_0$ are applied to the matrix, the only necessary information from that is the $x$ value, named $x_0$ here. Thus it is the only
calculated intermediate similarity value. From it, the complex
phase is required, and $c_2$, $s_2$ are the final result. Find $k_2$ to scale these for unitary action, and the result gives equal diagonal values. See the
python code at the end for more mundane details regarding cases where $c,s=0,0$ occur.

The algorithm tests well and works, here is a hopefully enlightening description.

First for reference, the unitary (complex Givens) rotation on a 2x2 matrix gives:
\begin{align}
&
   \pmatrix{c                & s \\
            -\bar{s}         & \bar{c} \\
   }
   \pmatrix{d_0              & x \\
            y                & d_1 \\
   }
   \pmatrix{\bar{c}          & -s \\
            \bar{s}         & c \\
   }
\\
=&
   \pmatrix{c                & s \\
            -\bar{s}         & \bar{c} \\
   }
   \pmatrix{d_0\bar{c}+x\bar{s}    & -d_0 s + x c \\
            y \bar{c} + d_1 \bar{s}              & -y s + d_1 c \\
   }
\\
=&
   \pmatrix{d_0c\bar{c} + d_1s\bar{s} + c\bar{s}x + \bar{c}sy            & (d_1 - d_0)cs + x c^2 - ys^2 \\
            (d_1 -d_0)\bar{c}\bar{s} -x\bar{s}^2 + y\bar{c}^2   & d_0s\bar{s} + d_1c\bar{c} - c\bar{s}x - \bar{c}sy
   }
\end{align}
To set $d_1 = d_0 $ we solve
$$d_0c\bar{c} + d_1s\bar{s} + c\bar{s}x + \bar{c}sy = d_0s\bar{s} + d_1c\bar{c} - c\bar{s}x - \bar{c}sy $$
or
$$\overbrace{(d_1 - d_0)}^{\delta}(s\bar{s} - c\bar{c}) + c\bar{s}(2x) +\bar{c}s(2y) = 0 \tag{1}$$
Without loss of generality, use real $c\in [0,1]$, and the fact $c\bar{c} + s\bar{s} =1$ (from the unitary form of the similarity) to have
$$s\bar{s} - c\bar{c} = (1 - c\bar{c}) - c\bar{c} = 1 - 2c\bar{c}$$
and
$$s = (1 - c^2)^{\frac{1}{2}}e^{j\beta}$$
Then (1) becomes
$$\delta(1-2c^2) + 2c(1 - c^2)^{\frac{1}{2}}\left[e^{-j\beta}x +e^{j\beta}y\right]=0 \tag{2}$$
From here it seems best to perform an
intermediate similarity to achieve $|x| = |y|$. When that is true, the ellipse term $e^{-j\beta}x +e^{j\beta}y$ is easier
to deal with. To achieve it, diagonalize the skew portion of the matrix, giving $y = \bar{x}$. This then causes the ellipse to "collapse" to a line on the real axis.
The phase
angle $\measuredangle xe^{j\beta} = \pm\frac{\pi}{2}$ will give the ellipse term as zero.
So with $\beta = \frac{\pi}{2} + \measuredangle x$ we have the phase of $s$. For the magnitude, use the now reduced
 equation (2)
$$\delta(1-2c^2)=0 \tag{3}$$
We see here that $c=\frac{1}{\sqrt{2}}$ solves it, or using the scaled forms, $|c| = |s|$:
\begin{align}
  c &= |x| \\
  s &= j|x|e^{j\beta}  = x\\
\end{align}

Here is the complete python code:


def hollow(d0, x, y, d1):
  ''' return c,s for a 2x2 unitary (complex Givens) with result of diagonals equal '''
  ys=y - x.conjugate()
  if ys==0: # if the matrix's skew portion is already diagonal
    c = complex(0, abs(x))
    s = x # s has the same phase as x and |c| = |s|
  else:
    d = d1 - d0
    ds =  d.imag
    ys2 = (ys*ys.conjugate()).real
    rad = ds*ds + ys2
    s0 = complex(0,ds + math.sqrt(rad)) # can do plus or minus square root here
    x0 = d*ys*s0 + ys*ys*x - s0*s0*y # the x after diagonalizing the skew portion
    if x0==0: c,s = ys, s0
    else:
      ax0 = abs(x0)
      s1 = complex(0,abs(x0))
      c = ys*s1 + x0*s0.conjugate()
      s = s0*s1 - x0*ys.conjugate()
  n = abs(complex(abs(s),abs(c))) # the square root of the sum of the norm squared of c and s
  c=c/n
  s=s/n
  return c,s



The function is called hollow since if the trace is zero, the result is a zero diagonal matrix, also known as a hollow matrix.
A: $ \left( \begin{array}{ccc}
p & q \\
r & s 
\end{array} \right) $
$ \left( \begin{array}{ccc}
a & b \\
c & d 
\end{array} \right) $
$ =\left( \begin{array}{ccc}
\frac{a+d}{2} & x \\
y & \frac{a+d}{2} 
\end{array} \right) $
$2pa+2qc=a+d$
$pb+qd=x$
$ra+sc=y$
$2rb+2sd=a+d$
$2bra+2sda=a^{2}+da$
$2bra+2bsc=2by$
$s=\frac{2by-a^{2}-da}{2bc-2da}$
$r=\frac{by-scb}{ab}$
$p=\frac{2ax-b^{2}-cb}{2ad-2bc}$
$q=\frac{x-pb}{d}$
