[proofplan]
We first extract from equality of transfer functions the equality of all Markov parameters $C_1A_1^kB_1=C_2A_2^kB_2$. Using these identities, we define a [linear map](/page/Linear%20Map) from the reachable state space of the first realization to the reachable state space of the second by matching finite input-generated sums. Observability proves that this assignment is independent of the chosen representation, reachability proves that it is onto, and observability of the first system proves that it is one-to-one. The resulting bijection intertwines the state, input, and output matrices; its inverse is the required similarity matrix.
[/proofplan]
[step:Extract equality of the Markov parameters from equality of transfer functions]
Let $n_i$ denote the state dimension of the $i$-th realization, so $A_i \in \mathbb{R}^{n_i \times n_i}$, $B_i \in \mathbb{R}^{n_i \times m}$, $C_i \in \mathbb{R}^{p \times n_i}$, and $D \in \mathbb{R}^{p \times m}$. For a matrix $M \in \mathbb{C}^{r \times r}$, let $\|M\|_{\mathrm{op}}$ denote the operator norm of the linear map $M:\mathbb{C}^r \to \mathbb{C}^r$ with respect to the Euclidean norm on $\mathbb{C}^r$.
For each $i \in \{1,2\}$, define the transfer function
\begin{align*}
G_i:\rho(A_i) \to \mathbb{C}^{p \times m}
\end{align*}
by
\begin{align*}
G_i(z)=C_i(zI_{n_i}-A_i)^{-1}B_i+D,
\end{align*}
where $\rho(A_i)=\{z \in \mathbb{C}: zI_{n_i}-A_i \text{ is invertible}\}$ is the resolvent set of $A_i$.
Choose $R>0$ such that $|z|>R$ implies $z \in \rho(A_1)\cap \rho(A_2)$ and
\begin{align*}
\left|\frac{1}{z}\right|\|A_i\|_{\mathrm{op}}<1
\end{align*}
for both $i=1,2$. For such $z$, the Neumann series gives
\begin{align*}
(zI_{n_i}-A_i)^{-1}
= z^{-1}\left(I_{n_i}-z^{-1}A_i\right)^{-1}
= \sum_{k=0}^{\infty} z^{-k-1}A_i^k.
\end{align*}
Hence
\begin{align*}
G_i(z)-D=\sum_{k=0}^{\infty} z^{-k-1}C_iA_i^kB_i.
\end{align*}
Since $G_1(z)=G_2(z)$ for all sufficiently large $|z|$, each scalar entry of the matrix-valued Laurent expansion at infinity agrees with the corresponding scalar entry of the other expansion. Uniqueness of Laurent coefficients applied entrywise gives
\begin{align*}
C_1A_1^kB_1=C_2A_2^kB_2
\end{align*}
for every integer $k\ge 0$.
[guided]
The transfer function contains the Markov parameters as the coefficients of its expansion at infinity. We now make that statement precise.
Let $n_i$ be the state dimension of the $i$-th realization. Thus $A_i \in \mathbb{R}^{n_i \times n_i}$, $B_i \in \mathbb{R}^{n_i \times m}$, $C_i \in \mathbb{R}^{p \times n_i}$, and $D \in \mathbb{R}^{p \times m}$. Whenever $M \in \mathbb{C}^{r \times r}$ appears inside $\|M\|_{\mathrm{op}}$, this denotes the operator norm of the linear map $M:\mathbb{C}^r \to \mathbb{C}^r$ with respect to the Euclidean norm on $\mathbb{C}^r$.
For $i \in \{1,2\}$, define
\begin{align*}
G_i:\rho(A_i) \to \mathbb{C}^{p \times m}
\end{align*}
by
\begin{align*}
G_i(z)=C_i(zI_{n_i}-A_i)^{-1}B_i+D.
\end{align*}
Here $\rho(A_i)$ denotes the set of complex numbers $z$ for which the matrix $zI_{n_i}-A_i$ is invertible.
Choose $R>0$ large enough that, whenever $|z|>R$, both matrices $zI_{n_1}-A_1$ and $zI_{n_2}-A_2$ are invertible and
\begin{align*}
\left|\frac{1}{z}\right|\|A_i\|_{\mathrm{op}}<1
\end{align*}
for $i=1,2$. This allows us to expand the resolvent by the Neumann series:
\begin{align*}
(zI_{n_i}-A_i)^{-1}
= z^{-1}\left(I_{n_i}-z^{-1}A_i\right)^{-1}
= \sum_{k=0}^{\infty} z^{-k-1}A_i^k.
\end{align*}
Multiplying on the left by $C_i$ and on the right by $B_i$, we obtain
\begin{align*}
G_i(z)-D=\sum_{k=0}^{\infty} z^{-k-1}C_iA_i^kB_i.
\end{align*}
The hypothesis says $G_1(z)=G_2(z)$ on the common resolvent set, hence in particular for all sufficiently large $|z|$. Therefore the two [Laurent series](/page/Laurent%20Series) at infinity agree as $p \times m$ matrix-valued series. This means that, for every row index $a \in \{1,\dots,p\}$ and column index $b \in \{1,\dots,m\}$, the scalar Laurent series given by the $(a,b)$ entries agree. By uniqueness of scalar Laurent coefficients applied to each entry, the matrix coefficients are equal:
\begin{align*}
C_1A_1^kB_1=C_2A_2^kB_2
\end{align*}
for every integer $k\ge 0$.
[/guided]
[/step]
[step:Define the state map on reachable sums and prove it is well-defined]
Let $\mathcal{R}_1 \subset \mathbb{R}^{n_1}$ denote the reachable subspace of the first realization:
\begin{align*}
\mathcal{R}_1
=
\operatorname{span}\{A_1^kB_1u : k\ge 0,\ u\in \mathbb{R}^m\}.
\end{align*}
Define a map
\begin{align*}
S:\mathcal{R}_1 \to \mathbb{R}^{n_2}
\end{align*}
as follows. If
\begin{align*}
x=\sum_{k=0}^{N}A_1^kB_1u_k
\end{align*}
for some $N\in\mathbb{N}\cup\{0\}$ and vectors $u_0,\dots,u_N\in\mathbb{R}^m$, set
\begin{align*}
Sx=\sum_{k=0}^{N}A_2^kB_2u_k.
\end{align*}
We prove that $S$ is well-defined. Suppose
\begin{align*}
\sum_{k=0}^{N}A_1^kB_1u_k=0.
\end{align*}
Set
\begin{align*}
y=\sum_{k=0}^{N}A_2^kB_2u_k \in \mathbb{R}^{n_2}.
\end{align*}
For every integer $\ell\ge 0$, using the Markov parameter identities from the previous step gives
\begin{align*}
C_2A_2^\ell y
=
\sum_{k=0}^{N}C_2A_2^{\ell+k}B_2u_k
=
\sum_{k=0}^{N}C_1A_1^{\ell+k}B_1u_k
=
C_1A_1^\ell\sum_{k=0}^{N}A_1^kB_1u_k
=
0.
\end{align*}
Since the second realization is observable, the only state $y\in\mathbb{R}^{n_2}$ satisfying $C_2A_2^\ell y=0$ for every $\ell\ge 0$ is $y=0$. Thus every linear relation among the vectors $A_1^kB_1u$ is sent to the corresponding zero relation among the vectors $A_2^kB_2u$, so $S$ is well-defined and linear.
[/step]
[step:Use reachability and observability to prove the state map is an isomorphism]
Since the first realization is reachable, $\mathcal{R}_1=\mathbb{R}^{n_1}$, so
\begin{align*}
S:\mathbb{R}^{n_1}\to\mathbb{R}^{n_2}
\end{align*}
is defined on the whole first state space. Since the second realization is reachable, every $y\in\mathbb{R}^{n_2}$ has the form
\begin{align*}
y=\sum_{k=0}^{N}A_2^kB_2u_k
\end{align*}
for some $N\in\mathbb{N}\cup\{0\}$ and $u_0,\dots,u_N\in\mathbb{R}^m$; hence $y=Sx$ for
\begin{align*}
x=\sum_{k=0}^{N}A_1^kB_1u_k.
\end{align*}
Thus $S$ is surjective.
To prove injectivity, let $x\in\mathbb{R}^{n_1}$ satisfy $Sx=0$. Since the first realization is reachable, write
\begin{align*}
x=\sum_{k=0}^{N}A_1^kB_1u_k.
\end{align*}
For every integer $\ell\ge 0$,
\begin{align*}
C_1A_1^\ell x
=
\sum_{k=0}^{N}C_1A_1^{\ell+k}B_1u_k
=
\sum_{k=0}^{N}C_2A_2^{\ell+k}B_2u_k
=
C_2A_2^\ell Sx
=
0.
\end{align*}
Since the first realization is observable, $x=0$. Therefore $S$ is injective. Since $S$ is both injective and surjective, it is a linear isomorphism. In particular $n_1=n_2$.
[/step]
[step:Verify the intertwining identities for $A$, $B$, and $C$]
For $u\in\mathbb{R}^m$, the definition of $S$ with $N=0$ gives
\begin{align*}
SB_1u=B_2u.
\end{align*}
Therefore
\begin{align*}
SB_1=B_2.
\end{align*}
Next let $x\in\mathbb{R}^{n_1}$. By reachability, choose $N\in\mathbb{N}\cup\{0\}$ and $u_0,\dots,u_N\in\mathbb{R}^m$ such that
\begin{align*}
x=\sum_{k=0}^{N}A_1^kB_1u_k.
\end{align*}
Then
\begin{align*}
SA_1x
=
S\left(\sum_{k=0}^{N}A_1^{k+1}B_1u_k\right)
=
\sum_{k=0}^{N}A_2^{k+1}B_2u_k
=
A_2Sx.
\end{align*}
Hence
\begin{align*}
SA_1=A_2S.
\end{align*}
Finally, for the same representation of $x$,
\begin{align*}
C_2Sx
=
\sum_{k=0}^{N}C_2A_2^kB_2u_k
=
\sum_{k=0}^{N}C_1A_1^kB_1u_k
=
C_1x.
\end{align*}
Thus
\begin{align*}
C_2S=C_1.
\end{align*}
[/step]
[step:Convert the state isomorphism into the required similarity matrix]
Define
\begin{align*}
T:\mathbb{R}^{n_2}\to\mathbb{R}^{n_1}
\end{align*}
by
\begin{align*}
T=S^{-1}.
\end{align*}
Since $S$ is a linear isomorphism, $T$ is invertible. From
\begin{align*}
SA_1=A_2S
\end{align*}
and $S=T^{-1}$, we obtain
\begin{align*}
T^{-1}A_1=A_2T^{-1}.
\end{align*}
Multiplying on the right by $T$ gives
\begin{align*}
A_2=T^{-1}A_1T.
\end{align*}
From $SB_1=B_2$, we obtain
\begin{align*}
B_2=T^{-1}B_1.
\end{align*}
From $C_2S=C_1$, we obtain
\begin{align*}
C_2=C_1T.
\end{align*}
These are precisely the required similarity identities.
[/step]