  For the Poisson problem in two dimensions, posed on a domain partitioned into
axis-aligned rectangles with up to one hanging node per edge, we envision an
efficient error reduction step in an instance-optimal $hp$-adaptive finite
element method. Central to this is the problem: Which increase in local
polynomial degree ensures $p$-robust contraction of the error in energy norm?
We reduce this problem to a small number of saturation problems on the
reference square, and provide strong numerical evidence for their solution.
