The Harris recurrence of a stepping-out slice-sampling-within-Gibbs MCMC I want to use a multistage version of the MCMC here. That is, I want to use a Gibbs sampler to draw from a general joint distribution $p(x_1, x_2, x_3, \ldots)$ with a Gibbs step for each full conditional $p(x_i | \ldots x_{i - 1}, x_{i + 1}, \ldots)$, where I use Neal's (2003) stepping-out slice sampler to draw from each full conditional. (For each index $i$, sample $p(x_i | \ldots, x_{i - 1}, x_{i + 1}, \ldots)$ by first sampling auxiliary variable $u_i$ from $p(u_i | \ldots, x_{i - 1}, x_i, x_{i + 1}, \ldots)$ and then $x_i$ from $p(x_i | u_i, x_1, \ldots x_{i - 1}, x_{i + 1}, \ldots)$). 
Edit 4/7/16
I'm assuming a variant of slice sampling that does not necessarily sample from $p(x_i | u_i, x_1, \ldots x_{i - 1}, x_{i + 1}, \ldots)$ in a direct fashion. This full conditional is a uniform density on some "slice" $S$ of the $x_i$ component of the state space. Rather than calculate $S$ using an optimization, Neal proposed stepping-out and shrinkage procedures to draw an approximate sample of $x_i$ (Figures 3 and 5 of the 2003 paper).
End edit 4/7/16
Question: when is this Markov chain Harris recurrent? 
Harris recurrence would in some sense ensure eventual convergence to the stationary distribution from any starting state, and it would also give us useful ergodic theorems such as the MCMC analogue of the strong law of large numbers (Robert and Casella (2004)). From Roberts and Rosenthal (2006), we know Metropolis-within-Gibbs algorithm is Harris recurrent if the chain will eventually move in every coordinate direction with probability 1. I want the same thing to be true for slice-samping-within-Gibbs, but I cannot find a discussion of this specific point.
 A: I prove that the stepping-out and shrinkage procedure satisfies detailed balance, but this is of course not enough to show irreducibility or ergodicity.  And it's easy to construct examples in which there are regions with zero probability density in which a slice sampler that looks only at the part of the slice found by stepping out won't be ergodic. (Assuming the step size for stepping out is fixed - one could recover ergodicity by randomly picking the stepsize from an unbounded distribution.) 
If the conditional distributions to which univariate slice sampling with stepping out are applied have non-zero probability density everywhere within their range, then I think one could show that the sampler is ergodic, but there are probably cases where it's still not geometrically ergodic.
By the way, thanks for looking at the stepping out / shrinkage methods.  I've been puzzled why so many people just think of slice sampling as sampling independently from the slice, when the cases where this isn't possible were my main motivation for looking at slice sampling. And the use of slice sampling as a method that can be applied automatically to many problems is crucially dependent on being able to do something like stepping out and shrinkage.
A: Slice sampling is a special case of Gibbs sampling, to the point that in Monte Carlo Statistical Methods, we started our chapters on Gibbs sampling with a first chapter on slice sampling (Chapter 8).
To quote verbatim from the book (p.326):

Slice sampling relies upon the decomposition of the density $f(x)$ as
  $$ f(x) \propto \prod_{i=1}^k f_i(x)\,, $$ where the $f_i$'s are
  positive functions, but not necessarily densities. For instance, in
  a Bayesian framework with a flat prior, the $f_i(x)$ may be chosen as
  the individual likelihoods. This decomposition can then be associated
  with $k$ auxiliary variables $\omega_i$, rather than one as in the
  fundamental theorem, in the sense that each $f_i(x)$ can be written as
  an integral $$ f_i(x) = \int \mathbb{I}_{0\le \omega_i\le
  f_i(x)}\,d\omega_i\,, $$ and that $f$ is the marginal distribution of
  the joint distribution $$(x,\omega_1,\ldots,\omega_k) \sim
  p(x,\omega_1,\ldots,\omega_k) \propto \prod_{i=1}^k \mathbb{I}_{0\le
  \omega_i\le f_i(x)}\,. $$This particular demarginalization of $f$
  introduces a larger dimensionality to the problem and induces a
  generalization of the random walk of Section 8.1 which is to have
  uniform proposals one direction at a time.

Now, why is slice sampling a particular case of Gibbs sampling? Simply because in the "augmented space" made of the original variables $X_i$'s and of the auxiliary slice variables $U_i$'s, the steps are sheer simulations from the full conditionals:

1.Generate from $p(u_1 | u_2,\ldots, u_n,x_{1},  \ldots,x_n)=\mathbb{I}_{0\le u_1\le p(x_1|x_2,\ldots,x_n)}$
2.Generate from $p(x_1| u_1,\ldots, u_n,x_{2},  \ldots,x_n)=\mathbb{I}_{0\le u_1\le p(x_1|x_2,\ldots,x_n)}$ $$\vdots$$
  2i-1.Generate from $p(u_i | u_1,\ldots,
  u_{i-1},u_{i+1},\ldots,u_n,x_{1},  \ldots,x_n)=\mathbb{I}_{0\le u_i\le
  p(x_i|x_1,\ldots,x_{i-1},x_{i+1},\ldots,x_n)}$
2i.Generate from $p(x_i| u_1,\ldots, u_n,x_{1},\ldots,x_{i-1},x_{i+1},
  \ldots,x_n)=\mathbb{I}_{0\le u_i\le
  p(x_i|x_1,\ldots,x_{i-1},x_{i+1},\ldots,x_n)}$ $$\vdots$$
  2n-1.Generate from $p(u_n | u_1,\ldots, u_{n-1}x_{1}, 
  \ldots,x_n)=\mathbb{I}_{0\le u_n\le p(x_i|x_1,\ldots,x_{n-1})}$
2n.Generate from $p(x_n| u_1,\ldots, u_n,x_{1}, 
  \ldots,x_{n-a})=\mathbb{I}_{0\le u_n\le p(x_i|x_1,\ldots,x_{n-1})}$

So this is Gibbs sampling 101, only using uniform draws. Hence a Metropolis-Hastings move with probability one of accepting. If the support satisfies the connectivity property set by Besag (1994), the chain is irreducible and hence Harris positive recurrent.
