The search for experimental signatures of the critical point (CP) of strongly interacting matter is one of the main objectives of the NA61/SHINE experiment at CERN SPS. In the course of the experiment, a beam momentum (13A – 150A GeV/$c$) and system size (p+p, p+Pb, Be+Be, Ar+Sc, Xe+La, Pb+Pb) scan is performed. Local proton density fluctuations in transverse momentum space represent an order parameter of the chiral phase transition and are expected to scale according to a universal power-law in the vicinity of the CP. They can be probed through an intermittency analysis of the proton second scaled factorial moments (SSFMs) in transverse momentum space. Previous such analyses [1] revealed power-law behavior in NA49 Si+Si collisions at 158A GeV/$c$, the fitted power-law exponent being consistent with the theoretically expected critical value, within errors. Probes of NA61/SHINE systems at the maximum SPS energy revealed no intermittency effect in Be+Be, whereas Ar+Sc analysis is inconclusive due to large uncertainties. The analysis has recently been extended to Pb+Pb collision data at lower energies.

We present a summary of the current status of NA61/SHINE intermittency analysis, and review novel techniques developed and employed to subtract non-critical background and estimate statistical and systematic uncertainties. Additionally, we use Monte Carlo simulations to assess the statistical significance of the observed intermittency effect.