{'fraction_non_alphanumeric': 0.038716229941092827, 'fraction_numerical': 0.014422100345317896, 'mean_word_length': 4.194344798480692, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Statistics of wave interactions in nonlinear disordered systems

We study the properties of mode-mode interactions for waves propagating in nonlinear disordered one-dimensional systems. We focus on i) the localization volume of a mode which defines the number of interacting partner modes, ii) the overlap integrals which determine the interaction strength, iii) the average spacing between eigenvalues of interacting modes, which sets a scale for the nonlinearity strength, and iv) resonance probabilities of interacting modes. Our results are discussed in the light of recent studies on spreading of wave packets in disordered nonlinear systems, and are related to the quantum many body problem in a random chain.


I. INTRODUCTION
In the absence of nonlinearity (or many-body interactions in quantum systems) all eigenstates in onedimensional random lattices with disorder are spatially localized. This is Anderson localization [1], which has been discovered fifty years ago in disordered crystals as a localization of electronic wavefunction. It can be interpreted as an interference effect between multiple scatterings of the electron on random defects of the potential. Recent experiments on the observation of Anderson localization were performed with light propagation in spatially random optical media [2,3], with noninteracting Bose-Einstein condensates expanding in random optical potentials [4,5], and with wave localization in a microwave cavity filled with randomly distributed scatterers [6].
In many situations nonlinear terms in the wave equations (respectively, many body interaction terms in quantum systems) have to be included. Thus, a fundamental question which has attracted the attention of many researchers is what happens to an initial excitation of arbitrary shape in a nonlinear disordered lattice. Nonlinearity renormalizes excitation frequencies, thereby inducing interaction between NMs. Numerical studies show that wave packets spread subdiffusively and Anderson localization is destroyed [7][8][9][10]. In the regime of strong nonlinearity, far from where it can be treated perturbatively, new localization effects of selftrapping occur [11]. A theoretical explanation of the subdiffusive spreading was offered in Refs. [7,9,12]. It is based on the fact that the considered models are in general nonintegrable. Therefore deterministic chaos will lead to an incoherent spreading. Estimates of the excitation transfer rate across the packet tail are obtained by calculating probabilities of mode-mode resonances inside the packet. Some predictions of this approach include the effect of different degrees of nonlinearity and were successfully tested in [13].
In this work we study the statistical properties of mode-mode interactions. We focus on i) the localization volume of a mode which defines the number of interacting partner modes, ii) the overlap integrals which determine the interaction strength, iii) the average spacing between eigenvalues of interacting modes which sets a scale for the nonlinearity strength, and iv) resonance probabilities of interacting modes. We discuss the results in the light of recent studies [7][8][9][11][12][13] on spreading of wave packets in disordered nonlinear systems, and relate our findings to the quantum two interacting particle problem in a random chain.

II. NONLINEAR SCHRÖDINGER CHAIN
We consider the disordered discrete nonlinear Schrödinger equation (DNLS) with the Hamiltonian Here ψ l are complex variables, l are lattice site indices and β ≥ 0 is the nonlinearity strength. The random on-site energies ǫ l are chosen uniformly from the interval − W 2 , W 2 , with W denoting the disorder strength. The equations of motion are generated byψ l = ∂H D /∂(iψ ⋆ l ): Eq. (2) conserves the energy (1) and the norm S = l |ψ l | 2 . Varying the norm of an initial wave packet is strictly equivalent to varying β, therefore we choose S = 1. Note that Eq. (2) is used to qualitatively describe the evolution of a dilute Bose-Einstein condensate trapped into a deep periodic potential [4], and also the evolution of a light wave in disordered one-dimensional waveguide lattices with cubic Kerr-type nonlinearity (under the paraxial approximation) [3].
For β = 0, Eq. (1) with ψ l = A l exp(−iλt) reduces to the eigenvalue problem The width of the eigenfrequency spectrum λ ν of (3) is The normalized eigenvectors A ν,l ( l A 2 ν,l = 1) are the NMs, and the eigenvalues λ ν are the frequencies of the NMs. We order the NMs in space by increasing value of the center-ofnorm coordinate X ν = l lA 2 ν,l . The equations of motion of (1) in normal mode space read with the overlap integrals The variables φ ν determine the complex time-dependent amplitudes of the NMs.

III. PROPERTIES OF NORMAL MODES
A. Localization length, volume and participation number The asymptotic spatial decay of an eigenvector is given by A ν,l ∼ e −l/ξν , where ξ ν is the localization length of a mode ν with the eigenvalue λ ν . We calculate the average ξ ν at a given energy using the standard transfer matrix approach [14] and show the results in Fig. 1. As expected, the most extended modes correspond to the bandwidth center with ξ(λ = 0, W ) ≈ 100/W 2 for W ≤ 4 [14] (see Fig. 2). In what follows we refer only to the localization length near the bandwidth center. We also observe a small peak at λ = 0 for W ≤ 1. The smaller the disorder strength, the more pronounced the peak is. For instance 1 10 100 W for W = 0.5, the additional peak height is about 8% of the total value. Its magnitude will not exceed roughly 10% for W → 0, as first discussed in Refs. [15]. The origin of this anomaly is the deviation from single parameter scaling due to the symmetry A l (λ) = (−1) l A l (−λ) at W = 0 [16].
Next we estimate the number of NMs which interact with a chosen mode ν. This is assumed to be equivalent to estimating the number of sites, where the norm density |A ν,l | 2 of the eigenvector is not exponentially small. This number is coined localization volume, V ν . V ν is related to the localization length ξ ν , though quantitatively the two quantities might differ. We consider two ways of estimating this number. A widely used quantity is the participation number p ν = 1/ l A 4 ν,l . It is a measure of the inhomogeneity of the distribution of eigenvector amplitudes in real space. Another quantity is the effective distance between the exponential tails of the eigenvector, which is given by 12m is the second moment of the norm density distribution. Both quantities yield the exact width of a flat and compactly distributed norm density distribution. However when fluctuations are included, p ν will be reduced and underestimates the correct volume, while the effective distance does not. Therefore, we use V ν = 12m (ν) 2 + 1 as a measure of the localization volume.
We calculate numerically the average localization volume V = V ν and participation number p = p ν of NMs as a function of W (see Fig. 2). For this purpose we fix the strength of disorder W and take a chain which is much longer than ξ(W ). We calculate V and p, taking into account only those modes whose eigenvalues are located near the bandwidth center and which are not close to the boundaries of a lattice. Then, we take another re-alization and repeat the procedure. Finally, we perform the averaging with respect to different disorder realizations. We find that the localization volume V scales on average as 3.3ξ for weak disorder (see the dashed lines on Fig. 2), and tends to V = 1 in the limit of strong disorder. We also note that the participation number is almost identical to the localization length for weak disorder, and therefore misses the localization volume by a factor of three.
Our numerical results indicate that multi-humped NMs are rare and do not affect the statistical results on V. These multi-humped NMs are localized over a number of lattice sites which can be at far distance from each other in real space (see e.g. [17]). If they were statistically relevant, the second moments would be overestimated. However, we find that in the limit of strong disorder W → ∞ the localization volume V → 1, therefore multi-humped states do not significantly contribute (although we admit that there might be a measurable contribution in this limit of strong disorder). In the limit of weak disorder we observe that V, P, ξ scale in the same way with the disorder strength, therefore we can exclude any significant statistical relevance of multi-humped NMs in this regime.

B. Overlap integrals
We study statistical properties of the absolute values of the overlap integrals (5) perturbatively for weak disorder, and numerically using two different methods. In particular we aim at estimating the average absolute value of these overlap integrals for NMs which are interacting with each other within the range of one localization volume, in order to exclude statistically irrelevant exponentially weak interactions of distant NMs. Note, that in following notations, the absolute value is omitted for the sake of simplicity. In order to avoid multiple repetitions in (5) we use ν 3 ≥ ν 2 ≥ ν 1 ≥ ν.
Perturbative calculations. Let us consider a chain with finite size N and fixed boundary conditions: yields eigenvectors A q,l = 2 N +1 s ql . Eq. (6) transforms to with the coupling which mixes standing waves with the eigenvalues λ q = 2 cos (πq/(N + 1)) in the presence of disorder. The overlap integral at W = 0 will be zero for all combinations of indices except if a selection rule is satisfied [18]. It is enough to replace this rule byq 4 = ±q 1 ± q 2 ± q 3 . In short we will denote byq 4 a mode number which satisfies the selection rule for a given triplet of mode numbers (q 1 , q 2 , q 3 ). The selection rule applies to N 3 overlap integrals I 0 ∼ 1/N . The other N 4 integrals I 1 = 0. Therefore, the average overlap integral becomes I (W = 0) ∼ 1/N 2 . Let us estimate the corrections to this average when disorder is added. We first consider integrals I 1 which were strictly zero at the limit W = 0. We perform a perturbation calculation (W small) for a mode q 4 such that Straightforward calculation gives (see also [19]) Assuming now a triplet of modes (q 1 , q 2 , q 3 ) is given, and that q 4 =q 4 , the first order nonzero correction to the corresponding overlap integral reads We started with q 4 =q 4 , but in the presence of disorder the mode with number q =q 4 will become excited. Therefore, after summation over l in (12) we find Note that the indices (q 1 , q 2 , q 3 ) are implictely hidden in the quantityq 4 . In order to estimate the average, we have to take the absolute value of (13), to sum over each index q i , i = 1, 2, 3, 4 and each time to divide by N . Let us perform the averaging over q 4 . The denominator λ q4 − λq 4 will become of the order of 1/N when q 4 is close toq 4 . Replacing the sum by an integral, we estimate Since the disorder average K p,q = 0 and its variance is finite (i.e. not depending on N ) the final averaging over q 1 , q 2 , q 3 yields The overlap integrals I 0 for q 4 =q 4 were of the order of 1/N for W = 0. It is straightforward to obtain that the disorder induced correction will be of the same order as (15), which is still smaller than the unperturbed value. Thus the average value of I up to the first order of perturbation in W is given by (remember that κ = 2W √ N +1 ) where a is some constant independent on the system's parameters. One can conclude from Eq. (16), that for small enough W such that W < ( √ N ln(N )) −1 , the first term prevails and the total average integral is I ∼ 1/N 2 . In the opposite case, when W > ( √ N ln(N )) −1 , the second term in Eq. (16) dominates and, as a result, we get I ∼ W N −3/2 ln(N ).
Note that the perturbed eigenvectors A ql = A For W → ∞ the NM eigenvectors become single site profile, and the overlap integrals tend to zero. Therefore, for a given size N , the average overlap integral will start to increase with W for W > ( √ N ln(N )) −1 , reach a maximum at W max , and decay down to zero for infinitely strong disorder. It is reasonable to assume that the localization volume ξ(W max ) ∼ N . In that case for small values of W we obtain which is an estimate of the interaction strength of NMs within the spatial range of one localization volume. Numerical calculations. Method I. We fix a chain size N and calculate the average value of the overlap integrals taking into account all integrals (5). Then, we take another realization and repeat the procedure. Finally, we perform the averaging with respect to different disorder realizations. Each averaged integral I is a function of W . As derived in the above perturbation approach, it has a maximum value I (W max ) at a certain W max , as is shown in the inset of Fig. 3 for N = 40. Now we vary the chain size N , and repeat the procedure. In Fig. 3 we plot the maximum values of I (W max ) as a function of W max (red curve). We find that ξ(W max )/N ≈ 8/3, as expected in the above perturbation approach. We also find that for large N the data can be fitted with the power law I ≈ 0.0034W α , with α = 3.40 ± 0.02 (see Fig. 3). The fit was done using different number of numerical points (from 3 to 14) starting from the smallest W. In all cases the RMS relative error was better than 10 −3 . We expect that this method will overestimate the corresponding prefactors. This is due to the fact that N ≈ 0.38ξ(W max ) and therefore states overlap more strongly than in an extended system, as seen in the next method.
Method II. We fix the strength of disorder W and choose a chain size N ≫ ξ(W ). We select a middle part of a smaller size (core) of the width L, and do not consider the edges in order to avoid boundary effects. We use only modes within the core (ν = 1..L). For each mode ν we calculate its localization volume V ν . Now we consider only NMs which happen to reside in a corresponding neighborhood, i.e. we select V ν /2 modes from the right and left (in case V ν is odd, one mode is randomly taken from left or right in addition). Therefore, we have defined a subset of NMs which interact with the νth NM. We calculate all overlap integrals for this subset. Then, we move on to the next reference mode from the core. This procedure is performed for all NMs from the core, for many realizations. For small W the data can be fitted with the power law I ≈ 3.84 · 10 −5 W 3.4 (see Fig. 3) . Note that both methods yield the same exponents. Note also that we lack more data to distinguish between the numerically found law W 3.4 and the perturbation result −W 4 ln(W ).

IV. FREQUENCY SCALES
There are two frequency scales set by the linear equations (3): the average spacing d of NMs within the range of a localization volume and the width of the spectrum ∆ [7,9,12]. The two scales d ≤ ∆ determine the packet evolution details in the presence of nonlinearity. In order to calculate the average spacing and its distributions numerically, we fix the strength of disorder W and take a chain which is much longer than V (W ). We select a middle part of a smaller size (core), and do not consider the edges. For each mode ν within a core we form its subspace which consists of those modes which live in its localization volume V ν [see Sec. III A,III B for details].
We take the eigenvalues of these modes (including the eigenvalue of ν-th mode), sort them and compute absolute values of spacings between them. Then, we proceed to the next reference mode from the core. This procedure is performed for all NMs from the core and for many realizations, such that we end up with a large number of spacings (usually of the order of 10 6 ) [20].
Typical probability density functions (PDF) of the spacings s are shown in Fig. 4. For strong disorder W ≫ 1 the relative contribution of small spacings to the PDF becomes smaller. The reason is that the localization volume tends to one and for each reference mode ν we take into account only a single neighboring mode. As a result, spacings between eigenvalues increase.
We also note that the computed PDFs are far from following a Wigner-Dyson distribution Especially for small W , such a distribution could be expected due to large localization lengths. However we find systematic deviations towards a Poisson distribution with an enhancement of the probability density at small spacings. This is due to the fact that NMs overlap in general only partially in real space. The Wigner-Dyson distribution is recovered only for very short chains, when N < ξ(W ) [see Fig. 4]. In this case, all eigenmodes occupy the same volume, and level repulsion is recovered as expected. Interestingly, similar level repulsion occurs in short resonators in one-dimensional random lasers. In that case the PDF of spacings between frequencies of the neighboring lasing modes tends also to the Wigner-Dyson distribution [21]. In Fig. 5 we plot the result for the average spacing d. In the limit of small W we estimate the average spacing as d = ∆/V ∝ W 2 . For large disorder strength W ≫ 1 the localization volume tends to one. Thus, only two modes form a subspace of each reference mode ν (one of which is the mode ν itself). Therefore, the spacing can be calculated by considering two numbers (emulating two corresponding eigenvalues) which are randomly distributed within the width of the spectrum ∆. The average distance between these numbers x and y, assuming that x ≥ y is It follows d = ∆/3. In Fig. 5 the two theoretical estimates are shown to be close to the numerical data.

V. RESONANCES
When a wave packet spreads, its size grows, and the norm density inside the packet drops. Previous studies indicate that this process of spreading is due to resonances in mode-mode interactions. These resonances lead to chaotic dynamics inside the packet, and to a subsequent incoherent spreading. Let us estimate the number of resonant modes in the packet. Excluding secular interactions, the amplitude of a NM with |φ ν | 2 = n ν is modified by a set of three other modes µ ≡ (µ 1 , µ 2 , µ 3 )  in first order in β as (4) (see [9]) The perturbation approach breaks down and resonances set in when √ n ν < |φ (1) ν |. Since all considered NMs belong to the packet, we assume their norms to be equal to n.
We perform a statistical numerical analysis by computing the PDF of R ν, µ . For a given NM ν we obtain R ν, µ0 = min µ R ν, µ . Collecting R ν, µ0 for many ν and many disorder realizations, we find the probability density distribution W(R ν, µ0 ). We also analyze separate contributions from three different types of interactions, namely from quadruplets (all four modes are different), triplets (only three of four modes are different) and pairs (only two different modes participate in the interaction). For quadruplets all indices in (21) should be different, i.e. µ 1 = ν, µ 2 = ν, µ 3 = ν, µ 1 = µ 2 , µ 1 = ν, µ 1 = µ 3 , µ 2 = µ 3 . For triplets either µ 1 = ν such that The remaining cases form the subset of pairs. The probability densities W(R ν, µ0 ) of NMs being resonant when taking into account all contributions are shown in Fig. 6. The main result is that W(R ν, µ0 → 0) → C(W ) = 0. The constant C drops with increasing disorder strength W (see Fig. 7). We also calculate C by taking into account only quadruplets and triplets (see Fig. 7). We find that for weak disorder the quadruplet contributions are the dominant ones, while for strong disorder their contribution diminishes as compared to the triplet contribution.
For small R the probability densities W(R) can be approximated as The probability P for a mode, which is excited to a norm n, to be resonant at a given value of the interaction parameter β is given by

VI. DISCUSSION
We have studied statistical properties of eigenvalues and eigenvectors of waves in disordered one-dimensional systems as a function of the disorder strength. We estimated the localization volume of a mode which defines the number of interacting partner modes. We obtained the dependence on the disorder strength of the overlap integrals which determine the interaction strength. We analyzed the statistics of level spacings of normal modes within one localization volume. Finally, we obtained distribution functions for resonance probabilities of normal modes interacting in the presence of nonlinearity. Let us discuss some of the consequences of our findings.

A. Overlap integrals
In order to estimate the absolute value of the overlap integral (5) for modes within one localization volume for weak disorder, Shepelyansky [22] and Imry [23] assumed that the sum extends roughly over the localization volume V , with each term in the sum A ν,l A ν1,l A ν2,l A ν3,l having a random sign. The absolute value of the eigenvector is of the order of 1/V 1/2 due to normalization. Then (5) can be evaluated using the central limit theorem, for which the average absolute total value I rp ∼ V −3/2 ∼ W 3 . Our numerical finding I ∼ W 3.4 clearly rules out the random sign resut W 3 . As shown in the perturbation calculation in section III, the reason for the random sign failure is that NMs are similar to plane waves with definite phases on each lattice site (inside the localization volume). These phases enforce selection rules, which become strict in the very limit W = 0.
While we can now exclude the random sign result W 3 , we can not tell whether the numerical estimate I ∼ W 3.4 is correct, or the perturbation result I ∼ −W 4 ln(W ) will set in for small enough W . Ponomarev and Silvestrov [24] have also stressed the importance of phase correlations in Eq. (5). A numerical calculation of the average of the squared overlap integral was performed by Frahm et al [25] for 1.4 < W < 4 yielding I ∼ W 3.3 , in a good agreement with our numerical data.
The random sign estimate W 3 was taken to predict a strong increase of the localization length of two interacting particles in a one-dimensional random quantum chain [22,23]. The two particle localization volume V 2 , within a renormalization group approach, is given by where V is the single particle localization volume. For the random phase result, this yields V 2 ∼ V 2 [22,23,26]. We can clearly rule out such an outcome. Instead, we expect either V 2 ∼ V 1.6 (numerical data) or V 2 ∼ V ln 2 V (perturbation approach), which give a much weaker effect. These controversies call for more detailed investigations.

B. Asymptotic spreading of wave packets in nonlinear chains
According to a recent analysis of the spreading scenaria of wave packets [12], the only scale which separates different dynamical spreading regimes is the average spacing d. Therefore, the constant C from the previous section is inversely proportional to the mean level spacing: Following the theory developed in [7,9,12] for the asymptotic spereading, an exterior mode φ µ which is heated up by the packet obey the following evolution equation in accordance with (4) iφ µ ≈ λ µ φ µ + β I V 3 P(βn)n 3/2 f (t), where f (t)f (t ′ ) = δ(t − t ′ ) ensures that f (t) has a continuous frequency spectrum. Note, that here we also introduce the contribution of the overlap integrals estimated as I V 3 . Repeating the previous derivations [7,9,12], we finally get the following expression for the asymptotic growth of the second moment of spreading wave packets in nolinear chains From our numerical data for weak disorder it follows m 2 ∼ W −3.07 β 4/3 t 1/3 , while the perturbation approach yields m 2 ∼ W −8/3 (− ln W ) 2/3 β 4/3 t 1/3 . The prefactor dependence of W is another intriguing test which awaits numerical verification.

VII. CONCLUSION
In conclusion, we performed a statistical analysis and calculated the average localization volume occupied by an eigenmode as a function of disorder strength which determines the average number of a nonexponentially interacting eigenmodes. Then, we calculated the frequency spacings of the normal modes which happen to interact in a nonexponentially weak way and their distributions and the average numerically. This result is very important for the classification of different regimes of wave packet spreading in the presence of nonlinearity. We also studied statistical properties of the overlap integrals which determine the coupling strength between the interacting modes and, thus, influence properties of spreading. Finally, we estimated the number of resonant modes in the packet and proved that the most significant contribution to the spreading comes from the quadruplet and triplet resonances for small to moderate values of disorder strengths, and from triplets for the case of large disorder.