Localization, fractality, and ergodicity in a monitored qubit

We study the statistical properties of a single two-level system (qubit) subject to repetitive ancilla-based measurements. This setup is a fundamental minimal model for exploring the intricate interplay between the unitary dynamics of the system and the nonunitary stochasticity introduced by quantum measurements, which is central to the phenomenon of measurement-induced phase transitions. We demonstrate that this “toy model” harbors remarkably rich dynamics, manifesting in the distribution function of the qubit’s quantum states in the long-time limit. We uncover a compelling analogy with the phenomenon of Anderson localization, albeit governed by distinct underlying mechanisms. Specifically, the state distribution function of the monitored qubit, parameterized by a single angle on the Bloch sphere, exhibits diverse types of behavior familiar from the theory of Anderson transitions, spanning from complete localization to almost uniform delocalization, with fractality occurring between the two limits. By combining analytical solutions for various special cases with two complementary numerical approaches, we achieve a comprehensive understanding of the structure delineating the “phase diagram” of the model. We categorize and quantify the emergent regimes and identify two distinct phases of the monitored qubit: ergodic and nonergodic. Furthermore, we identify a genuinely localized phase within the nonergodic phase, where the state distribution functions consist of delta peaks, as opposed to the delocalized phase characterized by extended distributions. Identification of these phases and demonstration of transitions between them in a monitored qubit are our main findings.


I. INTRODUCTION
Statistical properties of random systems constitute one of the most outstanding topics of theoretical physics, which keeps attracting a lot of attention, in spite of a long history.Such an interest results, in particular, from the beauty and fundamental importance of the physics of phase transitions.One seminal example is the theory of the Anderson transition between localized and delocalized phases [1,2].Its archetypal setup involves noninteracting particles in a random potential.When strong disorder is introduced in the system the wave functions of particles become spatially localized.Conversely, for weak disorder, the wave functions may spread uniformly throughout the system, exhibiting extended behavior.The transition between the two phases is characterized by critical fluctuations of wave functions, leading to the formation of a multifractal spectrum.
Recently, a similar type of phase transition governed by randomness has been discovered in hybrid random quantum circuits, subject to both unitary evolution and local measurements (monitoring) [3][4][5][6][7].In monitored systems, randomness is introduced by the quantum probabilistic nature of measurement outcomes rather than by the disordered potential.The stochastic nonunitary evolution of the system emerges from the inherently projective nature of measurements.The influence of measurements on the dynamics of quantum systems has attracted considerable attention, largely because of recent advances in the field of quantum information processing.Regardless of the particular quantum hardware, the nuisance of environmental noise is a formidable challenge [8][9][10].In this context, it is crucial that measurements can act both as a tool to monitor the properties of a quantum system and as a source of controllable disturbances.
A direct link between measurement-induced entanglement transitions and the phenomenon of Anderson localization was recently established in Refs.[47,48,[50][51][52] by deriving nonlinear sigma-models for monitored free fermions.These field theories share significant similarities [50] with those used in the context of Anderson localization.Particularly noteworthy is the striking resemblance between the predicted entanglement transition for free fermions in spatial dimensions larger than one (D > 1) and the Anderson transition observed in disordered systems of dimension D + 1 [51,52].In the field of measurement-induced dynamics, the entanglement transition can be considered as a "metal-insulator transition" for quantum information.In this analogy [51], mutual information serves as a counterpart to dimensionless conductance, shedding light on the evolving properties of the quantum system under the influence of measurements.
Importantly, both the Anderson transition and the measurement-induced entanglement transition are phenomena that take place in infinite systems (thermodynamic limit).The main question we are addressing here pertains to whether a macroscopic spatial size or a macroscopically large number of degrees of freedom is a necessary condition for the observation of complex behaviors driven by randomness, thereby leading to transitions between distinct phases.More specifically, our primary focus lies in ascertaining whether a transition akin to the Anderson localization-delocalization transition can be driven by repeated measurements in a microscopic quantum system (cf.Ref. [77]).
In the present paper [78], we come across such complexity even in a single monitored two-level system-a qubit (e.g., a Loss-DiVincenzo spin qubit [79]), i.e., in the smallest possible quantum system with nontrivial dynamics.The thermodynamic limit that is necessary for true phase transitions can be reached here in the limit of infinite observation time and, correspondingly, an infinite number of measurements.Our setup comprises two twolevel systems interacting with each other, one of them representing the qubit, while the other serves as the detector, see Sec.II A. This is arguably the simplest model implementing "ancilla-based" generalized measurements.Experimentally, monitoring and steering of qubit quantum trajectories has become feasible [80][81][82][83][84].
We demonstrate that this simple measurement model possesses extremely rich dynamics induced by the interplay of unitary evolution and stroboscopic measurements.This shows up, in particular, in the statistical properties of the distributions of quantum states of the monitored qubit, which correspond to different quantum trajectories.We find that, in the long-time limit, the probability of finding a qubit in a given state right after the measurement can be described by the time-averaged distribution function, W (θ), of a single angle variable parameterizing the state, see definitions in Sec.III A. This is because the measurement protocol attracts the quantum trajectories to a specific one-dimensional manifold on the Bloch sphere of the qubit.One relevant observable, where the state statistics manifests itself, is the expectation value of the occupation of the monitored level (site) of the qubit.Upon varying the strength of coupling between the system and the detector, as well as the measurement period, W (θ) exhibits a remarkably diverse and nontrivial behavior.It is worth mentioning that recent work [77] predicted a cascade of dynamical transitions in a similar model approaching the quantum Zeno limit of frequent measurements (continuous monitoring).Here, we are interested in the stationary phase diagram of the model in the full parameter space, which includes stroboscopic measurements with a finite period.
One could argue that such complexity is kind of similar to that arising during the evolution of chaotic systems [104].However, there is a crucial difference between the chaotic and monitored systems: The complexity in the Hamiltonian or dissipative chaotic systems results from the intrinsic nonlinearity of the evolution or nonintegrability induced by the system boundaries; an additional source of randomness is not needed there.On the con-trary, the complexity of the monitored system is unavoidably related to the quantum-mechanical uncertainty in the measurement outcomes (Born's rule) and disappears if one applies a deterministic post-selection procedure.Technically, the nonlinearity in our problem emerges at the level of a functional master equation governing the distribution functions for our stochastic maps.
The importance of measurement-induced randomness points towards properties of disordered systems and the Anderson localization transition rather than to chaotic systems in the present context.Having in mind the above-mentioned link between the measurement-induced transitions and Anderson localization in large systems, it is very natural to anticipate that, in the limit of the infinite observation time, statistical properties of W (θ) should be reminiscent of those of the wave functions in disordered systems.Indeed, we observe that W (θ) may look localized (sharply peaked, see a representative example in the upper panel of Fig. 1 and the first two panels of Fig. 10), delocalized (see a representative example in the middle panel of Fig. 1 and last three panels of Fig. 10) and even fractal (self-similar, see a representative example in the lower panel of Fig. 1 and Fig. 12).
However, this visual similarity between the wavefunction statistics and the statistics of the states of the monitored qubit is not conclusive.The classification of various types of behavior of the monitored qubit and identification of distinct phases is only possible by using a combination of several mutually complementary quantitative indicators.Inspired by the striking parallels to the statistics of wave functions and the local density of states in the theory of Anderson localization, we employ here the standard indicators of the Anderson transition: (A) the participation ratio and its scaling, Sec.V A; (B) the support of W (θ), Sec.V B; (C) the typical value of W (θ), Sec.V C; (D) the Hausdorff dimension of the curve W (θ), Sec.V D. In addition, we characterize the stochastic evolution resulting in the steady-state distribution W (θ) by the (E) ergodicity marker of the corresponding Markov process, Sec.V E, which is common in the studies of chaotic dynamical systems.
By combining the analytical solutions for several special cases, which are related to various commensurability conditions, with numerical approaches, we understand the overall structure of the phase diagram of the monitored qubit, see Fig. 15.The extensive numerical analysis of the indicators has allowed us to categorize and quantify the rich variety of regimes exemplified in Figs. 1 and  10.Using indicator (E), we have identified two distinct phases of the monitored qubit: ergodic and nonergodic.Furthermore, we have found a genuinely localized phase within the nonergodic phase, where the state distribution functions consist of delta peaks, as opposed to the delocalized phase characterized by W (θ) with nonzero support.Identification of distinct phases and demonstration of transitions between them in a monitored qubit are our main findings.
Our analysis of the qubit toy model signifies that measurement-induced transitions are characterized by an order parameter that is not a conventional scalar quantity but rather a distribution function.This distinguishing feature again resonates with the peculiarities of the Anderson transition, where the order parameter is represented by a distribution function of the local density of states [2].Although the "single-spin" system we study cannot demonstrate the genuine entanglement transition, the discovered complexity of the model manifested in the probability distribution W (θ) suggests that in macroscopic monitored systems additional "hidden" transitions are feasible, which could possibly be observed in various distribution functions rather than average quantities.The paper is organized as follows.In Sec.II, we introduce the model and measurement protocol.The evolution of the qubit state on the Bloch sphere is discussed in terms of measurement operators.We further define the angle distribution function (ADF) W (θ) and explain its manifestation in various averages (stationary solution of a Master equation and time average of a single typical quantum trajectory).In Sec.III, the measurement operators of our model are derived.Based on these operators, we explain why the time evolution can be described asymptotically in terms of a single angle variable.We obtain the ADF analytically for several special parameter choices in Sec.IV.The quantifiers of localization, fractality, and ergodicity are introduced in Sec.V, where we illustrate and explain their properties using generic ADFs obtained from numerical simulation.Section VI is dedicated to the systematic numerical calculation of these indicators in the parameter plane of the model.Diagrams of the indicators in the parameter plane are shown and their structure is explained in terms of the analytically understandable special cases.The emergence of different phases in the parameter plane is demonstrated.We discuss the relation of our results to previous works on similar topics, implications of our findings to experiment, as well as further directions, in Sec.VII and conclude in Sec.VIII.Technical details are relegated to Appendices.

A. Model and measurement protocol
The model we study is illustrated in Fig. 2. The entire setup (green, solid box) consists of a system (orange, dashed box) and a detector (blue, dotted box).The twolevel system is represented by two tunnel-coupled sites s 1 and s 2 , which can be occupied by a single spinless electron (one can equivalently consider any other realization of a qubit as the system).One of the two sites (labeled by 1) is coupled to the two-state detector that can be a single spin-1/2 (or, equivalently, another pair of sites occupied by a single electron).The particular form of the coupling between the system and the detector is chosen to conserve the number of electrons in the system; in our case, the detector monitors the occupation of site s 1 .2. Scheme of the model: particles can occupy two sites, s1,2, which are connected by tunneling with the strength γ.The site s1 interacts with a two-level system (detector D), with the coupling strength being M .In the text we refer to the the chain of sites s1,2 (orange, dashed box) as "the system", the two level system D (blue, dotted box) is called "the detector", and the entirety of system and detector (green, solid box) is "the setup".

↑ ↓
The specific model to be analyzed is described by the following Hamiltonian: Here, Ĥs is the system Hamiltonian, Ĥsd couples the system and the detector, â1,2 (â † 1,2 ) are the fermionic annihilation (creation) operators on sites s 1 and s 2 , the Pauli matrix τx acts in the detector space, and γ and M are the tunneling and interaction constants, respectively.The occupation of site s 1 facilitates transitions between the two levels of the detector.If the detector is realized, e.g., with two auxiliary sites (cf.Ref. [69]), the chosen form of Ĥsd would correspond to the hopping between the auxiliary sites modulated by the density on s 1 of the system.We assume that the system sites are unbiased, i.e., they have the same on-site energies.The value of the on-site energy determines the origin and we have set it to zero.Besides, we disregard an unimportant phase of γ, which can be gauged out.The dimension of the Hilbert space of the setup (two-level system plus two-state detector) is four.The orthonormal basis vectors of the system can be chosen as follows: where the state |1, 0⟩ (|0, 1⟩) shows whether the first (second) site is occupied.
Within our protocol, the detector is initiated in a given state |−⟩ at time t 0 = 0, such that the initial state of the entire setup (measured system plus detector) reads as Initial setup's state: where ψ 0 is the initial state of the system.Measurements are performed stroboscopically at time instants t j = jT, j = 1, 2, 3, . . .(T is the measurement period), when the detector is projected on one of its two states in the basis of its initial state.The probabilities of the detector readouts are given by the standard Born's rule.The projection to state |−⟩ (with probability P − ) will be called a "no-click" event; the projection to the flipped state |+⟩ (with probability P + ) will be referred to as a "click" event.The system's state ψ(t j ) after the projection of the detector depends on the outcome ± ("postmeasurement state"): After each projection, the detector is reinitialized in the |−⟩ state at post-measurement times t > j ≡ t j + 0: This corresponds to "resetting" in the problem, which removes memory effects [70] that would appear if the joint system-detector evolution between times t > j and t j+1 started with the fiducial state of the detector after projection at t j .
The dynamics of the entire setup between t 0 and t 1 , as well as between any two successive measurements, can be described by the unitary evolution operator, ÛT = exp −i ĤT .
This unitary evolution entangles the system and the detector, so that, generically, the state of the setup is not a product state for all times except for t j (and t > j ).However, after each measurement on the detector, the total wave function collapses into a separable product of the system and detector states.Importantly, the unitary evolution of the setup is governed by the full Hamiltonian (1) at all times between the successive measurements, combining both the system's own tunneling dynamics and the one induced by the system-detector coupling.This should be contrasted with other models known in the literature (see, e.g., Ref. [77]), where the periods of "measurement dynamics" (governed only by Ĥsd , with Ĥs switched off) follow the periods of free (decoupled from the detector, Ĥsd switched off) evolution of the system.In such models, the measurement Kraus operators are independent of the system Hamiltonian (would not involve γ for the system considered).Our model is thus more realistic for studying the interplay and competition of the measurement-induced dynamics with the system dynamics.It is this type of competition that results in the measurement-induced entanglement transitions, see the discussion of a related toy model in Ref. [41].

B. Electron states and post-measurement mapping
In what follows, we consider the system states at times t j , ψ(t j ), which are normalized two-component spinors that can be written as After omitting the overall phase of ψ, the corresponding amplitudes can be parameterized by two angle variables, Now, we can map the Hilbert space of the system onto the space of unit vectors, starting at the origin and pointing to points on the Bloch sphere in three dimensions, see Fig. 3 for examples of states on the Bloch sphere.Note that the resulting parameterization of the Bloch sphere differs from the one based on conventional Euler angles, where the azimuthal angle varies from 0 to π.The coefficients (or angles) parameterizing the state ψ take some values at the initial time, α(t 0 ) = α 0 , β(t 0 ) = β 0 [or, equivalently, θ(t 0 ) = θ 0 , φ(t 0 ) = φ 0 ] and change during the quantum evolution.Since between the measurement events, the system is entangled with the detectors, these coefficients are only defined at postmeasurement times t > j .Such evolution at discretized time instants is described by the mapping that depends on random measurement outcomes: where After the j-th measurement, there are 2 j different endpoints (in principle, some endpoints can describe the same state).Each generated state (including all endpoints) is labeled by a bitstring showing the sequence of measurement outcomes along the quantum trajectory leading to the state.
The probabilities of the measurement readout at t > j generically depend on the state of the system at time t > j−1 , i.e., on α j−1 and β j−1 .To reflect this property, we introduce notations with µ = ± for click-and no-click outcomes.The mapping ( 8) can be formulated in the matrix form: Singling out the normalization factor with P (µ) j in this matrix equation renders the matrices M± (Kraus operators) independent of the system's state.Note that α and β on the right-hand side of Eq. (11) are not marked by the upper index (µ) because the outcome at t > j−1 (whether it is "click" or "no-click") is not important for obtaining the state at t > j .This is a consequence of the ancilla's resetting after its projection, which also implies that the probabilities P (µ) j do not depend explicitly on the previous measurement outcome (hence, only one outcome label µ).
If we start with the state determined by α 0 and β 0 and explore all quantum trajectories of length j, we come across the tree-like graph with the branching number 2 and 2 j endpoints [105], see Fig. 4.Each endpoint of the graph represents a quantum state of the system obtained after acting on {α 0 , β 0 } T by a random product of 2 × 2 matrices M+ and M− .It corresponds to a given sequence of click and no-click outcomes -a "quantum trajectory" parameterized by the kth "bit-string" T k ≡ {µ 1 , µ 2 , . . ., µ j }, where k = 1, ..., 2 j labels one of the endpoints of the tree, see Fig. 4.Each quantum trajectory T k is weighted by the total Born's probability, Following Eq. ( 11), the normalization factor of the resulting state is given by the square root of the same total probability (13).
It is worth noting that random products of matrices [106] appear in various contexts ranging from Anderson localization in one-dimensional arrays of impurities [107][108][109][110][111] to biological evolution and computer science (see, e.g., Refs.[112][113][114][115][116] and references therein).However, in most of the applications, the focus is the Lyapunov exponent characterizing such random products.In particular, the maximum Lyapunov exponent determines the localization length for low-dimensional disordered systems described by the transfer-matrix techniques [117].In our case, the renormalization factor introduced by the total Born's probability does not allow the state vector to change its length, so that we are not interested in the Lyapunov exponent of the random matrix M. Instead, we are focusing on the statistics of the resulting system's states as parameterized by the angles θ j and φ j .Furthermore, the probabilities of applying the two matrices to the state-vector are generically state-dependent in our case, in contrast to most works on random products of matrices.Below, we will analyze the statistics of these Bloch-vector angles in the limit of infinite time, j → ∞.

C. Angle distribution function and long-time limit
Owing to the probabilistic nature of measurement outcomes, the model requires a statistical description.We characterize the system by the statistics of its pure states at times t > j , i.e., by the statistics of both Bloch-sphere angles.In most of the setups studied in the present paper, the statistics of quantum trajectories in the longtime limit turn out to be fully described by the statistics of angle θ.In particular, in the case of finite hopping, γ ̸ = 0, we encounter an attraction of almost all quantum trajectories to the circle formed at the intersection of the Bloch sphere with the Y Z-plane, see Sec.III B below.This plane is parameterized by φ = π/2, so that α = cos(θ/2) and β = i sin(θ/2) on this manifold, see Fig. 5.In what follows, we refer to it as the Grand Circle (GC).Importantly, the GC is an invariant manifold of evolution: starting from the state belonging to the GC, the system never escapes from it.
This important property of generic evolution in our setup allows us to use the distribution function of the single angle θ to describe the statistics of the system's states in the long-time limit.In view of this special role of the GC, we will, for simplicity, focus on quantum trajectories starting on the GC, without going into detail of the transient behavior of quantum trajectories approaching the GC.We will discuss those special fine-tuned cases, for which the evolution of the system does not have the GC as an attractor, separately.
The matrix mapping (11) for the GC can be equivalently rewritten in terms of the discrete mapping of angle θ j−1 → θ j in the course of quantum evolution: with probability P − (θ j−1 ), Θ + (θ j−1 ) with probability P + (θ j−1 ).( 14) Here, the functions Θ ± (θ) are derived from the form of the matrices M± .The probabilities of click and noclick outcomes, P ± , are functions of θ corresponding to Eq. ( 10) via Eq.( 6) with φ = π/2.The explicit forms of these functions are presented in Appendix C. It is also useful to introduce the functions for the inverse ("retrospective") mapping θ j → θ j−1 , describing the angles from which a given angle θ is obtained by application of M− or M+ (when these are invertible matrices): The following relations clearly hold for invertible mappings: To characterize all possible endpoints of a quantum trajectory after j steps, we consider the probability distribution of these states.Specifically, we take the angle θ k for each endpoint (labeled by k = 1, . . ., 2 j ) of the above-described tree-like evolution graph, thus accounting for all possible quantum trajectories T k of length j, and sum the corresponding δ-functions with their Born weights (13): Here, θ 0 in the argument marks the starting point of the graph, i.e., the θ-angle of the initial state on the GC.Both θ k (for each k) and P Born generically depend on θ 0 (exceptions will be discussed separately).The superscript (ε) on both sides of Eq. ( 17) denotes a regularization of the δ-function parametrized by ε.This procedure is important, because the distribution of pure states for any number of measurement steps j and also for j → ∞ is a set of delta-peaks.Taking the limit ε → 0 at the end of the calculation allows one to define a continuous distribution.Evidently, W j (θ|θ 0 ) is normalized: The structure of the angle distribution is similar to the pattern of the local density of states in disordered systems (see also Sec.V C below for more details).The local density of states captures the weight of states of a given energy at a given position, while the angle distribution captures the weight of quantum trajectories at a given angle, thus quantifying the probability to find the system in the state parameterized by this angle.In a finite closed disordered system, the local density of states as a function of energy is represented as a series of exact deltapeaks located at the eigenenergies of the system, similarly to Eq. (17).The weights of these peaks are determined by the amplitudes of the corresponding eigenfunctions at the point where the density of states is calculated [the counterpart of P Born (T k ) in Eq. ( 17)].With increasing size of the disordered system, the density of energy levels increases (similarly, the density of the endpoints θ k of a quantum trajectory generically increases with increasing the number of steps j).However, the distance between levels with high weight may remain finite as a result of the spatial localization of eigenstates.The inhomogeneity of a disordered system leads to a nontrivial structure of the density of states, while the complex structure of the angle distribution is determined by the M± maps and the corresponding Born probabilities.
The distribution function of the local density of states plays a role of the functional order parameter in the theory of Anderson transitions [2].The transition to the metallic state in the thermodynamic limit is characterized by changing the character of the locally defined spectrum from quasi-discrete to continuous.Importantly, in order to obtain the continuous spectrum from an infinite set of weighted δ-functions, one introduces infinitesimal broadening of the states, and sends the width of δ-functions to zero at the end of the calculation, after taking the thermodynamic limit.We adopt a similar procedure in our case; this once again underscores the relation of the present problem to Anderson localization.
We do not know a priori whether or not the distribution of states keeps a memory of the initial state in the limit j → ∞.As we are going to show, both situations (dependence on and independence of the initial angle) are possible, depending on the parameters of the setup.We introduce a time-independent steady-state distribu-tion in the long-time limit (when such a limiting function exists) as where It is worth emphasizing that the "thermodynamic limit" j → ∞ is taken in Eq. ( 21) first.Below, we will encounter situations (e.g., period-2 trajectories, Sec.IV C) when the limit j → ∞ does not exist, implying the absence of a unique steady state.In such a case, it is convenient to introduce the "time-averaged" distribution of states via This distribution is uniquely defined by the system parameters and the initial angle θ 0 .The above equation is, therefore, a suitable definition for a general angle distribution function (ADF) that describes the asymptotic state distribution of the system for the parameter tuple (M/γ, T γ, θ 0 ).If the limit j → ∞ in Eq. ( 21) exists, the stationary state distribution W (θ|θ 0 ) is equivalent to the ADF.In what follows, we will use the definition (22) when discussing the steady-state distributions.
For any given quantum trajectory on the GC, the probability of a single-step transition θ j → θ j+1 is determined by Born's rule and is equal to the corresponding probability P µ (θ j ).If the quantum trajectory visits the vicinity of every point θ with W (θ|θ 0 ) > 0 many times, all possible transitions between the discretized angle intervals are probed.The probability of any transition is then repeatedly sampled according to the Born rule, which means that the time-averaged distribution of states from a single quantum trajectory should converge to the ADF.If this is the case, the ADF has a simple interpretation in terms of quantum trajectories: If a single trajectory is observed for a sufficiently long time, the fraction of time it spends in a certain interval of the GC is determined by the integral of W (θ|θ 0 ) over that interval.The long-time behavior of almost any quantum trajectory is in this case completely described by the model parameters, being independent of outcome sequences.This time-averaging [cf.Eq. ( 22)] is, in particular, naturally implemented in the numerical simulations based on the Monte-Carlo procedure.
Following this reasoning, we can investigate a distribution for an individual quantum trajectory T of length m-without performing any explicit average over outcomes (or trajectories): where θ T (t > j ) is the θ angle at time t > j for the quantum trajectory and θ T (t 0 ) = θ 0 .Equation (23) allows one to approximate the distribution (22), if any typical path for sufficiently large m reproduces the distribution (22), as described above.This approach serves as a basis for the Monte-Carlo numerical simulations.
In any numerical analysis, representation of the ADF necessitates discretization of the angles.As a result, δfunctions in Eq. ( 17) are necessarily regularized.We will discuss the procedures of numerical evaluation and characterization of W (θ|θ 0 ) in Sec.V. Instead of taking the limit ε → 0 after introducing peak broadening of order ε, we will use the equivalent regularization where lim ε ′ →0 W (ε ′ ) is the set of unregularized δ-peaks from Eq. (17).At finite ε this defines a discretized distribution function, as for example in Eq. (23).

D. Master equation
The time-dependent distribution of states Eq. ( 18), can be obtained from an iterative integral master equation (ME) with the initial condition W 0 = δ(θ − θ 0 ): The distributions W (±) j here describe the endpoints of the evolution graph obtained after the click (no-click) readout at the last step.For brevity, we skip the initial-angle argument of W and encode θ 0 in the initial condition W 0 (θ) = δ(θ − θ 0 ); note that the ME has the same form for arbitrary W 0 .In fact, in most of the cases we consider, the limiting distribution will not depend on the initial angle.
In general, the stationary ME can also have solutions that do not correspond to an ADF, if the steady state is degenerate.As we will see below, the existence of a unique steady state depends on the setup's parameters.However, in most situations, the steady state is nondegenerate, allowing us to infer the ADF from it.In the following, we focus on such cases when the steady state of the ME is described by the ADF.In this case, the θ 0 argument in the ADF can be dropped and the ADF can be investigated with the help of the ME.We will point out special cases when such a treatment is not justified.
Performing the integration in Eq. ( 26) by using (for invertible maps) the identity we arrive at the functional equation for the steady-state distribution: This type of functional equations [118] was addressed in the literature devoted to random products of matrices (cf.Refs.[113][114][115]), where the probabilities of applying the matrices are typically state-independent.In the present problem, the dependence of P µ on θ is dictated by the quantum-mechanical Born's rule.

E. Characteristic features of the ADF, quantitative indicators, and classification of phases and regimes
The solution to the functional equation ( 27) cannot be obtained analytically in a closed form, except for some fine-tuned special cases.In particular, the reduction of the master equation to the Fokker-Planck form (cf. Refs.[77,98]) is possible in the limiting case of frequent measurements, when the system's own evolution governed by Ĥs is slow compared to the measurement rate, so that the change in θ after the no-click measurement is small.Another possibility is related to various kinds of commensurability in the setup's parameters, which leads to simplifications in the functions P µ (θ) and F µ (θ).In general, however, the solution can be obtained only numerically.Our strategy below is to identify the relevant special cases, where the analytical treatment is possible, and guess the overall "phase diagram" describing different types of ADF's behavior based on the exactly (or nearly exactly) solvable cases.The rest of the parameter space will be analyzed numerically.
Experience from related works (e.g., Refs.[118], [115], [77]) suggests that the steady-state angle distribution W (θ) and, thus, the ADF can show a rich phenomenology of characteristic features, being either smooth or singular.Specifically, it may exhibit, for instance (cf.Fig. 1), (i) isolated narrow peaks with vanishing background between them (akin to "localized phase" in the terminology of Anderson localization); (ii) a smooth background covering all the angles, with fluctuations on top (akin to "ergodic metallic phase"); (iii) coexisting regions of nonzero values separated by segments where the function vanishes (akin to "granular metal"); (iv) a fractal-like pattern of singularities distributed over zero or nonzero background.
Indeed, in our analysis, we encounter all these types of behavior, see Secs.IV and VI.
In what follows, we will establish a classification of the regimes existing in our setup and demonstrate the existence of finite areas (distinct phases) in the overall two-dimensional phase diagram.For this purpose, we will fix the value of the hopping matrix element γ in Ĥs and vary the parameters γT and M/γ.For each point in this two-dimensional parameter space, based on the steadystate angle distribution, we will calculate the following quantitative indicators: 1. Participation ratio R and its scaling exponent upon coarse-graining, Sec.V A; 2. Support of W (θ): the fraction of angles yielding a given fraction of the total probability, Sec.V B; 3. Position of the maximum in the histogram of heights for the discretized angle distribution, Sec.V C; 4. Hausdorff dimension of the curve W (θ) obtained by covering it with square boxes, Sec.V D; 5. Ergodicity marker of the Markov process derived from Eq. ( 27) for a given discretization scale, Sec.V E.
These indicators distinguish between localized and extended distributions and describe the degree of nonergodicity and fractality in a given setup.Although none of them can serve simultaneously as both the necessary and the sufficient condition of localization or delocalization, the combination thereof is foreseen to give a consistent description yielding a well-defined phase diagram of the system's state in the long-time limit.

III. SINGLE-STEP SOLUTION
In this Section, we analyze the quantum dynamics between the initial time t 0 = 0 and the first postmeasurement time t > 1 , or, equivalently, between two successive post-measurement times t > j and t > j+1 .In terms of the matrix formulation of the problem, we elaborate on a single-step solution for the mapping given by Eqs.(8,11), and (14).We start with a warm-up example of the system with decoupled sites, γ = 0, and study the state of the first site.We use this simple example to set the stage, illustrating the basic steps of the approach.Much more important is the general case with finite tunneling, γ ̸ = 0. We will explore it for fixed half-filling where a single electron tunnels between sites s 1 and s 2 .
A. Warm-up exercise: γ = 0 The model with γ = 0 was discussed extensively in Ref. [69] (see also Refs.[96] and [102] for a different form of Ĥsd ); here, we repeat the derivation, as it turns out to be instructive for compreheding the general case of γ ̸ = 0.If there is no tunnelling between the sites, we can consider only the first site-the one that is coupled to the detector.With a single electron residing in the two-site system, this site can be either empty, the state |0⟩, or occupied, the state |1⟩.The initial state of the system, given by Eq. ( 5), can be prepared by hybridizing the two sites at times t < 0 with the subsequent switching off the hopping between the sites at t = 0.The repeated measurements will then disclose whether the first site is occupied or not.Using the parameterization by angles θ and φ, Eq. ( 6), the initial expectation value of occupancy of the first site is equal to cos 2 (θ 0 /2).
Using the basis b(1) we construct the matrix form of the Hamiltonian and the evolution operator here T ≡ t k+1 − t k .The setup state between the initial time and the first measurement reads Right after the first measurement and reinitialization of the detector, the setup state becomes where denote the setup state after click, µ = +, and no-click, µ = −, outcomes, and P (±) 1 are probabilities of the click and no-click outcomes of the first measurement, This procedure is repeated for later post-measurement times and yields the mapping In terms of matrices Mµ , the mapping is given by Note that M+ is a projecting matrix.Equations ( 35) and (36) show that the stroboscopic values of the angle θ defined at the post-measurement times, θ j , change during the quantum evolution while those of φ are equal to its initial value, φ j = φ 0 .Besides, φ 0 has no effect on the evolution of θ.In particular, regardless of the value of φ 0 .If we are interested in an observable that is not sensitive to φ 0 , e.g., the occupation of the first site, we can choose any value of φ 0 and explore the evolution of the trajectory parameterized by a single angle θ.

B. Main model: γ ̸ = 0
In the previous Section, we have set γ = 0 such that the expectation value of the occupation of the measured site is changed only by the measurement backaction.Let us now take into account the quantum dynamics of the system (governed by Ĥs ) between two successive measurements, which is due to finite tunneling between the measured and non-measured sites.Finite tunneling introduces a new energy scale γ, or In this case, none of the angles parameterizing the system state remains constant in the post-measurement mapping.The states could be anywhere on the Bloch sphere and should be parameterized by both Bloch sphere angles, θ and φ.However, as we show below, almost all trajectories are attracted to the GC in the long-time limit and, in this limit, the quantum trajectories can be parameterized only by the angle θ: The first steps of the description of two tunnel-coupled sites interacting with the detector are very similar to those explained in Sec.III A. We use the basis b(2) and construct the matrix Hamiltonian Next, we choose the initial state in a full analogy with the γ = 0 case, and solve the evolution equation between measurements.The phases of amplitudes α and β in ψ now enter expressions for Ψ (µ) (t > j ) in a nontrivial way.Therefore, the system state should generically be parameterized by two angles: However, we show at the end of this section that φ j → π/2 for j ≫ 1, yielding the GC state, Eq. ( 42), in the long-time limit for almost all quantum trajectories.The discrete evolution of the system state is described by the matrix mapping in the form of Eq. (11).Expressions for the probabilities P (µ) j are rather cumbersome and, since they are of secondary importance for the current explanations, we present them in Appendix A 1. The "no-click" and "click" matrices have the following form (see Appendix A 1 for algebraic details): where we have introduced short-hand notations: Determinants of these matrices read Other properties of matrices Mµ are described in Appendix A 2. Let us briefly recapitulate here the most important ones which will be used for the analysis of the system dynamics.The matrices M± are symmetric, M± = M T ± , not Hermitian, but generically invertible.The latter property allows one to find the previous state by backward-time evolution of the current state with the account of the measurement outcomes.The exception includes those system parameters, at which det M± = 0, see Fig. 8, such that at least one of the matrices is a projector up to the normalization of its nonzero eigenvalue, and the evolution cannot be inverted.Depending on the parameters M/γ and T γ, eigenvalues of M− and i M− can be either complex-valued (and then complex-conjugated to each other), or purely real and generically not equal to each other.In the former case, both eigenvectors point to the equator of the Bloch sphere, θ = ±π/2, while, in the latter one, they point to the GC.As we emphasized above, the GC plays a special role in our consideration.Multiplying the matrices M± with the column vectors {cos(θ/2), i sin(θ/2)} T and restoring the normalization, one proves that the GC is an invariant manifold of the post-measurement mapping.Furthermore, a generic trajectory is attracted to the GC, see an example in Fig. 5. Arguments supporting the attraction to the GC are given in Appendix B. The main point of this consideration is that there always exist products of matrices, M(k, l) = i M+ k M− l , whose eigenvectors point to the GC and, hence, these products tend to project the electron state to the GC in the limit of an infinite number of measurements Figure 6 shows minimum powers k at which the product M k + M − starts projecting the state to the GC.More precisely, the invertibility of the maps together with the GC forming an invariant set for the inverse maps means that a trajectory can never truly arrive at the GC.However, the closer they get to the GC, the harder it is to escape from its vicinity, because of the smoothness of the map.This allows us to restrict ourselves mainly to the initial conditions located on the GC and to parameterize the quantum trajectories in terms of a single angle θ.It is worth mentioning that attraction to the GC is a consequence of the equivalence of the two on-site energies in the system.If the on-site energies are different and the symmetry between the sites is broken, attraction to the GC disappears.

IV. ANALYSIS OF SOLVABLE CASES
The definition of the ADF, Eq. ( 22) implies averaging over all quantum trajectories, which is generically very nontrivial at γ ̸ = 0 and often requires extensive numerical simulations.We have managed to develop the (mostly) analytical description for finite tunneling only in several solvable cases described in this Section.We focus mainly on the GC and point out extensions of our results beyond the GC wherever necessary.
The special cases that can be addressed analytically at γ ̸ = 0 are related to (i) commensurability effects, (ii) the existence of periodic orbits on the GC, or (iii) the projecting nature of matrices M± .We will show how localization in W (θ) emerges in these cases.The structure of the entire phase diagram of the system is essentially determined by these special cases.
A. Localization at γ = 0 Before delving into the analysis of γ ̸ = 0 cases, let us return to the simple model of Sec.III A and find the ADF, W (θ) for γ = 0.The initial state is assumed to be θ(t = 0) = θ 0 and φ 0 = 0 (the latter equality makes intermediate equations shorter).Dependence of the final state on φ 0 will be trivially restored at the end.Equation (36) suggests that, after the very first click event, the system state is projected to |1⟩: the measurement reveals that site s 1 is occupied with probability one and it remains occupied for all later post-measurement times (since there no hopping to site s 2 ).In this case, the system is in the state |1⟩ at j → ∞ and the ADF becomes δ(θ).If the quantum trajectory contains only no-click outcomes, the final state at j → ∞ is |0⟩, see Eq. (35), which yields another contribution to the ADF: δ(θ − π).
Let us calculate the probability of this specific ("null-measurement") outcome, where L is the length of the no-click-trajectory.Recalling that we have chosen Im[α j ] = Im[β j ] = 0, the following equality holds true: The normalization of the post-measurement system state requires the left-hand side of Eq. ( 55)-and, hence, the right-hand side-to be unity.Combining Eqs. ( 35) and ( 55), we find: This shows that the ADF is given by The weights of delta functions coincide with the initial probabilities for the first site to be occupied, cos 2 (θ 0 /2), and empty, sin 2 (θ 0 /2).Thus, repeated generalized measurements of one site yield, in the limit of an infinite number of measurement steps, the same distribution of final states as a single projective measurement.An exception is a trivial case where the probability of the click event vanishes and there is only one quantum trajectory consisting of no-click events: with l being integer.In this special case, θ j = θ 0 for all j and W (θ|θ 0 ) = δ(θ − θ 0 ) coincides with the initial ADF.Another special ("commensurate") case is realized for with integer l.According to Eq. ( 32), a single no-click outcome in this case also immediately projects the system state (now, to |0⟩), thus, both matrices M ± are projectors then.This implies that the generalized measurement becomes a strong (projective) measurement (cf.Ref. [69]).
If instead of a single initial state with fixed θ 0 one prepares a set of states characterized by an initial distribution of θ 0 , W (θ|θ 0 ) in Eq. ( 57) should be averaged over θ 0 : The distribution functions (57,58) fully describe the statistics of the occupation of the measured site: since φ does not influence evolution of θ, the full ADF is factorized into a product W (θ) × δ(φ − φ 0 ), and similarly for W being averaged over the initial angles.
To conclude this section, let us reiterate that the qubit is always found within one of two states, {θ = 0, φ = φ 0 } and {θ = π, φ = φ 0 }, at γ = 0 in the limit j → ∞.The initial condition for θ 0 can only change the relative weights of these two states but is unable to modify the structure of the ADF containing only two δ-peaks.Such a distribution falls into our operational definition of the localized phase.The localization of the angle θ can also be considered as a measurement-induced steering [96,100] of an arbitrary state of the qubit to the states with θ = 0, π.This example also allowed us to pinpoint two special cases: M T = πl and M T = πl/2, corresponding to "frozen dynamics" and strong-measurement dynamics, respectively.In what follows, commensurability conditions of this sort will enable us to identify special cases of dynamics that will determine the phase diagram of the model at γ ̸ = 0.

B. Effects of (in)commensurability and commensurability transitions on GC and beyond
We now return to the general case of finite hopping γ ̸ = 0.

Frozen case
Let us start with a simple case, where Y T = 2πq, with q being integer.
Using equations from Appendix A 1, it is straightforward to show that the outcome probabilities are state-independent and the electron post-measurement states differ from previous ones only by a total phase: Thus, θ is not changed by the post-measurement mapping and its ADF coincides with the initial distribution at t = 0: Moreover, similar to the γ = 0 case, we can trivially restore the entire distribution function of states even beyond the GC: both angles θ and φ are frozen at their initial values.Hence, the initially localized distribution remains unchanged at any j = 1, 2, . ... One can refer to this commensurate case as the "frozen" case.

Shift case
The outcome probabilities also do not depend on the system state if M T = qπ with q being integer: The choice of the upper index of P (±) depends on the parity of q.To be more specific, choosing odd q = 2l + 1, we arrive at the following expressions for the matrices M± : where σ0,1,3 denote Pauli matrices acting in the system state.These matrices satisfy the relations: where p is integer.Equation ( 65) suggests that any product of matrices M± can be reduced to one out of four possible forms: This is another rare case where there is no attraction to the GC.Nevertheless, the evolution of a GC state, Eq. ( 42), is remarkable.Matrix σ3 trivially inverts the sign of the angle θ.Matrices M p + shift this angle: where ϕ is the phase of the eigenvalue υ + (the upper index of υ + is not important), see Appendix A 2. If the phase ϕ is commensurate, ϕ = p 1 π/p 2 , the matrices M+ yield p 2 -periodic (or 2p 2 -periodic) trajectories along the GC.Hence, the initially localized distribution of θ remains localized at infinite time.Examples of these scenarios are presented in Fig. 7.
The set of the phases ϕ that correspond to localization is somewhat similar to the set of parameters that yield the structure of the well-known Hofstadter butterfly [119], where the density of states is also determined by certain commensurability.Specifically, the structure of the energy levels in this model is governed by the ratio of magnetic flux through a lattice cell and the flux quantum.For rational values of this parameter entering the Harper equation [120], the density of states is represented by a finite set of peaks, similar to the structure of W (θ) for commensurate values of the shift angle (upper panel in Fig. 7).If the phase ϕ is not commensurate, the 0.0 0.5 1.0 1.5 shifts pϕ fill the entire GC at infinite time.This breaks any initial localization of the θ distribution (see an example in lower panel of Fig. 7).Whether or not the GC becomes filled homogeneously or the ADF becomes fractal similar to the spectrum of the Harper equation in the incommensurate regime depends on the combinatorics of various quantum trajectories and can be checked numerically.The case of M T = 2lπ is analogous to the above case of M T = (2l + 1)π up to reshuffling M+ ↔ M− .
We note in passing that the post-measurement states in quantum trajectories at M T = qπ are located on two 1D circles even beyond the GC.For example, in the example of odd q, the trajectories are governed by the intersection of the Bloch sphere with two planes being parallel to the GC plane and located symmetrically with respect to the GC.The matrix M− reflects the state with respect to the GC plane, while M+ rotates the state around the x-axis.The expression for a single step of the rotation is rather cumbersome and we omit it for the sake of brevity.Let us emphasize that such rotations can also be commensurate or incommensurate.

C. Period-2 trajectory
Above, we have considered two cases where the localized nature of W (θ) at j → ∞ is related to localized initial conditions -quantum dynamics and measurements are unable to destroy localization due to commensurability.Let us now analyze other solvable examples where localization of the ADF emerges at long times regardless of the nature of the initial conditions.
In this Section, we focus on the setup with where the post-measurement map again simplifies: For the GC, the probabilities of outcomes take the form: Noting that we conclude that the states {1/ √ 2, ±i/ √ 2} T form a period-2 trajectory of the matrices m± .Moreover, this is a limit cycle that generically attracts quantum trajectories.In fact, it is a kind of weak attraction: a typical trajectory will always return to any previously visited point on the GC; it is just that the most time is spent around the attractive points.The origin of attraction can be explained as follows: since m2 µ ∝ σ0 , any product of the matrices Mµ reduces to one of four possible forms: The eigenvalues of the products mµ m−µ are real numbers 2ηM γ − Y 2 sin(M T )/2.Importantly, absolute values of the eigenvalues are different, with the maximum being The eigenvectors of these products are {1/ √ 2, ±i/ √ 2} T and point to the intersections (two opposite points) of the equator with the GC.In the limit of large time, typical quantum trajectories are characterized by k 2 that is of the order of a large number of measurements for a long-time quantum trajectory, k 2 ∼ j ≫ 1.This follows from the analogy with classical random walks in 1D, where, following the diffusion law, the mean-square displacement (the counterpart of k 2 ) is proportional to time (the counterpart of the number of measurements), see Appendix D for details.
If the typical value of k is large, one can keep only the "main" eigenvector in the expansion of the matrix products: Thus, the matrices Mµ M−µ k project any state onto one of the two states {1/ √ 2, ±i/ √ 2} T .Finally, we note that which means that, in the long-time limit, all combinations in Eqs.(73)(74)(75) with k ≫ 1 project any state onto one of two states: {1/ √ 2, ±i/ √ 2} T .Thus, we conclude that, if Y T = (2l + 1)π, the typical quantum trajectories generate two universal (independent of the initial state) peaks in the ADF W (θ) at θ = ±π/2.This can be regarded as steering an arbitrary quantum state to these specific states [96,100].Numerical simulations confirm the above semi-phenomenological explanation, see Fig. Lines on the plane of dimensionless parameters {M/γ, T γ}, where det M± = 0.

D. Projecting matrices
There is a special case of real eigenvalues of the matrices M± , where one of them is zero and, therefore, the corresponding matrix is a projector.This holds true if the determinant of one of the matrices vanishes: In agreement with our prediction (Sec.IV C), the distribution consists of two peaks of (approximately) equal height, located at angles θ = ±π/2.Middle panel : Projective case, where one of the matrices M ± has a vanishing eigenvalue, projecting to the "main eigenvalue" (the leftmost peak).The blue curve was obtained numerically by performing a time average over a single random Monte-Carlo post-measurement trajectory.The dotted curve is the analytical prediction, Eq. ( 81), truncated at 20 terms.Lower panel: Almost projective case in the limit M/γ → ∞ and T γ → 0, where the matrix M+ projects to θ = 0, while matrix M− introduces a small shift.The inset shows the same distribution in the interval θ ∈ [−0.1, 0.1].
Let, for instance, M− be the projector.After the very first no-click event along the quantum trajectory, the system is projected to the eigenstate Υ and the process starts over.The resulting distribution is the limiting case of the pattern predicted in Ref. [115].
Such a structure is insensitive to the initial distribution and, if there is a gap between the main peak and the satellite, as well as between the satellites, W (θ) appears localized at j → ∞.Numerical simulations confirm this scenario, see the middle panel of Fig. 9.If the determinant of one of the matrices is close to zero but not equal to zero, the structure of localized peaks becomes more and more smeared with increasing deviation from conditions (79)(80).
The ME (27) for M− projecting to Υ (−) p is formally solved by the following ADF: where Here, the normalization factor N contains information about all of the satellite peaks, which are projected back onto the mean peak by an application of matrix M− and develop again by sequences of M+ .In the stationary state, the weight of the main peak is determined by the projections from all satellites, weighted by the respective probabilities.The weight of the satellites is, in turn, balanced against the main peak.It can be seen that this is achieved by first generating all satellites and then renormalizing this pattern.A good approximation to Eq. ( 81) can be obtained by truncating the sum, since higherorder terms are exponentially suppressed, see Fig. 9.
An interesting case of projective matrices is realized in the frequent-measurements limit M/γ → ∞ (so that Y ≃ M ) and T γ → 0, with M 2 T /γ kept constant.In this limit, the matrix M+ becomes projecting, whereas M− induces very small shifts: |Θ − (θ) − θ| ≪ 1.The ME (27) then reduces to a Fokker-Planck equation with resetting, see Refs.[77,98] for details.The limiting ADF contains an extended region in a finite domain of angles; beyond this domain W (θ) = 0 (see the lower panel of Fig. 9).Importantly, however, the smooth curve in this distribution appears only in the limit T γ → 0; for any finite T , the extended region comprises well-defined isolated peaks (see inset in Fig. 9).
In the special case where ⇒ det M+ = det M− = 0; both matrices are projectors and W (θ) consists of two peaks generated by these projectors: Here θ ± are the angles that correspond to eigenstates with nonzero eigenvalues of the projectors.Using equations from Appendix A 2, one can show that the probabilities of applying the matrix M− at angle θ + (and vice versa) are determined by choice of the solution of Eqs. ( 84) and ( 85) and are generically nonzero.This means that the "double-projecting limit" is not equivalent to the strongmeasurement one (cf.the case of M T = πl/2 for γ = 0 in Sec.IV A).
The weights A ± depend on the dwelling time of a given peak, i.e., on the typical length of sequences M− M− . . .M− and M+ M+ . . .M+ .Since the escape probabilities are the same for both peaks, Eq. ( 87), we conclude that Equation (86) gives an example of localization.Since the peak positions do not depend on the initial state, we come across the measurement-induced steering [96].

V. NUMERICAL SIMULATIONS: CHARACTERIZATION OF THE ANGLE DISTRIBUTION
In order to explore the phase diagram beyond the analytically solvable special cases, we have performed numerical studies of the model.A combination of two complementary methods-(i) solution of a discretized ME and (ii) Monte-Carlo simulations (see Appendixes E and F)has allowed us to systematically explore the ADF in the (M, T ) parameter plane (fixing γ = 1).The approach based on the iterative solution of the discretized ME is our main working tool, see Appendix E 1.The convergence of the numerically obtained distributions to the ADF as defined and motivated in Sec.II is extensively discussed in Appendix F. In essence, excellent agreement between results from the two completely different approaches proves their validity.We consider the domain M, T ∈ (0, 5], which turns out to include a broad range of regimes.Our understanding of generic distributions is based on the consideration of four cases (frozen, shift, period-2, and projective) in Sec.IV, which shape the rich "phase diagram" of the model.
To give an impression of various angle distributions for generic parameters, we fix (arbitrarily) M = 2.92 and present six different distributions, corresponding to different values of T in Fig. 10.Comparing several distributions, there is an immediate observation: For T = 2.5 and 2.7, W (θ) has heavy peaks around a few points and is close to zero at most angles.At T = 3.0, there are some high peaks, but the distribution has small finite values at all angles.For T = 3.0554, 3.1, and 3.722, the distribution is close to uniform (note the scale of the y-axis) but features an intricate structure on smaller scales.Based on the immediate visual difference between the distributions, it is tempting to refer to them as localized and delocalized.We now proceed with a more detailed description of the indicators that allow one to characterize the localized and delocalized angle distributions.

A. Participation ratio and its scaling
An ADF could be referred to as localized if there is a large chance to find the angle within a small interval on the GC [heavy peaks in W (θ)].If no such subset exists, the distribution is delocalized.To quantify this, we calculate a participation ratio (PR) for discretized angles where are the probabilities obtained from integrating the ADF over N discretization cells c i of equal size ∆θ = 2π/N (see Appendix E for details).A perfectly localized distribution Accordingly, large (small) values correspond to more delocalized (localized) distributions.
In addition to the PR, we also calculate its scaling with the number of coarse-graining cells, N g .Namely, we take the ADF at the highest available resolution (characterized by N ), superimpose a broader grid (characterized by N g ), and sum up the terms Pr 2 i in each broader cell-as opposed to numerically solving the ME at every considered discretization level with N cells separately.The scaling can be described approximately by a power-law with the PR exponent ζ.
In the preceding Sections, we referred to (de)localization as a property of wave functions in space in disordered systems.In a one-dimensional Anderson-localized system, the probability amplitude |ψ i | 2 corresponding to an arbitrary eigenstate of the disordered Hamiltonian falls off exponentially with the distance to the center site.In this case, a participation ratio is calculated as where L is the system size.In a localized system, the PR becomes independent of the system size if L exceeds the localization length ξ, whereas delocalization is defined by RL → L in the thermodynamic limit.Our model has a fixed size (the angle θ is compact), so that the role of the system size L/ξ → ∞ is played by the number of the GC discretization grid cells N .In the present problem, the counterpart of the modulus-squared amplitude of the wave function is the stationary probability Pr i on the GC.In a sense, increasing the number of "bins" in our model is similar to decreasing ξ (by increasing the strength of disorder) in the Anderson-localization problem.
The scaling of the PR with N does not generally have an equivalent meaning as the scaling of the PR with L in the context of Anderson localization.For example, for a box distribution which can be arbitrarily narrow, |I 1 | ≪ 1, we get R N (W ) ∝ N , as the fraction of the GC covered by the distribution is independent of N .The situation is different, if the distribution is given by a sum of delta-peaks: in this case, the discretized distribution W (∆θ) (θ) becomes narrower as ∆θ → 0, and the PR is constant as a function of N .In this sense, a distribution would only be localized, if its support on the GC decreased with increasing resolution.There can be situations where distributions are localized in the sense of a small PR value (narrow peaks), but delocalized in the sense that the support is independent of the discretization (if ∆θ is sufficiently small to resolve the distribution).Such distributions can be regarded as "metallic grains" (coexisting delocalized and localized regions), which are not generally expected in the problem of Anderson localization in homogeneously disordered systems.At the same time, we know of two special cases (period-2 trajectories and the projective case), where the distribution W (θ) is localized in the strict sense: the observed peak width scales to zero with increasing the resolution of discretization.Without additional analytical arguments, based only on the analysis of PR, we are limited by the minimum resolution ∆θ when distinguishing the true localization from the localization in the sense of metallic grains.

B. Support of the angle distribution
To further characterize the "localized-looking" ADF, we introduce another observable S C N , which captures the minimum support N C /N needed to cover a fraction C ≤ 1 of the total probability: where sorted(Pr) are descendingly sorted probabilities {Pr i } and N C is the smallest integer, such that the inequality is fulfilled.If the support is one, the distribution can be regarded as extended; otherwise, it can be localized or "granular".It is worth emphasizing, however, that this indicator explicitly depends on the cut-off scale C, which should be compared with the scale introduced by the discretization resolution in numerical simulations.

C. The typical height of the distribution
In the theory of Anderson transitions [2], the transition between insulating and metallic phases manifests itself in the local density of states [2,[121][122][123], On the localized side, there is just a small number of wave functions contributing to the sum (95) at any given site i (because most wave functions are localized away from i and their contribution to the local density of states at a given site is exponentially suppressed).As a result, the typical value of the local density of states (essentially the value at the maximum of the distribution) vanishes [2,121,123].On the metallic side close to the transition, many states contribute at any site, because the eigenfunctions ⟨i|ϵ n ⟩ are extended.The typical value becomes finite and the distribution is spread around the typical value [122,123].The order parameter of the Anderson transition is then a functional one: the distribution function of the local density of states.Characterizing this distribution function by the corresponding typical value, one can invent a simple-minded scalar "order parameter" that could be of practical use for a variety of questions in the context of Anderson transitions.
Inspired by this experience, we introduce a third observable to distinguish localized ADFs from delocalized.Making use of the similarity between the patterns of the local density of states and the angle distribution in our case, we introduce the height distribution H(h j ) of where i ∈ [0, N h − 1] and N h ≪ N is the number of distinct heights, see Fig. 11.This is a discretized version of the distribution of the "local density of states" (the local density of post-measurement trajectories on the GC): We count the number of bins, where the probability distribution at a given resolution lies within a given window of values (a "histogram of heights").We analyze the typical value of this distribution by considering the position of its maximum h max , such that H(h max ) = max(H).Based on numerical results, we distinguish three categories: Here h 0 = min(W (∆θ) ) + ∆h is the height corresponding to the leftmost bin in the histogram of heights with H(h i ) > 0, see Eq. (96).The first category implies a vanishing typical value in analogy to the insulating phase of a disordered system.The third category means a finite typical value in analogy to the metallic phase of a disordered system.The second, intermediate case corresponds to a nonvanishing typical value, however, at the left boundary of the distribution.This category describes, for instance, isolated peaks on top of a nonzero background in the ADF.

D. Fractal dimension
Having introduced observables to quantify "localization", we now take a closer look at the apparent substructure in some of the distributions.As an example, we consider M = 2.92, T ≈ 3.729 in Fig. 12 at high grid resolution N = 10 7 .The upper-left panel shows the entire distribution W (θ) with θ ∈ [−π, π).The other panels show sections of the distribution taken from progressively smaller intervals on the GC.The blue shaded areas indicate the intervals that are displayed in the next (zoom-in) panels.Remarkably, these four sections look similar to each other, suggesting that the distribution "repeats itself" on different scales, with the interval considered in the lower panel corresponding to 3 • 10 −4 /(2π) ≈ 5 • 10 −5 fractions of the GC.Numerically, we cannot further resolve this pattern without going to larger N .A heuristic argument suggests that this self-similarity can exist on any scale, rendering the distribution fractal.
We quantify fractality of the distribution by calculating its Hausdorff fractal dimension d: Overlaying W (θ) with a uniform grid of m −1 × m −1 cells, we count the number of cells C(m) required to fully cover the curve [124].The relation defines the box counting dimension d [125].If the structure can be fully resolved at finite m, we get d = 1.The dimension 1 < d < 2 corresponds to a fractal structure.Numerically calculating the fractal dimension, we cannot increase m above the number of grid cells N without trivializing the box-counting dimension.Thus, any curve with d > 1 should turn out to scale trivially when the box size is reduced beyond our numerical resolution.The emergence of fractality of the distributions can be exemplified by a heuristic consideration of the vicinity of the projective cases (cf.Ref. [115]).Let Θ + (θ) map a large fraction of the GC to a narrow interval around its main eigenangle, resulting in a slightly broadened peak.The matrix M− translates this peak to another angle interval, slightly "distorting" the peak shape (because F ′ − and P − are not constants).If we start from a peaked distribution around the main eigenangle, the translating map generates a set of decaying "peak clones" on the GC.Many of these peak clones are reflected back onto the main peak by the almost projecting map.Selfconsistency requires those modulations to be translated to the secondary peaks as well [cf.Eq. ( 81)].Recursively applying this argument suggests that the stationary limit is given by a fractal.
Generally, this mechanism is not limited to the projective cases.As soon as there is some back-and-forth copying between two points, fractality can emerge.Quantifying when exactly this breaks down towards the uniform cases requires a more careful consideration, which we do not provide here.A related analysis was performed in Ref. [115] for correlated random products of two matrices, where the appearance of singular peaks at "strong" eigenangles, as well as their "cloning", was found at the condition analogous to P ± (θ ± )|F ′ ± (θ ± )| = 1.The emergence of fractality is another parallel to the theory of Anderson transitions: At an Anderson transition, the wave function of the system becomes multifractal [2], which means that its self-similarity can be characterized by a whole set of nontrivial fractal dimensions by attributing a fractal dimension f (α) to a subset of points of the wave function that is characterized by scaling as L −α with the system size [2,[126][127][128].The function f (α) is called the singularity spectrum and can be extracted from the scaling of moments of the wave function (like the inverse participation ratio) with the system size [2,126].The box-counting dimension of the entire wave function is closely related to the singularity spectrum but contains less information [128].

E. Ergodicity indicator for the Markov process
Finally, we address ergodicity of the Markov chain defined by the matrix MN , Here, f µ (c i ) is the image of c i under the action of map F µ (θ ∈ c i ), Eq. ( 15), which is continuous on the periodic interval [−π, π) and invertible for all sets of parameters except for the projective limit; | . . .| denotes the length of the corresponding angle interval.The sparse matrix MN relates the cell probabilities Pr i , in the discretized ME, see Appendix E for details.There exist several different notions of ergodicity for Markov processes in the literature [129,130].In the following, we call the system ergodic, if the Markov process is irreducible [131].Irreducibility means that any state i (bin c i ) can be reached from every state k [129].If this is the case, the ADF is probed by any typical post-measurement trajectory implying ergodicity of the dynamical system [132].In our notation, for any i 0 , k 0 there exists a natural number s such that For the coarse-grained ADF, ergodicity implies: 1. Support of the ADF on the entire GC.
2. A unique stationary state of the ME.
3. Equivalence of the ADF, stationary state [133], and time average of a typical quantum trajectory.

Independence of the ADF of the initial angle.
To numerically check whether the process is irreducible for given parameters, we can consider its Markov matrix MN as the transition matrix of a directed graph where the N nodes from the set V N represent the grid cells c i and the edges from the set E N represent transitions k → i between the cells with finite (nonzero) probability in the matrix MN .Ergodicity of the Markov process with transition matrix MN is equivalent to G N having a single strongly connected component (SCC) containing all N nodes.The SCC is a maximal set of nodes v ⊆ V N , such that for every pair of nodes s, s ′ ∈ v there is a path s → s ′ (and also s ′ → s) that only traverses edges between nodes in v. Here, "maximal" means that no further nodes can be added to the set without breaking the latter property.
The SCCs of a graph G N can be calculated efficiently, within O(|V N | + |E N |) operations [134, 135] [136].Our ergodicity indicator has two values: 1 (ergodic, G N has a single SCC) and 0 (nonergodic, G N has more than one SCC).It is however not obvious whether or not this finite-N indicator can be used to classify the continuous process.Indeed, ergodicity can be an artifact of discretization.As an extreme example, suppose we discretized the entire GC into a single cell c 1 = [−π, π).This cell forms a single SCC, thus corresponding to an ergodic process.This dismisses any sub-intervals of the GC that may prove unreachable at higher discretization, which would render the process nonergodic.
Importantly, nonergodicity in M N does have implications for the continuous process.To see this, a simple example of a discrete, nonergodic process is helpful: let us assume that the set of nodes V N of the graph G N splits into two SCCs v 1 and v 2 .These SCCs can be either (i) disconnected, or (ii) connected only by edges in one direction, e.g.v 1 → v 2 .Edges in both directions v 1 → v 2 and v 2 → v 1 are not possible because this would mean that G N consists only of a single SCC.
In both situations (i) and (ii), the GC subset represented by the nodes from v 2 , I 2 ≡ i∈v2 c i , forms an invariant subset.This subset is also invariant for the continuous post-measurement evolution on the GC regardless of the value of N for which it was found, see Appendix G for details.The above argument can easily be generalized to situations where G N splits into more than two SCCs, the key insight being that, by construction, there always exists at least one SCC that corresponds to an invariant subset of the dynamical process.
We stress that the finite value of N in this argument does not imply any discretization of the dynamical process.Indeed, the mapping used in the construction of a finite-N graph G N is the mapping of continuous intervals of the angles (rather than a mapping of discrete interval labels), which preserves the continuous nature of the stochastic process.Discretization with finite N refers here only to the resolution of the instrument (a "microscope") employed to explore the continuous process by constructing a corresponding graph.
With the above prerequisites, we define nonergodicity of the continuous process in the following way: The process is nonergodic, if an invariant subspace of the dynamical process can be found that does not include the entire GC.Otherwise, the process is ergodic (for any discretization, paths between any two bins exist).Importantly, once our "microscope" detects such an invariant subset in the continuous process at some N , there is no need to further enhance its resolution: nonergodicity of the continuous process is already demonstrated.Finally, it is worth mentioning that nonergodicity does not imply degeneracy of the stationary state.In the above example with two SCCs, both cases (i) and (ii) are nonergodic, but only case (i) corresponds to a degenerate stationary state.

F. Summary of regimes and respective indicators
We have introduced four indicators which are suitable for identifying the behaviour of the ADF: 1. Localization and delocalization are reflected by the participation ratio and its scaling, the support measure, and the typical value of the "histogram of heights".These indicators have been used based on a loose analogy of the shape of the asymptotic angle distributions (possibly at different discretizations) with localized and delocalized wave functions in a disordered system.In the truly localized regime, we may expect R → 1, ζ → 0 regardless of the smallness of the bin size.
2. Fractality and self-similarity of the ADF curves are described by their box-counting Hausdorff dimension.
3. Ergodicity and nonergodicity of a discrete Markov process corresponding to the transition matrix MN is described by the connectivity of the graph G N .

VI. NUMERICAL RESULTS: (DE)LOCALIZATION, (NON)ERGODICITY, AND FRACTALITY
In this Section, we present the results of our numerical study of the structure of the "phase diagram" in the M -T plane for γ = 1.We obtain the stationary angle distribution W (θ) from the ME for a high-resolution grid in the domain 0 < {M, T } < 5, see Appendix E 1 for details.We classify the distribution into (non)ergodic, (de)localized, and fractal types by means of the characterization schemes described in Sec.V.

A. Cross-section of the phase diagram
We first investigate a single cross-section of the "phase diagram" through the parameter space at fixed M = 2.92, considering 640 equally spaced values of T in the interval T ∈ [10 −3 , 5], see Figs. 13 and 14.The distributions were obtained using the ME method with N = 10 5 grid cells, starting from uniform distributions and iterating for up to 10 4 steps.Different special conditions are indicated by vertical lines in Figs. 13 and 14, see captions of these Figures.Additionally, distributions for projective T -values were calculated from Eq. ( 81).

Delocalization and localization regions
The upper and middle panels of Fig. 13 show, respectively, the PR, R N , and the PR scaling exponent, ζ, calculated for N = 10 5 .First, we note that the PR exponent essentially follows the behavior of the PR and, therefore, we can focus on the dependence R N (T ).All regions with the values of R N that are attributed to localization or delocalization, can be understood based on the properties of the special cases.
In the limit T → 0, as well as at the values of T that correspond exactly to the frozen (dashed lines) and shift (dotted lines) cases, the PR indicates delocalization: R N ≈ N .However, this is an artefact of the uniform initial condition.At T → 0, the detector state is Zenofrozen (cf.Ref. [77]), since there is no (joint) unitary time evolution between the measurements.Therefore, the detector state never changes, and the measurement outcome is always no-click.At T = 0, the matrix M− acts trivially on the system state, such that it remains frozen as well (though not in an eigenstate of the projective density measurement).
Similarly, exactly along the frozen cases, both matrices M± act trivially and the system state never changes, while the analysis of the ME indicates a delocalized distribution.In these cases, any state is an eigenstate of the Markov matrix, and these degenerate stationary states do not correspond to the actual (frozen) ADF, which is given by the initial state.
In the shift case, the ME also always has a uniformly delocalized stationary state, as every bin has exactly the same incoming contributions.The actual ADF depends on the shift angle ϕ and the initial condition: If the shift angle is commensurate with π, nϕ = 2π for some n ∈ N, the ADF only has support on a finite set of points that depends on the initial condition.If the angle is not commensurate with 2π, any point on the invariant manifold can be approached arbitrarily closely and the ADF is delocalized.In any (shift) case, generic postmeasurement quantum trajectories do not converge to the GC.
Thus, the indicators based on the solution of the discretized ME with a uniform initial condition may fail in the special cases.In particular, they may overlook the localized points: "fake delocalization" at special points can emerge as a result of averaging over localized distributions.In this regard, the structure of the matrix MN can provide additional valuable insights into the expected behavior of the ADF, see Appendix E.
The situation is different already in the immediate vicinity of the frozen and shift cases: Freezing of state and Bloch-angle is clearly broken away from the commensurability points.The deviation lifts the exact de-generacy of eigenvalues in the frozen case, and the ME approach can be used-agreeing well with the MC time average (Appendix F).Note that the maximum (for given N ) value of the PR is achieved for finite segments of Tvalues around special commensurability points, suggesting the existence of delocalized phases.Specifically, the PR indicates regions of true delocalization of width ∆T ≈ 1/2 with almost saturated PR values R N ≈ N = 10 5 , where the distributions are close to uniform.
The localized behavior, R N ≪ N , is correlated with period-2 trajectories (solid lines) and the projective limit (dash-dotted lines).The behavior of the indicators shows, however, that localization at the period-2 trajectory is destroyed by slight deviations from these commensurability points.Perturbatively, we might expect that a product M± M∓ still has eigenstates close to the period-2 peaks.However, the attraction to the period-2 peaks relies on long chains of products M+ M− , which emerge from the contraction argument M 2 ± = 1, see Sec.IV C and Appendix D. If this contraction is not fulfilled exactly, the deviations are accumulated in long chains.This may lead to the broadening of the peaks in the ADF.
In the projective limit, the solution of the ME is given by a discrete set of delta-peaks, whose strength decays exponentially with the number of necessary transitions from the main projective peak, resulting in apparent localization.Indicators in Fig. 13 demonstrate that the stability of this projective localization with respect to variation in T depends on the value of T , in particular, through the "interaction" with neighboring commensurability conditions.Specifically, for projective lines at T ≈ 0.6, 2.5, 2.8, and 4.8, we come across relatively wide regions of small PR values.Regions around T ≈ 1.5 and 3.9 lie within narrow "valleys", while those at T ≈ 1.4 and 3.3 are enclosed by the shift and frozen lines and only show a sharp dip.At the same time, the combined effect of projective and period-2 cases favors localization in a relatively broad range of T , see the region around T = 2.7.Interestingly, all broad regions of apparent localization observed in Fig. 13 correspond to cases with a strong hierarchy in both post-measurement matrices.
The support measure S C N with C = 0.99 shows behavior similar to that of the PR and the PR exponent.Since nonergodic regions turn out to correlate with low values of this measure this establishes a link between (de)locatization and ergodic properties of the steady-state distributions.Finally, yet another indicator of (de)localization-the typical value of the "histogram of heights" shown in Fig. 14-also exemplifies the correlations between the indicators oriented on the (de)localization and (non)ergodicity.
In summary, such defined "localization" manifests in analogs of several common localization measures, all of them yielding similar regions in the cross-section, where the distribution appears localized.Whether these regions correspond to genuine localization (with the width of the ADF peaks tending to zero in the continuous limit, as it does exactly at the projecting limit or for period-2 trajec-tories) or to granularity (with the peak width saturating with increasing the resolution of discretization) requires additional analysis, which is beyond the scope of this paper.The consideration of (non)ergodicity and fractality below sheds more light on this question.

Ergodicity and fractality
We move on to consider the ergodicity indicator shown in the upper panel of Fig. 14.The plot shows whether or not the Markov process corresponding to MN with N = 10 5 is ergodic [whether or not G N defined in Eq. ( 103) has a single SCC].Importantly, nonergodicity does not necessarily imply localization.An example is given by a granulated case where delocalized grains are embedded in a localized background, which is characterized by the ADF not having full support.Such a distribution could look similar to the curve in the lower panel of Fig. 9, where a gap of very low density opens between intervals where the distribution is finite.As argued above, nonergodicity of the discrete process does imply nonergodicity of the continuous process-there are invariant disconnected subspaces on the GC.However, an analogous statement does not directly hold true for ergodic regions: ergocity of a discrete process does not necessarily mean ergodicity in the continuous limit.Nevertheless, the extension of ergodic regions at finite discretization observed in Fig. 14 suggests that genuine ergodic phases should exist in our parameter space.Indeed, these regions are remarkably stable with respect to the number of discretization cells, implying their stability and finiteness in the continuous limit.This is quite natural, in fact, once the existence of finite nonergodic regions has been established for continuous processes: all the rest should be then regarded as ergodic.
Most of the considered values of T in the cross-section correspond to ergodic distributions.However, we find nonergodic intervals of T in the vicinity of those projective (dash-dotted) lines, where both post-measurement matrices have strong eigenvalue hierarchies, i.e., the vicinity of the double projective case.To understand this, consider two small but finite intervals around the two "strong" eigenangles: A strong eigenvalue hierarchy means that a large fraction of the GC is mapped into the vicinity of the eigenangle by the corresponding map (see Fig. 18, upper-right panel).In particular, where Ξ µ (I µ ) is the image of the interval I µ under the action of Θ µ (θ ∈ I µ ).If these "attractive regions" of the intervals I µ are sufficiently large to include the respective other eigenangle, the condition facilitate escape from the intervals I µ , which thus form an invariant subset for the post-measurement state [137].
The fractal (box-counting) dimension displayed in the middle panel of Fig. 14 is trivial, d ≈ 1, around the frozen and shift cases.This is expected, since we learned from the PR values that these cases correspond to almost uniform distributions.For all projective cases, the fractal dimension shows dips to d ≈ 1.Based on the structure of the ADF that shows a series of exponentially shrinking peaks in this limit [see middle panel of Fig. 9 and Eq. ( 81)], we indeed expect trivial scaling C(m) ∝ 1/m in the exact projective limit.Some of the observed dips are extremely narrow, which suggests an increased susceptibility of the vicinity of the projective case towards forming a fractal pattern, in agreement with the argument in Sec.V D. Some of the period-2 cases also correspond to dips in the fractal dimension.This agrees with our expectation that, exactly in this case, there are only two peaks in the ADF and, therefore, d = 1.All other cases have a nontrivial fractal dimension-this appears to be the generic case in our system.
Although (non)ergodicity and (de)localization are not necessarily correlated types of behavior, a connection is established by the behavior of the support and the maximum in the histogram of heights.Interestingly, regions of nonergodicity are correlated with the first category of the height indicator (analogous to the manifestation of the insulating phase in the local density of states).Around all of the nonergodic projective cases, the maximum ventures into the first category.Additionally, around the projective case at T ≈ 4, there is a small first-category dip, which is not present in the ergodicity cross-section.Further, the behavior of h max fluctuates, which can be due to numerical limitations [138].The surrounding parameter regions of category-one behavior are "transitional" category-two regions (this does not have an analog in the theory of Anderson transitions).Finally, ergodic regions correspond to categorythree behavior of heights (delocalization in the Anderson picture) with occasional fluctuations into the second category around special lines.

B. Phase diagram of the monitored qubit
Having investigated a generic cross-section in the M -T -parameter space, we proceed with studying (de)localization, fractality, and (non)ergodicity in the whole parameter plane.Figure 15 shows the same quantities as in Figs. 13 and 14, calculated for a 160 × 160 grid in the domain M, T ∈ [10 −2 , 5] at N = 10 4 .The analysis of the entire diagrams is pretty much similar to that of the cross-sections, which we have elaborated in great detail.
The upper-left panel shows the PR values R N .Again, the localized and delocalized behavior can be distinguished based on the special cases, in analogy to what was discussed before.Delocalized phases are found around frozen and shift cases (dashed and dotted lines): as discussed above, the ME with a uniform initial condition yields delocalization along these special lines due to averaging over localized trajectories.In the vicinity of these lines, genuine delocalization quickly sets in for all quantum trajectories.Almost uniform delocalization around frozen and projective cases is remarkably stable with respect to changes in parameters, manifesting in broad yellow "bands" around commensurate lines.Localization is observed around the projective limit (orange and purple lines) and period-2-trajectories (solid lines).The period-2 trajectory crosses the delocalized bands through narrow "bridges" of localization (for example at T = 0.3, M ≈ 4.7).The projective lines are surrounded by a narrow or broad region of localization, depending on the parameters.The case where both GC maps M± are almost projective is special in that such regions always correspond to apparent localization in finite bands.Support measure and PR scaling phase diagrams are visually very similar to the PR diagram.
Comparing ergodicity of the discrete process (middleleft panel of Fig. 15) with the h max indicator (middleright panel of Fig. 15) we observe that the locations of the nonergodic regions are completely correlated with the locations of the most localized regions.There is an astonishing agreement between category three of the height distribution and uniform delocalization, as well as between category one and nonergodicity.Categories one and three are separated by the "transient" category two.
The fractal dimension as a function of M and T (lowerright panel of Fig. 15) confirms our conclusions from the analysis of the cross-section.Extended regions of almostuniform distributions around frozen and shift cases correspond to a trivial dimension d = 1.Away from these regions, the fractal dimension is nontrivial (except for the vicinity of the projective lines, which are not drawn to avoid covering fine lines of d ≈ 1).
As discussed above, the regions where both matrices Mµ have a strong eigenvalue hierarchy (double-projecting limit) are also special for ergodicity.As expected, in the vicinity of double-projective points (intersections of orange and blue curves in Fig. 8) we find finite regions of nonergodicity embedded into the mostly ergodic phase diagrams.To estimate the expected extension of these regions, we consider the following conditions for a nonergodic region: 1.Both matrices Mµ have eigenvectors on the GC, establishing the existence of a region according to relation (103); 2. The eigenangles θ −µ lie within the attractive region of the "partner" maps Mµ , fulfilling the relation (104).This corresponds to [139] Θ The nonergodicity criterion covers all regions which we find to be nonergodic at finite discretization (and which belong to category one of the height indicator), see the left panel of Fig. 16.Furthermore, it predicts that the regions of nonergodicity actually extend further in the parameter space (which would only be visible at higher discretization in our numerical procedures).
Under certain conditions, a direct connection between nonergodicity and localization can be demonstrated.Given a nonergodic process with images Ξ ± (I) ⊆ I for the invariant subset I, the inequality guarantees a localized distribution, starting from an initial state that is nonzero within the invariant subset, e.g., To see this, we consider the support S 1 of the distribution after applying the protocol once: This argument can be iterated indefinitely with the support shrinking exponentially with the number of iterations.Thus, the limiting distribution function consists of a set of δ-peaks on the invariant subset and is, therefore, localized.
The above argument provides a sufficient condition for localization of the ADF within the chosen invariant subset I. Generically, there can be several disjoint invariant sets, I j , in the GC and a similar analysis should be performed for each set.If the localization condition ( 106) holds true only for some I j , a "granular" structure, where the ADF is localized in some invariant subsets and delocalized over other angles, is not forbidden.If the initial angle θ 0 belongs to a "localized subset" of this structure, the state of the system is localized in the long-time limit.
More interesting is the truly localized regime where Eq. ( 106) is satisfied in all I j .At this point, we note that, besides for the double-projective points, the condition (106) can never be satisfied for the entire GC, because the maps Θ ± are bijections from the GC to the GC: The trivial map Θ(θ) = θ has slope one everywhere, and any map that has a slope smaller than one necessarily has the slope greater than one somewhere, such that the image of the GC is equal to the GC.However, this does not contradict the existence of true localization, since regions with µ |Θ ′ µ (θ)| ≥ 1 can be located in the complement Ū ≡ [−π, π)\U of the union U of all I j .Crucially, Ū is by construction unstable because of "leakage" into invariant subsets [140] and, thus, cannot be part of the support of the converged distribution.We conclude that the corresponding steady-state distributions are fully localized, since localization (in the sense of absence of states) within the complement Ū is established in the limit N → ∞.
We numerically confirmed that the condition (106) holds for most of the regions of nonergodic parameters, implying that the associated ADFs feature δ-peaks on the constructed invariant subsets.To evaluate this condition numerically, we calculate the invariant subspaces for the nonergodic parameters identified in Fig. 15 as the leaf nodes of the corresponding graphs (see also Sec.II and Appendix G).Parameters for which the condition is fulfilled for all invariant subspaces are indicated in the right panel of Fig. 16 as dark regions, indicating the localized phases within the nonergodic phase [141].105).Purple and orange lines correspond to projecting M− and M+ matrix, respectively.The M+ matrix has a strong eigenvalue hierarchy in the lower purple region as well, see Fig. 6.The color code is the same as in the middle left panel of Fig. 15.Right panel: Localization condition ( 106) evaluated for all the invariant subspaces obtained as the leaf nodes in nonergodic graphs GN , see Sec.V and Appendix G.The condition for localization is fulfilled in the dark purple regions.

C. Summary of the numerical analysis
Combining analytical solutions for the special cases with extensive numerical analysis, we have managed to identify various distinct regimes and phases in the monitored qubit's statistics and describe its phase diagram.The analytical approach has been used to build the skeleton of the diagrams, the numerical study has allowed us to understand the properties of the generic cases.
Analysis of the PR has allowed us to distinguish localized and delocalized regimes.The behavior of the system with respect to (de)localization can be qualitatively understood in terms of the special cases (Sec.IV).From the number of SCCs in the graph of the discretized Markov process, we have found ergodic and nonergodic phases.We have established the existence of finite nonergodic phases in the continuous process; hence, ergodic phases also survive in the continuous case.Thus, we predict a transition between the ergodic and nonergodic phases.The behavior of h max establishes a connection between (de)localization and (non)ergodicity [142].Further, we have identified the regions of a truly localized phase inside the nonergodic phase.Whether the nonergodicityergodicity and localization-delocalization transitions are equivalent in the monitored qubit is a challenging question, which we relegate for future work.Finally, we have found that fractal distributions are rather generic for solutions of a continuous functional equation governing the ADF.
The phase diagrams that are presented in this section, which reveal fractality and transitions in the statistics of a single monitored qubit, are our main results.

VII. DISCUSSION
A. Relation to previous work: Continuous measurements limit In the context of weak measurements, the concept of continuous measurements [86][87][88] of a qubit has been extensively discussed in the literature, both from the theoretical (see, e.g., Refs.[77,92,95,98] and references therein) and experimental [80-82, 143, 144] perspective.In our notations, the corresponding limit is M/γ → ∞ and T γ → 0, with the product M 2 T /γ being fixed.In this limit, after averaging over quantum trajectories, the evolution of a qubit is described by the famous Lindblad equation.
A system-detector setup similar to ours was theoretically considered in the continuous limit in Ref. [77] (see also Refs.[95,98] and the related earlier work [92]).Fo-cusing on the onset of the quantum Zeno effect, the authors of Ref. [77] have demonstrated that this onset is characterized by a series of dynamical transitions, which also manifested in the steady-state probability distribution of states (analogous to our ADF).It has been discovered that one of these transitions is seen as an opening of a region of forbidden states in the ADF.
Let us interpret this finding in terms of our classification of regimes and phases.The relevant example of the ADF is shown in the lower panel of Fig. 9.It implies establishing nonergodicity: the stationary ADF splits into two disconnected sectors.Specifically, a delocalized segment of angles is embedded in the background of zero height.This is not genuine localization in our language, as the delocalized region does not shrink to zero in the continuous limit.Following the list of expected types of behavior given in Sect.one can refer to such a regime as a "granulated metal".
We note in passing that there is another transition predicted in Ref. [77] which reveals itself as the emergence of a singularity in the ADF.The vicinity of the projective cases is associated with the emergence of power-laws in the ADF (cf.Ref. [115,116]), which further develop into the fractal pattern by the mechanism of "cloning."This type of possible transitions was not addressed in our analysis of phase diagram, as we focused on the indicators directly related to localization and ergodicity.
The consideration in Ref. [77] has focused on the far lower-right corner of the parameter plane, which is far beyond the range presented in our phase diagrams.Our study of the whole parameter plane reveals extended areas that correspond to ergodic and nonergodic phases of the dynamical process.Physically, the strength of the system-detector coupling and a finite time interval of their joint evolution conspire nontrivially to give rise to regimes that are much different from the quantum Zeno effect.
For instance, in the double-projective case, the system state can be deduced from the outcome of a single detector measurement, without knowing the state before the measurement.However, the probabilities of the detector readouts are independent of the previous measurements-in contrast to the Zeno regime.Instead, the system mimics a "quantum coin" with just two different states which are realized at random.In the frozen cases, the system state is completely unaffected by the detector readouts.Clearly, deviating from the continuous limit (stroboscopic measurements) adds complexity to the behavior and physical interpretation of the stationary probability distribution.This point is further underlined by astonishing parallels between the ADF and wave functions in the theory of Anderson transitions.
On a more technical level, allowing for arbitrary M and T can lead to a nonlinear implicit master equation describing the dynamical process on the GC.At any point on the GC, both measurement operators can induce large jumps of the state, which give rise to extremely complex typical quantum trajectories.This new (compared to the continuous limit) ingredient is reflected, for example, in fractal peak structures that are associated with back-andforth cloning of peaks in the ADF, or in the relevance of commensurability conditions.The entire parameter space is structured by special cases, where the associated Markov process reduces to a one-dimensional random walk (period-2, shifting, or projective cases), complete standstill (freezing), or displays mechanisms of projection and translation associated with the continuous limit [77]-but different in that both measurement operators still correspond to finite-angle steps.Despite these differences, the above-mentioned works, as well as the present paper, demonstrate the importance of the study of the distribution functions for identifying hidden dynamical transitions that are hardly observed in the averaged properties of a monitored qubit.This key message is anticipated to be relevant also for measurement-induced dynamics in large systems with a macroscopic number of degrees of freedom.

B. Implications of our findings
Our results have several immediate implications.We have already pointed out that the attraction of the quantum trajectories to the period-2 limiting cycle or to two main angles in the double-projective case realizes passive steering of the monitored qubit.The attraction of the quantum trajectories to the GC, where all states have the form {cos(θ/2), i sin(θ/2)} T , results in an effective finetuning of the tunneling amplitude to a purely imaginary quantity, γ → i|γ|.
The intriguing question is whether regimes and phases visible in the diagrams of Fig. 15 are experimentally detectable.The important argument in favor of observability is the broad range of system parameters where the regimes and phases show themselves.However, the phases can not be detected by monitoring conventional observables: Consider an observable that depends on the angle θ and introduce its mean value in terms of the ADF Two common examples are the expectation value of the site occupation numbers, N 1 = cos 2 (θ/2) and N 2 = sin 2 (θ/2), and the formally introduced entropy of entanglement between the sites: S 2 = −Tr[ρ 2 log(ρ 2 )], where ρ2 ≡ Tr 1 ρ, and ρ is the density matrix of the qubit.Note that S 2 does not describe any entanglement of particles, since there is only one electron in the system; nevertheless, this quantity carries some useful information about the monitored qubit.Using the results of Appendix A, it is straightforward to show that S 2 (θ) at post-measurement times can be expressed via N 1,2 (θ) as S 2 = −[N 1 log(N 1 ) + N 2 log(N 2 )] and is a smooth function of θ.
Note that ⟨N 1,2 ⟩ θ and ⟨S 2 ⟩ θ do not necessarily distinguish between the localized and strongly delocalized regimes.Indeed, the same expectation value N 1 = 1/2 can be obtained for both the uniform ADF and the localized ADF with two equal peaks at θ = ±π/2.For the same reason-the smooth dependence of these observables on θ-their averaged values are also unable to reflect the ergodicity-to-nonergodicity transition.
Instead, as emphasized at the end of Sec.VII A, the transition can be observed in "tomographic" measurements of individual quantum trajectories [77,[80][81][82], which yields the full state distribution function.It would also be interesting to relate the transition to the statistics of sequences of measurement outcomes (cf.Refs.[92,98]).A search for convenient observables that are sensitive to the transitions in the statistical properties of the monitored microscopic qubit remains an open question.

C. Outlook
We have demonstrated that the study of the smallest nontrivial models of ancilla-based measurements uncovers a lot of interesting properties of dynamical processes, physical implications of indirect measurements, and even the possibility of phase transitions in the dynamical behavior.Let us list some questions, which we leave for future studies, including possible extensions of our model and methods.
Introducing variations to specifics of our model, for example, imbalanced energy levels of the qubit, or a different kind of coupling to the ancilla, may generate a family of interesting quantum systems.It is a priori not clear whether or not all models can be described by a single angular variable.Another natural modification of the protocol would be to randomize the times between consecutive measurements, such that the protocol would be characterized by a single period T only on average.One can also omit the reinitialization of the detector, allowing the system to keep a memory of the measurement outcome (which corresponds to considering "correlated random products" of post-measurement matrices [115] in our approach).Robustness of the various regimes and phases as well as of the predicted phase transition to these modifications is worth investigating.
Regardless of the model details, it is interesting to understand better the mechanism of fractality, if it appears.
To this end, one can use the locator expansion, e.g., close to the double-projective case of our setup, see modern implementations of the method in Refs.[145][146][147].Exploring the possibility of multifractality in the category-two regions of the height indicator and analyzing fractality in terms of the singularity spectrum (see also the comment at the end of Sec.V D) may establish an even closer connection of transitions in the monitored qubits to the theory of Anderson localization.Regarding the numerical analysis of multifractality, a challenge of calculating scaling exponents for higher moments is the necessity to obtain the distribution at higher resolution N than used for the participation ratio.Besides fractal exponents of the limiting N → ∞ distribution, it could also be interesting to investigate the transient "dynamical" scaling of fractal exponents with the number of time steps of the protocol.
The similarities between Anderson (de)localized wave functions and the state distribution of the monitored qubit gives a hint to their common mathematical origin.A search for such an origin may start from a field-theoretical connection between Anderson and measurement-induced transitions which has been recently reported in Refs.[50,51].Understanding (de)localization in the continuous case based on the Master equation could rigorously establish the localizationdelocalization transition in the continuous model.
Especially interesting further studies are related to increasing the number of particles in the system which can be achieved by considering either a two-site chain with spinful fermions or bosons, or by slightly enlarging the number of sites that could then accommodate more spinless fermions.This would allow one to explore the interrelations between the ergodicity/nonergodicity (or possible localization) transitions in the monitored qubit with the entanglement between particles, thus going towards the field of measurement-induced entanglement transitions.Finally, it would be interesting to understand the influence of interparticle interactions in larger setups on the phases and regimes discovered in the present work.

VIII. SUMMARY AND CONCLUSIONS
We have studied statistical properties of a single qubit under the influence of stroboscopic ancilla-based measurements.The detector (ancilla) is represented by another two-level system.The qubit and the detector are coupled and evolve unitarily between the measurements.The detector is projectively measured at the end of each interval of the unitary evolution and, after this, is reinitialized in a given state.This protocol defines the two measurement (Kraus) operators (represented by 2 × 2 matrices in our model), which map a qubit state onto one of two possible new post-measurement states, depending on the measurement outcome.Any quantum trajectory of the discrete-time evolution of the qubit's post-measurement state is described by a random product of these measurement matrices, with the probabilities of applying them dictated by Born's rule.
The qubit state can be parameterized by two angles on the Bloch sphere.We showed that, in the long time limit, the quantum trajectories in our model are generically attracted to a one-dimensional circle on the Bloch sphere (Grand Circle -GC), which is described by a single angle, θ.We have characterized statistics of the qubit states in terms of the ADF W (θ) -the angle distribution function on the GC, Eq. (22).We have investigated the ADF by combining the analytical approach, suitable for some special cases, with two complementary numerical methods (solution of a discretized master equation and Monte-Carlo simulations of individual quantum trajectories).We have revealed the richness of the dynamics of the monitored qubit and demonstrated remarkable similarities to Anderson localization in disordered systems.Namely, as the system parameters are varied, the ADF of a monitored qubit exhibits a number of regimes reminiscent of those in the theory of Anderson transitions, appearing either localized or delocalized or even fractal.
Motivated by the analogy with Anderson localization, we have used several standard indicators from the localization theory and supplemented them with those from the theory of random evolution.The indicators encompass the participation ratio [Eq.(88)] and its scaling exponent [Eq.( 90)], the support measure [Eq.( 93)], the typical value of the state distribution (Fig. 11), the boxcounting dimension of the curve representing the ADF [Eq.( 98)], and the ergodicity indicator that characterizes the structure of Markov matrices [Eq.( 103)] for the qubit evolution.Analysis of their properties has allowed us to quantitatively describe W (θ), despite fundamental distinctions between the mechanisms underlying localization in a monitored qubit and in disordered systems.
Our combined approach has resulted in a solid classification of various emergent regimes and phases, which are reflected in the phase diagrams, Fig. 15.The analytical consideration of special cases (in particular, related to commensurability of the stroboscopic measurements) has been used to build the skeleton of the diagrams, while the numerical study has allowed us to understand the properties of the generic cases.Concretely, we have established the existence of ergodic and nonergodic phases and identified regimes of localization and delocalization behavior of the ADF.In a large portion of the parameter space, W (θ) exhibits fractality, which ranges from strong to weak, cf.discussion in Refs.[148][149][150] and references therein.We have predicted a transition between the ergodic and nonergodic phases and demonstrated genuine localization inside the nonergodic phase.This is our main result that may pave the way to the theory of various transitions in monitored qubits.Our work highlights the importance of studying the distribution functions of quantum states, as they can reveal concealed transitions that remain unnoticed when focusing on the averaged properties of systems subject to quantum measurements.The explicit forms of the mapping matrices M± are given by and (not normalized) eigenvectors Υ − of M± read as follows: 1.Two eigenvectors pointing to the GC and the other two pointing to the equator Let the eigenvectors of the matrix Mµ ≡ MEQ belong to the equator and eigenvectors of M−µ ≡ MGC point to the GC.The matrix MEQ itself is unable to provide the attraction to the GC.This is clear from the properties of a sequence of these matrices, as follows from Eq. (A23): The kth power of the "equator matrix," MEQ k , is a matrix with the eigenvectors located on the equator and generically complex eigenvalues [except for commensurate cases where arg v There is no reason to assume that such a matrix is able to project an arbitrary state onto the GC.Numerical results confirm the absence of attraction in postselected realizations, where there is no matrix MGC .
Next, consider the product where the bar means that the eigenvectors are normalized.Since the GC is the invariant manifold, the vectors Υ GC also belong to the GC, such that M k EQ M q GC with q > 1 generically projects a state |ψ⟩ onto the vector Υ (η,k) GC that corresponds to the larger eigenvalue v are close to intersections of the GC with the equator.The typical value of q can be estimated as q typ ∼ ⟨1/P EQ ⟩ (the typical length of the continuous chain of GC matrices).The position of the attractor depends on the typical value of the exponent k, which can be estimated as ⟨1/P GC ⟩ (the typical length of the continuous chain of equatorial matrices).
possible products of M+ and M− -corresponding to all possible quantum trajectories-simplify to be proportional to after N t time evolution steps.Here k is an integer with k ∈ [0, ⌊N t /2⌋].
From Eq. (D1), we already know all states with a finite probability after N t steps of time evolution.All that is left to do is to determine the probability of ending up in a specific state, taking into account all possible quantum trajectories that lead to this point.For this purpose it is instructive to draw the possible matrix products as a tree: Generally, for angle dependent transition probabilities P ± (θ), it is still difficult to determine the weight of a given node in above tree as there are still exponentially many possible outcome sequences that traverse the tree.Therefore, as a first step, consider a simplified situation with P + (θ) = P = 1 − P − (θ).We limit our consideration to even N t ; the generalization to odd N t is achieved by performing one more step of time evolution.If P = 1 − P = 1/2 counting the number of paths leading to one node in the tree suffices to determine the weight of the corresponding state, because each path of length N t has a probability of 2 −Nt .The counting generates the binomial coefficients: We thus find a distribution of the following form (note that this is not an ADF): where B[n, m] denotes the binomial distribution, ψ 0 is the initial state of the system (characterized by angle θ 0 ) represented as column vector [cf.Eq. ( 11)], and we introduced the short-hand notation angle[ψ] for the angle θ of state ψ on the GC.It was shown in Sec.IV that the matrix product ( M−µ M−µ ) n attracts the state to either θ = π/2 or to θ = −π/2 (depending on µ) for large n.This means that we can find a number n 1 , such that angle for any small ε and l ≥ n 1 .We compare the probability within these ε-regions to the probability anywhere else on the GC: This proves that the distribution converges to arbitrarily narrow peaks around the angles θ = ±π/2 for N t → ∞.
Since the tree structure is completely symmetric between the two attractive points for P + = P − , the time-averaged ADF reads after regularizing the δ-peaks of individual trajectories In order to generalize this result to imbalanced and nonconstant P ± , it is useful think of the tree structure as a one-dimensional random walk.The states of the random walk represent the different states that can be reached by forming combinations of M+ and M− matrices: Clearly, for constant P + ̸ = P − the random walk on this chain still follows a diffusion law without a preferred direction.The distribution of states at any given time can still be found exactly by diagonalizing the transition matrix where indices label the unit cells (dotted boxes) and the matrix elements are written in the A − B space of orange and blue nodes.Without calculating the state distribution explicitly, we can see from the symmetry of the random walk that the time-averaged ADF in this case still converges to the form (D4).
In the general case of angle-dependent probabilities P ± (θ), the calculation of the distribution of states more tedious.However, it stands to reason that the probability for the random walker to be found within the central region of the chain (around the node labeled 1) still vanishes as N t → ∞.If the random walker is sufficiently far from the central region, all states are close to the two special points ±π/2.In this situation, both matrices mediate transitions to the vicinity of the respectively opposite special point, such that the time evolution is again described by a simple random walk in state space, corresponding to the limiting time-averaged ADF (D4).We thus conclude that this distribution describes the period-2 case, as long as there is no other special mechanism confining the state to the central region.A stationary solution of the discretized ME corresponds to an eigenvector of MN with eigenvalue 1.As can be seen from Eq. ( 99), all entries of MN are positive.Furthermore, probability conservation is manifestly built into the ME,

N i=1
[ MN ] ik = 1. (E7) For this reason, MN is a positive Markov matrix and, as such, necessarily has an eigenvector corresponding to the eigenvalue 1 [152]: At least one stationary state of the ME always exists, regardless parameters and discretization.Such a state can be efficiently found for sparse matrices MN by repeatedly applying the matrix to an arbitrary initial state Pr 0 until convergence is reached ("power iteration")-provided that the eigenvalue λ = 1 is not degenerate and there is no other eigenvalue of modulus one [153,154].Assuming a hierarchy between the eigenvalues of MN , convergence to the stationary state happens at an exponential rate, which is proportional to the difference ∆λ = |λ 1 | − |λ 2 | between the dominant and the second-dominant eigenvalues (largest and second-largest eigenvalues by modulus), where n is the iteration number and a 0 is a positive number.There are two cases where the ME method can run into problems.(i) If the eigenvalue λ = 1 is degenerate, the converged state depends on the initial state.In this case, we cannot make an immediate connection to the asymptotic behavior of quantum trajectories.(ii) If there is only one eigenvector corresponding to λ = 1, but other eigenvectors have |λ ′ | = 1, the power iteration can converge to a state that does not correspond to a stationary solution of the ME.
In practice, we are able to find nondegenerate stationary states by power iteration efficiently for most parameter sets.There are two exceptions, corresponding to lines of special parameters: • Frozen trajectories are maximally degenerate, as every basis state (a probability vector comprising N − 1 zeros and one unity) solves the discretized ME.
• For certain "periodic" trajectories an iterative solution does not necessarily converge at all, as it can get stuck switching periodically between different angle bins (say, if there is a period-2 trajectory between two bins, there is a vector Pr, such that MN Pr ̸ = Pr but ( MN ) 2 Pr = Pr, corresponding to λ = −1).The nondegenerate solution to the ME in this case corresponds to balanced occupation of the two peaks.
The structure of the matrix MN allows us to understand important features of the resulting ADFs regarding localization or delocalization, see Sec.VI.For example, in the lower-left panel of Fig. 18 we show the structure of the Markov matrix close to the frozen case: The degenerate diagonal "band" of matrix elements splits into two separate bands which are slightly shifted from the main diagonal and bent towards opposite sites (crossing close to the middle and the ends).On the discretized level, it can be understood that a slight shift of the bands away from the main diagonal should lead to delocalization, since it essentially introduces transitions between neighboring grid cells.If the cells are smaller than the distance of the band to the diagonal at some point, then there is a next-nearestneighbor transition, but the other band can facilitate a back-transition to the cell in between, and so on.At the same time, perturbing the frozen case, the outcome probabilities P ± are almost independent of the state, such that the stationary state is almost translationally invariant.From the trajectory point of view, we can think of the state performing short-distance hops on the GC in a random direction, eventually covering many points on the GC in a diffusive fashion.Similarly, in the shift case, existing commensurability is also broken by a perturbation, leading to delocalization around the corresponding curve.
As a consequence of the above mechanism, the PR would reveal a sharp jump moving onto or off a frozen line, if we had chosen a localized initial state.Dynamically, this jump would appear as a crossover, since the distance from the commensurate line controls the "diffusion coefficient" in a post-measurement trajectory.Close to the frozen case, the time average of a post-measurement trajectory converges slowly.Similarly, the difference between the two dominant eigenvalues of the Markov matrix is controlled by the distance to the frozen line, such that the convergence rate of power iteration goes to zero as the commensurability is approached.

Monte-Carlo simulation
The idea of the Monte-Carlo (MC) approach here is to simulate (at most) a few post-measurement trajectories for a given initial state θ 0 , by randomly drawing measurement outcomes according to the Born rule, and performing time evolution according to the maps M± corresponding to the chosen outcomes.Specifically, for every measurement time t j , we draw a random number q(t j ) ∈ [0, 1] (uniformly distributed).Considering the probability P (+) j to apply the (+) map at that time (depending on the current position on the GC), we evolve the state with M+ , if q(t j ) ≤ P (+) j , and with M− otherwise.Importantly, one instance of the simulation follows one individual quantum trajectory (corresponding to a pure state).The resulting post-measurement trajectories are then time averaged; a histogram of visited angles is obtained.
This procedure provides information about the coarse-grained distribution W (ε) (θ|θ 0 ), if the outcome average (containing exponentially many terms in the number of time steps) is well described in terms of typical post measurement trajectories, see Eq. ( 23) of the main text.As this numerical procedure is based on randomly sampling possible outcomes according to their quantum mechanical probabilities, we refer to this method as MC simulation.The stability of the obtained histogram with respect to the number of initial time evolution steps, as well as the number of successive time points after the initial phase, indicates that the chosen m is sufficient for convergence.
We use the MC method to 1. verify that quantum trajectories generically converge to the GC (see Figs.While the MC method is useful owing to its simplicity, it converges rather slowly compared to the Markov method discussed below (with decreasing ε).
In Appendix F, we numerically confirm the convergence of generic post-measurement trajectories to the GC and find good agreement between the time-averaged distribution obtained from the MC simulations and the power iterated solution of the ME.The stationary state of the discretized ME is nondegenerate and its eigenvalue λ = 1 is always well separated from the other eigenstates of the Markov matrix (99), except for the vicinity of frozen and shift cases.This means that the GC distribution is generically independent of the initial state and can be obtained from the time average or stationary state equivalently.The general independence of the GC distribution from the initial state allows us to drop the initial state argument θ 0 from W (θ|θ 0 ) and investigate the distributions W (θ) within the chosen parameter range.All following distributions are obtained with the verified ME approach, allowing us to go to higher N and parameter grid resolutions, see Fig. 10. which should give a good qualitative picture of the distributions.We compare the distributions obtained from the MC simulation (time average of random quantum trajectories) to the stationary states of the discretized ME in the right panel of Fig. 19.We obtain the ME results using 10 5 grid cells and up to 10 4 iteration steps.Using a uniform distribution as an initial state for power iteration, the iterated solution is guaranteed to contain all stationary states, in case the stationary state is degenerate.If this is the case, the iterated state cannot be the same as the time-averaged state from a single post-measurement trajectory.The iterated distribution is coarse-grained to N = 10 3 cells.The distributions are compared by calculating the χ 2 -distance between the probability vectors from ME and MC calculations, We observe good agreement, d χ 2 < 10 −2 ≪ 1, for most values of parameters.Some exceptions are related to the frozen and shift cases (note bright markers exactly on dashed and dotted lines).If the freezing condition is perfectly fulfilled, both ME and MC methods preserve the initial state, which is localized at one point of the GC for MC evolution, and a uniform distribution on the GC for the ME calculation.In the numerical comparison, we ignore the fact that the MC initial condition does not necessarily lie on the GC and just compare the distributions over θ, keeping in mind that the frozen cases should be excluded from further analysis.This gives a large deviation d χ 2 N →∞ → 1 between the distributions.In the shift cases, there is again no convergence of the MC trajectory to the GC, while the ME calculation takes place entirely on the GC.Therefore, these cases must also be excluded from the comparison between MC and ME distributions.
Other points where the difference between MC and ME distributions is large are related to a strong eigenvalue hierarchy for M+ or M− .Vanishing λ µ min for µ ∈ {+, −} corresponds to the orange lines in the plot, and for small M λ + min /λ + max becomes small (see the dark blue region in Fig. 19, right panel).In the projective limit, the GC distribution can be calculated to high accuracy from a small number of terms in Eq. ( 81) and is independent of the initial angle.This distribution is also the unique stationary state of the ME.However, writing the ME according to Eq. ( 99) would no longer be valid, as the projective GC map is not invertible.This can also lead to inaccuracies close to the projective case.
Other inaccuracies with d χ 2 ∼ 10 −2 do not correspond to systematic problems with the methods but only to insufficient convergence (number of time points for the MC method, number of grid cells for the ME method).To demonstrate this, we show in the left panel of Fig. 20 a comparison between MC and ME distributions at one of the yellow parameter points with d χ 2 ∼ 10 −2 , but increase the number of time steps in the MC simulation to 10 8 , and the grid size for the ME to 10 6 , coarse graining both results to ∆θ = 2π/10 3 to approximate the corresponding W (∆θ) .As the inset shows, the distributions agree perfectly.Note that d χ 2 ∼ 10 −2 is still good agreement between the distributions, sufficient for the following numerical investigations.In the discussed example distribution it corresponds to a slightly different weight distribution in the heavy right peak between the two methods.
We have thus verified the convergence of our MC distribution and confirmed the validity of obtaining the GC distribution from the ME approach.A comparison between MC and ME calculations shows agreement between the distributions in all but the frozen and shift cases, demonstrating that generic trajectories initiated at random states on the Bloch sphere eventually evolve according to a unique distribution on the GC.The ME-matrix M10 3 has for all but the frozen, shift, and period-2 cases a single eigenvalue of modulus one.The nondegeneracy of the stationary probability vector confirms the uniqueness of the GC distribution at the chosen grid size.In the special frozen and shift cases, the stationary distribution is not unique and depends on the initial state.In the period-2 case, we can analytically calculate the stationary distribution, see Sec.IV C. In this case, the distribution W (θ|θ 0 ) is unique and independent of the initial state, but it cannot be reliably found using power iteration.It should be instead obtained either from the MC simulation or by explicitly finding the eigenvector of MN corresponding to λ = 1.
Appendix G: Nonergodicity in the discrete and continous process In this Appendix, we explain in detail how finding nonergodicity in the process described by a finite matrix MN at any number of discretization cells N has immediate implications for the continuous process.To see this, a simple example of a nonergodic process is helpful: Let us say that the set of nodes V N of the graph G N splits into two strongly connected components (SCCs) v 1 and v 2 Consider the directed edges that lead either from a node in v 1 to a node in v 2 or vice versa.In general, two situations are possible: (i) No such edges exist.
(ii) Edges exist only in one direction, say v 1 → v 2 without loss of generality.
Edges in both directions v 1 → v 2 and v 2 → v 1 are not possible, because this would mean that G N consists only of a single SCC.In both cases (i) and (ii), the subset of the nodes in v 2   There is an edge from v1 to v2 (but not vice-versa), because there are edges in G4 leading from the orange set of nodes to the blue one (but not vice-versa).The condensation is indeed acyclic as there are no circular paths.v2 is a leaf node of the condensation, because no transitions out of v2 are possible.Thus, it defines an invariant subset for the continuous process.This is because a matrix element [ MN ] i,k is non-zero ((k, i) ∈ E N ) if and only if one of the maps permits a transition k → i with nonzero probability, (G4) Calculating the overlap |f µ (c i )∩c k | does not involve any discretization-induced approximation for the angles, since the intervals c i and f µ (c i ) are continuous.In other words, this continuous mapping between the intervals is not equivalent to a discrete mapping between the labels i of the cells.The discretization of the GC into a finite number of intervals establishes the resolution of an "instrument" exploring the mapping between the continuous angles.Therefore, our nonergodicity construction leads to a statement about the continuous process and is independent of the fact that the subset v 2 was found by discretizing the GC.Examples of the finite-N representation of an ergodic process and a non-ergodic one described by situation (ii) are shown in Figs.21 and 22, respectively.
The argument is readily generalized to situations where G N splits into more than two SCCs.Formally, by using the language of the graph theory, one can define one "supernode" for each SCC of G N and introduce a directed edge between two supernodes, if there exists a corresponding transition between the two represented sets of nodes [155].This defines the acyclic condensation graph of G N [156].The construction of the condensation is illustrated in Fig. 22 and also for a more complex example with more than two SCCs in Fig. 23.The leaf nodes (nodes without outgoing edges) in this graph correspond to invariant subsets of the GC.In the steady state, all of the weight of the distribution is accumulated in the invariant subset corresponding to the leaf nodes, with zeros everywhere else.

FIG. 4 .
FIG.4.Exponential branching of states generated from the initial state (gray) by different combinations of postmeasurement matrices M± (the j = 2 trajectories are shown).After the j-th measurement, there are 2 j different endpoints (in principle, some endpoints can describe the same state).Each generated state (including all endpoints) is labeled by a bitstring showing the sequence of measurement outcomes along the quantum trajectory leading to the state.

FIG. 5 .
FIG. 5. Example of a trajectory (colorful markers) on the Bloch-sphere (light blue) converging to the GC.The parameters are M = 2.92 and T = 1 and the trajectory is initialized at (θ, φ) = (1.3,2.5).Scattered dots correspond to the positions of the quantum trajectory immediately after measurements.The time instances are color-coded, with light (dark) markers corresponding to early (late) times.The (YZ)-plane (GC-plane) is indicated in gray.

FIG. 6 .
FIG. 6. Upper left panel: Minimum values of the power k, at which both eigenvectors of M k + M− point to GC up to errors of the order of the numerical precision.We have analyzed k ∈ [0, 10 3 ].White regions correspond to k = 0. Black dashed lines correspond to the special (frozen) case Y T = 2lπ with l being integer.Upper right panel : Configurations of eigenvectors of M±. (0): all eigenvectors off GC; (1): only eigenvectors of M− on GC; (2): only eigenvectors of M+ on GC; (3): eigenvectors of both matrices on GC.Lower panels: Ratios of eigenvalues of the matrices.

FIG. 7 .
FIG. 7. Time histogram of a post-measurement trajectory in the shift case M T = π over 10 5 time instances.By tuning the shift angle to a multiple of π, a localized distribution emerges (upper panel ).If an incommensurate shift angle is chosen, all points on the GC are eventually covered (lower panel ).

FIG. 9 .
FIG. 9. Examples of GC-angle distributions W (θ) for the period-2 and projective cases (see Secs.IV C and IV D, respectively).Both panels show the numerical results (blue) obtained for γ = 1 with the different M and T -values shown in the plot titles.Upper panel : Period-2-case, with T Y = π.In agreement with our prediction (Sec.IV C), the distribution consists of two peaks of (approximately) equal height, located at angles θ = ±π/2.Middle panel : Projective case, where one of the matrices M ± has a vanishing eigenvalue, projecting to the "main eigenvalue" (the leftmost peak).The blue curve was obtained numerically by performing a time average over a single random Monte-Carlo post-measurement trajectory.The dotted curve is the analytical prediction, Eq. (81), truncated at 20 terms.Lower panel: Almost projective case in the limit M/γ → ∞ and T γ → 0, where the matrix M+ projects to θ = 0, while matrix M− introduces a small shift.The inset shows the same distribution in the interval θ ∈ [−0.1, 0.1].
corresponding to the nonzero eigenvalue and "eigenangel" θ − .The subsequent click event occurs with the probability P + = 1−(M/Y ) 2 sin 2 (Y T ) and drives the system out of Υ (−) p .The asymptotic angle distribution will consist of the highest peak governed by Υ P , its smaller satellite resulting from M+ Υ (−) p , even smaller satellite of the second generation M 2 + Υ (−) p , etc. Application of M− to any of the satellites moves the state back to Υ (−) p

FIG. 10 .
FIG. 10.Different GC distributions obtained for M = 2.92 by solving the discretized Master equation (N = 10 5 grid cells) numerically for different values of T , starting from a homogeneous initial condition.

10 FIG. 11 .
FIG. 11.Examples of histograms of heights falling into the three different categories introduced in Sec.V C. Upper panels: Distributions obtained from the master equation with N = 10 5 , at (M = 2.263, T = 3.498), (M = 0.990, T = 1.811), (M = 4.052, T = 3.768) from left to right.Lower panels:The corresponding histograms of heights.From left to right categories 1 (localized, maximum at the "numerical zero," h = ∆h, see inset), 2 (maximum at the left boundary h0 of the histogram but at h > ∆h, see inset), and 3 (delocalized, maximum not at the left boundary).Here, h0 = min(W (∆θ) ) + ∆h is the height corresponding to the leftmost bin in the histogram of heights with H(hi) > 0.

FIG. 12 .
FIG.12.Example of a fractal ADF.Angle distribution for M = 2.92, T ≈ 3.729 is calculated from the ME with N = 10 7 grid cells, starting from a uniform initial condition.The upper panel shows W (θ) on the entire GC, θ ∈ [−π, π).The lower panels show progressively smaller segments of the GC, with the respective interval indicated by blue shading in the preceding panel.

FIG. 14 .
FIG. 14. (Non)ergodicity and fractality vs. localization in the cross-section of the M -T "phase diagram" with T ∈ [10 −3 , 5].All curves were obtained from discretized ME, parameters are the same as in Fig. 13.Upper panel: Ergodicity indicator of the Markov process corresponding to MN : 0 (1) indicates nonergodicity (ergodicity) in terms of connection of graph GN , Eq. (103).Middle panel: Fractal dimension d from box counting [Eq.(98)]; d = 1 corresponds to a one-dimensional curve, 1 < d < 2 gives a nontrivial Hausdorff dimension of the curve.Lower panel: Category of the position of the maximum of the height histogram [Eq.(96)]; 1 and 3 indicate localization and delocalization, respectively; 2 indicates peaks on top of the extended background.These categories establish a connection between localization and nonergodicity as well as delocalization and ergodicity, see discussion in the text.The vertical lines denote the same cases as in Fig. 13.

FIG. 15 .
FIG. 15.Diagrams illustrating localization, fractality, and ergodicity of the distributions obtained by using the ME method.All panels were generated with N = 10 4 grid cells and a maximum of 10 4 iteration steps starting from a uniform initial distribution.Upper-left panel: participation ratio, Eq. (88); large (small) values correspond to delocalization (localization).Upper-right panel: Scaling exponent ζ of the PR [Eq.(90)]; 0 (1) indicate localization (delocalization).Middle left panel: Ergodicity marker of the Markov process corresponding to the transition matrix MN , indicating whether the graph defined in Eq. (103) has a single strongly connected component; yellow (purple) regions are ergodic (nonergodic).Middle right panel: Category of the position of the maximum of the height histogram, Eq. (96); 1 and 3 indicate localization and delocalization, respectively; 2 indicates peaks on top of the extended background.Lower-left panel: Support S C N , Eq. (93), with C = 0.99.Lower-right panel: Fractal dimension d of the distributions, Eq. (98); d = 1 corresponds to a one-dimensional curve, 1 < d < 2 describes a fractal pattern.Lines show special cases described in Sec.IV.Solid lines mark period-2 cases: Y T = (2k + 1)π with k ∈ N0.Dashed lines denote frozen cases: Y T = 2kπ.Dotted lines denote shift cases: M T = kπ.Dash-dotted lines correspond to the projective limit for M− [purple, Eq. (79)] and for M+ [orange, Eq. (80)].

FIG. 16 .
FIG.16.Evaluation of nonergodicity and localization conditions.Left panel: Ergodic and nonergodic regions obtained by using the phenomenological nonergodicity condition in the continuous Markov process, Eq. (105).Purple and orange lines correspond to projecting M− and M+ matrix, respectively.The M+ matrix has a strong eigenvalue hierarchy in the lower purple region as well, see Fig.6.The color code is the same as in the middle left panel of Fig.15.Right panel: Localization condition (106) evaluated for all the invariant subspaces obtained as the leaf nodes in nonergodic graphs GN , see Sec.V and Appendix G.The condition for localization is fulfilled in the dark purple regions.
GC .Convergence (and the speed) of the attraction is determined by the ratio max v(η) GC /min v (η)GC and by the typical value of the exponent q.The former approaches unity if the eigenvectors Ῡ(η) GC T Due to the reduction rule M 2 µ = 1 the exponentially growing number of states in the generic case (cf.Fig.4) reduces to linear growth for the period-2 case.

5 and 17); 2 .
compare the time average over a single trajectory to the ME result (see Fig.19in Appendix F).

20 FIG. 20 .
FIG. 20.Examples for ADFs(22) as the stationary state of the discrete master equation (E3) (solid blue line) and the time average of a single Monte-Carlo post-measurement trajectory (23) (dotted orange line).

I 2 ≡
i∈v2 c i c i = [θ i − ∆θ/2, θ i + ∆θ/2] (G2)forms an invariant subset for the continuous post-measurement evolution on the GC: A trajectory with θ 0 ∈ I 2 can never transition into the set v 1 ,Ξ µ (I 2 ) ⊆ I 2 ⇒ θ j / ∈ i∈v1 c i ∀j ∈ N (G3)as none of the maps M± will allow for such a transition by construction of G N via MN .

4 FIG. 21 .
FIG.21.Example of a graph G5 corresponding to an ergodic Markov process described by M5.It can be verified that every node can be reached from every node, such that all nodes belong to the same SCC.

2 FIG. 22 .
FIG.22.Left panel : Example of a graph G4 that corresponds to a nonergodic Markov process with M4 and thus to a nonergodic continuous process.Nodes belonging to different SCCs are drawn in different colors.Transitions from orange to blue are possible, but not vice-versa.Right panel : The condensation of the graph shown on the left.It contains two supernodes v1 and v2, corresponding to the two SCCs of G4.There is an edge from v1 to v2 (but not vice-versa), because there are edges in G4 leading from the orange set of nodes to the blue one (but not vice-versa).The condensation is indeed acyclic as there are no circular paths.v2 is a leaf node of the condensation, because no transitions out of v2 are possible.Thus, it defines an invariant subset for the continuous process.

4 FIG. 23 .
FIG.23.Left panel : More complicated example for a graph G12 that corresponds to a nonergodic Markov process with M12 and thus to a nonergodic continuous process.Nodes belonging to different SCCs are drawn in different colors.Right panel : The condensation of the graph shown on the left.It contains four supernodes, corresponding to the four SCCs of G12.Edges of the condensation correspond to possible transitions between SCCs.The condensation is acyclic, as there are no circular paths.The condensation has a single leaf node (blue, v4) which does not have any outgoing edges.The leaf node again defines an invariant subset of the continuous process.