Mott glass and criticality in an <mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML"> <mml:mrow> <mml:mi>S</mml:mi> <mml:mo>=</mml:mo> <mml:mfrac> <mml:mn>1</mml:mn> <mml:mn>2</mml:mn> </mml:mfrac> </mml:mrow> </mml:math> bilayer Heisenberg mode

We employ the stochastic series expansion quantum Monte Carlo method to investigate the $S=1/2$ antiferromagnetic Heisenberg model on a bilayer square lattice with diluted interlayer couplings. Both regular and random dilution patterns are considered. In systems with regular dilution, tuning the interlayer interaction drives a quantum phase transition from a N\'eel-ordered phase to a quantum disordered phase, consistent with the $O(3)$ universality class. In contrast, random dilution gives rise to a two-step transition: from the N\'eel phase to an intermediate Mott glass (MG) phase, followed by a transition to the quantum disordered phase. Within the MG phase, the uniform magnetic susceptibility exhibits a stretched-exponential temperature dependence ${\ensuremath{\chi}}_{u}\ensuremath{\sim}exp(\ensuremath{-}b/{T}^{\ensuremath{\alpha}}), 0&lt;\ensuremath{\alpha}&lt;1$. At the N\'eel-to-glass transition, the quenched disorder modifies the critical exponents in a manner consistent with the Harris criterion. These findings provide new insights into disorder-driven quantum phase transitions and the emergence of glassy phases in diluted bilayer quantum magnets.