Extended symmetric quantum phase in a honeycomb Heisenberg model with sublattice-selective interactions

Motivated by two-dimensional bilayer systems we numerically study an anti-ferromagnetic spin-$1/2$ Heisenberg model on the honeycomb lattice with a nearest-neighbour exchange coupling $J_1$, and a next-nearest-neighbour exchange coupling $J_2^{\prime }$ for the B sublattice sites only. Using infinite Projected Entangled-Pair States (iPEPS), variational uniform matrix-product states on cylinders, and exact diagonalization, we find an extended symmetric regime $0.4\lesssim J_2^{\prime }/J_1\lesssim 0.6$, in which local observables show no magnetic, valence-bond, or chiral spin order. At $J_2^{\prime }/J_1=0.5$, the PEPS correlation length grows systematically with bond dimension, and an inverse-correlation-length extrapolation favors a vanishing limit. We also find that various observables scale algebraically with the finite bond-dimension-induced correlation length, which points to a gapless spin liquid ground state. We propose a ${\mathbb{Z}}_2$ Dirac spin liquid parton state, with Dirac points that are protected by translation, time-reversal and three-fold rotation symmetry, as a promising candidate state to explain the numerical results. We also discuss the possibility that the symmetric ground state is a featureless, short-range entangled state with a small gap.