[proofplan]
We work in the incidence algebra of a finite poset and use the defining identity that the Möbius function is the convolution inverse of the zeta function. Multiplying $\mu_P$ by the zeta function gives a convolution sum over the interval $[x,y]$, and comparing with the delta function gives the diagonal value and the zero-sum relation. Solving the zero-sum relation for the final term gives the displayed recursion. Finally, induction on interval cardinality shows that the recursion only uses the induced order structure of the interval, so $\mu_P(x,y)$ is invariant under interval isomorphism.
[/proofplan]
[step:Extract the recursion from the convolution inverse identity]
We use the incidence algebra construction for finite posets, as in [citetheorem:8100]. Let $I(P;\mathbb Z)$ denote the incidence algebra of $P$ over $\mathbb Z$: its elements are functions
\begin{align*}
f:\{(a,b)\in P\times P:a\le b\}\to \mathbb Z.
\end{align*}
For $f,g\in I(P;\mathbb Z)$, their convolution product $f*g\in I(P;\mathbb Z)$ is defined by
\begin{align*}
(f*g)(a,b)=\sum_{a\le c\le b}f(a,c)g(c,b).
\end{align*}
Define the zeta function $\zeta_P\in I(P;\mathbb Z)$ by $\zeta_P(a,b)=1$ for every $a\le b$. Define the delta function $\delta_P\in I(P;\mathbb Z)$ by requiring $\delta_P(a,a)=1$ for every $a\in P$ and $\delta_P(a,b)=0$ for every $a,b\in P$ with $a<b$.
By definition, the Möbius function $\mu_P\in I(P;\mathbb Z)$ is the convolution inverse of $\zeta_P$, so
\begin{align*}
\mu_P*\zeta_P=\delta_P.
\end{align*}
Evaluating this identity at a comparable pair $x\le y$ gives
\begin{align*}
\sum_{x\le z\le y}\mu_P(x,z)\zeta_P(z,y)=\delta_P(x,y).
\end{align*}
Since $\zeta_P(z,y)=1$ for every $z$ with $x\le z\le y$, this becomes
\begin{align*}
\sum_{x\le z\le y}\mu_P(x,z)=\delta_P(x,y).
\end{align*}
If $x=y$, the sum has the single term $z=x$, and $\delta_P(x,x)=1$, hence
\begin{align*}
\mu_P(x,x)=1.
\end{align*}
If $x<y$, then $\delta_P(x,y)=0$, hence
\begin{align*}
\sum_{x\le z\le y}\mu_P(x,z)=0.
\end{align*}
Separating the term $z=y$ from this finite sum gives
\begin{align*}
\sum_{x\le z<y}\mu_P(x,z)+\mu_P(x,y)=0.
\end{align*}
Therefore
\begin{align*}
\mu_P(x,y)=-\sum_{x\le z<y}\mu_P(x,z).
\end{align*}
[guided]
The incidence algebra packages the recursion into one multiplication identity. Its elements are functions
\begin{align*}
f:\{(a,b)\in P\times P:a\le b\}\to \mathbb Z,
\end{align*}
and the convolution product is
\begin{align*}
(f*g)(a,b)=\sum_{a\le c\le b}f(a,c)g(c,b).
\end{align*}
The finiteness of $P$ ensures that every interval $\{c\in P:a\le c\le b\}$ is finite, so the displayed sum is a finite integer sum.
Define $\zeta_P\in I(P;\mathbb Z)$ by $\zeta_P(a,b)=1$ for every comparable pair $a\le b$. Define $\delta_P\in I(P;\mathbb Z)$ by requiring $\delta_P(a,a)=1$ for every $a\in P$ and $\delta_P(a,b)=0$ for every $a,b\in P$ with $a<b$.
The Möbius function $\mu_P$ is defined as the convolution inverse of $\zeta_P$. Thus
\begin{align*}
\mu_P*\zeta_P=\delta_P.
\end{align*}
Now fix $x,y\in P$ with $x\le y$. Evaluating the convolution identity at $(x,y)$ gives
\begin{align*}
(\mu_P*\zeta_P)(x,y)=\sum_{x\le z\le y}\mu_P(x,z)\zeta_P(z,y).
\end{align*}
Because $z\le y$ for every index $z$ in this sum, the definition of $\zeta_P$ gives $\zeta_P(z,y)=1$. Therefore the identity $\mu_P*\zeta_P=\delta_P$ becomes
\begin{align*}
\sum_{x\le z\le y}\mu_P(x,z)=\delta_P(x,y).
\end{align*}
There are two cases. If $x=y$, then the only element $z$ satisfying $x\le z\le x$ is $z=x$, by antisymmetry of the partial order. Hence the sum reduces to
\begin{align*}
\mu_P(x,x)=\delta_P(x,x)=1.
\end{align*}
If $x<y$, then $\delta_P(x,y)=0$, so
\begin{align*}
\sum_{x\le z\le y}\mu_P(x,z)=0.
\end{align*}
The term with $z=y$ is exactly $\mu_P(x,y)$. Splitting it off from the finite sum gives
\begin{align*}
\sum_{x\le z<y}\mu_P(x,z)+\mu_P(x,y)=0.
\end{align*}
Solving for $\mu_P(x,y)$ yields
\begin{align*}
\mu_P(x,y)=-\sum_{x\le z<y}\mu_P(x,z).
\end{align*}
This is the desired interval recursion.
[/guided]
[/step]
[step:Prove interval dependence by induction on interval cardinality]
Let $P$ and $Q$ be finite posets, let $x\le y$ in $P$, and let $u\le v$ in $Q$. Suppose
\begin{align*}
\varphi:[x,y]\to [u,v]
\end{align*}
is a poset isomorphism. We prove $\mu_P(x,y)=\mu_Q(u,v)$ by induction on $N=|[x,y]|$.
If $N=1$, then $x=y$. Since $\varphi$ is an isomorphism, $u=v$, and the diagonal formula gives
\begin{align*}
\mu_P(x,y)=\mu_P(x,x)=1=\mu_Q(u,u)=\mu_Q(u,v).
\end{align*}
Assume $N>1$, so $x<y$, and assume the result holds for all isomorphic intervals of cardinality smaller than $N$. A poset isomorphism preserves the unique minimum and maximum of an interval, so $\varphi(x)=u$ and $\varphi(y)=v$. The recursion in $P$ gives
\begin{align*}
\mu_P(x,y)=-\sum_{x\le z<y}\mu_P(x,z).
\end{align*}
For each $z\in [x,y)$, the restriction of $\varphi$ defines a poset isomorphism
\begin{align*}
[x,z]\to [u,\varphi(z)].
\end{align*}
The interval $[x,z]$ has smaller cardinality than $[x,y]$ because $y\notin [x,z]$. By the induction hypothesis,
\begin{align*}
\mu_P(x,z)=\mu_Q(u,\varphi(z)).
\end{align*}
Since $\varphi$ bijects the set $\{z\in P:x\le z<y\}$ with the set $\{w\in Q:u\le w<v\}$, substitution in the finite sum gives
\begin{align*}
\mu_P(x,y)=-\sum_{u\le w<v}\mu_Q(u,w).
\end{align*}
Applying the same recursion in $Q$ gives
\begin{align*}
-\sum_{u\le w<v}\mu_Q(u,w)=\mu_Q(u,v).
\end{align*}
Thus $\mu_P(x,y)=\mu_Q(u,v)$, completing the induction. Therefore the value of $\mu_P(x,y)$ depends only on the isomorphism type of the interval $[x,y]$.
[/step]