We numerically study level statistics of disordered interacting quantum many-body systems. A two-parameter plasma model which controls level repulsion exponent β and range h of interactions between eigenvalues is shown to reproduce accurately features of level statistics across the transition from ergodic to many-body localized phase. Analysis of higher order spacing ratios indicates that the considered β-h model accounts even for long range spectral correlations and allows to obtain a clear picture of the flow of level statistics across the transition. Comparing spectral form factors of β-h model and of a system in the ergodic-MBL crossover, we show that the range of effective interactions between eigenvalues h is related to the Thouless time which marks the onset of quantum chaotic behavior of the system. Analysis of level statistics of random quantum circuit which hosts chaotic and localized phases supports the claim that β-h model grasps universal features of level statistics in transition between ergodic and many-body localized phases also for systems breaking time-reversal invariance.Many-body localization (MBL) [1, 2] is a robust phenomenon of ergodicity breaking in disordered interacting quantum many-body systems [3][4][5][6]. It has attracted considerable attention over the last decade, notable findings include an emergent integrability of MBL phase due to the existence of local integrals of motion (LIOMs) [3, 7-10] and an associated unbounded logarithmic growth of the bipartite entanglement entropy after a quench from a separable state [11,12]. A wide regime of subdiffusive transport on the ergodic side of the transition was found [13][14][15]. Signatures of MBL have been observed experimentally in 1D [16,17] and in 2D system [18], however, the stability of MBL in 2D is still debated [19].Spectral statistics of ergodic systems with (without) time reversal invariance follow predictions of Gaussian orthogonal (unitary) ensemble (GOE, GUE, respectively) of random matrices [20,21] while eigenvalues of localized systems are uncorrelated resulting in Poisson statistics (PS). A ratio of consecutive spacings between energy levelswas proposed as a simple probe of the level statistics in [22] with n = 1 and employed in investigation of ergodicity breaking in various settings [23][24][25][26][27][28][29][30]. Higher order spacing ratios (n > 1), studied in [31][32][33][34], are valuable tools to assess properties of level statistics. In contrast to standard measures such as level spacing distribution or number variance they do not require the so called unfolding, i.e. the procedure of setting density of energy levels ρ(E) to unity which can lead to misleading results [35].Recently, an analytical understanding of an appearance of random matrix theory statistics in systems without a clear semiclassical limit have been developed in a periodically driven Ising models [36,37] or in random Floquet circuits [38]. Variants of such systems have been argued to undergo ergodic-MBL transition [39,40].In this work we consider ...