  Let $X$ be the constrained random walk on ${\mathbb Z}_+^d$ representing the
queue lengths of a stable Jackson network and $x$ its initial position. Let
$\tau_n$ be the first time the sum of the components of $X$ equals $n$. $p_n
\doteq P_x(\tau_n < \tau_0)$ is a key performance measure for the queueing
system represented by $X$, stability implies $p_n\rightarrow 0$ exponentially.
Currently the only analytic method available to approximate $p_n$ is large
deviations analysis, which gives the exponential decay rate of $p_n$. Finer
results are available via rare event simulation. The present article develops a
new method to approximate $p_n$ and related expectations. The method has two
steps: 1) with an affine transformation, move the origin onto the exit boundary
of $\tau_n$, take limits to remove some of the constraints on the dynamics,
this yields a limit unstable constrained walk $Y$ 2) Construct a basis of
harmonic functions of $Y$ and use them to apply the classical superposition
principle of linear analysis. The basis functions are linear combinations of
$\log$-linear functions and come from solutions of "harmonic systems," which
are graphs whose vertices represent points on the "characteristic surface" of
$Y$, the edges between the vertices represent conjugacy relations between the
points, the loops represent membership in "the boundary characteristic
surfaces." Using our method we derive explicit, simple and almost exact
formulas for $P_x(\tau_n < \tau_0)$ for $d$-tandem queues, similar to the
product form formulas for the stationary distribution of $X$. The same method
allows us to approximate the Balayage operator mapping $f$ to $x \rightarrow
{\mathbb E}_x \left[ f(X_{\tau_n}) 1_{\{\tau_n < \tau_0\}} \right]$ for a range
of stable constrained random walks in $2$ dimensions. We indicate how the ideas
of the paper relate to more general processes and exit boundaries.
