The sandpile model, polyominoes, and a $q,t$-Narayana polynomial

We give a polyomino characterisation of recurrent configurations of the sandpile model on the complete bipartite graph $K_{m,n}$ in which one designated vertex is the sink. We present a bijection from these recurrent configurations to decorated parallelogram polyominoes whose bounding box is a $m×n$ rectangle. Other combinatorial structures appear in special cases of this correspondence: for example bicomposition matrices (a matrix analogue of set partitions), and (2+2)-free posets. A canonical toppling process for recurrent configurations gives rise to a path within the associated parallelogram polyominoes. We define a collection of polynomials that we call $q,t$-Narayana polynomials, the generating functions of the bistatistic $(\mathsf{area ,parabounce} )$ on the set of parallelogram polyominoes, akin to Haglund's $(\mathsf{area ,hagbounce} )$ bistatistic on Dyck paths. In doing so, we have extended a bistatistic of Egge et al. to the set of parallelogram polyominoes. This is one answer to their question concerning extensions to other combinatorial objects. We conjecture the $q,t$-Narayana polynomials to be symmetric and discuss the proofs for numerous special cases. We also show a relationship between the $q,t$-Catalan polynomials and our bistatistic $(\mathsf{area ,parabounce}) $on a subset of parallelogram polyominoes.