Absorption to Fluctuating Bunching States in Non-Unitary Boson Dynamics

We show that noisy nonunitary dynamics of bosons drives arbitrary initial states into a novel fluctuating bunching state, where all bosons occupy one time-dependent mode. We propose a concept of the noisy spectral gap, a generalization of the spectral gap in noiseless systems, and demonstrate that the exponentially fast absorption to the fluctuating bunching state takes place asymptotically. The fluctuating bunching state is unique to noisy nonunitary dynamics, with no counterpart in any unitary dynamics and nonunitary dynamics described by a time-independent generator. We also argue that the times of relaxation to the fluctuating bunching state obey a universal power law as functions of the noise parameter in generic noisy nonunitary dynamics.


I. INTRODUCTION
Long-time relaxation dynamics of a system is a central issue in non-equilibrium physics.Recent studies show that rich long-time behaviors arise in isolated, Floquet, and open systems.
Spectral analysis has been an inevitable procedure to understand noiseless linear dynamics described by a time-independent generator.Indeed, various intriguing phenomena have been uncovered through the eigenvalues and eigenmodes of the generator.
Examples of such phenomena include thermalization dynamics and its breakdown in isolated and Floquet systems in light of the eigenstate thermalization hypothesis  and drastic changes of steady states and relaxation times through spectral transitions in open systems .
In stark contrast with noiseless cases, however, the direct spectral analysis for noisy dynamics has not been fully developed.
This is because diagonalizing one generator does not bring full information for dynamics in noisy cases.While the spectral analysis for the averaged dynamics over noise is possible [56][57][58][59][60][61][62], this procedure may discard some important information, such as how the state after long time is affected by the temporal noise.
In this work, we show that noisy non-unitary dynamics of free bosons results in a novel fluctuating bunching state, where all bosons are absorbed to one timedependent mode from any initial states after long times, as schematically shown in Fig. 1.For this purpose, we propose a new analysis based on a noisy spectral gap, a generalization of the spectral gap in noiseless systems.Using non-unitary free-boson dynamics subject to random losses and postselections, we demonstrate that the noisy spectral gap takes a non-zero value.This suggests the exponentially fast absorption to the fluctuating bunching state, which does not arise in any unitary dynamics and non-unitary dynamics by time-independent generators.We also argue that the relaxation times for the absorption are universally proportional to the inverse square of the noise strength.
The rest of this paper is organized as follows.In Sec.II, we explain the general framework of the fluctuating bunching state, where the noisy spectral gap is also introduced.In Sec.III, we show that fluctuating bunching states emerge in a concrete model of an optical network with photon loss effect and postselection.In Sec.IV, we find the non-trivial power law of relaxation times and discuss that a wide range of non-unitary dynamics universally exhibit the power law.Section V is devoted to the conclusion.

II. FLUCTUATING BUNCHING STATES IN THE GENERAL FRAMEWORK
We consider general non-unitary dynamics of free bosons where pure states are mapped to pure states and the number of bosons is conserved; if we set an initial state of n bosons as |ψ t=0 ⟩ ∝ n p=1 b † x in p |0⟩, the quantum state after t-step discrete dynamics described by V t becomes [48] Here x in p is the initial position of the pth boson, |0⟩ is the vacuum state, and √ N t is the normalization factor to ensure ⟨ψ t |ψ t ⟩ = 1.We consider the non-unitary dynamical matrix V t decomposed as with {Q t } being non-unitary matrices for one step.We show that the spectrum of V t plays an important role in noisy dynamics with i.i.d.random nonunitary matrices {Q t }.To see this, assuming the diagonalizability, we consider the eigenequations of V t , Here, ϕ i t ( φi t ) is the ith right (left) eigenmode for V t with the bi-orthonormality φi † t ϕ j t = δ ij , and eigenvalues {λ t,i } are ordered as x is a creation operator for an eigenmode ϕ i t with ϕ i † t ϕ i t = 1 and Ñt = ⟨0| bn t,1 b †n t,1 |0⟩ is the normalization factor.We refer to Eq. ( 4) as the bunching state because all bosons occupy a single mode ϕ 1 t , as schematically shown in Fig. 1.We note that, if the bunching state emerges, classical computers can easily sample the probability distribution of bosons by sampling that of distinguishable particles [48].In the following, we show that bosons become distinguishable even in noisy non-unitary dynamics, while Ref. [48] explores noiseless dynamics.The classically computable distribution is contrasted to the unitary dynamics where the boson sampling problem can exhibit quantum supremacy [79][80][81][82][83][84][85][86][87][88][89][90][91][92][93][94].
Equation (4) independent of {x in p } indicates the memory-loss effect, where all initial states are absorbed into the bunching state.Such memory loss of quantum states is absent for unitary dynamics, where |λ t,i | are the same for all eigenmodes and thus |ψ t ⟩ always includes information about the initial state.
To understand the absorption into the bunching state in noisy non-unitary dynamics, we introduce the noisy spectral gap as where the overline denotes the average over random noises.Here, we assume that the sample fluctuation of ∆ t is small, ∆ | are eigenvalues of the one-step time-independent evolution matrix Q t = Q.However, there are qualitative differences between noisy and noiseless systems.One is that asymptotic convergence of ∆ t is non-trivial in noisy dynamics while it is explicitly obtained in noiseless dynamics.This is because the spectrum of the timedependent generator Q t does not give the gap in the noisy case, whereas that of the time-independent generator Q does in the noiseless case.We later demonstrate that ∆ exists and takes non-zero values for a concrete noisy model.Another feature unique to noisy non-unitary dynamics, which is distinct from the noiseless dynamics through a time-independent generator, is that the system does not reach any stationary (i.e.time-independent) state but leads to the bunching state largely fluctuating in time.This is because the dominant eigenmode ϕ 1 t depends on time, unlike the noiseless case where {ϕ i t } corresponding to eigenmodes of Q are independent of t.We note that such an absorption into the fluctuating state in quantum dynamics has not been known before.
We also explore relaxation times τ that bosonic systems take to reach the fluctuating bunching states.In noisy dynamics, where Q t involves a random variable R satisfying R = 0, we argue that relaxation times typically exhibit a power law for small β with β ∝ (R 2 ) 1 2 being the strength of the noise.This power law is discussed toward the end of the manuscript in detail.While we can evaluate relaxation times through various quantities, one natural choice is the definition from the noisy spectral gap, τ ∆ ∝ 1/∆.
Before closing this section, we briefly mention the gap discrepancy problem.In the non-unitary dynamics by a time-independent generator, due to small overlaps of left and right eigenmodes the spectral gap obtained from the generator may not give plausible predictions of relaxation times [62,[95][96][97], which is referred to as the gap discrepancy problem.In the present work, where the spectral gap is generalized from noiseless dynamics to noisy dynamics, we consider cases where there is no gap discrepancy problem and thus the relaxation time can be evaluated through the noisy spectral gap.

III. FLUCTUATING BUNCHING STATES IN A CONCRETE MODEL
We demonstrate that the fluctuating bunching state emerges in noisy non-unitary dynamics whose timeevolution matrix for one step is given by Here, the time-independent unitary matrix U is where {θ j } with j = 1, 2, 3 are real parameters, ζ represents a two-site unit-cell covered by a blue rectangle in Fig. 2 (a), and the periodic boundary condition is imposed.The time-dependent non-unitary matrix G t is where z η,t takes random values uniformly sampled from the box distribution, z η,t ∈ [−1/2, +1/2].Here, η represents a two-site unit-cell covered by an orange rectangle in Fig. 2 (a).The real parameter β determines the strength of noise leading to non-unitary dynamics, and βz η,t corresponds to R in the previous section.
Equation (1) describes the dynamics of photons in an optical network with photon loss effect and postselection.The creation operator b † x is transformed by U when photons pass through linear optical elements, such as beam splitters, phase shifters, and wave plates.Nonunitary dynamics by G t occurs when we postselect the cases where all photons remain in the system, despite the effect of photon loss.The photon loss can be artificially introduced through optical elements coupled to the environment [48].We consider situations where the optical elements have imperfections leading to noise, that is, the strength of loss randomly deviates from a mean value depending on η and t.This type of noise makes our setting distinct from conventional regimes such as fluctuational electrodynamics treating thermal noises and quantum fluctuations [98].We note that the dynamics with n conserved is distinct from the dynamics averaged over all outcomes where the number of photons can be decreased, like the dynamics governed by the Gorini-Kossakowski-Sudarshan-Lindblad equation [99].The long-time limit of the averaged dynamics becomes the trivial vacuum state with n = 0, which is contrasted to our setting where the fluctuating bunching states emerge.
Figure 2 (b) shows that ∆ t converges to a constant value ∆ in time, where the overline denotes the average over random realizations of {z η,t }.The non-zero gap ∆ implies that exponentially fast absorption to the fluctuating bunching state asymptotically occurs for this model.We note that ∆ is proportional to 1/X, which suggests that the fluctuating bunching state emerges with the timescale τ ∆ = O(X), as detailed in Appendix A. We note that the non-zero noisy spectral gap in this paper is defined for each finite-size system, and we do not consider the thermodynamic limit X → ∞.
We find evidence that ∆ gives the decay rate to the fluctuating bunching state even for each trajectory.To see this, we evaluate the long-time behavior of e t,1 −e t,2 as an approximation of ∆ at t = 10 8 .Here, e t,1 and e t,2 are the first and second instantaneous Lyapunov exponents defined as e t,i = ln(Λ t,i )/t (10) without ensemble average, where {Λ t,i } are singular values of V t , We can also capture relaxation dynamics through {Λ t,i } as well as {λ t,i }.Indeed, both |λ t,1 | ≫ |λ t,2 | and Λ t,1 ≫ Λ t,2 indicate that V t can be appoximated by the rank-1 matrix.As shown in Fig. 3, Λ t,2 /Λ t,1 exponentially decays when t is increased.Since limitations on the numerical precision and the exponential decay of Λ t,2 /Λ t,1 prevent us from directly computing Λ t,2 /Λ t,1 for large t, we obtain e t,1 − e t,2 through evaluating growth rates of randomly chosen vectors, as detailed in Appendix B. Figure 2 (c) shows that e t,1 − e t,2 converges to the same value as ∆ for all samples, which indicates |λ t,2 /λ t,1 | ∼ e −∆t for each random realization.Note that the spectral gap and Lyapunov exponents of noisy non-unitary quantum dynamics have seldom been explored, while the first Lyapunov exponents for products of transfer matrices have been studied for obtaining localization lengths in time-independent systems with random potentials [101].
As a consequence of the absorption to the fluctuating bunching state, physical observables exhibit fluctuating dynamics independent of initial states in the long run.Figure 4 shows that expectation values of an observable, x bx |ψ t ⟩ /n, approach the same fluctuating trajectory with large t for all initial states.We note that the dominant mode ϕ 1 t is spatially localized as detailed in Appendix C, where the localization position depends on time as indicated by Fig. 4. The fluctuation and localization of the bunching state mean that there is no equilibration and thermalization.

IV. POWER LAW OF RELAXATION TIMES
We argue that relaxation times at which the bosonic system reaches the fluctuating bunching state are universally given by τ ∝ β −2 for small noise strength β.Our argument follows from two steps.First, we introduce four versions of relaxation times, τ ∆ , τ λ , τ Λ , and τ x , and numerically show that they exhibit similar behavior, in particular, the power law in Eq. (6).Second, we analytically evaluate τ λ and τ Λ instead of τ ∆ , for which the analytical estimation is difficult.In addition, we make plausible arguments for the universality of the power law τ ∝ β −2 .
As the first step, we define relaxation times as the smallest time steps at which f t ≤ ln(c) is satisfied, where c ≪ 1 is a threshold and a function f t is defined through a quantity related to |λ t,2 /λ t,1 |, Λ t,2 /Λ t,1 , or observables.On the basis of the noisy spectral gap, we define the relaxation time τ ∆ with f t = −∆t, which leads to τ ∆ = | ln(c)/∆|.Figure 3 shows that other two quantities, − ln[|λ t,1 /λ t,2 | 2 ]/2 and ln (Λ t,2 /Λ t,1 ), also decay in time.In particular, behaviors of −∆t and ln(Λ t,2 /Λ t,2 ) are quite similar.Then, we define relaxation times τ λ and τ Λ through f t = − ln |λ t,1 /λ t,2 | 2 /2 and f t = ln(Λ t,2 /Λ t,1 ), respectively.In Fig. 5, we discover that τ λ , τ Λ , and τ ∆ all show the same β dependence, τ ∝ β −2 , for sufficiently small β.In particular, τ Λ and τ ∆ are almost overlapped.We also find that the same power law appears for an additional timescale τ x , which is obtained from Here, the subscripts for x 2 t represent the different initial states.Relaxation times as functions of the noise parameter, which exhibit the power law in Eq. ( 6).Green, blue, red, and purple symbols represent τ∆, τΛ, τ λ , and τx, respectively, with c = 10 −6 and X = 20.Here, the average in ft is taken over 10 2 samples for obtaining one τ , and then we iterate it 10 2 times and obtain the averaged τ over them.Expectation values are computed with n = 2, where initial states are (x in 1 , x in 2 ) = (−5, 5) for x 2 t,a and (−1, 0) for x 2 t,b .The black symbols and orange broken line correspond to τ λ Ω and τ Λ Ω , in Eqs. ( 13) and ( 16), respectively, where the latter is exactly proportional to β −2 .
As the second step, we semi-analytically give an upper bound of τ λ .To this end, we focus on As detailed in Appendix E, we find a lower bound which is satisfied asymptotically.Here, µ and ν are respectively the largest eigenvalues of  (12). Figure 5 indicates the power law of τ λ Ω , This is because ) for small β with U and F (U and F) being some complicated but β-independent matrices whose explicit forms are written in Appendix E. Since ln(| √ ν/µ|) should be zero for β = 0, we believe that the subleading terms β 2 of Q and Q typically lead to ln(| √ ν/µ|) ∝ β 2 .While τ λ Ω only gives the upper bound of τ λ , Fig. 5 shows that both timescales are proportional to β −2 in this model.
We also analytically obtain an approximated value of τ Λ using a perturbative analysis.To this end, we focus on satisfied.Thus, τ Λ can be evaluated through ln(Ω Λ t ) ≤ ln(2c 2 ).When β is small, we can neglect higher-order terms of β, which results in G t ≃ η (σ 0 + βz η,t σ 1 + β 2 z 2 η,t σ 0 /2), where σ 0 and σ 1 are the Pauli matrices.Thus, we can approximate Q t as where As detailed in Appendix F, this expansion leads to ln which is valid when β 2 t/X is small.The orange broken line in Fig. 5 corresponds to defined by is parallel to τ Λ , τ λ , and τ x , which is consistent with our argument that all relaxation times satisfy τ ∝ β −2 .We argue that the power law in Eq. ( 6) for small β emerges in a wide range of noisy non-unitary dynamics, on the basis of the following two assumptions satisfied in our model.We first assume that ln(| √ ν/µ|) ∝ β 2 and thus τ λ Ω ∝ β −2 , which is plausible because subleading order terms of Q and Q become β 2 when the average of noise is zero, as detailed in Appendix F. We also assume that various relaxation times such as τ λ , τ λ Ω , τ ∆ , and τ Λ Ω exhibit the same dependence on the noise parameter β, as confirmed in the concrete model.In addition, details of models have no effect on the power law τ Λ Ω ∝ β −2 , as long as the average of the noise is zero.We indeed discover the same power law τ ∝ β −2 for a model proposed in Ref.
[102], which is different from that in the present work, as shown in Appendix G.

V. CONCLUSION
We have shown that the novel fluctuating bunching state, where all bosons occupy one time-dependent mode independent of initial states, appears in noisy non-unitary dynamics.This is analyzed through the spectrum of noisy non-unitary time-evolution matrices.In particular, the noisy spectral gap ∆, the extension of the spectral gap in noiseless systems, plays a key role.We have also argued that times of relaxation to the fluctuating bunching state universally exhibit the power law as functions of the noise parameter, τ ∝ β −2 .
The spectral analysis can be applied to general noisy dynamics, not restricted to the dynamics of bosons explored in this work.For example, our analysis should be applicable to the dynamics of fermions, while fermions do not occupy one dominant mode due to the Pauli exclusion rule.In addition, it may be possible to analyze transitions of relaxation times, such as purification transitions in monitored quantum many-body systems [70], through the noisy spectral gap ∆ in the present work. where , and I X is the X × X identity matrix.Here, Π s = [w 1 s ][w 1 s ] † is the projection matrix onto the direction of w 1 s , and thus w 2 s is orthogonal to w 1 s .The Lyapunov exponents are obtained through 2 as a function of the system size X.The IPR is almost independent of X, which indicates that ϕ 1 t is localized.The average is taken over 40 samples and 100 time steps after a time at which |λt,2/λt,1| < c = 10 −6 is satisfied.(b) The probability distribution Pt({x out q }) of a bunching state with n = 3, (x in 1 , x in 2 , x in 3 ) = (−6, 1, 8), β = 0.3, and t = 42418.The value of Pt({x out q }) is represented as the size and color of the symbols: the big (tiny) and dense green (shallow gray) symbols correspond to large (small) Pt({x out q }).
the largest and second-largest Lyapunov exponents.This is because increasing s makes w 1 s (w 2 s ) approach the first (second) Lyapunov vector, the eigenmode of lim t→∞ Ṽt corresponding to the largest (second-largest) eigenvalue.Since the vectors w 1 s and w 2 s are normalized for every M step, we can make t much larger than a time step at which Λ t,2 /Λ t,1 becomes smaller than the numerical precision.Indeed, numerical results with M = 10 4 and T = 10 4 in Fig. 2 (c) give a plausible value of the noisy spectral gap for t = 10 8 , while the directly computed Λ t,2 /Λ t,1 becomes smaller than our numerical accuracy 10 −16 in a range 10 4 ≤ t ≤ 10 5 , as indicated by Fig. 3.

Appendix C: Localization of the bunching state
We find that the occupied state ϕ 1 t is spatially localized in the model explored in the main text; the average of the inverse participation ratio (IPR) does not decrease monotonically when we increase the system size X, as shown in Fig. 7 (a).Figure 7 (b) shows the probability distribution [48]  = n is the number of bosons at x regarding the input (output) configuration of bosons, and Per[W t ] is the permanent of an n-by-n matrix W t defined as the submatrix of V t , [W t ] qp = [V t ] x out q x in p .Here, x in(out) q with q = 1, 2, • • • , n is the position of the qth input (output) boson.If the fluctuating bunching state is realized, initial states have no effect on the probability distribution, and thus P t ({x out q }) is independent of {x in p }.We can see that the three bosons gather around one position due to the localization of ϕ 1 t .Note that the localization position depends on t, as indicated by Fig. 4 in the main text.

Appendix D: Power law of various relaxation times
As discussed in the main text, we conjecture that various types of possible relaxation times for the absorption exhibit the same power-law dependence on β.In this section, we numerically confirm this conjecture for the model introduced in the main text.
We define a relaxation time at which the fluctuating bunching state is realized as the smallest time at which f τ ≤ ln(c) is satisfied with c ≪ 1.In the main text, we explore four relaxation times τ λ , τ Λ , τ ∆ , and τ x respectively , where the number of bosons is n = 2 and subscripts a and b represent different initial states.In Fig. 8, these relaxation times correspond to red, blue, green, and purple symbols, respectively.In addition to these relaxation times explored in the main text, Fig. 8  respectively correspond to cyan and grey symbols.From Fig. 8, we can understand that all relaxation times exhibit the same power law τ ∝ β −2 ; the orange broken line τ Λ Ω = 3X|2 ln(c) + ln(2)|β −2 , derived in Sec.F, is exactly proportional to β −2 , and all relaxation times are parallel to this line.From Fig. 8, we can also understand that behaviors of the eigenvalues and singular values are quite similar; relaxation times based on −∆t (green symbols) and ln (Λ t,2 /Λ t,1 ) (blue symbols) are almost identical.and ft = − ln |Λt,1/Λt,2| 2 /2.The average included in ft is taken over 10 2 samples for obtaining one τ .Then, we reiterate obtaining τ for 10 2 times and take the average over them.The subscripts a and b for x 2 t correspond to two initial states of bosons, (x in 1 , x in 2 ) = (−5, 5) and (−1, 0).Regarding τ∆ with ft = −∆t, ∆ is obtained through the ensemble average over 10 4 samples.The orange broken line corresponds to the relaxation time τ Λ Ω defined through the perturbative analysis on ln (Λt,2/Λt,1).Since this relaxation time is obtained as −β 2 τ Λ Ω /6X − ln( √ 2) = ln(c), the orange broken line is exactly proportional to β −2 .The black symbols represent τ λ Ω = ln(c)/ ln (| √ ν/µ|), which is above the red symbols for large t, due to the inequality obtained in Eq. (E4).The parameters are (θ1, θ2, θ3) = (0.37π, 0.19π, 0.25π) with the system size being X = 20.

Appendix E: analysis through the upper bound on τ λ
We evaluate the relaxation time τ λ based on eigenvalues of V t , where | is smaller than the threshold value c ≪ 1, the bunching state is realized within the precision c.
As discussed in the main text, the timescale τ λ , which is obtained from , has an analytical upper bound τ λ Ω , which typically shows β −2 decay.In this section, we show this fact in detail.

General bound
In this subsection, we give discussions applicable to a wide range of models that exhibit noisy non-unitary dynamics, not restricted to the model explored in the main text.Then, in the next subsection, we demonstrate the concrete application of our general discussions for the model in the main text.
To discuss the decay of |λ t,2 /λ t,1 |, we focus on Ω λ t defined as Below, we evaluate the left-hand side of Eq. (E4).To this end, we use tr(A)tr(B) = tr(A ⊗ B), where the trace in the right-hand side is carried out in H ⊗2 with H ⊗m being the extended space whose dimension is X m .Then, |tr(V t )| 2 can be written as where t is a time-independent matrix that acts on the extended space H ⊗2 .The largest eigenvalue of Q, which we write as µ, determines the growth rate of |tr(V t )| 2 .In the same way, |tr where ) is a time-independent X 4 × X 4 matrix and the trace in the right-hand side is taken in H ⊗4 .For evaluating the other terms in |tr( where S is the swap operator acting on H ⊗2 .Thus, we obtain where where We notice that S and Q commute each other, which means that they have the same eigenmodes.The eigenvalues of S, denoted as ρ, take 0 or 4 since eigenvalues of S are ±1 and all ingredients of S such as S ⊗ S and I X 2 ⊗ S commute with each other.Therefore, we focus only on eigenmodes with ρ = 4 and neglect the other eigenmodes with ρ = 0, which do not contribute to dynamics of |tr(V t ) 2 − tr(V 2 t )| 2 .Then, the largest eigenvalue of QS/4, which we write as ν, determines the growth rate of |tr( t )| 2 is determined through |µ 2 /ν| t /4 in the long-time limit.We assume that non-unitary matrices {Q t } are determined through a random variable R with R = 0 and that the system is characterized by the noise parameter β ∝ R 2 , where the overline represents the average over the random realizations of R. Here, β is the strength of randomness, and the second moment of R/β is independent of β.We note that R corresponds to βz η,t in the case of the model studied in the main text.
For small β, we can expand Q t as where B t and C t are random matrices.Here, B t = 0 is satisfied because B t only includes the first-order terms of R/β.We also assume that A t is independent of R and A † t = A −1 t such that the system exhibits unitary dynamics if β = 0. We note that other randomness independent of R can be included in dynamics.For example, {A t } in Eq. (E11) can be random unitary matrices determined through a random variable R ′ that is independent of R. The average represented by the overline is over all kinds of randomness in noisy dynamics.Due to B t = 0, Eq. (E11) leads to where We argue that these Taylor expansions of averaged matrices Q and Q in Eq. (E12) typically make subleading order terms of µ, ν be β 2 in a wide range of noisy non-unitary dynamics.In addition, we assume that the leading terms of |ν| and |µ| are 1, resulting in ln(| √ ν/µ|) = 0 with β = 0.This assumption is valid since [|tr(V t )| 2 ] 2 and |tr(V t ) 2 − tr(V 2 t )| 2 should not decay in unitary dynamics.Thus, Eq. (E12) suggests that the decay rate of [|tr(V t )| 2 ] 2 /|tr(V t ) 2 − tr(V 2 t )| 2 typically exhibits the β 2 dependence with respect to the noise parameter, which gives the lower bound of |λ t,1 /λ t,2 | 2 in the long run through the Cauchy-Schwarz inequality in Eq. (E4).We also stress that the analysis explained above is applicable to the long-time regime as well as the short-time regime.Thus, τ λ Ω defined through f t = t ln(| √ ν/µ|) gives the upper bound of τ λ .While there is a possibility that τ λ decays faster than τ λ Ω , we believe that the τ λ Ω and τ λ typically exhibit the same power law τ ∝ β −2 , which is confirmed in the model studied in the main text.

Example
We can actually show that Q and Q are written in the form of Eq. (E12) regarding the model explored in the main text.In the model, the time-evolution matrix for one step is defined as Here, the matrices U and G t are defined as where z η,t is chosen from the box distribution See Fig. 2 (a) in the main text.From Eqs. (E15)-(E18), we can understand that this model satisfies the assumptions that A t = U is unitary and B t = 0 with B t = ( η z η,t σ 1 )U .Equations (E16) and (E17) lead to where In this model, we can write down Q and Q through calculating G and G explicitly, not relying on the perturbative expansion of Q t in Eq. (E11).To this end, we decompose G t as where G d t is a diagonal matrix and H is the Hadamard matrix The non-unitary matrix in Eq. (E21) describes the dynamics of photons in optical networks exposed to loss effects and postselections [48], and such dynamics of a photon has been observed experimentally [38].Equations (E20) and (E21) result in with G d and G d being diagonal matrices where H = H ⊗ H and H = H ⊗ H.The diagonal matrix G m athrmd is written as where where g(β) = sinh(β/2) β/2 and I 4 is the 4 × 4 identity matrix.When The diagonal matrix G d is written as where G(η 1 , η 2 , η 3 , η 4 ) is a 16 × 16 matrix.When all of η i with i where I 16 is the 16 × 16 identity matrix.When where G(η 1 , η 1 ) and G(η 3 , η 4 ) respectively correspond to Eqs. (E26) and (E25).We can obtain G(η 1 , η 2 , η 3 , η 4 ) regarding other combinations such as η 2 = η 4 , through the corresponding permutation for the matrix elements of Eq. (E29).When η 1 = η 2 , η 2 ̸ = η 3 , and where G(η 1 , η 1 ) corresponds to Eq. (E26).We can obtain G(η 1 , η 2 , η 3 , η 4 ) regarding the other combinations, such as η 1 = η 3 and η 2 = η 4 , through the permutation for the matrix elements of Eq. (E30).When where I 2 and 0 2 are the 2 × 2 identity matrix and the null matrix, respectively.We can obtain G(η 1 , η 2 , η 3 , η 4 ) regarding the other combinations, such as η 2 = η 3 = η 4 , through the corresponding permutation for the matrix elements of Eq. (E31).When Since g(β) ≃ 1 + β 2 /24 for small β, Eqs.(E24)-(E32) lead to Eq. ( E12) where E = U and E = U.Here, F (F) is a complicated matrix obtained from the Taylor expansion of G (G) and the multiplication of H and U (H and U).

Appendix F: Perturbation analysis on the ratio of singular values
As discussed in the main text, we can also evaluate the relaxation time τ Λ through the singular values {Λ t,i } of V t in Eq. (11), where Λ t,1 ≥ Λ t,2 • • • ≥ Λ t,X .Here, we discuss the perturbation analysis for evaluating τ Λ .
If Λ t,1 ≫ Λ t,2 is satisfied, the fluctuating bunching state emerges.We perturbatively calculate which approaches 2(Λ t,2 /Λ t,1 ) 2 in the long run provided that Λ t,1 ≫ Λ t,2 ≫ Λ t,3 is satisfied.We again expand Q s as for small β, where B s and C s are random matrices.Remember that B s = 0 and A † s = A −1 s are assumed.We note that the analysis explained below is not restricted to the model explored in the main text and thus is applicable to a wide range of noisy non-unitary dynamics, as long as the assumptions above are satisfied.Up to the second-order terms of β, V t can be expanded as where A s2,s1 = A s2 A s2−1 • • • A s1+1 A s1 .For consistency of the notation, we define A 0,1 = A t,t+1 = I X with I X being the X × X identity matrix.In Eq. (F3), second order terms of β like β 2 A t,s2+1 B s2 A s2−1,s1+1 B s1 A s1−1,1 with s 2 − 1 > s 1 + 1 are ignored.This is because we carry out the average over randomness at the end of the analysis, and such terms vanish due to B s = 0 with B s1 and B s2 being independently distributed for s We note that (A † s B s ) 2 , (B † s A s ) 2 , and B † s B s do not depend on s.Equation (F8) is valid in the short-time regime for which β 2 t/X ≪ 1 is satisfied, since tr 2B † s B s + (A † s B s ) 2 + (B † s A s ) 2 should be proportional to the system size X.However, the power law of the relaxation time τ Λ Ω ∝ β −2 obtained from Eq. (F8) is consistent with actual relaxation time τ Λ obtained through numerical simulations, which can be confirmed regarding two models explored in the main text and in Sec.G.For the model explored in the main text, we have A s = U , B s = Z s U , and C s = Z 2 s U/2, where Z s = η z η,s σ 1 with σ 1 being one of the Pauli matrices.Then, the trace in Eq. (F8) becomes tr 2B † s B s + (A † s B s ) 2 + (B † s A s ) 2 = 4tr Z 2 s = X/3, which leads to Eq. ( 15).

FIG. 1 .
FIG. 1. Schematic picture of the fluctuating bunching state.Various initial states, |ψ a t=0 ⟩ , |ψ b t=0 ⟩ , • • • , are absorbed into the fluctuating bunching state, where all bosons described by yellow circles occupy one dominant eigenmode of Vt.We can define the time-dependent gap ∆t from the eigenvalues of ln(Vt)/t.

FIG. 2 .
FIG. 2. (a) Schematic picture of the non-unitary dynamics of photons in an optical network.Photons evolve via the unitary matrix U (blue rectangles) and then experience non-unitary dynamics by Gt (orange rectangles) due to the photon loss effect and postselection.(b) The convergence of ∆t into the noisy spectral gap ∆ ≃ 6 × 10 −4 , where the number of samples is 10 5 .(c) The difference between the first and second Lyapunov exponents.Symbols with different colors correspond to different realizations of noises {zη,t}.For large t, et,1 − et,2 for each trajectory converges to the sample-independent value same as ∆.In (b) and (c), β = 0.3 and X = 20.
FIG.5.Relaxation times as functions of the noise parameter, which exhibit the power law in Eq. (6).Green, blue, red, and purple symbols represent τ∆, τΛ, τ λ , and τx, respectively, with c = 10 −6 and X = 20.Here, the average in ft is taken over 10 2 samples for obtaining one τ , and then we iterate it 10 2 times and obtain the averaged τ over them.Expectation values are computed with n = 2, where initial states are (x in 1 , x in 2 ) = (−5, 5) for x 2 t,a and (−1, 0) for x 2 t,b .The black symbols and orange broken line correspond to τ λ Ω and τ Λ Ω , in Eqs.(13) and (16), respectively, where the latter is exactly proportional to β −2 .

B3) e 1 +FIG. 7 .
FIG. 7. (a) The averaged inverse participation ratio (IPR) of the dominant mode x |ϕ 1t (x)| 4 /[ x |ϕ 1 t (x)| 2 ]2 as a function of the system size X.The IPR is almost independent of X, which indicates that ϕ 1 t is localized.The average is taken over 40 samples and 100 time steps after a time at which |λt,2/λt,1| < c = 10 −6 is satisfied.(b) The probability distribution Pt({x out q }) of a bunching state with n = 3, (x in 1 , x in 2 , x in 3 ) = (−6, 1, 8), β = 0.3, and t = 42418.The value of Pt({x out q }) is represented as the size and color of the symbols: the big (tiny) and dense green (shallow gray) symbols correspond to large (small) Pt({x out q }).