Quantum chaos in a harmonic waveguide with scatterers

A set of zero-range scatterers along its axis lifts the integrability of a harmonic waveguide. Effective solution of the Schr\"odinger equation for this model is possible due to the separable nature of the scatterers and millions of eigenstates can be calculated using modest computational resources. Integrability-chaos transition can be explored as the model chaoticity increases with the number of scatterers and their strengths. The regime of complete quantum chaos and eigenstate thermalization can be approached with 32 scatterers. This is confirmed by properties of energy spectra, the inverse participation ratio, and fluctuations of observable expectation values.

However, a generic system is not completely chaotic nor integrable (see examples in ).Certain incompletely-chaotic systems -the systems with no selection rules -relax to a state whose properties are governed by the inverse participation ratio (IPR) [28,31].Inverse of this parameter estimates the number of integrable system eigenstates comprising the nonintegrable one.IPR ranges from 0 for completely-chaotic systems to 1 for integrable ones.Then it can serve as a measure of the system's chaoticity [53].IPR also governs fluctuations of eigenstate expectation values [32].The energy-spectrum statistics of incompletely-chaotic systems lie between the Wigner-Dyson and Poisson ones.Certain systems demonstrate the Šeba statistics [22].
The most obvious objects of chaotic property simulation are lattice systems.However, they have a finite Hilbert space and its dimension is restricted due to computational difficulties (complexity of lattice system simulations increases as a high power of the lattice site number and exponentially with the number of particles).Then, on increase of the system chaoticity, each eigenstate can fill the full Hilbert space.A system with infinite Hilbert space -the Sinai type billiard -was analyzed in [54], where ∼ 3 × 10 5 eigenstates were calculated.However, chaoticity of such billiard cannot be tuned.
The present model -a particle in a harmonic waveguide with zero-range scatterers along its axis -has an infinite Hilbert space.As the scatterers are a particular case of independent perturbations [55], IPR should be inversely proportional to the number of scatterers.The model chaoticity can also be tuned by the scatterer strengths.This model was already used in [55] for numerical confirmation of the general relations between properties of wavefunctions and the number of scatterers.The present paper is devoted exclusively to the harmonic waveguide with scatterers and analyzes properties of wavefunctions for weak perturbations and for additional models, as well as properties of energy spectra.
Since a zero-range scatterer is a particular case of separable interactions, the present model belongs to systems with high-rank separable perturbations [56].Energy spectra of several physical systems of such type have already been considered.They are the flat rectangular billiards -generalization of the Šeba billiard [21] -with 1-3 [23], 6 [24], and 2 [46] scatterers.Theoretical predictions for a single scatterer in a harmonic potential were compared to experiments [57].Series of separable interactions can also approximate the dipole-dipole ones [58,59].Energy spectra of two dipolar particles in a harmonic trap were calculated [60] using such expansion.An advantage of systems with separable rank-s interactions is that calculations require diagonalization of a s × s matrix, (cf. to α × α matrix in the direct diagonalization method for α eigenstates).In addition, the present model allows an analytical summation over axial states.Then the system properties are calculated here for millions of eigenstates.
The paper has the following organization.The model is described in Sec.(2).Section (3) analyzes the energy spectra statistics.Properties of wavefunctions, including expectation value fluctuations and IPR, are presented in Sec. 4. Appendices provide derivation details and additional technical information.
A system of units in which Planck's constant is ħ h = 1 is used below.

The model
The Hamiltonian of a particle with the mass m in an axially-symmetric harmonic waveguide with the transverse frequency ω ⊥ contains the kinetic and potential energies, Here z is the axial coordinate, ρ = x 2 + y 2 is the transverse radius, △ ρ is the transverse Laplacian, and A is the vector potential (its role will be discussed below).Integrability of the perturbed Hamiltonian is lifted by the zero-range scatterers where δ r e g is the Fermi-Huang pseudopotential and the scatterers are located along the waveguide axis, i.e., their positions R s ′ = (0, 0, z s ′ ) have zero transverse coordinates.The scatterers are numbered from left to right (z s ′ > z s ′′ if s ′ > s ′′ ).The model is restricted in the sector of the axially-symmetric states, as other states vanish at the waveguide axis, and, therefore, are not affected by the scatterers.Then the eigenstates of Ĥ0 , labeled by the axial l and radial n ≥ 0 quantum numbers, are 〈ρ, z|nl〉 = 〈ρ|n〉 〈z|l〉 with the radial wavefunctions Here a ⊥ = (mω ⊥ ) −1/2 is the transverse oscillator range and L (0) n are the Laguerrre polynomials (see [61]).The discrete energy spectrum is provided either by the periodic boundary conditions (PBC) 〈z + L|l〉 = 〈z|l〉, or by the hard-wall box (HWB) 〈z = L|l〉 = 〈z = 0|l〉 = 0. Then the axial wavefunctions are either with 1 ≤ l < ∞ for HWB.
For PBC, the eigenstate |nl〉 of Ĥ0 has the eigenenergy where λ = (L/a ⊥ ) 2 characterizes the aspect ratio and l 0 = LA/(2π) is the scaled vector potential.If A = 0, the inversion (P) invariance of the Hamiltonian Ĥ0 leads to the degeneracy of the energies E nl and E n−l .This degeneracy can be lifted by any P-noninvariant perturbation.
The vector potential lifts it as well, with no effect on the simple wavefunctions (5), though the Hamiltonian losses the time-reversal (T) invariance.
Four kinds of the model are considered here.The first three kinds correspond to PBC.The first, non-symmetric, model has A ̸ = 0 and is T-noninvariant.The scatterer positions form irregular sequence due to random shifts −0.25 ≤ δ s ′ < 0.25.The shifts are calculated once for each number of scatterers and there is no average over the shifts.In the second, symmetric, model with is non-degenerate as l is positive.Together, the four kinds of the model cover different symmetries of the Hamiltonian (Tinvariant, PT-invariant, and non-symmetric), as well as different boundary conditions (PBC and HWB).
The eigenstates of the non-integrable system |α〉, solutions to the Schrödinger equation Ĥ |α〉 = E α |α〉, are labeled in the increasing order of the eigenenergies E α .Expansion over the integrable system eigenstates |nl〉 transforms the Schrödinger equation to the form According to (3) where the value of the regular part of |α〉 at R s ′ is The last equality above follows from the spherical symmetry of 〈r|α〉 in the vicinity of R s ′ [62].
As a result, we get the following system of linear equations for 〈R s ′ |α〉 r eg For the wavefunctions (4) and ( 5) or (6) and energies (7) or (9) the sum over l above can be calculated analytically (see Appendix A).Then the system (13) attains the form with Here [] denote the integer part, ) is the Hurwitz zeta function (see [61]), The system ( 14) has a non-trivial solution at ϵ = ϵ α ≡ mL 2 (E α − ω ⊥ )/2 where an eigenvalue of its matrix has a root as a function of ϵ.The matrix S s ′ s ′′ (ϵ) has poles at ϵ = ϵ nl , as it is seen from ( 13).Then the eigenvalues can have poles at ϵ = ϵ nl as well.Between these poles each eigenvalue is a monotonic function of ϵ, as demonstrated by direct calculations (see Appendix B).Than all eigenenergies ϵ α in each interval between neighboring ϵ nl can be calculated as roots of s eigenvalues.Although the eigenvalue monotonicity was not proved exactly, this algorithm provides the number of eigenenergies ϵ α which differs from the number of ϵ nl in the same energy interval by not more than s.It is an evidence that no eigenenergies ϵ α are lost.
There seems to be no fundamental obstacle for experimental realization of the present model.In the case of cold trapped atoms, atoms of other kind in optical tweezers might play the role of scatterers, and the interaction strength might be tuned by a Feshbach resonance.T-noninvariant models might be realized with trapped ions in a magnetic field.In optics, optical defects might work as scatterers [63] for photons in an optical cavity or waveguide (see also [64,65] and the references therein).The PBC models might be realized with circular atomic or optical waveguides.

Statistics of energy spectra
The differences in energy spectra between integrable and chaotic systems were the first distinctive properties of quantum chaos (see [9][10][11]).These properties are defined in terms of the unfolded energy ᾱ(ϵ α ) -the smooth part of the dependence α(ϵ α ).For the present model, the unfolding function is the same as for the underlying integrable system.The number of states below the scaled energy ϵ is the staircase function For PBC, using Eq. ( 7) for ϵ nl , we have The smooth part is extracted by replacing the integer part [x] with x − 1/2.The limits of the sum over l, [l 0 ± ϵ/π], are replaced in the same way.As a result, we get The ϵ-independent terms are dropped here, since only differences between ᾱ(ϵ α ) appear in the following expressions.Similar expression is obtained for the HWB model The first property of the energy spectrum considered here is the nearest-neighbor distribution (NND) -the density of probability to have the given value of the unfolded energy difference ᾱ = ᾱ(ϵ α ) − ᾱ(ϵ α−1 ) between the neighboring energy levels [9][10][11].Integrable systems have the Poisson NND, while completely-chaotic ones have the Wigner-Dyson distributions for Gaussian ensembles of random orthogonal matrices (GOE) and unitary matrices (GUE) in the cases of T-invariant and T-noninvariant systems, respectively.The Šeba NND with A Se ba ≈ 2.1266 and B Se ba ≈ 0.3481 was obtained [22] for certain incompletely-chaotic systems.
All states of a non-integrable system correspond to the same symmetry and then their energies demonstrate repulsion.Then NND ( 22), (23), and (24) of non-integrable systems vanish at the zero level spacing and decrease approaching this point.The integrable system states of different symmetry can be energy degenerate, and then NND (21) decreases exponentially with the level spacing.For non-integrable systems NND decreases at large spacing too, although completely-chaotic systems are characterized by Gaussian decrease [see Eqs. ( 22) and ( 23)], while the Šeba NND (24) decreases exponentially.
Another property of energy spectra is the spectral rigidity ∆ 3 (∆ ᾱ) -the least-square deviation of the staircase function α(ϵ) from the best fit to a straight line on a given interval of the unfolded energy ∆ ᾱ.The spectral rigidity for integrable and completely-chaotic (T-invariant and T-noninvariant) systems are given, respectively, by [9-11] where γ Eul ≈ 0.5772 is the Euler's constant [61].
The energy spectrum properties are calculated below for the four kinds of the models.The parameters l 0 = 0.25−e −4 ≈ 0.232 (for the T-noninvariant models) and λ = π 3 (1+ 5) ≈ 100 are expressed in terms of transcendent numbers [(1 + 5)/2 is the golden ratio].Most of the results are obtained for 10 6 eigenstates in the unitary regime, V s ′ = 10 6 V 0 for all scatterers.
Figure 1(a) shows NND calculated for the non-symmetric model with different numbers of scatterers in the unitary regime.For s = 2, NND follows the Šeba plot, as well as for the case of s = 1 considered in [27].When the number of scatterers increases, NND tends to the GUE prediction and approaches it at s = 32.GUE is approached as the model is T-noninvariant.The calculated spectral rigidity [see Fig. 1(b)] demonstrates the same tendency.
For the symmetric model (see Fig. 2), NND for s = 2 again follows the Šeba predictions (indeed, the case of two scatterers of the same strength is always P-invariant).However, at s = 32, NND and spectral rigidity approach the GOE predictions, although the system is Tnoninvariant.It is a consequence of the real matrix of the interaction with scatterers obtained when the z coordinate origin is shifted to  This assumption is confirmed by the NND and spectral rigidity for the HWB model (see Fig. 4).This model with non-degenerate energy spectrum is more chaotic than the PBC T-invariant one, as now Šeba and GOE predictions are approaching at s = 16 and s = 32, respectively.Then, this model is less chaotic than the T-noninvariant PBC ones.There is also a noticeable difference between these models in the statistics of integrable system energy spectra -for the HWB model NND at small spacings and spectral rigidity are below the Poisson predictions.
Thus, for all kinds of the model the statistics tend to the Wigner-Dyson predictions on increase of the number of scatterers.This agrees with the behavior of spectral rigidity of flat 2D billiards [23].However, the present model does not demonstrate another property of the 2D flat billiards -the shifting toward Poisson statistics at higher energy [23].It is clearly shown in Fig. 5, where the plots for different energy regions are close together and do not demonstrate a systematic dependence on the energy.This difference is related to the nature of the logarithmic asymptotic freedom revealed in [23].This effect is caused by the decreased effective interaction strength v eff ∼ 1/ ln ϵ (see Eq. ( 19) in [23]), while the characteristic energy level separation ∂ ϵ α /∂ α is independent of the energy for 2D billiards with ϵ α ∝ α.In contrast, if ϵ α ∝ α γ (γ ̸ = 1), the derivation [23] would lead to v eff ∝ ϵ 1−1/γ , while ∂ ϵ α /∂ α ∝ ϵ 1−1/γ has the same energy dependence and ratio of the effective interaction strength to energy level separation is independent of energy.Therefore, the logarithmic asymptotic freedom does not  appear in the present model with ϵ α ∝ α 2/3 as well as in generic systems with ϵ α ∝ α γ (γ ̸ = 1), being a specific property of 2D billiards.
The transition between the T-invariant and non-symmetric models due to the change of the vector potential is demonstrated in Fig. 6.The GUE and Šeba statistics take place at l 0 < 10 −4 and l 0 > 10 −2 , respectively.
The system chaoticity depends also on the scatterer strength V s ′ .NND approaches this unitary regime already at V s ′ = 10 −1 V 0 , as Fig. 7(a) shows.For V s ′ = 10 −4 V 0 NND almost coincides with the integrable system one.Spectral rigidity demonstrates the same behavior (see Fig. 7(b)).Thus, the system's chaotic properties depend on two parameters: the number of scatterers and their strengths.Interaction of these parameters is illustrated by Fig. 8, which demonstrates that the NND and spectral rigidity dependencies in the unitary regime for 4 scatterers are approached at V s ′ = 10 −2 V 0 and V s ′ = 5 × 10 −3 V 0 for 8 and 32 scatterers, respectively.
Figure 9 shows dependence of the system statistics on the scatterer locations.All nonsymmetric cases (1 and 2, corresponding to different sets of the random shifts δ s ′ in (8), and 3, where ζ 1 = 0 and ζ s ′ with s ′ > 1 are chosen randomly from the interval [0, 1] and sorted) provide close results approaching the GUE predictions.The plots for the symmetric distribution are clearly different and approach GOE predictions (see the discussion above).
If the scatterer positions form a periodic sequence, ζ s ′ = (s ′ − 1)/s and V s ′ is constant, the picture is completely different.In this case, according to the Bloch's theorem, the eigenstate can be expressed as 〈ρ, z|α〉 = ρ, z|α p exp(ipz).The L-periodicity plays the role of the Bornvon Karman boundary conditions, leading to the discrete spectrum of the quasimomentum p = 2πk p /L with integer k p .The function ρ, z|α p has the period L/s and satisfies the Schrödinger equation with single scatterer Here the integrable Hamiltonian Ĥ0 (A− p) of the form (1) contains the vector potential A− p.Therefore, the total energy spectrum is a superposition of s spectra of the one-scatterer systems with scaled vector potentials l 0 − k p (k p + s gives the same result as k p ).This is the reason (see [9]) why NND for the periodic case does not have a dip at small spacings and both NND and spectral rigidity are close to the Poisson predictions.
The random matrix theory [9][10][11] predicts universal spectral rigidity plots corresponding to the Poisson, GOE, and GUE NNDs.However, the Šeba NND can correspond to various spectral rigidity plots, as it is shown in Fig. 10.It is worth noting that the plots for the T-invariant PBC and HWB models are close together, while the one for the T-noninvariant model is completely different.
Statistics of energy spectra can be also characterized by average level spacing ratio [66,67]which increases with the system chaoticity.In the present case (see Appendix C) , this monotonic increase takes place only in the vicinity of the Poisson statistics.However, the average level spacing ratio becomes almost the same for the Šeba and GOE statistics and has strong fluctuations on the transition between them.This may be related to the small number of the degrees of freedom in the present models compared to many-body models, where the level spacing ratio is generally used.An additional advantage of the level spacing ratio is that unfolding the spectrum is not required.However, this advantage is not essential for the present models as the unfolding functions are well defined.For these reasons, the level spacing ratio is not used here.

Properties of wavefunctions
Possibility of statistical description of quantum-chaotic systems is based, through ETH, on properties of their wavefunctions.The number of integrable system eigenstates comprising the non-integrable one is characterized by the number of principal components (NPC) η −1 , where η = nl |〈n, l|α〉| 4 is IPR.Equations ( 10) and (11) allow us to express the expansion coefficients here in the form where 〈R s ′ |α〉 r e g are solutions to the system ( 14) and the normalization factor N α is determined by the normalization condition nl |〈n, l|α〉| 2 = 1.For the energies (7) the sums over n here and in IPR can be expressed in terms of the Hurwitz zeta functions (see [61]) where Then the normalization condition takes the form ∞ l=−∞ P l = 1, where is the occupation of the states with the given axial quantum number l and Similarly, for IPR we have The expressions above are used for T-noninvariant models (non-symmetric and symmetric), where A ̸ = 0 and 〈R s ′ |α〉 r e g are complex.In the T-invariant models (A = 0) 〈R s ′ |α〉 r eg are real.
For PBC the normalization condition can be expressed as ∞ l=0 P T l = 1, where and Respectively, IPR can be expressed as For HWB we have the normalization condition and where and q B l = (π 2 l 2 /4 − ϵ α )/λ.A recurrence relation η −1 s = η −1 s−1 + ν was derived [55] for NPC η −1 s of the system with s scatterers.This means that NPC increases and, respectively, IPR decreases with the number of scatterers.In the case of weak interaction the dependence of NPC on the number of scatterers is nonlinear (see Fig. 11).This is a consequence of the strong dependence of the system's chaotic properties on the number of scatterers.NPC also increases with the eigenstate energy due to increase of the energy level density.In the case of the statistics of energy spectra, this increase was compensated by decrease of the effective interaction strength (see Fig. 5 and the related discussion above).Here we see that the wavefunction properties are determined by the interaction strength V s ′ rather than the effective one.The NPC dependence on the number of scatterers can be approximated by η −1 s = 1 + ν ′ s γ .For a weak interaction V s ′ = 10 −3 V 0 the power γ ≈ 0.48 becomes independent of the eigenstate energy (see Fig. 11(a)).For stronger interaction V s ′ /V 0 = 10 −2 (see Fig. 11(b)) the power γ increases with the eigenstate energy  and the dependence of NPC on s tends to the linear one.This means that ν is independent of s since the system's chaotic properties are independent of the number of scatterers.In the unitary regime, this dependence is confirmed by IPR calculated for all kinds of the model (see Fig. 12(a)) which is approximated by inverse-linear functions with a good accuracy.We can see that for each number of scatterers the non-symmetric model has the minimal IPR and, therefore, demonstrates the highest chaoticity, the T-invariant PBC model has the highest IPR, and the HWB one lies between them.This order agrees with the NND and spectral rigidity of energy spectra for these models discussed in Sec. 3 above.However, the symmetric T-noninvariant model has substantially higher IPR than the non-symmetric one, although properties of energy spectra of these models demonstrate similar chaoticity.This difference can be related to properties of real and complex random Gaussian variables [55].As well as any characteristic of chaos, IPR depends also on the interaction strength (see Fig. 12(b)).This figure also demonstrates that the systems with 4, 8, and 32 scatterers have approximately the same IPR (η ≈ 0.2) at V s ′ /V 0 = 10 −1 (and in the unitary regime), 10 −2 , and 5 × 10 −3 ,respectively, in agreement with the energy spectra statistics (see Fig. 8).Chaotic properties of physical systems are also characterized by fluctuations of observable expectation values.Expectation value of the observable Ô in eigenstates of the non-integrable system is related to ones in integrable system eigenstates where the expansion coefficients 〈n, l|α〉 are given by (28).Four observables are considered here.The transverse potential energy mω 2 ⊥ ρ 2 /2 is nondiagonal in the integrable system eigenstates As the potential energy increases with the total energy, the part U ⊥ of the transverse potential energy in the total energy is considered here.Its expectation value in the non-integrable system eigenstates can be expressed as (see Appendix D) for T-noninvariant models.In the T-invariant PBC case we have In the last case, HWB, the expectation value takes the form Other observables, diagonal in integrable system eigenstates, are the axial momentum where θ (l) = 0 for l < 0, 1/2 for l = 0, and 1 for l > 0, and the occupation of the odd axial modes nl Podd n ′ l ′ = δ n ′ n δ l ′ l δ lmod2,1 , where lmod2 is the reminder of the division of l by 2. For T-noninvariant models their expectation values are expressed in terms of the occupations P l (31), 2l−1 are expressed in terms of probabilities ( 34) and (37), respectively.For an observable Ô, the variance of its expectation value fluctuations between nonintegrable system eigenstates is defined as According to [32], this variance is proportional to IPR and the variance between the integrable system eigenstates Var α ( Ô) = ηVar nl ( Ô) .
The latter variances are calculated in Appendix for the four observables presented above.The variance of the axial momentum depends on the averaging interval [ϵ min , ϵ max ] boundaries.The variances of other observables are independent of the interval, Var nl ( Ppos ) = 1/4, Var nl ( Podd ) = 1/4, and Var nl ( Û⊥ ) = 1/45.Figure 12(b) confirms the rule (47) for the integrability-chaos transition on variation of the scatterer strength in the non-symmetric model, both for 4, 8, and 32 scatterers.This rule is also confirmed when the number of scatterers is changed for all four models considered here (see Fig. ( 13)).

Conclusion
An effective method of numerical solution, based on properties of high-rank separable perturbations, is developed for a harmonic waveguide with a vector potential and either PBC or HWB in the axial direction, perturbed by zero-range scatterers along the waveguide axis.The energy-degeneracy of the unperturbed system can be lifted by the vector potential which also lifts T-invariance.The energy spectra properties -near-neighbor distribution and spectral rigidity, as well as IPR and fluctuation variance of observable expectation values, are calculated for 10 6 eigenstates.The chaoticity measures of the model increase with the number of scatterers and their strengths.This allows exploring the integrability-chaos transition.decreases exponentially with n and, due to the translational invariance, it is independent of ζ.In the limit of ζ → 0, the sum of the first terms in T n (ζ, 0) was calculated in [62] ∞ in terms of the Hurwitz zeta function (see [61]).The first, proportional to ζ −1 , term here is removed by the derivative in (13).Then we get Eqs.( 14) and ( 15).
For T-invariant models, when A = 0, we have real T n (ζ, 0).In the case of PBC, we can just set l 0 = 0 in Eqs.(A.5), (A.7), and (A.8) and get 1 − e −2|p n | (λn > ϵ, ζ > 0) , (A.10) In the case of HWB, substitution of Eqs. ( 6) and ( 9) to (A.1) leads to 2), it is not a function of z − z ′ only, since HWB is not translational invariant.Using partial fraction decomposition and the real part of the summation formula (A.4), we get for For λn > ϵ α and ζ > ζ ′ we have The term causing the divergence is separated in the same way as in Eq. (A.8), providing where pk = p n k sign(l k − l 0 ).In the T-invariant case, when l k ̸ = 0, they can be expressed as For HWB, when l k ̸ = 0, we have and has a form of the matrix with i max eigenvalues B i .When ϵ approaches ϵ k , the singular part dominates and the eigenvalues tend to ±∞.Then in the T-invariant PBC case with l k ̸ = 0 we have i ma x = 2 and two eigenvalues of the matrix S s ′ s ′′ (ϵ) have singularities at ϵ → ϵ k , there are no singular eigenvalues (i ma x = 0) in the case of HWB with l k = 0, and single eigenvalue has a singularity in other cases when i max = 1.Results of numerical calculations in Fig. 14 demonstrate these properties.They also show that the eigenvalues decrease monotonically with ϵ.Then each eigenvalue can have single root in the interval [ϵ k , ϵ k+1 ].In Fig. 14, the number of eigenvalues with roots increases from 0 to 4 in parts (a)-(e).
In the close vicinity of ϵ k direct numerical diagonalization of the matrix S s ′ s ′′ (ϵ) becomes inaccurate if i ma x > 0. However, in this vicinity i max eigenvalues are approximated by B i with good accuracy.In order to calculate other eigenvalues, the matrix S cont s ′ s ′′ (ϵ) is projected out of the envelope of the vectors b

D Expectation values
Substituting Eqs. ( 28) and (41) into Eq.( 40) we can get the following expression for the expectation value of the transverse potential energy in the non-integrable system eigenstates , (D.1) where the last transformation uses the normalization condition, Eq. ( 7) for ϵ nl , and Eq. ( 30) for q i .The sum over n here can be transformed as = −q l ζ(2, q l ) + (q l − 1) 2) where the summation over n with Eq. ( 29) for the first term in the square brackets and partial fraction decomposition for the second term are used.The last sum over n is reduced to 1/(q l − 1) due to cancellation of the terms.This leads to Eq. ( 42).The derivation above is related to the T-noninvariant models.The same transformation of the sum over n leads to the expectation values for the T-invariant PBC (43) and HWB (44) models.
The variance between the integrable system eigenstates can be evaluated analytically.The product ᾱ(ϵ)U ⊥ can be approximated by the sum It has the same ϵ dependence as ᾱ(ϵ) taken with the same accuracy [the first term in Eq. ( 19)].Then the average U ⊥ = 1/3 is independent of the averaging interval (this value agrees to the virial theorem).In the same way we find U 2 ⊥ = 2/15 and, therefore, Var nl ( Û⊥ ) = 1/45.Although in the HWB model l ≥ 0, we get the same results due to the distinction between Eqs. (7) and (9).
For the average axial momentum, we approximately evaluate the sum

Figure 1 :
Figure 1: Near-neighbor distribution (a) and the spectral rigidity (b) for the nonsymmetric model with different numbers of scatterers in the unitary regime .

2 .Figure 2 :
Figure 2: Near-neighbor distribution (a) and the spectral rigidity (b) for the symmetric model with different numbers of scatterers in the unitary regime .

Figure 3 :
Figure 3: Near-neighbor distribution (a) and the spectral rigidity (b) for the Tinvariant PBC model with different numbers of scatterers in the unitary regime .

Figure 4 :
Figure 4: Near-neighbor distribution (a) and the spectral rigidity (b) for the HWB model with different numbers of scatterers in the unitary regime .

Figure 5 :Figure 6 :
Figure 5: Near-neighbor distribution (a) and the spectral rigidity (b) for the nonsymmetric model with 4 and 8 scatterers in the unitary regime at various regions of the non-integrable system eigenstate labels α .

Figure 7 :Figure 8 :
Figure 7: Near-neighbor distribution (a) and the spectral rigidity (b) for the nonsymmetric model with 32 scatterers for various scatterer strengths.

Figure 9 :
Figure 9: Near-neighbor distribution (a) and the spectral rigidity (b) for the nonsymmetric model with 8 scatterers for various scatterer locations.

Figure 10 :
Figure 10: Near-neighbor distributions (a) which are close to the Šeba one for various models and the corresponding spectral rigidity (b).

Figure 11 :
Figure 11: Inverse participation ratio as a function of the number of scatterers calculated for the non-symmetric model with V s ′ /V 0 = 10 −3 (a) and V s ′ /V 0 = 10 −2 (b) at various regions of the non-integrable system eigenstate labels α.The lines show the dependencies η s = 1/(1 + ν ′ s 0.48 ) with the ν ′ values presented in the legend in the part (a) and η s = 1/(1 + ν ′ s γ ) with the ν ′ and γ values presented in the legend in the part (b) .

Figure 12 :
Figure 12: (a) Inverse participation ratio as a function of the number of scatterers calculated in the unitary regime for four kinds of the model.The lines show the best fit by inversely linear functions.This part uses the same data as Fig. 3(a) in [55].(b) Ratio of fluctuation variances between the non-integrable and integrable systems eigenstates as a function of the scatterer strength for the non-symmetric model with 4, 8, and 32 scatterers.The points correspond to the four observables, the lines connect the calculated IPR values.The data for 32 scatterers are the same as in Fig. 3(b) of [55].

Figure 13 :
Figure 13: Ratio of fluctuation variances between the non-integrable and integrable systems eigenstates as a function of the number of scatterers for T-noninvariant (a) and T-invariant (b) models in the unitary regime.The points correspond to the four observables, the lines connect the calculated IPR values.The data for the nonsymmetric model in the part (a) are the same as in Fig. 3(c) of [55].

Figure 14 :
Figure 14: Examples of eigenvalue dependence on the energy between two neighboring eigenenergies of the integrable system for the non-symmetric model with 4 scatterers (a-e) and the T-invariant model with 8 scatterers (f).

Figure 15 :
Figure 15: The level spacing ratio averaged over different eigenstate label intervals (a) for the non-symmetric model with 32 scatterers as a function of the interaction strength and (b) for the symmetric model in the unitary regime as a function of the number of scatteres.
2, the scatterer positions are invariant over inversion under (z 1 +z s )/2.This inversion changes the sign of the term (i/m)A∂ /∂ z in the Hamiltonian Ĥ0 .This sign is also changed by the time-reversal (complex conjugation).Then the symmetric model with equal V s ′ is PT-invariant.The third, T-invariant, model has A = 0 and the same scatterer positions as the non-symmetric one.Only this model has a degenerate energy spectrum of the integrable Hamiltonian.The fourth, box, model corresponds to HWB.The scatterer positions are z s ′ = (s ′ + δ s ′ )L/(s + 1).Although A = 0, the energy spectrum is the scale of the interaction strength, and the summandsT n (ζ s ′ , ζ s ′′ ) and T (ζ s ′ , ζ s ′′ ) with ζ s ′ ≥ ζ s ′′ have to be calculated.T n (ζ s ′ , ζ s ′ ) and T T n (ζ s ′ , ζ s ′′ ), as well as the matrix S s ′ s ′′ (ϵ), is real, and S s ′ s ′′ (ϵ) is symmetric.