[state] (t2) [right of = b2, node distance = 2cm, fill = orange!30] $t_2$;
[state] (tn-1) [right of = bn-1, node distance = 2cm, inner sep = 0cm, fill = orange!30] $t_{n-1}$;
[state] (tn) [right of = bn, node distance = 2cm, fill = orange!30] $t_{n}$;

[state] (y) [fill = gray!40, right of = rhsdotstemp, node distance = 3cm] $y$;

[state, draw, rectangle] (init) [left of = rhsdots, node distance = 2.1cm, minimum width = 2cm, minimum height = 1cm] $\mathbf{D(n)}$;
(initr) [above of = init, node distance = 1.2cm] ;

(initr) edge (init)
(init) edge [bend left = 10] node $\frac{1}{n}$ (b1)
(init) edge [bend left = 10] node [below, pos = 0.7] $\frac{1}{n}$ (b2)
(init) edge [bend right = 10] node [above, pos = 0.7] $\frac{1}{n}$ (bn-1)
(init) edge [bend right = 10] node  [left]$\frac{1}{n}$ (bn)
(b1) edge node [above,sloped] $1$ (t1)
(b2) edge node [below,sloped] $1$ (t2)
(bn-1) edge node [above,sloped] $1$ (tn-1)
(bn) edge node [above,sloped] $1$ (tn)

(t1) edge node [above,sloped] $0.5 {-} \varepsilon$ (y)
(t2) edge node [above,sloped] $0.5 {-} 3\varepsilon$ (y)
(tn-1) edge node [above,sloped] $0.5 {-} (2n{-}3)\varepsilon$ (y)
(tn) edge node [above,sloped] $0.5 {-} (2n{-}1)\varepsilon$ (y)

(t1) edge [loop right, looseness = 4] node $0.5 {+} \varepsilon$ (t1)
(t2) edge [loop above, looseness = 4] node $0.5 {+} 3\varepsilon$  (t2)
(tn-1) edge [loop below, looseness = 4] node [xshift = -0.5cm] $0.5 {+} (2n{-}3)\varepsilon$ (tn-1)
(tn) edge [loop right, looseness = 4] node $0.5 {+} (2n{-}1)\varepsilon$ (tn)

(y) edge [loop above, looseness = 4] node $1$ (y)

[right of = s, node distance = 1.2cm] $\emptyset$;
[above of = a1, node distance = 0.7cm] $\{l_1\}$;
[above of = a2, node distance = 0.7cm] $\{l_2\}$;
[below of = an-1, node distance = 0.7cm] $\{l_{n-1}\}$;
[below of = an, node distance = 0.7cm] $\{l_n\}$;
[left of = s1, node distance = 0.75cm, yshift = 0.2cm] $\{a\}$;
[right of = s2, node distance = 0.75cm] $\{a\}$;
[below right of = s3, node distance = 0.75cm] $\{a\}$;
[above left of = sn-1, node distance = 0.75cm, yshift = 0.1cm, xshift = 0.2cm] $\{a\}$;
[right of = sn, node distance = 0.75cm] $\{a\}$;
[left of = sn+1, node distance = 0.75cm, yshift = -0.2cm] $\{a\}$;
[right of = x, node distance = 0.75cm] $\{b\}$;

[right of = init, node distance = 1.2cm] $\emptyset$;
[above of = b1, node distance = 0.7cm] $\{l_1\}$;
[above of = b2, node distance = 0.7cm] $\{l_2\}$;
[below of = bn-1, node distance = 0.7cm] $\{l_{n-1}\}$;
[below of = bn, node distance = 0.7cm] $\{l_n\}$;
[above of = t1, node distance = 0.7cm] $\{a\}$;
[below of = t2, node distance = 0.7cm] $\{a\}$;
[above of = tn-1, node distance = 0.7cm] $\{a\}$;
[below of = tn, node distance = 0.75cm] $\{a\}$;
[right of = y, node distance = 0.75cm] $\{b\}$;

Visualization of LMCs $\mathcal{M}_{n}$ and $\mathcal{N}_{n}$ used in the proof of <Ref>.

This theorem is a consequence of the following, more general result.

There is a family $\{(\mathcal{M}_{n}, \mathcal{N}_{n}) \mid n \in \mathbb{N}_{\geq 1} \}$ of pairs of finite LMCs such that for all $n \in \mathbb{N}_{\geq 1}$ we have $\mathcal{M}_{n} \sim_\varepsilon \mathcal{N}_{n}$ for $0 < \varepsilon \leq \frac{1}{4n}$, but $\mathcal{M}_n \not \simeq_\delta \mathcal{N}_n$ for any $\delta < 2n  \varepsilon$.

Let $n \in \mathbb{N}_{\geq 1}$ and $\varepsilon \in \left (0, \frac{1}{4n}\right]$. We consider the family $\{(\mathcal{M}_{n}, \mathcal{N}_{n}) \mid n \in \mathbb{N}_{\geq 1}\}$ consisting of pairs of LMCs constructed as those depicted in <Ref>.

In such a pair, both LMCs $\mathcal{M}_n$ and $\mathcal{N}_n$ consist of two parts. The first one, abstractly denoted $\mathbf{D}(n)$ in the figure, contains the initial state of the respective LMC and consists of a Markov chain that, by solely doing coin flips, generates a uniform distribution over the values $\{1,\dots,n \}$. The case for $n = 6$ is the well-known Knuth-Yao-Dice <cit.>, which can be generalized to uniformly sample over arbitrary nonnegative integers $n \in \mathbb{N}_{\geq 1}$ <cit.>. We assume that none of the states in $\mathbf{D}(n)$ is labeled by any atomic proposition and that every inner node of the LMC contains at most $2$ outgoing edges, both with probability exactly $\frac{1}{2}$. In the terminology of trees, the model has precisely $n$ leaves, and from the $i$-th leaf a transition from $\mathbf{D}(n)$ to the state $a_i$ and $b_i$, respectively, occurs with probability $1$.

The second part of the respective chains is given as follows: for $\mathcal{M}_{n}$, it has states $\{a_1, \dots, a_n, s_1, \dots, s_{n+1}, x\}$ with transition probabilities

(a_i)(s_i)    = 1/2           and     (a_i)(s_i+1)       = 1/2

(s_i)(s_i)    = 1/2 + 2  (i-1) ε           and     (s_i)(x)       = 1/2 - 2(i-1)ε

(x)(x)    = 1,

while its state space in $\mathcal{N}_{n}$ is $\{b_1, \dots, b_n, t_1, \dots, t_n, y\}$ and the transition probabilities are given as

\begin{align*}
\prob(b_i)(t_i) &= 1\\
\prob(t_i)(t_i) &= \frac{1}{2} + (2i -1) \varepsilon \qquad \text{and} \qquad \prob(t_i)(y) = \frac{1}{2} - (2i - 1) \varepsilon\\
\prob(y)(y) &= 1.
\end{align*}

Furthermore, $AP$ consists of the pairwise distinct atomic propositions $l_1, \dots, l_n, a, b$. The state labels are $l(a_i) = l(b_i) = \{l_i\}$ and $l(s_i) = l(t_i) = \{a\}$ for $1 \leq i \leq n$ as well as $l(x) = l(y) = \{b\}$. Note that the states in $\mathbf{D}(n)$ all have empty labels, so they cannot be related to any of the states $a_i, b_i, s_i, t_i, x$ or $y$ by any relation that preserves the state labels.

The relation $\sim_\varepsilon$ in $\mathcal{M}_n \oplus \mathcal{N}_n$ is the symmetric and reflexive closure of the relation that contains all pairs of corresponding states in the respective copies of $\mathbf{D}(n)$, as well as

\begin{align*}
\{(a_i, b_i) \mid 1 \leq i \leq n \} \cup \{(x, y)\} \cup \{(s_i, t_i), (s_{i+1}, t_i) \mid 1 \leq i \leq n \}.
\end{align*}

In particular, $s_{init}^{\mathcal{M}_n} \sim_\varepsilon s_{init}^{\mathcal{N}_n}$.

Let $\delta \geq 0$ such that $\mathcal{M}_n \simeq_\delta \mathcal{N}_n$, and let $\mathcal{M}_{n}', \mathcal{N}_{n}'$ be the corresponding $\delta$-perturbations, i.e., $\mathcal{M}_n' \sim \mathcal{N}_n'$. We now show $\delta \geq 2n\varepsilon$.

To this end, we first observe that every two states $a_i, a_j$ and $a_i, b_j$ and $b_i, b_j$ with $i \neq j$ have different labels, while $l(a_i) = l(b_i)$ for every $1 \leq i \leq n$. Furthermore, $\mathcal{M}_n' \sim \mathcal{N}_n'$ iff ${s_{init}^{\mathcal{M}_n}}' \sim {s_{init}^{\mathcal{N}_n}}'$.
However, this can only be the case if $a_i'\sim b_i'$ for all $i$. Otherwise, we would have to perturb the states in $\mathbf{D}(n)$ in such a way that, for some $i$, the states $a_i$ and $b_i$ are not reachable anymore when traversing $\mathbf{D}(n)$. But, as all transition probabilities from one state in $\mathbf{D}(n)$ to another one in $\mathbf{D}(n)$ are equal to $\frac{1}{2}$, and the probability to transition from the $i$-th leaf of $\mathbf{D}(n)$ to $a_i$ or $b_i$ equals $1$ for all $i$, this requires a perturbation by at least $1 > 2n\varepsilon$. Hence, the only possible way to ensure $\mathcal{M}_n' \sim \mathcal{N}_n'$ for some $\delta$-perturbations with $\delta < 2n\varepsilon$ is to enforce $a_i' \sim b_i'$ for all $i$.

Let $1 \leq i \leq n$. First assume that, in $\mathcal{M}_n'$ and $\mathcal{N}_n'$, we have $s_i' \nsim t_i'$ and $s_{i+1}' \nsim t_i'$. Then it holds that $\prob(a_i)([t_i']_{\sim}) = 0$ while $\prob(b_i)([t_i']_{\sim}) = 1$, and it is easy to see that a perturbation of at least $1 \geq 2n\varepsilon$ is required to make the states $a_i'$ and $b_i'$ bisimilar. Now, assume that $s_i' \sim t_i'$ and $s_{i+1}' \nsim t_i'$. Here, $\prob(a_i)([t_i']_{\sim}) = \frac{1}{2}$ while $\prob(b_i)([t_i']_{\sim}) = 1$, and thus a perturbation of at least $\frac{1}{2} \geq 2n\varepsilon$ is necessary to obtain $a_i' \sim b_i'$ (by assigning, for example, a probability of $\frac{3}{4}$ to the transition $a_i' \to s_i'$ and $b_i' \to t_i'$, as well as a total probability of $\frac{1}{4}$ to transitions $a_i' \to [s_{i+1}']_\sim$ and $b_i' \to [s_{i+1}']_\sim$, respectively). By an analogous argument we can show that if $s_i' \nsim t_i'$ and $s_{i+1}' \sim t_i'$, a perturbation of at least $2n \varepsilon$ is required as well.

Thus, we can w.l.o.g. assume that, if $\delta < 2n\varepsilon$, it has to holds that $s_i' \sim t_i' \sim s_{i+1}'$ for every $1 \leq i \leq n$. But then the transitivity of $\sim$ implies $s_1' \sim s_{n+1}'$, which by construction is only possible when perturbing by at least $2n\varepsilon$: in the resulting model, probabilities of $\frac{1}{2} + n \varepsilon$ are assigned to the transitions $s_1' \to s_1'$ and $s_{n+1}' \to s_{n+1}'$, while the transitions $s_1' \to x'$ and $s_{n+1}' \to x'$ are taken with probability $\frac{1}{2} -n\varepsilon$.

All in all, this shows that for $\mathcal{M}_n \simeq_\delta \mathcal{N}_n$ a tolerance of at least $\delta = 2n\varepsilon$ is required, or equivalently, that $\mathcal{M}_n \not \simeq_\delta \mathcal{N}_n$ for any $\delta < 2n\varepsilon$.

Let $\varepsilon \in \left( 0, \frac{1}{4} \right]$ and set $n = n(\varepsilon) = \left \lfloor \frac{1}{4 \varepsilon} \right \rfloor$. As $\left \lfloor \frac{1}{4\varepsilon} \right \rfloor \leq \frac{1}{4 \varepsilon}$, this choice of $n$ implies both $\varepsilon \leq \frac{1}{4n}$ and $\frac{1}{4\varepsilon} < n+1$.

From the first inequality, it follows by <Ref> that $\mathcal{M}_n \sim_{\varepsilon} \mathcal{N}_n$ and $\mathcal{M}_n \not \simeq_\delta \mathcal{N}_n$ for any $\delta < 2 n \varepsilon$, where $\mathcal{M}_n$ and $\mathcal{N}_n$ are again as depicted in <Ref>.
Furthermore, $\frac{1}{4 \varepsilon} < n+1$ implies $\frac{1}{2} \cdot \frac{n}{n+1} < \frac{1}{2} \cdot \frac{n}{1/(4\varepsilon)} = 2 n \varepsilon$, so in particular $\mathcal{M}_n \not \simeq_\delta \mathcal{N}_n$ for any $\delta \leq \frac{1}{2} \cdot \frac{n}{n+1}$. As $f(x) = \frac{1}{2} \cdot \frac{x}{x+1}$ is strictly increasing in $x > 0$, and because $n = \lfloor \frac{1}{4\varepsilon} \rfloor \geq 1$ for any $\varepsilon \in (0, \frac{1}{4}]$, it follows that for any such $\varepsilon$ and $n$ we have $\frac{1}{4} = \frac{1}{2} \cdot \frac{1}{1+1} = f(1) \leq f(n) = \frac{1}{2} \cdot \frac{n}{n+1} < 2n\varepsilon$. In combination, this yields $\mathcal{M}_n \not \simeq_\delta \mathcal{N}_n$ for any $\delta \leq \frac{1}{4}$.

[->,>=stealth',shorten >=1pt,auto, semithick]
every state = [text = black]

[state] (s) [fill = yellow] $s$;
(stemp) [right of = s, node distance = 2.8cm] ;
(stempl) [left of = stemp, node distance = 2cm, yshift = 0.1cm] $\vdots$;
(stempr) [right of = stemp, node distance = 2cm, yshift = 0.1cm] $\vdots$;
[state] (s1) [above of = stemp, node distance = 3cm, fill = green] $s_1$;
[state] (s2) [above of = stemp, node distance = 1.5cm, fill = green] $s_2$;
[state] (sn) [below of = stemp, node distance = 1.5cm, fill = green] $s_{n}$;
[state] (sn+1) [below of = stemp, node distance = 3cm, fill = green, inner sep = 0cm] $s_{n+1}$;
[state] (x) [right of = stemp, node distance = 2.8cm] $x$;

(sinit) [above of = s, node distance = 1cm] ;

(sinit) edge (s)
(s) edge[bend left = 20] node [above,sloped] $\frac{n+2}{(n+1)^2}$ (s1)
(s) edge node [above,sloped] $\frac{n+2}{(n+1)^2}$ (s2)
(s) edge node [below,sloped] $\frac{n+2}{(n+1)^2}$ (sn)
(s) edge [bend right = 20] node [below,sloped] $\frac{1}{(n+1)^2}$ (sn+1)

(s1) edge [bend left = 20] node [above,sloped] $0.5 {+} \varepsilon(n{-}1)$ (x)
(s2) edge node [above,sloped] $0.5 {+} \varepsilon(n{-}3)$ (x)
(sn) edge node [below,sloped] $0.5 {+} \varepsilon({-}n {+} 1)$ (x)
(sn+1) edge [bend right = 20] node [below,sloped, pos = 0.6] $0.5 {+} \varepsilon({-}n{-}1)$ (x)

(s1) edge [loop right] node $0.5 {-} \varepsilon(n{-}1)$ (s1)
(s2) edge [loop below] node $0.5 {-} \varepsilon(n{-}3)$ (s2)
(sn) edge [loop above] node $0.5 {-} \varepsilon({-}n{+}1)$ (sn)
(sn+1) edge [loop right] node $0.5 {-} \varepsilon({-}n{-}1)$ (sn+1)

(x) edge [loop above] node [above] $1$ (x)

[state] (t) [fill = yellow, right of = x, node distance = 2cm] $t$;
(ttemp) [right of = t, node distance = 2.8cm] ;
(ttempr) [right of = ttemp, node distance = 2cm, yshift = 0.1cm] $\vdots$;
(ttempl) [left of = ttemp, node distance = 2cm, yshift = 0.1cm] $\vdots$;
[state] (t1) [above of = ttemp, node distance = 3cm, fill = green] $t_1$;
[state] (t2) [above of = ttemp, node distance = 1.5cm, fill = green] $t_2$;
[state] (tn) [below of = ttemp, node distance = 1.5cm, fill = green] $t_{n}$;
[state] (tn+1) [below of = ttemp, node distance = 3cm, fill = green, inner sep = 0cm] $t_{n+1}$;
[state] (y) [right of = ttemp, node distance = 2.8cm] $y$;
(tinit) [above of = t, node distance = 1cm] ;

(tinit) edge (t)
(t) edge [bend left = 20] node [above,sloped] $\frac{1}{(n+1)^2}$ (t1)
(t) edge node [above,sloped] $\frac{n+2}{(n+1)^2}$ (t2)
(t) edge node [below,sloped] $\frac{n+2}{(n+1)^2}$ (tn)
(t) edge [bend right = 20] node [below,sloped] $\frac{n+2}{(n+1)^2}$ (tn+1)

(t1) edge [bend left = 20] node [above,sloped] $0.5 {+} n\varepsilon$ (y)
(t2) edge node [above,sloped] $0.5 {+} \varepsilon(n{-}2)$ (y)
(tn) edge node [below,sloped] $0.5 {+} \varepsilon({-}n {+} 2)$ (y)
(tn+1) edge [bend right = 20] node [below,sloped] $0.5 {+} \varepsilon({-}n)$ (y)

(t1) edge [loop right] node $0.5 {-} n\varepsilon$ (t1)
(t2) edge [loop below] node $0.5 {-} \varepsilon(n{-}2)$ (t2)
(tn) edge [loop above] node $0.5 {-} \varepsilon({-}n{+}2)$ (tn)
(tn+1) edge [loop right] node $0.5 {-} \varepsilon({-}n)$ (tn+1)

(y) edge [loop above] node[above] $1$ (y)

[below of = s, node distance = 0.7cm] $\{a\}$;
[left of = s1, node distance = 0.75cm] $\{b\}$;
[above of = s2, node distance = 0.7cm] $\{b\}$;
[below of = sn, node distance = 0.7cm] $\{b\}$;
[left of = sn+1, node distance = 0.75cm] $\{b\}$;
[right of = x, node distance = 0.75cm] $\{c\}$;

[below of = t, node distance = 0.7cm] $\{a\}$;
[left of = t1, node distance = 0.75cm] $\{b\}$;
[above of = t2, node distance = 0.7cm] $\{b\}$;
[below of = tn, node distance = 0.7cm] $\{b\}$;
[left of = tn+1, node distance = 0.75cm] $\{b\}$;
[right of = y, node distance = 0.75cm] $\{c\}$;

The LMCs $\mathcal{M}_n$ (left) and $\mathcal{N}_n$ (right) for $n \in \mathbb{N}_{\geq 1}$ used in the proof of <Ref>.

Let $n \in \mathbb{N}$ with $n \geq 1$. The family $\mathcal{F} = \{(\mathcal{M}_n, \mathcal{N}_n) \mid n \geq 1\}$ of LMCs that we consider is constructed as in <Ref>. All of these models are well-defined, as every transition probability is in $[0,1]$ by the choice of $\varepsilon$, and

\begin{align*}
\frac{1}{(n+1)^2} + n \cdot \frac{n+2}{(n+1)^2} = \frac{1 + n(n+2)}{(n+1)^2} = 1.
\end{align*}

To illustrate our construction, the LMCs $\mathcal{M}_2$ and $\mathcal{N}_2$ can be found in <Ref>. By checking the conditions of $\varepsilon$-bisimilarity, it follows that always $\mathcal{M}_n \sim_\varepsilon \mathcal{N}_n$, and that the relations $\sim_{\varepsilon}$ are the reflexive and symmetric closures of

\begin{align*}
\textcolor{white}{\cup}&\{(s,t), (x,y)\}
\cup \{(s_i, t_i) \mid i = 1,\dots,n + 1 \} \cup \{ s_i, t_{i+1}) \mid i = 1, \dots, n\}.
\end{align*}

In particular, we have $\vert {\sim_\varepsilon}(\{s_{n+1}\}) \vert = \vert {\sim_\varepsilon}(\{t_{1}\}) \vert = 2$, while $\vert {\sim_\varepsilon}(\{s_{i}\}) \vert = \vert {\sim_\varepsilon}(\{t_{j}\}) \vert = 3$ for all remaining choices of $i, j$.

Now assume that $\mathcal{M}_n \simeq_\delta \mathcal{N}_n$, and let $\mathcal{M}_n', \mathcal{N}_n'$ be the corresponding $\delta$-perturbations of $\mathcal{M}_n$ and $\mathcal{N}_n$, i.e., let $\mathcal{M}_n' \sim \mathcal{N}_n'$.

For $\mathcal{M}_n' \sim \mathcal{N}_n'$ to hold, it follows by construction that we need to have $s' \sim t'$ and $x' \sim y'$ in $\mathcal{M}_n' \oplus \mathcal{N}_n'$. While the latter is obviously possible for all $\delta \geq 0$, the former requires for every $\sim$-class $C$ that $\prob'(s')(C') = \prob'(t')(C')$, and in particular that

\begin{align}
\prob'(s')([s_{n+1}']_\sim) = \prob'(t')([s_{n+1}']_\sim). \label{Proof: Theorem: No common nepsilon quotient with graph isomorphism - Star 1}
\end{align}

By the choice of $\varepsilon$ it follows that if we perturb $\mathcal{M}_n$ such that the probability to move to $s_{n+1}' \in [s_{n+1}']$ becomes $0$, $\delta$ already has to be $\geq n \varepsilon$, which shows the claim. Hence, we can w.l.o.g. assume that $\prob'(s')([s_{n+1}']_\sim) > 0$, and therefore also that $\prob'(t')([s_{n+1}']_\sim) > 0$.
The equivalence class in question is thus of the form

\begin{align*}
[s_{n+1}']_{\sim} = \{m_1', \dots, m_k', n_1', \dots, n_l' \mid m_1', \dots, m_k' \in S^{\mathcal{M}_n'}, n_1', \dots, n_l' \in S^{\mathcal{N}_n'}\},
\end{align*}

for $k, l \geq 1$. Further, the construction of $\mathcal{M}_n$ yields that the probability to move from $s$ to any set of states $A \subseteq S$ in $\mathcal{M}_n$ that contains $s_{n+1}$ is precisely

\begin{align*}
\prob(s)(A) = \frac{1}{(n+1)^2} + (\vert A \vert -1) \frac{n+2}{(n+1)^2} = \frac{1 + (\vert A \vert - 1)(n+2)}{(n+1)^2}.
\end{align*}

Because we assume $\mathcal{M}_n'$ to be a $\delta$-perturbation of $\mathcal{M}_n$, it now follows from <Ref> that in $\mathcal{M}_n'$ we have $\prob'(s')([s_{n+1}']_\sim)  \in [lb(s'), ub(s')]$ for

\begin{align*}
lb(s') &= \frac{1+(k-1)(n+2)}{(n+1)^2} - \frac{\delta}{2} \quad \text{ and } \quad
ub(s') = \frac{1+(k-1)(n+2)}{(n+1)^2} + \frac{\delta}{2}.
\end{align*}

For the probability $\prob'(t')([s_{n+1}']_\sim)$, on the other hand, we have to consider two cases. If $t_1' \in [s_{n+1}']_\sim$ then $t_1' \sim s_{n+1}'$. But, by construction, making these two states bisimilar requires a perturbation by at least $\delta = 2n \varepsilon$ for both $\mathcal{M}_n$ and $\mathcal{N}_n$. Hence, in this case, $\delta \geq n \varepsilon$, which proves the claim.

Otherwise, if $t_1' \notin [s_{n+1}']_\sim$, then it follows again from the construction and <Ref> that $\prob'(t')([s_{n+1}']_\sim) \in [lb(t'), ub(t')]$, where the bounds are given as

\begin{align*}
lb(t') = \frac{l(n+2)}{(n+1)^2} - \frac{\delta}{2} \quad \text{ and } \quad ub(t') = \frac{l(n+2)}{(n+1)^2} + \frac{\delta}{2}.
\end{align*}

In particular, we <Ref> can only hold if there is a $p \in [lb(s'), ub(s')] \cap [lb(t'), ub(t')]$ with $p = \prob'(s')([s_{n+1}']_\sim) = \prob'(t')([s_{n+1}']_\sim)$.
We utilize this fact to show $\delta \geq n\varepsilon$, by doing a case distinction on how $[lb(s'), ub(s')]$ and $[lb(t'), ub(t')]$ can intersect.

Case 1. $lb(s') \leq lb(t') \leq ub(s') \leq ub(t')$.
The first and third inequality are equivalent to $1 \leq (l - k + 1)(n+2)$, and since $n+2 \geq 3$ and $l,k \in \mathbb{N}_{\geq 1}$ this implies $l - k + 1 \geq 1$. Moreover, the inequality $lb(t') \leq ub(s')$ is equivalent to $\frac{(n+2)(l-k+1) -1}{(n+1)^2} \leq \delta$. As $\varepsilon \in \left(0, \frac{1}{n(n+1)^2} \right]$, we can combine all of this to get

\begin{align*}
n \varepsilon \leq \frac{1}{(n+1)^2} \leq \frac{n+2 - 1}{(n+1)^2} \leq \frac{(n+2)(l-k+1) - 1}{(n+1)^2} \leq \delta.
\end{align*}

Case 2. $lb(s') \leq lb(t') \leq ub(t') \leq ub(s')$.
Similar to Case 1, $lb(s') \leq lb(t')$ implies $l - k + 1\geq 1$, and thus $l \geq k$. On the other hand, $ub(t') \leq ub(s')$ yields $(l - k + 1)(n+2) \leq 1$, which, because $n +2 \geq 3$ and $l,k \in \mathbb{N}_{\geq 1}$, gives us $l-k+1 \leq 0$. But then $l + 1 \leq k$, so $l < k$, a contradiction to $l \geq k$. This case can therefore not occur.

Case 3. $lb(t') \leq lb(s') \leq ub(s') \leq ub(t')$.
As for Case 2 we can reach a contradiction by observing that $lb(t') \leq lb(s')$ implies $l < k$ while $ub(s') \leq ub(t')$ yields $l \geq k$.

Case 4. $lb(t') \leq lb(s') \leq ub(t') \leq ub(s')$.
The first and third inequality are equivalent to $(n+2)(l-k+1) \leq 0$. Together with $n + 2 \geq 3$ and $l,k \in \mathbb{N}_{\geq 1}$ this yields $l-k+1 \leq 0$, or, equivalently, $k - l - 1 \geq 0$. Furthermore, $lb(s') \leq ub(t')$ can be transformed into $\frac{(n+2)(k-l-1) +1}{(n+1)^2} \leq \delta$. Similar to Case 1 we conclude

\begin{align*}
n \varepsilon \leq \frac{1}{(n+1)^2} \leq \frac{(n+2)(k-l-1) + 1}{(n+1)^2} \leq \delta.
\end{align*}

All in all, this shows that $\mathcal{M}_n \simeq_\delta \mathcal{N}_n$ only if $\delta \geq n \varepsilon$, as otherwise <Ref> is violated.

We set $S = S^\mathcal{M} \oplus S^\mathcal{N}$ and $\prob = \prob^{\mathcal{M} \oplus \mathcal{N}}$.

(i) $\Leftrightarrow$ (ii). The equivalence is clear due to the $1$-to-$1$-correspondence of $\varepsilon$-perturbations of $\mathcal{M} \oplus \mathcal{N}$ and pairs of $\varepsilon$-perturbations of $\mathcal{M}$ and $\mathcal{N}$.

(ii) $\Rightarrow$ (iii). Let $\mathcal{L}$ be an $\varepsilon$-perturbation of $\mathcal{M} \oplus \mathcal{N}$ in which $s_{init}^\mathcal{M} \sim s_{init}^\mathcal{N}$. Following an argument similar to the one in the proof of <Ref>, $\sim_\mathcal{L}$ on $\mathcal{L}$ is an $\varepsilon$-bisimulation on $\mathcal{M} \oplus \mathcal{N}$. Further, as $s_{init}^\mathcal{M} \sim_\mathcal{L} s_{init}^\mathcal{N}$ by assumption and $\sim_\mathcal{L}$ is clearly transitive, this relation satisfies all the requirements posed in (iii).

Now let $A \in \xfrac{S}{\sim_\mathcal{L}}$. The distribution $\prob_A^* \in Distr(\xfrac{S}{\sim_\mathcal{L}})$ with $\prob_A^*(B) = \prob_\mathcal{L}(u')(B')$ for some (or, equivalently, all) $u \in A$ satisfies the requirements of (iii), as by <Ref> it holds for all $\sim_\mathcal{L}$-closed $C \subseteq S$ that

\begin{align*}
\vert \prob(s)(C) - \prob_A^*(C) \vert = \vert \prob(s)(C) - \prob_\mathcal{L}(s')(C') \vert \leq \frac{\varepsilon}{2}.
\end{align*}

(iii) $\Rightarrow$ (ii). We set $S = S_\mathcal{M} \oplus S_\mathcal{N}$. Let $R$ be a transitive $\varepsilon$-bisimulation on $\mathcal{M} \oplus \mathcal{N}$ such that $(s_{init}^\mathcal{M}, s_{init}^\mathcal{N}) \in R$, and such that for every $A \in \xfrac{S}{R}$ there is a distribution $\prob_A^*\in Distr(\xfrac{S}{R})$ which satisfies for all $R$-closed sets $C \subseteq S$ and all $s \in A$ that $\vert \prob(s)(C) - \prob_A^*(C) \vert \leq \frac{\varepsilon}{2}$.

For $B \in \xfrac{S}{R}$, we define for fixed $A \in \xfrac{S}{R}$ and $s \in A$ the value $\Delta(s, B) = \prob(s)(B) - \prob_A^*(B)$, and set

\begin{align*}
B^+_s = \{B \in \xfrac{S}{R} \mid \Delta(s, B) > 0\} \qquad \text{ and } \qquad B^-_s = \{B \in \xfrac{S}{R} \mid \Delta(s, B) < 0\}.
\end{align*}

Note that, if $B \in B^+_s$ we can pick, for all $t \in B$, values $0 \leq \delta^+(s,t) \leq \prob(s)(t)$ such that $\sum_{t \in B} \delta^+(s,t) = \Delta(s,B)$. Similarly, if $B \in B^-_s$, there are $0 \leq \delta^-(s,t) \leq 1-\prob(s)(t)$ for all $t \in B$ such that $\sum_{t \in B} \delta^-(s,t) = \vert \Delta(s,B) \vert$.

Define a new LMC $\mathcal{L}$ on $S$, whose labeling is just like that of $\mathcal{M} \oplus \mathcal{N}$ and whose transition distributions $\prob_\mathcal{L}(s) \in Distr(S)$ for any $s \in S$ are, for $B \in \xfrac{S}{R}$ and $t \in B$, given by

\begin{align*}
\prob_\mathcal{L}(s)(t) = \begin{cases}
\prob(s)(t), & \text{if } \Delta(s,B) = 0 \\
\prob(s)(t) - \delta^+(s, t) & \text{if } B \in B^+_s \\
\prob(s)(t) + \delta^-(s,t) & \text{if} B \in B^-_s
\end{cases}.
\end{align*}

Then $\prob_\mathcal{L}(s)$ is a distribution for all $s \in S$, as

\begin{align*}
&\sum_{t \in S} \prob_\mathcal{L}(s)(t) \\&= \sum_{B \in S/R} \prob_\mathcal{L}(s)(B) \\
&= \sum_{\substack{B \in S/R,\\ \Delta(s,B) = 0}} \prob_\mathcal{L}(s)(B) + \sum_{B \in B^+_s} \prob_\mathcal{L}(s)(B) + \sum_{B \in B^-_s} \prob_\mathcal{L}(s)(B) \\
&= \sum_{\substack{B \in S/R,\\ \Delta(s,B) = 0}} \prob(s)(B) + \sum_{B \in B^+_s} (\prob(s)(B) - \Delta(s,B)) + \sum_{B \in B^-_s} (\prob(s)(B) + \vert \Delta(s,B) \vert) \\
&= \sum_{\substack{B \in S/R,\\ \Delta(s,B) = 0}} \prob_A^*(B) + \sum_{B \in B^+_s} \prob_A^*(B) + \sum_{B \in B^-_s} (\prob(s)(B) + \vert \prob(s)(B) - \prob_A^*(B) \vert) \\
&= \sum_{\substack{B \in S/R,\\ \Delta(s,B) \geq 0}} \prob_A^*(B) + \sum_{B \in B^-_s} (\prob(s)(B) - \prob(s)(B) + \prob_A^*(B)) \\
&=  \sum_{B \in S/R} \prob_A^*(B) = 1.
\end{align*}

Furthermore, $R$ is a probabilistic bisimulation on $\mathcal{L}$, as for all $s,t \in A$ and $B \in \xfrac{S}{R}$ the construction yields $\prob_\mathcal{L}(s)(B) = \prob_A^*(B) = \prob_\mathcal{L}(t)(B)$.

Now fix $s \in A$. Then

\begin{align*}
1&= \sum_{B \in S/R} \prob(s)(B) = \sum_{B \in S/R} (\prob(s)(B) + \prob_A^*(B) - \prob_A^*(B)) \\&= \sum_{B \in S/R} \prob_A^*(B) + \sum_{B \in S/R} \Delta(s, B) = 1 + \sum_{B \in S/R} \Delta(s, B),
\end{align*}

so $\sum_{B \in S/R} \Delta(s, B) = 0$
and hence $\sum_{B \in B^+_s} \Delta(s, B) = - \sum_{B \in B^-_s} \Delta(s, B)$, or, equivalently,

\begin{align}
\sum_{B\in B^+_s} \Delta(s,B) = \sum_{B \in B^-_s} \vert \Delta(s,B) \vert \label{Proof: CharacterizationPertubedEpsBisimTransitiveEpsBisim - Eq1}.
\end{align}

\begin{align*}
\Vert \prob(s) - \prob_\mathcal{L}(s) \Vert_1 &= \sum_{t \in S} \vert \prob(s)(t) - \prob_\mathcal{L}(s)(t) \vert \\
&= \sum_{B \in S/R} \sum_{t \in B} \vert \prob(s)(t) - \prob_\mathcal{L}(s)(t) \vert \\
&= \underbrace{\sum_{\substack{B \in S/R, \\ \Delta(s, B) = 0}} \sum_{t \in B} \vert \prob(s)(t) - \prob_\mathcal{L}(s)(t) \vert}_{=0} + \sum_{\substack{B \in S/R, \\ \Delta(s, B) \neq 0}} \sum_{t \in B} \vert \prob(s)(t) - \prob_\mathcal{L}(s)(t) \vert \\
&= \underbrace{\sum_{B \in B^+_s} \sum_{t \in B} \vert \prob(s)(t) - \prob_\mathcal{L}(s)(t) \vert}_{= \sum_{B \in B^+_s} \Delta(s, B)}  + \underbrace{\sum_{B \in B^-_s} \sum_{t \in B} \vert \prob(s)(t) - \prob_\mathcal{L}(s)(t) \vert}_{= \sum_{B \in B^-_s} \vert \Delta(s, B) \vert \overset{(\ref{Proof: CharacterizationPertubedEpsBisimTransitiveEpsBisim - Eq1})}{=} \sum_{B \in B^+_s} \Delta(s,B)}\\
&= 2 \cdot \sum_{B \in B^+_s} \Delta(s,B) \\
&= 2 \cdot \sum_{B \in B^+_s} (\prob(s)(B) - \prob_A^*(B)) \\
&= 2 \cdot \left \vert \prob(s)\left(\bigcup_{B \in B^+_s}B\right) - \prob_A^*\left(\bigcup_{B \in B^+_s}B\right) \right \vert \\
&\leq 2 \cdot \frac{\varepsilon}{2} = \varepsilon,
\end{align*}

where the inequality follows from the properties of $\prob_A^*$ and the fact that $\bigcup_{B \in B^+_s} B$ is, as a union of $R$-equivalence classes, $R$-closed. All in all, this shows that $\mathcal{L}$ is indeed an $\varepsilon$-perturbation of $\mathcal{M} \oplus \mathcal{N}$ in which, as $R$ is a bisimulation on this model and since $(s_{init}^\mathcal{M}, s_{init}^\mathcal{N}) \in R$ by assumption, it holds that $s_{init}^\mathcal{M} \sim s_{init}^\mathcal{N}$.

(i) $\Rightarrow$ (ii). Let $\mu^* \in Distr(X)$ be as in (i), and let $l \in \{1, \dots, k\}$. We define $B_{l}^+ = \{i \in X \mid \mu_l(i) \geq \mu^*(i)\}$ and $B_{l}^- = \{i \in X \mid \mu_l(i) < \mu^*(i)\}$. Then $X = B_{l}^+ \uplus B_l^{-}$ and

\begin{align*}
\Vert \mu_l - \mu^* \Vert_1 &= \sum_{i \in X} \vert \mu_l(i) - \mu^*(i) \vert \\
&= \sum_{i \in B_{l}^+} \vert \mu_l(i) - \mu^*(i) \vert + \sum_{i \in B_l^-} \vert \mu_l(i) - \mu^*(i) \vert \\
&= \sum_{i \in B_{l}^+} (\mu_l(i) - \mu^*(i)) + \sum_{i \in B_l^-} (\mu^*(i) - \mu_l(i)) \\
&= \underbrace{\mu_l(B_{l}^+) - \mu^*(B_{l}^+)}_{\geq 0} + \underbrace{\mu^*(B_{l}^-) - \mu_l(B_l^-)}_{\geq 0} \\
&= \vert \mu_l(B_{l}^+) - \mu^*(B_{l}^+) \vert + \vert \mu_l(B_{l}^-) - \mu^*(B_l^-) \vert \\
&\leq \frac{\varepsilon}{2}+ \frac{\varepsilon}{2} = \varepsilon
\end{align*}

where the last inequality follows from the fact that $B_{l}^+, B_l^- \subseteq X$ and the properties of $\mu^*$. As we have chosen $l \in \{1, \dots, k\}$ arbitrarily it follows that $\mu = \mu^*$ is a suitable choice for a distribution that satisfies the conditions in (ii).

(ii) $\Rightarrow$ (i). This is just a special case of <Ref> and can be proved analogously.

(ii) $\Rightarrow$ (iii). Let $\mu \in Distr(X)$ be a distribution as in (ii). We have to show that there is a non-negative solution for the linear constraint system of (iii). To this end, set $x_i = \mu(i)$ and $\delta_{l, i} = \vert \mu_l(i) - \mu(i) \vert$ for all $i \in X, l \in \{1, \dots, k\}$. Then $x_i, \delta_{l, i} \geq 0$ for all $i, l,$ and:

* $\sum_{i \in X} x_i = \sum_{i \in X} \mu(i) = 1$ as $\mu \in Distr(X)$, so the first constraint is satisfied.

* For all $i, l,$ we have $\delta_{l, i} = \vert \mu_l(i) - \mu(i) \vert$ which implies $- \delta_{l, i} \leq \mu_l(i) - \mu(i) \leq \delta_{l, i}$. As $x_i = \mu(i)$ for all i it follows that $- \delta_{l, i} \leq \mu_l(i) - x_i \leq \delta_{l, i}$, so the second and third constraint are satisfied.

* $\sum_{i \in X} \delta_{l, i} = \sum_{i \in X} \vert \mu_l(i) - \mu(i) \vert = \Vert \mu_l - \mu \Vert_1 \leq \varepsilon$ for all $i, l,$ by the assumption on $\mu$, so the fourth constraint is satisfied.

All in all this shows that the chosen assignment of values to the variables $x_i$ and $\delta_{l, i}$ does indeed constitute a non-negative solution of the linear constraint system.

(iii) $\Rightarrow$ (ii). Let $\{x_i \mid i \in X\} \cup \{\delta_{l, i} \mid l \in \{1, \dots, k\}, i \in X\}$ be a non-negative solution of the linear constraint system of (iii). Define $\mu: X \to [0,1]$ via $\mu(i) = x_i$ for all $i \in X$. Then $\mu$ is a distribution as $\sum_{i \in X} \mu(i) = \sum_{i \in X} x_i = 1$ by the first constraint, and for $l \in \{1, \dots, k\}$,

\begin{align*}
\Vert \mu_l - \mu \Vert_1 = \sum_{i \in X} \vert \mu_l(i) - \mu(i) \vert = \sum_{i \in X} \vert \mu_l(i) - x_i \vert \overset{\text{Constr. } 2 + 3} = \sum_{i \in X} \delta_{l, i} \overset{\text{Constr. } 4}{\leq} \varepsilon.
\end{align*}

Thus, $\mu$ is a distribution satisfying the requirements of (ii).

We start with the first claim, i.e., we start by showing that it is -complete to decide if $\mathcal{M} \simeq_\varepsilon \mathcal{N}$ for given $\mathcal{M}, \mathcal{N}$ and $\varepsilon > 0$.
The proof ideas take inspiration from the proof of <cit.>, which can be found in the full version <cit.> of <cit.>.

Let $\mathcal{M}$ and $\mathcal{N}$ be two LMCs, and set $S = S^\mathcal{M} \oplus S^\mathcal{N}$. We consider the following nondeterministic polynomial time algorithm:

* Nondeterministically guess a partition $X$ of $S$ such that

* $X$ only relates states with the same label

* For each block $B \in X$ it holds that $B \cap S_1 \neq \emptyset \neq B \cap S_2$

* There is a block $B_{init} \in X$ such that $s_{init}^{\mathcal{M}_1}, s_{init}^{\mathcal{M}_2} \in B_{init}$

* For every $B \in X$, apply <Ref> to $X$ and $\{\prob(s) \mid s \in B\}$, where the distributions are lifted to $X$ in the natural way, to check if there is a $\prob_B \in Distr(X)$ that satisfies condition (ii) of the lemma. If not, go back to the first step.

This algorithm has a solution iff $\mathcal{M} \simeq_\varepsilon \mathcal{N}$. The direction from right to left is clear, so we focus on the implication from left to right. Let $X, \{\prob_B \mid B \in X\}$ be a solution of the algorithm, and let $s \in B$ for some $B \in X$. By Lemma 13 of the full version <cit.> of <cit.> we can compute, in polynomial time, a distribution $\widehat{\prob}(s) \in Distr(S)$ such that $\widehat{\prob}(s)(C) = \prob_B(C)$ for all $C \in X$ and $\Vert \prob(s) - \widehat{\prob}(s) \Vert_1 = \sum_{C \in X} \vert \prob(s)(C) - \prob_B(C) \vert$. Define new LMCs $\mathcal{M}', \mathcal{N}'$ that differ from $\mathcal{M}$ and $\mathcal{N}$ only in the transition distribution functions, which are defined as $\prob'(s)(t) = \widehat{\prob}(s)(t)$. Then $\mathcal{M}'$ and $\mathcal{N}'$ are $\varepsilon$-perturbations of $\mathcal{M}$ and $\mathcal{N}$, respectively, as for any $s \in S$ it holds that

\begin{align*}
\Vert \prob(s) - \prob'(s') \Vert_1 = \Vert \prob(s) - \widehat{\prob}(s) \Vert_1 = \sum_{C \in X} \vert \prob(s)(C) - \prob_B(C) \vert \leq \varepsilon
\end{align*}

Furthermore, for $B, C \in X$, and $s,t \in B$,

\begin{align*}
\prob'(s')(C) = \widehat{\prob}(s)(C) = \prob_B(C) = \widehat{\prob}(t)(C) = \prob'(t')(C),
\end{align*}

so the equivalence induced by $X$ is a probabilistic bisimulation on $\mathcal{M}' \oplus \mathcal{N}'$. Because there is a dedicated block $B_{init} \in X$ with $s_{init}^{\mathcal{M}}, s_{init}^{\mathcal{N}} \in B_{init}$ we have $\mathcal{M}' \sim \mathcal{N}'$, and hence $\mathcal{M} \simeq_\varepsilon \mathcal{N}$.

[->,>=stealth',shorten >=1pt,auto, semithick]
every state = [text = black]

[state] (s) [] $s$;
(stemp) [below of = s, node distance = 1.2cm] …;
[state] (s1) [left of = stemp, node distance = 2cm] $s_1$;
[state] (sn) [right of = stemp, node distance = 2cm] $s_n$;
[state] (sa) [below of = s1, node distance = 1.5cm] $s_a$;
[state] (sb) [below of = sn, node distance = 1.5cm, fill = green] $s_b$;

(sinit) [left of = s, node distance = 1.2cm] ;

(sinit) edge (s)
(s) edge node [pos = 0.4, left, yshift = 0.1cm] $\frac{p_1}{T}$ (s1)
(s) edge node [pos = 0.4, right, yshift = 0.1cm] $\frac{p_n}{T}$ (sn)
(s1) edge node [left, pos = 0.4] $\frac{1}{2}$ (sa)
(s1) edge node [above, pos = 0.25] $\frac{1}{2}$ (sb)
(sn) edge node [right, pos = 0.4] $\frac{1}{2}$ (sb)
(sn) edge node[above, pos = 0.25] $\frac{1}{2}$ (sa)
(sa) edge [loop left] node $1$ (sa)
(sb) edge [loop right] node $1$ (sb)

[state] (t) [right of = s, node distance = 8cm] $t$;
(ttemp) [below of = t, node distance = 1.2cm] ;
[state] (t1) [left of = ttemp, node distance = 2cm] $t_y$;
[state] (t2) [right of = ttemp, node distance = 2cm] $t_n$;
[state] (ta) [below of = t1, node distance = 1.5cm] $t_a$;
[state] (tb) [below of = t2, node distance = 1.5cm, fill = green] $t_b$;

(tinit) [left of = t, node distance = 1.2cm] ;

(tinit) edge (t)
(t) edge node [pos = 0.4, left, yshift = 0.1cm] $\frac{N}{T}$ (t1)
(t) edge node [pos = 0.4, right, yshift = 0.1cm] $1 - \frac{N}{T}$ (t2)
(t1) edge node [left, pos = 0.4] $\frac{1}{2} - \varepsilon$ (ta)
(t1) edge node [above, pos = 0.25, xshift = 0.2cm] $\frac{1}{2}+ \varepsilon$ (tb)
(t2) edge node [right, pos = 0.4] $\frac{1}{2}- \varepsilon$ (tb)
(t2) edge node[above, pos = 0.25] $\frac{1}{2}  + \varepsilon$ (ta)
(ta) edge [loop left] node $1$ (ta)
(tb) edge [loop right] node $1$ (tb)

[right of = s, node distance = 0.75cm] $\{a\}$;
[left of = s1, node distance = 0.75cm] $\{a\}$;
[right of = sn, node distance = 0.75cm] $\{a\}$;
[right of = sa, node distance = 0.75cm, yshift = -0.2cm] $\{a\}$;
[left of = sb, node distance = 0.75cm, yshift = -0.2cm] $\{b\}$;

[right of = t, node distance = 0.75cm] $\{a\}$;
[left of = t1, node distance = 0.75cm] $\{a\}$;
[right of = t2, node distance = 0.75cm] $\{a\}$;
[right of = ta, node distance = 0.75cm, yshift = -0.2cm] $\{a\}$;
[left of = tb, node distance = 0.75cm, yshift = -0.2cm] $\{b\}$;

The LMCs $\mathcal{M}$ (left) and $\mathcal{N}$ (right) used in the -hardness proof of <Ref> (adapted from <cit.>).

We adapt the -hardness proof of <cit.> (cf. the full version <cit.> of <cit.>) that uses a reduction from the well-known SubsetSum-problem.
An instance $(P, N)$ of SubsetSum consists of a nonempty set $P = \{p_1, \dots, p_n\} \subseteq \mathbb{N}$ together with a value $N \in \mathbb{N}$, and poses the question if there is an index set $I \subseteq \{1,\dots,n\}$ such that the sum of all elements associated to the indices in $I$ equals $N$, i.e., such that $\sum_{i \in I} p_i = N$.

For a SubsetSum-instance $(P,N)$, let $T = \sum_{i=1, \dots, n} p_i$ be the sum of all elements in $P$, and consider the LMCs $\mathcal{M}$ and $\mathcal{N}$ depicted in <Ref>, which can be constructed from $(P,N)$ in polynomial time.

In <cit.>, it is shown that the SubsetSum-instance $(P, N)$ has a solution iff the LMC $\mathcal{M} \oplus \mathcal{N}$ has an $\varepsilon=\frac{1}{2T}$-quotient with $k=5$ states.
We use this result to show hardness of the question whether $\mathcal{M} \simeq_\varepsilon \mathcal{N}$.

Let $\varepsilon=\frac{1}{2T}$.
First observe that there is no $\varepsilon$-pertubation $\mathcal{N}'$ of $\mathcal{N}$ in which two distinct states are bisimilar. Now assume that $\mathcal{M} \simeq_\varepsilon \mathcal{N}$, and let $\mathcal{M}', \mathcal{N}'$ be the corresponding $\varepsilon$-perturbations with $\mathcal{M}' \sim \mathcal{N}'$.
For $i \in \{ 1, \dots, n \}$, the perturbed states satisfy either $s_i' \sim t_y'$ or $s_i' \sim t_n'$.
Due to $s' \sim t'$, it can be shown that $I = \{ i \in \{1, \dots, n \} \mid s_i' \sim t_y'\}$ is a solution to $(P,N)$.
On the other hand, any index set $I \subseteq \{1, \dots, n \}$ induces $\varepsilon$-perturbations such that $s_i' \sim t_y'$ if $i \in I$  and $s_i' \sim t_n'$ if $i \notin I$.
If $I$ is a solution to $(P,N)$, this yields $s' \sim t'$ and hence that $\mathcal{M}' \sim \mathcal{N}'$.
For more elaborative calculations, we refer to the proof of <cit.> that can be found in the full version <cit.> of <cit.>.

Now we turn our attention to the second claim of the theorem, i.e., to the -completeness of deciding $\mathcal{M} \sim_\varepsilon^* \mathcal{N}$. As, by definition, this holds iff $s_{init}^\mathcal{M} \sim_\varepsilon^* s_{init}^\mathcal{N}$, we instead show -completeness of the problem to decide if, for given states $s,t$ of a finite LMC $\mathcal{M}$ and a given $\varepsilon \in (0,1]$, $s \sim_\varepsilon^* t$.

The following nondeterministic, polynomial time guess-and-check algorithm computes a transitive $\varepsilon$-bisimulation $R$ on $\mathcal{M}$ with $(s,t) \in R$:

* Nondeterministically guess a partition $X$ of $S$ that only relates states with the same label and such that there is a block $B \in X$ with $s,t \in B$.

* Check if the induced equivalence $R$ with $X = \xfrac{S}{R}$ is an $\varepsilon$-bisimulation on $\mathcal{M}$ by deciding for all $(u,v) \in R$ if $\prob(u)(A(u,v)) \leq \prob(v)(A(u,v)) + \varepsilon$, where $A(u,v)$ is the union of all $R$-equivalence classes $B \in X$ with $\prob(u)(B) > \prob(v)(B)$.

Obviously, if $s \sim_\varepsilon^* t$, the algorithm finds a solution as the partition induced by $\sim_\varepsilon^*$ satisfies all the requirements for $X$ and $R$, respectively. On the other hand, if $X$ is a  solution of the algorithm, then the equivalence $R$ induced by $X$, i.e., such that $X = \xfrac{S}{R}$, is an $\varepsilon$-bisimulation as the condition on the probabilities posed in the second part of the algorithm ensures that, for all $(p,q) \in R$, the maximal difference in transition probabilities to subsets of $S$ is bounded from above by $\varepsilon$. As, by the first part of the algorithm $(s,t) \in X$, it follows that $s \sim_\varepsilon^* t$.

-hardness:  The hardness follows similar to that of the first claim, as in $\mathcal{M}$ and $\mathcal{N}$ as in <Ref> it holds that $s \sim_\varepsilon^* t$ iff $\mathcal{M} \simeq_\varepsilon \mathcal{N}$. This is due to the fact that all states $s_1, \dots, s_n$ in $\mathcal{M}$ have the same next-state distribution, and so the only way for them to become bisimilar to either $t_y$ or $t_n$ is to be perturbed uniformly, making the induced $\varepsilon$-bisimulation transitive.

Let $R$ be an equivalence. As described in the -membership proof of (ii) in <Ref>, we can check in polynomial time if $R$ is a (transitive) $\varepsilon$-bisimulation by deciding, for all $(s,t) \in R$, if $\prob(s)(A(s,t)) \leq \prob(t)(A(s,t)) + \varepsilon$, where $A(s,t)$ is the union of all $R$-equivalence classes $C$ with $\prob(s)(C) > \prob(t)(C)$. This shows (i). Furthermore, we know from <Ref> that, given $R$, we can decide in polynomial time if, for each $C \in \xfrac{S}{R}$ the centroid property (see <Ref>) holds, i.e., if $R$ satisfies (iii) of <Ref>. As it is decidable in polynomial time if $R$ is a transitive $\varepsilon$-bisimulation, deciding if $R$ is a $\varepsilon$-perturbed bisimulation is also possible in polynomial time.

§ PROOFS OF SECTION 5

Let $s \approx_\varepsilon^w t$ and $t \approx_\delta^w u$. Define the relation

\begin{align*}
R= \{(p,q) \mid \exists \, a \in S \text{ with} p \approx_\varepsilon^w a \text{ and} a \approx_\delta^w q\}.
\end{align*}

Then $(s,u) \in R$, so the claim follows if $R$ is a weak $(\varepsilon + \delta)$-bisimulation. To this end, we first observe that for all $(p,q) \in R$ there is a $b \in 2^{AP}$ with $b = l(p) = l(q)$ as both $\approx_\varepsilon^w$ and $\approx_\delta^w$ only relate states with the same label. Furthermore, for any $A \subseteq S$ it holds that

\begin{align*}
\mathrm{Pr}_p(b \Until A) \overset{p \approx_\varepsilon^w a}{\leq} \mathrm{Pr}_a(b \Until {\approx_\varepsilon^w}(A)) + \varepsilon \overset{a \approx_\delta^w q}{\leq} \mathrm{Pr}_q(b \Until {\approx_\delta^w}({\approx_\varepsilon^w}(A))) + \varepsilon + \delta
\end{align*}

and, as ${\approx_\delta^w}({\approx_\varepsilon^w}(A)) \subseteq R(A)$, this implies

\begin{align*}
\mathrm{Pr}_p(b \Until A) \leq \mathrm{Pr}_q(b \Until R(A)) + \varepsilon + \delta.
\end{align*}

Therefore, $R$ is indeed a weak $(\varepsilon + \delta)$-bisimulation, so $s \approx_{\varepsilon + \delta}^w u$.

For the first claim, consider the LMC $\mathcal{M}$ depicted on the left of <Ref>. The equivalence $R$ given via the equivalence classes $R = \{\{s, t\}, \{s_1, t_1\}, \{x\}, \{y\}\}$ is the largest branching $\varepsilon$-bisimulation on $\mathcal{M}$. In particular, $s \not \approx^b_\varepsilon s_1$ because otherwise we would have

\begin{align*}
\vert \mathrm{Pr}_{s}([s]_R \Until \{y\}) - \mathrm{Pr}_{s_1}([s_1]_R \Until \{y\}) \vert = \left \vert \frac{3}{8} - \frac{3}{4} \right \vert = \frac{3}{8} > \varepsilon,
\end{align*}

and $t \not \approx^b_\varepsilon t_1$ as, otherwise,

\begin{align*}
\vert \mathrm{Pr}_{t}([t]_R \Until \{x\}) - \mathrm{Pr}_{t_1}([t_1]_R \Until \{x\}) \vert = \left \vert \frac{5}{8} + \frac{5}{4} \varepsilon - \varepsilon^2 - \left (\frac{1}{4} + \varepsilon \right) \right \vert = \left \vert \frac{3}{8} + \frac{1}{4}  \varepsilon - \varepsilon^2 \right \vert > \varepsilon,
\end{align*}

where the inequality holds for all $\varepsilon < \frac{1}{8} \cdot (\sqrt{33} - 3) \approx 0.34307$, so in particular for all $\varepsilon < \frac{1}{4}$.

As $(s,t) \in R$, $s \approx^b_\varepsilon t$ in the LMC. However, $s \not \approx^w_\varepsilon t$ since

\begin{align*}
\mathrm{Pr}_t(L(t) \Until \{x\}) = \frac{5}{8} + \frac{5}{4} \varepsilon - \varepsilon^2 > \frac{5}{8} + \varepsilon = \mathrm{Pr}_s(L(s) \Until \{x\})
\end{align*}

for all $0 < \varepsilon < \frac{1}{4}$.
Hence, $s \approx^b_\varepsilon t$ does, in general, not imply $s \approx^w_\varepsilon t$.

Regarding the second claim, consider the LMC on the right of <Ref>. There, $\approx_\varepsilon^w$ is the symmetric and reflexive closure of $\{(s,t), (s,u), (s,v), (t,v), (t,w), (u,v), (v,w)\}$. In particular, $s \approx_\varepsilon^w t$ and $u \not \approx_\varepsilon^w w$.

However, $(s,t) \notin R$ for any branching $\varepsilon$-bisimulation $R$ on $\mathcal{M}$. To see this, first observe that, due to the state labeling, $x$ and $y$ cannot be related to any other state. It follows that $(u,w) \notin R$ for any such $R$ as, e.g., $\vert \mathrm{Pr}_u([u]_R \Until [x]_R) - \mathrm{Pr}_w([w]_R \Until [x]_R) \vert = 2 \varepsilon > \varepsilon$.

Now, assume that $(s,u) \in R$ for a branching $\varepsilon$-bisimulation $R$. Then

\begin{align*}
\vert \mathrm{Pr}_u([u]_R \Until [x]_R) - \mathrm{Pr}_s([s]_R \Until [x]_R) \vert = \begin{cases}
\frac{1}{4} + \frac{\varepsilon}{2}, &\text{if } (s, v) \notin R \\
\frac{\varepsilon}{2}, & \text{if } (s,v) \in R
\end{cases}
\end{align*}

and since in the first case we have $\frac{1}{4} + \frac{\varepsilon}{2} > \varepsilon$ because $0 < \varepsilon < \frac{1}{4}$, $(s,u) \in R$ only if $(s,v) \in R$. But then

\begin{align*}
\vert \mathrm{Pr}_t([t]_R \Until [x]_R) - \mathrm{Pr}_s([s]_R \Until [x]_R) \vert = \begin{cases}
\frac{1}{2} + \frac{\varepsilon}{2}, &\text{if } (t, v), (t, w) \notin R \\
\frac{1}{4} + \frac{\varepsilon}{2}, &\text{if } (t,v) \in R, (t, w) \notin R \\
\frac{1}{4} + \varepsilon, &\text{if } (t,v) \notin R, (t, w) \in R \\
\varepsilon, & \text{if } (t,v), (t,w) \in R
\end{cases}
\end{align*}

and as in the first three cases the differences are $> \varepsilon$, $(s,t) \in R$ requires $(t,v), (t,w) \in R$. However, in this case, the transitivity of $R$ implies $(u,w) \in R$, which is a contradiction.

Otherwise, if $(s,u) \notin R$, then $[s]_R \neq [u]_R$. Thus, $(s,t) \in R$ requires in particular that $\vert \mathrm{Pr}_s([s]_R \Until [u]_R) - \mathrm{Pr}_t([t]_R \Until [u]_R) \vert < \varepsilon$. If $(u,v) \notin R$ we have $\mathrm{Pr}_t([t]_R \Until [u]_R) = 0$, so the above difference becomes $\frac{1}{2} > \varepsilon$. Otherwise, if $(u,v) \in R$ then $\vert \mathrm{Pr}_s([s]_R \Until [u]_R) - \mathrm{Pr}_t([t]_R \Until [u]_R) \vert = \vert 1 - \mathrm{Pr}_t([t]_R \Until [u]_R) \vert$ and, since we know that $(u,w) \notin R$, we have  $\mathrm{Pr}_t([t]_R \Until [u]_R) \leq \frac{1}{2}$, so $\vert \mathrm{Pr}_s([s]_R \Until [u]_R) - \mathrm{Pr}_t([t]_R \Until [u]_R) \vert \geq \frac{1}{2} > \varepsilon$. Hence, $(s,t) \notin R$.

All in all, this shows that there can be no branching $\varepsilon$-bisimulation $R$ on $\mathcal{M}$ that contains $(s,t)$, and therefore it holds that $s \approx_\varepsilon^w t$ while $s \not \approx_\varepsilon^b t$.

Assume that $\prob(s)(L(s)) = 0$ for all $s \in S$. Then ${\sim_\varepsilon^*} = {\approx_\varepsilon^b}$ and ${\sim_\varepsilon} = {\approx_\varepsilon^w}$.

We start by showing the first claim, i.e. that, if there are no stutter steps ${\sim_\varepsilon^*} = {\approx_\varepsilon^b}$.

Let $R$ be a branching $\varepsilon$-bisimulation, and let $A \subseteq S$ be $R$-closed. Then, for all $(s,t) \in R$, we have $l(s) = l(t)$ and further

\begin{align}
\prob(s)(A) = \mathrm{Pr}_s([s]_R \Until A) \qquad \text{ and } \qquad \prob(t)(A) = \mathrm{Pr}_t([t]_R \Until A) \label{eq-1}
\end{align}

as no successor of either $s$ or $t$ is in $[s]_R = [t]_R$. But then it follows that

\begin{align*}
\vert \prob(s)(A) - \prob(t)(A) \vert = \vert \mathrm{Pr}_s([s]_R \Until A) - \mathrm{Pr}_t([t]_R \Until A) \vert \overset{(s,t) \in R}{\leq} \varepsilon.
\end{align*}

Hence, $R$ is an equivalence and an $\varepsilon$-APB, and thus a transitive $\varepsilon$-bisimulation by <Ref>.

Now let $R$ be a transitive $\varepsilon$-bisimulation, and let $(s,t) \in R$. Then $R$ is also an $\varepsilon$-APB by <Ref>, and it holds for all $R$-closed sets $A \subseteq S$ that

\begin{align*}
\vert \mathrm{Pr}_s([s]_R \Until A) - \mathrm{Pr}_t([t]_R \Until A) \vert = \vert \prob(s)(A) - \prob(t)(A) \vert \leq \varepsilon,
\end{align*}

where the equality follows from <Ref>. Therefore, $R$ is a branching $\varepsilon$-bisimulation.

Regarding the second claim, the result that ${\sim_\varepsilon} = {\approx_\varepsilon^w}$ follows directly from the fact that, if $\prob(s)(L(s)) = 0$ for all $s \in S$, we have $\mathrm{Pr}_s(L(s) \Until A) = \prob(s)(A)$ for every $A \subseteq S$, so it follows that for any relation $R$

\begin{align*}
\prob(s)(A) &\leq \prob(t)(R(A)) + \varepsilon \qquad \text{ iff } \qquad \mathrm{Pr}_s(L(s) \Until A) \leq \mathrm{Pr}_t(L(t) \Until R(A)) + \varepsilon.
\end{align*}

Note that we cannot replace ${\sim_\varepsilon}$ by $\sim_\varepsilon^*$. To see this consider, for example, a modification of the LMCs of <Ref>, where we replace the self-loops at $v, w, [v]$ and $[w]$ with the new transition probabilities $\prob(v)(w) = \prob(w)(v) = \prob([v])([w]) = \prob([w])([v]) = 1$. Then, in the direct sum $\mathcal{M}$ of these LMCs, we have $\prob(q)(L(q)) = 0$ for all states $q$. In particular, the largest weak $\varepsilon$-bisimulation $\approx_\varepsilon^w$ on $\mathcal{M}$ contains the pair $(t, [s])$, but we have already argued in <Ref> that there can be no transitive $\varepsilon$-bisimulation that relates these two states.

!0.07[->,>=stealth',shorten >=1pt,auto, semithick]
every state = [text = black]

[state] (s0) [fill = yellow] $s$;
[state] (s1) [right of = s0, node distance = 3cm, fill = yellow] $t$;
[state] (x) [right of = s1, fill = green, node distance = 3cm] $x$;
(sinit) [left of = s0, node distance = 1.2cm] ;

(sinit) edge (s)
(s0) edge node [above] $1$ (s1)
(s1) edge node [above] $1$ (x)
(x) edge [loop above] node [left, pos = 0.15] $1$ (x)

[above of = s0, node distance = 0.75cm] $\{a\}$;
[above of = s1, node distance = 0.75cm] $\{a\}$;
[right of = x, node distance = 0.75cm] $\{b\}$;

LMC in which $s \approx^b t$ and $s \approx^w t$, but neither $s \sim_\varepsilon t$ nor $s \equiv_\varepsilon t$ for any $\varepsilon \in [0,1)$.

We first deal with (ii). Consider the LMC depicted in <Ref>, and let $\varepsilon \in [0,1)$. The equivalence $R$ induced by the classes $C_1 = \{s, t\}$ and $C_2 = \{x\}$ is both a branching bisimulation and a weak bisimulation, as $x$ has a unique label and every state in $C_1$ reaches $C_2$ via a $C_1$-path, resp. an $\{a\}$-labeled path, with probability $1$. Hence, $s \approx^b t$ and $s \approx^w t$, so in particular also $s \approx^b_\varepsilon t$ and $s \approx^w_\varepsilon t$ for any $\varepsilon$.
However, $\vert \prob(s)(\{x\}) - \prob(t)(\{x\}) \vert = 1 > \varepsilon$ and $\prob(s)(\{t\}) = 1 > \varepsilon = \prob(t)(R(\{t\}) + \varepsilon$ for any $R$ that only relates states with the same label, so neither $s \equiv_\varepsilon t$ nor $s \sim_\varepsilon t$.

To show (i) we do a case distinction on $\approx_\varepsilon$.
If ${\approx_\varepsilon} = {\approx^b_\varepsilon}$, consider the LMC depicted on the left-hand side of <Ref> where $\varepsilon_1, \varepsilon_2 \in (0,1)$, $\varepsilon_1 \neq \varepsilon_2$, $\varepsilon_1 + \varepsilon_2 < 1$, and $\varepsilon = \vert \varepsilon_1 - \varepsilon_2 \vert$.
In this LMC, both $s \sim_\varepsilon t$ and $s \equiv_\varepsilon t$. However, $\mathrm{Pr}_s([s]_R \Until \{x_1\}) = \frac{\varepsilon_1}{\varepsilon_1 + \varepsilon_2} = \mathrm{Pr}_t([t]_R\Until \{x_2\})$ and $\mathrm{Pr}_s([s]_R \Until \{x_2\}) = \frac{\varepsilon_2}{\varepsilon_1 + \varepsilon_2} = \mathrm{Pr}_t([t]_R \Until \{x_1\})$ for any equivalence $R$ that only relates states with the same label. Hence,

\begin{align*}
\vert \mathrm{Pr}_s([s]_R \Until \{x_1\}) - \mathrm{Pr}_t([t]_R \Until \{x_1\}) \vert = \frac{\vert \varepsilon_1 - \varepsilon_2\vert }{\varepsilon_1 + \varepsilon_2} \overset{\varepsilon_1 + \varepsilon_2 < 1}{>} \vert \varepsilon_1 - \varepsilon_2 \vert = \varepsilon.
\end{align*}

Thus, there can be no branching bisimulation $R$ with $(s,t) \in R$, i.e., $s \not \approx_\varepsilon^b t$.

On the other hand, if ${\approx_\varepsilon} = {\approx^w_\varepsilon}$, consider the LMC on the right-hand side of <Ref> with $\varepsilon \in (0,1)$. In this LMC, $\sim_\varepsilon$ is given as the symmetric and reflexive closure of $\{(s, t), (s_1, t_1)\}$. As $\sim_\varepsilon$ is an equivalence, it is also an $\varepsilon$-APB by <Ref>, which yields $s \sim_\varepsilon t$ and $s \equiv_\varepsilon t$. Furthermore, $\mathrm{Pr}_t(L(t) \Until \{x\}) = 1$ while $\mathrm{Pr}_s(L(t) \Until \{x\}) = \frac{4(1-\varepsilon)^2}{(2-\varepsilon)^2}$. But then the inequality

\begin{align*}
\mathrm{Pr}_t(L(t) \Until \{x\}) = 1 > \mathrm{Pr}_s(L(s) \Until \{x\}) + \varepsilon = \frac{4(1-\varepsilon)^2}{(2-\varepsilon)^2} + \varepsilon
\end{align*}

holds for all $\varepsilon \in (0,1)$, so $s \not \approx_\varepsilon^w t$.

As $\sim_\varepsilon$ and $\equiv_\varepsilon$ above are transitive, the claim follows analogously for $\sim_\varepsilon^*$ and $\equiv_\varepsilon^*$.

We start with the first claim. Let $R$ be a branching $\varepsilon$-bisimulation and $(s,t) \in R$. We write $C$ for $[s]_R = [t]_R$ and get, since $S \setminus C$ is $R$-closed,

\begin{align*}
\vert \mathrm{Pr}_s(\Box C) - \mathrm{Pr}_t(\Box C) \vert &= \vert 1 - \mathrm{Pr}_s(C \Until (S \setminus C)) - (1 - \mathrm{Pr}_t(C \Until (S \setminus C))) \vert \\
&= \vert \mathrm{Pr}_s(C \Until (S \setminus C)) - \mathrm{Pr}_t(C \Until (S \setminus C)) \vert \overset{(s,t) \in R}{\leq} \varepsilon.
\end{align*}

Regarding the second claim assume that, for $C \in \xfrac{S}{R}$, there is an $s \in C$ with $\mathrm{Pr}_s(\Box C) > 0$. Then $\mathcal{M}$ must contain a non-empty bottom strongly connected component $B \subseteq S$ such that $B \subseteq C$. In particular, for any $t \in B$ it holds that $(s,t) \in R$ and $\mathrm{Pr}_t(\Box C) = 1$. Hence, $\mathrm{Pr}_t(C \Until (S \setminus C)) = 0$ and it follows that, since again $S \setminus C$ is $R$-closed,

\begin{align*}
\vert \mathrm{Pr}_s(C \Until (S \setminus C)) - \mathrm{Pr}_t(C \Until (S \setminus C)) \vert = \vert \mathrm{Pr}_s(C \Until (S \setminus C)) \vert \overset{(s,t) \in R}{\leq} \varepsilon.
\end{align*}

Thus, the probability to reach any state outside of $C$ from $s$ is at most $\varepsilon$, and we get

\begin{align*}
\mathrm{Pr}_s(\Box C) = 1 - \mathrm{Pr}_s(C \Until (S\setminus C)) \geq 1 - \varepsilon.
\end{align*}

Further, $B \neq \emptyset$ implies that there can be no $s \in C$ with $\mathrm{Pr}_s(\Box C) = 0$, as otherwise $1 = \mathrm{Pr}_s(C \Until (S \setminus C)) \leq \varepsilon$, which contradicts $\varepsilon \in [0,1)$. So, $\mathrm{Pr}_s(\Box C) \geq 1 - \varepsilon$ for all $s \in C$.

To show the third claim, let $R$ be a weak $\varepsilon$-bisimulation on $\mathcal{M}$ with $(s,t) \in R$, set $b = l(s) = l(t)$ and let $A  = S \setminus L(b)$. Then $A$ is $R$-closed, and $R(A) = A$. Therefore,

\begin{align*}
\vert \mathrm{Pr}_s(\Box b) - \mathrm{Pr}_t (\Box b) \vert = \vert 1 - \mathrm{Pr}_s(b \Until A) - (1 - \mathrm{Pr}_t(b \Until A)) \vert = \vert \mathrm{Pr}_s(b \Until A) - \mathrm{Pr}_t (b \Until A) \vert \leq \varepsilon,
\end{align*}

where the inequality follows form the fact that $(s,t) \in R$, and hence

\begin{align*}
\mathrm{Pr}_s(b \Until A) \leq \mathrm{Pr}_t(b \Until R(A)) + \varepsilon = \mathrm{Pr}_t(b \Until A) + \varepsilon
\end{align*}

as well as $\mathrm{Pr}_t(b \Until A) \leq \mathrm{Pr}_s(b \Until A) + \varepsilon$ by a similar argument.

We start with the direction from left to right. Let $R$ be a branching $\varepsilon$-bisimulation on $\mathcal{M}$. By construction, $\xfrac{S_R}{R^b} = \{C \cup \{s_C\} \mid C \in div_R\} \cup \{C \mid C \in \xfrac{S}{R} \setminus div_R\}$, so the only difference between $\xfrac{S}{R}$ and $\xfrac{S_R}{R^b}$ is the addition of the “divergence states” $s_C$ to all equivalence classes in $div_R$.

As $R^b$ is an equivalence that only relates states with the same label, the only thing left to show is that $R^b$ satisfies the condition on the probabilities. Following <Ref>, however, it is sufficient to show that $R^*$ is an $\varepsilon$-APB on $S_R$, i.e., that for all $(s,t) \in R^b$ and all $R^b$-closed sets $A' \subseteq S_R$ it holds that

\begin{align*}
\vert \prob_R(s)(A') - \prob_R(t)(A') \vert \leq \varepsilon.
\end{align*}

Let $C' \in \xfrac{S_R}{R^b}$ and define $C = C' \cap S$. Then $C = C'$ if $C \notin div_R$ and $C = C' \setminus \{s_C\}$ otherwise. Furthermore, let $(s,t) \in R^b$ with $s,t \in C'$, let $A'$ as above and define $A = A' \cap S$. Then, as $A'$ is $R^b$-closed, it is clear that $A$ must be $R$-closed.

If $s_C \notin C'$, i.e., if $C \notin div_R$, then $\mathrm{Pr}_u(\Box C) = 0$ for all $u \in C$. Thus, neither $s$ nor $t$ can have a transition in $\mathcal{M}_R$ to any state in $A' \setminus A$. But then

\begin{align}
\vert \prob_R(s)(A') - \prob_R(t)(A') \vert = \vert \prob_R(s)(A) - \prob_R(t)(A) \vert = \vert \mathrm{Pr}_s(C \Until A) - \mathrm{Pr}_t(C \Until A) \vert \leq \varepsilon \label{Proof: EpsilonBisimFromBranchingBisimOnTranfsformation - Eq1}
\end{align}

where the inequality follows from $(s,t) \in R$ and $A$ being $R$-closed.

Now assume that $C \in div_R$, i.e., that $s_C \in C'$. If both $s = t = s_C$ then the claim is clear. Thus, we can w.l.o.g. assume that $s \in C$ and $t \in C'$ (the case $s \in C'$ and $t \in C$ works analogously). We do a case distinction on the inclusion of $C'$ in $A'$.

Case 1: $C'$ is not contained in $A'$. Note that, since $A'$ is $R^b$-closed, the assumption implies that  $u \notin A'$ for all $u \in C'$, and in particular that $s_C \notin A'$.

Now, if both $s,t \in S$, then again $\prob_R(s)(A') = \prob_R(s)(A)$ and $\prob_R(t)(A') = \prob_R(t)(A)$, so the claim follows as in <Ref>. Otherwise, $s \in S$ and $t = s_C$. But then $\mathrm{Pr}_s(\Box C) \geq 1-\varepsilon$ by <Ref> and $\prob_R(t)(s_C) = 1$. Hence, as $s_C \notin A'$, $\prob_R(t)(A') = 0$ and

\begin{align*}
\prob_R(s)(A') = \prob_R(s)(A) = \mathrm{Pr}_s(C \Until A) \leq \mathrm{Pr}_s(C \Until (S \setminus C)) = 1 - \mathrm{Pr}_s(\Box C) \leq \varepsilon,
\end{align*}

which implies $\vert \prob_R(s)(A') - \prob_R(t)(A') \vert = \vert \prob_R(s)(A') \vert \leq \varepsilon$.

Case 2: $C' \subseteq A'$.
By <Ref> we know that, as $C \in div_R$, both $\mathrm{Pr}_s(\Box C) \geq 1 - \varepsilon$ and $\mathrm{Pr}_t(\Box C) \geq 1- \varepsilon$, so it follows that $\prob_R(s)(s_C) \geq 1- \varepsilon$ and $\prob_R(t)(s_C) \geq 1 - \varepsilon$. Hence, it holds for $u \in \{s,t\}$ that

\begin{align*}
1-\varepsilon \leq \prob_R(u)(s_C) \overset{s_C \in A'}{\leq} \prob_R(u)(A') \leq 1
\end{align*}

and thus $\prob_R(s)(A') - \prob_R(t)(A') \leq \varepsilon \text{ and } \prob_R(s)(A') - \prob_R(t)(A') \geq - \varepsilon$ or, equivalently, $\vert \prob_R(s)(A') - \prob_R(t)(A') \vert \leq \varepsilon$.
Therefore, $R^b$ is indeed an $\varepsilon$-APB.

For the other direction, assume that $R^b$ is a transitive $\varepsilon$-bisimulation on $\mathcal{M}_R$, and let $(s,t) \in R$. Furthermore, let $C = [s]_R = [t]_R$ and $B \subseteq S$ be $R$-closed with $B \neq C$. In particular, $B$ is a disjoint union of $R$ equivalence classes $C_1, \dots, C_n$, i.e., $B = C_1 \uplus \dots \uplus C_n$. We can w.l.o.g. assume that $C \neq C_i$ for all $i$, as otherwise $\vert \mathrm{Pr}_s(C \Until B) - \mathrm{Pr}_t(C \Until B) \vert = 0$. But then it follows that $\mathrm{Pr}_s(C \Until B) = \prob_R(s)(B)$ and $\mathrm{Pr}_t(C \Until B) = \prob_R(s)(B)$ as $B \subseteq S \setminus C$. Moreover, since $\prob_R(s)(s_D) = 0$ for all $D \in \xfrac{S}{R}$ with $D \neq C$, we have $\mathrm{Pr}_s(C \Until B) = \prob_R(s)(B_R)$ and $\mathrm{Pr}_t(C \Until B) = \prob_R(s)(B_R)$, where $B_R = B \cup \bigcup_{i = 1, \dots, n} s_{C_i}$ is $R^b$-closed. Since $R^b$ is a transitive $\varepsilon$-bisimulation on $\mathcal{M}_R$, which by <Ref> is the case iff it is an $\varepsilon$-APB, it follows that

\begin{align*}
\vert \mathrm{Pr}_s(C \Until B) - \mathrm{Pr}_t(C \Until B) \vert = \vert \prob_R(s)(B_R) - \prob_R(t)(B_R) \vert \leq \varepsilon
\end{align*}

which, as $R$ only relates states with the same label, proves the claim.

Let $\mathcal{M}$ be a finite LMC. Let $\mathcal{L} = \{b \in 2^{AP} \mid \exists \, s \in S \colon \mathrm{Pr}_s(\Box b) > 0\}$ be the set of “divergent labels” in $\mathcal{M}$.
As a first step, we construct from $\mathcal{M}$ a transformed LMC $\mathcal{M}^w = (S^w, \prob^w, s_{init}, l^w)$ with

* $S^w = S \cup \{s_b \mid b \in \mathcal{L}\}$ where the $s_b$ are fresh, pairwise different states

* $l^w(s) = l(s)$ if $s \in S$ and $l(s_b) = b$ for all $b \in \mathcal{L}$

* For $s, t \in S^w$ the distribution $\prob^w(s)$ is defined by

\begin{align*}
\prob^w(s)(t) = \begin{cases}
\mathrm{Pr}_s(L(s) \Until t), & \text{if } s,t \in S \text{ and} l(s) \neq l(t) \\
\mathrm{Pr}_s(\Box b), & \text{if } s\neq s_b, t = s_b, b = l(s) \text{ and} b \in \mathcal{L} \\
1, &\text{if} s = t = s_b \text{ for some } b \in \mathcal{L}\\
0, &\text{otherwise}
\end{cases}.
\end{align*}

Then $\mathcal{M}^w$ is again a finite LMC, and it is clear that the transformation from $\mathcal{M}$ to $\mathcal{M}^w$ preserves the probabilities of any given state $s \in S$ to reach a differently labeled state $t \in S$.

Now let, for a given $\varepsilon$-bisimulation $R$ on $\mathcal{M}$, $R^w$ be the finest reflexive and symmetric relation on $S^w$ with $R \subseteq R^w$ and $(s, s_b) \in R^w$ iff $s \in S$, $b = l(s)$ and $\mathrm{Pr}_s(\Box b) \geq 1 - \varepsilon$. As $R \subseteq R^w$, it follows for all $(s,t) \in R$ that $(s,t) \in R^w$.

The main part of the proof is now to show that $R^w$ is actually an $\varepsilon$-bisimulation on $\mathcal{M}^w$. To this end, let $(s,t) \in R^w$, $b = l(s) = l(t)$ and $A \subseteq S^w$. We must show that

\begin{align*}
\prob^w(s)(A) \leq \prob^w(t)(R^w(A)) + \varepsilon
\end{align*}

which we do by a case distinction on $s$ and $t$.

Case 1: $s = s_b$.
As $(s,t) = (s_b, t)\in R^w$ it holds by definition of $R^w$ that $\mathrm{Pr}_t(\Box b) \geq 1 - \varepsilon$. If now $s_b \in A$ then

\begin{align*}
\prob^w(s)(A) = \prob^w(s_b)(A) = 1 = 1 - \varepsilon + \varepsilon \leq \prob^w(t)(s_B) + \varepsilon \overset{s_b \in A}{\leq} \prob^w(t)(R^w(A)) + \varepsilon,
\end{align*}

while if $s_B \notin A$ we have

\begin{align*}
\prob^w(s)(A) = \prob^w(s_b)(A) = 0 \leq  \varepsilon \leq \prob^w(t)(R^w(A)) +  \varepsilon.
\end{align*}

Case 2: $t = s_b$.
Again, as $(s,t) = (s, s_b) \in R^w$ we have $\mathrm{Pr}_s(\Box b) \geq 1 - \varepsilon$. If $s_b \in R^w(A)$ then $\prob^w(t)(R^w(A)) \geq \prob^w(t)(s_b)  = 1$, so in particular

\begin{align*}
\prob^w(t)(R^w(A)) +  \varepsilon \geq 1 +  \varepsilon \geq 1 \geq \prob^w(s)(A).
\end{align*}

Otherwise, if $s_b \notin R^w(A)$, then also $s_b \notin A$. Thus, $\prob^w(t)(R^w(A)) = \prob^w(t)(A) = 0$ and hence

\begin{align*}
\prob^w(s)(A) \leq \prob^w(s)(S \setminus \{s_b\}) = 1 - \underbrace{\prob^w(s)(\{s_b\})}_{\geq 1 - \varepsilon \text{ as} (s, s_b) \in R^w} \leq \varepsilon \leq \prob^w(t)(R^w(A)) + \varepsilon.
\end{align*}

Case 3: $s,t \in S$ and $s_b \notin A$.
Be definition, neither $s$ nor $t$ can transition to any state $s_{b'} \in S^w \setminus S$ with $s_b \neq s_b'$. Thus, as $s_b \notin A$ by assumption, we get

\begin{align}
\prob^w(s)(A) = \prob^w(s)(A \cap S) = \sum_{\substack{q \in A \cap S \\ l(q) \neq b}} \mathrm{Pr}_s(b \Until \{q\})  = \mathrm{Pr}_s(b \Until ((A \cap S) \setminus L(s))). \label{M*-eq4}
\end{align}

Next, we show that
\begin{align}
R((A \cap S) \setminus L(t)) = R(A \cap S) \setminus L(t). \label{M*-eq5}
\end{align}

To see this, let $q \in R((A \cap S) \setminus L(t))$. Then there is a $p \in (A \cap S) \setminus L(t)$ with $(p, q) \in R$. As $R$ only relates states with the same label, $q \notin L(t)$. Thus, as $p \in (A \cap S)$, we get $q \in R(A \cap S)$, and as $q \notin L(t)$ it follows that $q \in R(A \cap S) \setminus L(t)$. Hence, $R((A \cap S) \setminus L(t)) \subseteq R(A \cap S) \setminus L(t)$. Regarding the other direction, let $q \in R(A \cap S) \setminus L(t)$. Then $l(q) \neq l(t)$ and there is a $p \in A \cap S$ with $(p, q) \in R$. Again, as $R$ only relates states with the same label, we have $l(p) = l(q)$ and, in particular, $l(p) \neq l(t)$, so $p \in (A \cap S) \setminus L(t)$. Therefore $q \in R((A \cap S) \setminus L(t))$, yielding $R(A \cap S) \setminus L(t) \subseteq R((A \cap S) \setminus L(t))$.

Now, as $R \subseteq R^w$ by definition and $A \subseteq A \cap S$, it follows that

\begin{align}
\prob^w(t)(R^w(A)) &\geq \prob^w(t)(R^w(A \cap S)) \geq \prob^w(t)(R(A \cap S))  = \sum_{\substack{q \in R(A \cap S) \\ L(q) \neq b}} \mathrm{Pr}_t(b \Until \{q\}) \nonumber \\&= \mathrm{Pr}_t(b \Until (R(A \cap S) \setminus L(t))) \overset{(\ref{M*-eq5})}{=} \mathrm{Pr}_t(b \Until R((A \cap S) \setminus L(t))). \label{M*-eq6}
\end{align}

By combining <Ref> we finally get

\begin{align*}
\prob^w(s)(A)& \overset{(\ref{M*-eq4})}{=} \mathrm{Pr}_s(b \Until ((A \cap S) \setminus L(s))) = \mathrm{Pr}_s(b \Until ((A \cap S) \setminus L(t))) \\
&\overset{(s,t) \in R}{\leq} \mathrm{Pr}_t(b \Until R((A \cap S) \setminus L(t))) + \varepsilon \overset{(\ref{M*-eq6})}{\leq} \prob^w(t)(R^w(A)) + \varepsilon.
\end{align*}

Case 4: $s,t \in S$ and $s_b \in A$.
First, we show that

\begin{align}
R^w(S^w \setminus R^w(A)) \subseteq S^w \setminus A. \label{M*-eq2}
\end{align}

To see this, let $q \in R^w(S^w \setminus R^w(A))$ which means that there is $p \in S^w \setminus R^w(A)$ with $(p,q) \in R^w$. But then we have $p \notin R^w(A)$ and hence it must hold that $q \notin A$. Thus, $q \in S^w \setminus A$.

Next, we argue that

\begin{align}
\prob^w(t)(S^w \setminus R^w(A)) \le \prob^w(s)(R^w(S^w \setminus R^w(A))) + \varepsilon. \label{M*-eq3}
\end{align}

This follows as in Case 3 by swapping the roles of states $s$ and $t$ and by using $A' = S^w \setminus R^w(A)$. In particular we have $s_b \notin A'$ since $s_b \in A \subseteq R^w(A)$.

Together, this yields

\begin{align*}
\prob^w(s)(A)
~&=~  1- \prob^w(s)(S^w \setminus A)\\
~&\overset{(\ref{M*-eq2})}{\leq}~ 1 - \prob^w(s)(R^w(S^w \setminus R^w(A)))\\
~&\overset{(\ref{M*-eq3})}{\leq}~ 1 - \Big(\prob^w(t)(S^w \setminus R^w(A)) - \varepsilon\Big)\\
~&=~ \Big(1 - \prob^w(t)(S^w \setminus R^w(A))\Big) + \varepsilon\\
~&=~ \prob^w(t)( R^w(A)) + \varepsilon.
\end{align*}

All in all, this shows that $R^w$ is actually an $\varepsilon$-bisimulation on $\mathcal{M}^w$. The claim now follows from <Ref> and the fact that the transformation from $\mathcal{M}$ to $\mathcal{M}^w$ preserves the probabilities for any state to reach a differently labeled state.

The following nondeterministic, polynomial time guess-and-check algorithm computes a branching $\varepsilon$-bisimulation $R$ on $\mathcal{M}$ with $(s,t) \in R$:

* Nondeterministically guess a partition $X$ of $S$ that only relates states with the same label and such that there is a block $B \in X$ with $s,t \in B$.

* Check if, for the induced equivalence $R$ with $X = \xfrac{S}{R}$, the equivalence $R^b$ from <Ref> is a transitive $\varepsilon$-bisimulation in the LMC $\mathcal{M}_R$ defined as in <Ref>, by deciding for all $(u,v) \in R^b$ if $\prob_R(u)(A(u,v)) \leq \prob_R(v)(A(u,v)) + \varepsilon$, where $A$ is the union of all $R^b$-equivalence classes $B \in X$ with $\prob_R(u)(B) > \prob_R(v)(B)$.

If $s \approx_\varepsilon^b t$, the algorithm finds a solution as the partition induced by $\approx_\varepsilon^b$ satisfies all the requirements for $X$ and $R$, respectively. On the other hand, any solution $X$ of the algorithm induces a transitive $\varepsilon$-bisimulation $R^b$ on $\mathcal{M}_R$ that contains $(s,t)$, so $s \approx_\varepsilon^b t$ follows by <Ref>.

-hardness:  The hardness follows by polynomial reduction from the decision problem if, for given $s,t \in S$ and $\varepsilon \in (0,1]$, $s \sim_\varepsilon^* t$, which is -complete by (the proof of) <Ref>. Given LMC $\mathcal{M}$, we construct in polynomial time a new LMC $\mathcal{N} = (S', \prob', s_{init}, l')$ as follows:

* $S' = S \uplus \{s' \mid s \in S\}$, i.e., we add for each state $s \in S$ a fresh duplicate $s'$ to $S'$

* $\prob'(s)(s') = 1$ and $\prob'(s')(t) = \prob(s)(t)$ for all $s,t$

* $l'(s) = l(s)$ if $s \in S$ and $l'(s') = (l(s))' = l(s)'$ if $s \in S' \setminus S$, i.e., we duplicate the labels in $\mathcal{M}$ and label each state in $S' \setminus S$ by a primed version of the label of its corresponding state in $S$.

By construction, there are no stutter steps in $\mathcal{N}$, i.e., $\prob'(q)(L(q)) = 0$ for all $q \in S'$. Using <Ref> it follows that ${\sim_\varepsilon^*} = {\approx_\varepsilon^b}$ in $\mathcal{N}$. Furthermore, there is a $1$-to-$1$-correspondence between $\varepsilon$-bisimulations in $\mathcal{M}$ and $\mathcal{N}$: from a given $\varepsilon$-bisimulation on $\mathcal{M}$ we can construct an $\varepsilon$-bisimulation on $\mathcal{N}$ by adding the primed version of all pairs of related states to the relation, and we obtain an $\varepsilon$-bisimulation on $\mathcal{M}$ from one on $\mathcal{N}$ by removing all of these pairs. All in all, this shows that $s \sim_\varepsilon^* t$ in $\mathcal{M}$ iff $s \sim_\varepsilon^* t$ in $\mathcal{N}$ iff $s \approx_\varepsilon^b t$ in $\mathcal{N}$.