On $p$-robust saturation on quadrangulations

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.