Berezinskii approach to disordered spin systems with asymmetric scattering and application to the quantum boomerang effect

We extend the Berezinskii diagrammatic technique to one-dimensional disordered spin systems, in which time-reversal invariance is broken due to a spin-orbit coupling term inducing left-right asymmetric scattering. We then use this formalism to theoretically describe the dynamics of the quantum boomerang effect, a recently discovered manifestation of Anderson localization. The theoretical results are confirmed by exact numerical simulations of wave-packet dynamics in a random potential.


I. INTRODUCTION
In the presence of a spatially disordered potential, quantum wave packets may experience, after an transient temporal spreading, a complete freezing of their density distribution due to the proliferation of destructive interference in the multiple scattering process.This phenomenon, which generically occurs in low dimension, is one of the most representative manifestations of Anderson localization [1].As such, it has been primarily exploited in the experimental quest for the localization of cold atoms in random potentials [2][3][4][5].Recently, however, a variety of alternative signatures of Anderson localization has been identified.Those include the temporal freezing of the coherent backscattering effect in reciprocal space [6][7][8] or the universal growth of narrow peak structures in the density profile [9,10] and momentum distribution [11][12][13] of spreading wave packets (see [14] for a review).
Recently, yet another unexpected manifestation of Anderson localization, dubbed quantum boomerang effect (QBE), has been discovered [15].The QBE corresponds to a back-and-forth motion of the mean position of a quantum wave packet launched with a finite velocity in a given direction in a random potential.In one dimension, for instance, if the quantum particle is launched to the right, it will first move to the right over a distance of the order of the mean free path, then make a U-turn and eventually return to its starting point at long time.This phenomenon was also shown to exist in higher-dimensional random or pseudo-random systems [15,16], as well as in kicked-rotor models [16], where it was recently demonstrated experimentally [17].While originally described in time-reversal-invariant (TRI) systems, recently the QBE was also shown to exist in systems without time-reversal symmetry [18][19][20].In [18], in addition, the QBE was characterized in the presence of a spin-orbit coupling mechanism inducing left-right asymmetric scattering between different spin states.This is also the scenario addressed in the present paper.
At a theoretical level, describing the temporal dynamics of quantum wave packets in the presence of disorder is a challenging task [21][22][23].In one dimension, however, a very powerful analytical approach known as the Berezinskii diagrammatic technique has been developed [24].Originally, this method was successfully used for calculating the ac conductivity of electronic conductors in the localization regime or the long-time density distribution of spreading wave packets [24][25][26], the predictions being exact in the limit of weak disorder.More recently, it also allowed to describe the QBE in TRI systems [15] and, in the context of electron scattering, was extended to account for the presence of spin-orbit coupling [27].
In this paper, we extend the Berezinskii diagrammatic technique to TRI-broken spin-dependent systems in which a spin-orbit coupling term induces asymmetric scattering, as recently realized experimentally with cold atoms [28,29].This formalism is developed in sections II and III.In Sec.IV we then apply the method to the calculation of a specific observable, the mean position of a quantum-mechanical wave packet launched in a random potential with finite velocity.This provides a thorough theoretical description of the QBE in spinorbit coupled systems with asymmetric scattering, complementing results obtained in the recent work [18].Generally speaking, the formalism presented in this paper provides a practical analytical tool to characterize the dynamics of spinor wave packets in disordered systems with TRI-broken symmetry.

II. PRINCIPLES OF THE BEREZINSKII TECHNIQUE
We start by recalling the main ideas of the original Berezinskii diagrammatic technique used to compute the time-dependent transport properties of one-dimensional disordered systems.The starting point is the singleparticle Hamiltonian arXiv:2307.15354v1 [cond-mat.dis-nn]28 Jul 2023 where V (x) is a random (disorder) potential and H 0 is the disorder-free part of the Hamiltonian (e.g., H 0 = p 2 /2m).
We suppose that the random potential has a vanishing mean, V (x) = 0, and follows a Gaussian statistics characterized by the two-point correlation function , where η is called the disorder strength.Symbol (. ..) here denotes averaging over different disorder realizations.The function C(x ′ −x) quantifies the range of the spatial correlation of the disorder.
In the whole paper, we restrict ourselves to a delta correlated potential, i.e., C(x In this paper, we aim at describing the time evolution of quantum-mechanical wave packets governed by an Hamiltonian of the type of Eq. ( 1).In the localization problem, this evolution is characterized by considering the disorder average of observables that depend quadratically on the wave function, such as the density n(x, t) or the mean position ⟨x(t)⟩ = dx x n(x, t) of the wave packet.These observables, by definition, can be expressed in terms of the disorder-averaged correlator G R (x, x ′ , ϵ)G A (x ′′ , x, ϵ − ℏω) [30,31], where −1 |x ′ ⟩ are the singlerealization, retarded and advanced Green's functions at energy ϵ associated with Hamiltonian (1).The energy difference ℏω introduced in the correlator allows us to capture the time dependence of observables after an inverse Fourier transform.The precise connection between ⟨x(t)⟩ and the Green's function correlator, for instance, will be given in Sec.III.Both Green's functions G R/A that appear in the correlator may be computed in a perturbative fashion using the Born expansion [30]: where are the retarded and advanced Green's functions associated with the free part of the Hamiltonian.Physically, the expansion (2) describes a multiple scattering sequence involving scattering events on the random potential at points x 1 , x 2 , . . .The average product G R (x, x ′ , ϵ)G A (x ′′ , x, ϵ − ℏω) includes all possible correlations between two multiple scattering paths starting at the initial points x ′ and x ′′ , respectively, and both ending at the final point x.The starting point of the Berezinskii technique is to take advantage of the one-dimensional geometry, which enables us to order the scattering events on a line: Thanks to this ordering, each contribution to the product G R G A may be represented by a diagram, like the one shown in Fig. 1(a) [24], which combines a retarded (solid lines) and an advanced (dashed lines) multiple scattering sequence, respectively unfolded in the upper and lower parts of the diagram.The scattering events occur at the points x i .In the example of Fig. 1(a), the upper sequence involves 8 scattering events (twice at point x 3 ), and the lower one 7 scatterings events.
The diagrams effectively contributing to the Green's function correlator G R G A do not have arbitrary shapes.Indeed, because of the assumed Gaussian statistics and the corresponding Wick's theorem, only diagrams whose scattering events can all be paired appear.For instance, the diagram in Fig. 1(a) vanishes upon averaging because some scattering events cannot be paired.Pairing scattering events at different points would occur in the case of a weakly correlated disorder (for which the correlation length is smaller than the mean free path), a problem previously addressed in [25,26,32].
The second important approximation of the Berezinskii technique is to assume that among all possible diagrams contributing to the correlator, only those for which the phase factors induced by free-particle Green functions exactly compensate each other when ω → 0 matter.This approximation, which holds true in the regime of "weak disorder" (see Sec. III C), amounts to imposing that there is exactly the same number of retarded and advanced Green's function in between any two successive scattering events x i and x i+1 .In turn, this yields restrictions on the possible scattering vertices, which are building blocks of the diagrams: aside the trivial constraint that the solid/dashed lines have to be continuous, all vertices have to be phaseless.
With these conditions implemented, the Berezinskii diagrammatic technique provides a strategy to exactly sum all possible diagrams with nonzero contributions, as we now detail for the case of a spin-orbit Hamiltonian H 0 .

III. DIAGRAMMATIC APPROACH WITHOUT TIME-REVERSAL SYMMETRY
A. Free Hamiltonian and Green's function In this work, we extend the standard Berezinskii technique to a one-dimensional spin system with spin-orbit coupling and Zeeman splitting breaking all anti-unitary time-reversal symmetries [18,28,29].The corresponding disorder-free Hamiltonian reads: where σ i are the usual Pauli matrices.The Hilbert space is spanned by two-dimensional complex-valued spinors.γ is the strength of the spin-orbit coupling, Ω is the Rabi frequency and δ the detuning.Diagonalization of the Hamiltonian H 0 yields two energy bands denoted by ± with corresponding energies Due to this band structure, for a given energy ϵ the Hamiltonian hosts either 2 or 4 possible eigenstates.
From now on, we focus on the case where only two eigenstates are involved, which corresponds to a dynamics operating at energies belonging to the lower band only [18].We denote by k ± the momenta of these two states, and by v ± = 1 ℏ |dE − (k ± )/dk| the associated velocities.As compared to the standard single-particle Hamiltonian H0 = p 2 /2m, it should be noted that the two involved momenta are not just of opposite sign, i.e., k − ̸ = −k + (and, correspondingly, v − ̸ = v + ).The left-right symmetry is therefore broken which, as will be seen below, constitutes the most significant difference as compared to the usual Berezinskii approach.In [27], a much simpler situation was studied, where only the spin-orbit interaction is present (i.e.δ = Ω = 0); in such a case, the dispersion relation is symmetric with respect to k → −k, so that v − = v + , and the extension of the Berezinskii technique is rather easy.In contrast, the calculations presented in the present paper are more general and valid when (generalized) TRI is broken.
In the diagrammatic treatment of disorder scattering introduced in the previous section, a fundamental ingredient is the free Green's function (3), which we need to evaluate for the Hamiltonian (5).To this aim, we use the definition A careful calculation of the integral in momentum space provides us with where, in particular, the diagonal value G R 0 (x, x, ϵ) is obtained by properly accounting for all the real and complex poles in the denominator of Eq. ( 7).Note that, strictly speaking, when x − x ′ ̸ = 0 these expressions only hold at distances |x − x ′ | larger than the de Broglie wavelength 2π/|k ± |.This knowledge, however, is sufficient within the weak disorder limit (15) where the Berezinskii approach operates.Because of translation invariance, the free Green's function G R 0 (x, x ′ , ϵ) = G R 0 (x − x ′ , ϵ).Its Fourier transform is therefore diagonal in momentum space, with the diagonal value defined as G R 0 (k, ϵ) = dr e −ikr G R 0 (r, ϵ), where r = x − x ′ .Note that with this definition, Eq. ( 8) implies that G R/A (k, ϵ) ̸ = G R/A (−k, ϵ), contrary to TRI systems.From Eq. ( 8), finally, the advanced Green's function follows from Hermitian conjugation, In the following, we will also need the energy-shifted Green's function G A 0 (x ′ , x, ϵ − ℏω), where ω ≪ ϵ/ℏ.To evaluate this object, we use the Taylor expansions k ± (ϵ− ℏω) ≈ k ± ∓ ω/v ± , so that:

B. Mean free times
Before constructing the diagrammatic approach based on the Hamiltonian H 0 + V (x), let us introduce a few important scattering parameters that will be used in the List of all possible initial vertices.In the limit ω → 0 they correspond to the factors (a) (ℏv+ following.The central one is the concept of scattering mean free time, which gives the average time scale between two consecutive scattering events.In the present case, however, two different mean free paths can be defined due to the left-right asymmetry. where ρ(E − (k ∓ )) is the density of states evaluated at the energy E − (k ∓ ) of the final state.Using that the disorder is uncorrelated, i.e., V (x ′ )V (x) = ηδ(x ′ −x) (see Sec. II), we infer: where κ ≡ |⟨χ + |χ − ⟩| 2 is the overlap factor of the two spin states.In the following we will be also led to use the mean free time associated with the weighted sum of the two scattering processes: which turns out to be the relevant time scale governing the boomerang effect, as will be shown in Sec.IV.Note that the validity of the Fermi golden rule used above is only guaranteed in the weak disorder limit described by Eq. ( 15) below.

C. Vertices
At the core of the Berezinskii diagrammatic technique is the idea to transfer the propagating factors from the free Green's functions to the scattering events, called vertices.For example, assuming x i > x j , the free Green's function can be split as where we formally associate the weights and exponential factors to the vertices at points x i and x j .The difference between TRI system and the TRI-broken case is that these factors depend on the direction of propagation.For example, in the TRI system the opposite case x j > x i would result in just a change of sign of the phase factors in Eq. ( 13), whereas in the system with broken TRI also the velocities change.
Initial vertices.We start by selecting the relevant initial vertices effectively contributing to the correlator G R (x, x ′ , ϵ)G A (x ′′ , x, ϵ − ℏω).In general, scattering paths may start from any of the 4 vertices shown with their weights in Fig. 2. The vertices with advanced and retarded lines starting into opposite directions, i.e., vertices (c) and (d), carry exponential factors with phases i(k ± x ′′ −k ∓ x ′ ).Upon integration over the starting points x ′ and x ′′ [cf.Eq. ( 29) in Sec.IV], they typically yield negligible contributions.Thus, we can restrict the analysis to only two classes of initial vertices: (a) and (b).These classes, in turn, correspond to two different types of initial states for the dynamics, (a) with positive (v + ) and (b) with negative (v − ) initial velocity.
A second simplification is based on the assumption that no scattering happens between the initial points x ′ and x ′′ [15,33].This invites us to introduce the Wigner variables r = x ′ − x ′′ and x = (x ′ + x ′′ )/2.In the limit ω → 0, vertices (a) and (b) are thus approximated by their counterparts starting from a single point x.At the level of Green's functions, this simplification reads [26]: where k ϵ is the wave number satisfying the dispersion relation ϵ = E − (k ϵ ).
Phaseless scattering vertices Our ultimate goal is to sum all significant contributions to the product of Green's functions G R G A .This formidable task is, in general, out of reach except in the so-called weak-disorder limit where ℓ = τ + v + = τ − v − is the scattering mean free path.Under this condition, only a limited set of scattering vertices that do not accumulate any phase and, as such, are not vanishingly small upon disorder averaging, should be considered when constructing correlation diagrams.The procedure to identify this set is detailed in Appendix A for clarity.It yields four families of vertices that are listed in Fig. 3.One can easily check that the phase associated with each vertex is zero.For instance, the vertex a 1 originates from a factor of the type ) in the disorder average of the Born expansion (2).With the help of the splitting procedure (13) and of Eq. ( 8), this corresponds to a vertex weight whose phase is indeed zero.In turn, the weights of all phaseless scattering vertices are Notice that among all diagrams in Fig. 3, the vertices families c, e, and f involve a "backscattering event" in both the retarded and advanced parts.This implies that, in the spin system described by Hamiltonian (5), the associated weights include the spin-state overlap factor κ = |⟨χ + |χ − ⟩| 2 .

D. Correlation diagrams
Knowing all possible phaseless scattering vertices relevant to our problem, we now wish to write down the equations describing the diagrams contributing to the correlator G R G A .An example of such a correlation diagram is shown in Fig. 1(b).Its generic structure can be divided into three left, right and central blocks denoted by L, R, and Z, as illustrated in Fig. 1(b).These different blocks are characterized by their total number of incoming and outgoing solid (retarded) and dashed (advanced) lines.
We first consider the left blocks L. Because the scattering vertices change the number of lines by at most 2, these blocks always have the same even number 2m ′ (with m ′ an integer) of retarded and advanced lines attached.For instance, L in Fig. 1(b) has m ′ = 1.With this property in mind, let us denote by L m ′ (x) the sum of contributions from all L blocks that have their right boundary at point x with 2m ′ lines.To calculate L m ′ (x), we consider how it changes with an infinitesimal change of the boundary position, x → x + δx, by counting all possibilities of adding new scattering vertices to L m ′ (x).This counting in detailed in Appendix B for clarity.Taking the limit δx → 0, it yields the following differential equation [34]: This equation is solved by an Ansatz , which leads to an iterative equation for L m ′ : where s = 2 − 2/κ + iν with ν = ω(v + + v − )ℏ 2 /κη.The explicit solution of Eq. ( 19) is with Ψ(a, b; z) the confluent hypergeometric function of the second kind.Note that in the usual case of spinless TRI systems, L m satisfies a similar equation as Eq. ( 19), but with s TRI = 2iωv/η, where v is the velocity of the state at energy ϵ [24].The main difference is that s TRI is fully imaginary, while in our case s has a finite real part.The treatment of the right block R is fully analogous.Denoting by R m (x) the sum of all right-hand blocks which have their left boundary at point x with 2m lines (with m an integer), we find that R m (x) = L m (−x) and, with a similar Ansatz, R m = L m .
Let us finally consider the central block Z.As compared to L and R, this block has one additional line which connects points x and x, i.e., for left and right block having 2m ′ and 2m retarded and advanced lines attached, the central block Z m ′ ,m connecting them has 2m ′ +1 lines at its left boundary and 2m+1 lines at its right boundary.For instance, the diagram in Fig. 1(b) has 2m ′ + 1 = 3 and 2m + 1 = 1.To derive a differential equation for Z m ′ ,m (x, x), we have to make an assumption on the direction of the extra line.Its type depends on the sign of x − x and introduces a kind of asymmetry because our problem differentiates left and right directions.Here, we assume x − x < 0, that is the additional line is going from left to right, like in Fig. 1(b).The total derivative of Z m ′ ,m (x, x) with respect to x includes the contributions from scattering vertices but also it has to include the derivative of the final vertex.These final vertices are analyzed analogously to the initial vertices.Out of four possibilities only two are phaseless, and thus contribute to the final sum of diagrams.They correspond to vertices with lines incoming only from a single direction, i.e., both from left or both from right.The list of all phaseless initial and final vertices is summarized in Fig. 4, together with their corresponding weights, denoted by Γ ±,. and Γ .,±for initial and final vertices, respectively.
Computing the total derivative of the central block at the final point x, assuming x < x, we find: The sign of the first term in the right-hand side depends on the final vertex type, i.e., Γ •,+ or Γ •,− .It turns out, on the other hand, that this expression does not depend on the sign of x − x.Note that when v + = v − and κ = 1, Eqs. ( 18) and ( 21) reduce to the known spinless TRI case [24].While we are not aware of any analytic solution for the differentio-recursive equation (21), in general the direct knowledge of the full function Z m ′ ,m (x, x) is not required for the computation of observables.An example of this will be given in the next section, when discussing the quantum boomerang effect.We conclude this section by expressing the Green's function correlator in Eq. ( 14) in terms of the blocks L, R and Z described above.For x < x, and if we sup-pose that the initial wave function only populates the state with initial velocity v + (this is the practical case that will be considered in Sec.IV), the correlator is the sum of two contributions corresponding to the two possible final vertices c and d in Fig. 4, with In the opposite case x > x, finally, Eq. ( 22) still holds but with Γ x<x +,± changed to Together with the solution of Eq. ( 21), Eqs.(22)(23)(24)(25)(26) constitute the final solution of the localization problem.
In the next section, we apply this formalism to access the time evolution of a particular observable, the mean position of wave packets, featuring the quantum boomerang effect.

IV. QUANTUM BOOMERANG EFFECT WITHOUT TIME-REVERSAL SYMMETRY
We now apply the above formalism to the theoretical description of a concrete problem, the quantum boomerang effect (QBE).We recall that the QBE describes a back-and-forth motion of the mean position of a quantum particle launched with nonzero initial velocity in a disordered potential.Here we describe this phenomenon based on the TRI-broken Hamiltonian H 0 + V , with the free part H 0 defined by Eq. (5).

A. Mean position
To describe the QBE within the Berezinskii technique, we consider for definiteness a wave packet initially launched in a disordered potential with the mean eigen wave number k + of the Hamiltonian (5) in the corresponding spin state |χ + ⟩.We denote by ϵ 0 = E − (k + ) the associated energy.We thus write the initial wave function as where σ is the wave-packet width.As explained in the previous sections, the ensuing dynamics of this state in the disorder gives rise to a coupling with the backwardpropagating state of wave number k − and spin component |χ − ⟩.
By definition, the disorder-average mean position is Using that ψ(x, t) = dx ′ G R (x, x ′ , t)Ψ 0 (x ′ ), we can relate its Fourier transform ⟨x(ω)⟩ = dte iωt ⟨x(t)⟩ to the Green's function correlator as where we expressed the retarded and advanced Green's functions in the Fourier domain.To simplify this expression, we make use of Eq. ( 14), which leads to with W the Wigner distribution of the initial state: For an initial wave function (27) of spatial width σ much smaller than the mean free path, we find where for convenience we replaced the integral over x by an integral over ∆x ≡ x − x, using that the integrand depends only on x − x due to statistical translational invariance.Equation (31) directly connects the average mean position to the Green's function correlator, which we now compute using the results of the previous section.

B. Time evolution of the boomerang effect
Inserting the general Berezinskii result (22) into Eq.( 31), we infer where is the contribution of velocities v + [technically, of the final vertex c in Fig. 4], and is the contribution of velocities v − [final vertex d in Fig. 4].Notice that we here introduced for convenience the mean free path ℓ = τ + v + = τ − v − .In Eqs. ( 33) and ( 34), the two terms in the right-hand side are the contributions of x < x and x > x, respectively, with the coefficients L m defined by Eq. (20).The quantities S i m , on the other hand, are given by spatial integrals of the block functions R m and Z m ′ ,m .For instance, we have To compute the coefficients S 0 m ′ , we perform a partial integration in the right-hand side and use Eq. ( 21) to express the derivative dZ m ′ ,m ′ (x, x)/dx in terms of Z m ′ ,m .This provides us with the iterative equation where β m ≡ (2κm 2 + 2m + 1)/(ℏ 2 v + v − ) and we remind that ν ≡ ω(v and is deduced from an iterative equation similar to Eq. (36): The coupled system of equations ( 36) and ( 38) is closed, so that at a formal level it can in principle be solved to find the coefficients S 0 m and, in turn, to compute the first sum in the right-hand side of Eq. (33).The calculation of the coefficients S 1 m , S 2 m and S 3 m that appear in Eqs.(33) and (34) follows the same lines.We provide the corresponding iterative equations they obey in Appendix C for the sake of completeness.
Instead of exactly computing all the S i m coefficients, however, the mean position can be conveniently evaluated from its short-time expansion using a Padé approximant [33,34].The short-time expansion of ⟨x(t)⟩ + and ⟨x(t)⟩ − is systematically obtained by inserting the series S i m (ν) = n s i m,n /(iν) n and Q i m (ν) = n q i m,n /(iν) n in the iteratives relations (36), (38) and (C1)-(C6), and computing the s i m,n and q i m,n coefficients at arbitrary order in 1/ν.This procedure eventually yields the following short-time expansion for the mean position: ) for a disordered system whose free Hamiltonian part is given by Eq. ( 5).The numerical simulations consist in a temporal propagation of the initial state (27) with the Shrödinger equation.For these simulations we choose σ = 50 for the wave-packet width, γ = 0.4, δ = Ω = 0.4 for the Hamiltonian parameters, and an energy ϵ0 = 0. We take v−/v+ = 0.66 for the ratio of velocities, and a disorder strength η = 0.0049.The simulations are done using a system of length L = 10000 with a small discretization ∆x = 0.2, and numerical results are averaged over 40960 disorder realizations.
where ∆ ≡ v − /v + and τ is defined by Eq. ( 12).In Appendix D, we also provide the corresponding expansions for the partial components ⟨x(t)⟩ + and ⟨x(t)⟩ − that respectively describe right and left-moving particles after the last scattering event.In Fig. 5 we show a comparison between an exact numerical calculation of ⟨x(t)⟩ ± based on a temporal wave-packet propagation with the disordered Schrödinger equation (details on the numerical simulations are given in the figure caption), and the short-time expansion up to order 11 obtained by solving the Berezinskii equations as explained above.Numerical and theoretical results are in very good agreement without any fit parameter up to t/τ ≈ 3.This corresponds to a finite radius of convergence in time, which is also present in the TRI version of the quantum boomerang effect [15].This radius seems to be dependent of the ratio of velocities ∆.Most importantly, as indicated by Eq. (39), the expression of ⟨x(t)⟩ is no longer universal because it depends on the velocities' ratio starting from the 4th order.This is a significant difference with the TRI quantum boomerang effect, which solely depends on the dimensionless time scale t/τ at all times.The TRI solution is fully recovered when It is also instructive to compare the exact, quantummechanical short-time expansion (39) with the classical prediction of the Boltzmann equation, which discards any interference in the multiple scattering process.At a classical level, the mean position is simply given by ⟨x(t)⟩ class.= τ v + (1−e −t/τ ), which is essentially the same expression as in TRI systems (see the supplemental material of [18] for details on the classical calculation).This  The initial state and parameters of the system are the same as in Fig. 5.
classical result has the short-time expansion which starts to deviate from the quantum-mechanical prediction (39) starting from the 4th order.For completeness, in Appendix D we also provides the short-time expansions for the classical components ⟨x(t)⟩ class.
+ and ⟨x(t)⟩ class.− .With the short-time expansion (39) at hand, we can infer the ⟨x(t)⟩ using a Padé approximant of the full Taylor series [35].To this aim we use that, at long time, ⟨x(t)⟩ ∝ 1/t 2 (see below).With this knowledge, we compute the mean position at any time using where A n (t) is a diagonal Padé approximant [35] whose coefficients are computed from the Taylor expansion at a desired order n [e.g., Eq. (39) for n = 4].In Fig. 6 we compare the exact numerical simulations for ⟨x(t)⟩ + , ⟨x(t)⟩ − and ⟨x(t)⟩ to the corresponding Padé approximants contructed from the Berezinskii technique.The plots reveal the QBE: after a few mean free times, the mean position exhibits a maximum and eventually decays to zero.For all quantities, we find a very good agreement between the simulations and the Berezinskii approach up to long times.Let us finally come back to the long-time behavior of ⟨x(t)⟩.The latter is most easily visualized in the inset of Fig. 6, which shows the mean position obtained from numerical simulations of the Schrödinger equation in log-log scale.We find that its long-time asymptotics is well approximated by a function scaling as α log(βt/τ )/(t/τ ) 2 , which is of the same form as in spinless TRI Hamiltonians [33].In the present case of Hamiltonian (5), however, a direct derivation of this asymptotic limit appears to be much more involved, and is left for future work.

V. CONCLUSION
In this paper, we have extended the Berezinskii diagrammatic technique describing the dynamics of Anderson localization in one dimension to TRI-broken disordered Hamiltonians including a spin-orbit coupling term that induces an asymmetry between right and left scattering processes.As an application of the formalism, we have computed the time evolution of the mean position of a wave-packet launched in a given direction, and recovered the quantum boomerang effect discussed in [18].As an extension of this work, it would be interesting to extract analytical long-time, asymptotic expansions for the mean position in this system, or to characterize the dynamics of other observables such as the mean square width ⟨x 2 (t)⟩ or the full density distribution of the wave packet.
Using the same procedure as for S 0 m and explained in the main text, we find the following coupled iterative equations for the S i m , Q i m (i = 1, 2, 3): Appendix D: Partial components ⟨x(t)⟩± We finally provide the short-time expansions for ⟨x(t)⟩ + and ⟨x(t)⟩ − , and the exact expressions (valid at any time) of their classical counterparts ⟨x(t)⟩

FIG. 1 .
FIG. 1.(a) Example of diagram involved in the product G R (x, x ′ , ϵ)G A (x ′′ , x, ϵ − ℏω) computed from the Born expansion (2).Solid lines represent free retarded Green's functions G R 0 , and dashed lines free advanced Green's functions G A 0 .Scattering events occur at points xi.(b) Example of non-vanishing contribution to the correlator G R (x, x, ϵ)G A (x, x, ϵ − ℏω).The wavy lines refer to a twopoint correlation function of the random potential.At each scattering event, one finds a vertex belonging to the set presented in Fig 3.The correlation diagram can be divided into three blocks L, Z, and R, separated by the points x and x.

FIG. 3 .
FIG. 3.All possible phaseless scattering vertices to be considered in the Berezinskii technique.Vertices from a, b and c families have dashed-lines counterparts.The weights associated to the vertices are indicated in Eq. (17).

FIG. 7 .
FIG. 7. Schematic representation of Eq. (B1).L m ′ (x + δx) can be constructed with L m ′ (x) and all possible combinations of the scattering vertices.The figure shows one example for each scattering vertex.
To find them, let us denote by |±⟩ = |k ± ⟩ ⊗ |χ ± ⟩ the two eigenstes of H 0 , where |χ + ⟩ and |χ − ⟩ are the spin state components associated with the wave numbers k + and k − , respectively.This leads us to define τ + and τ − , the scattering mean free times for the processes |+⟩ → |−⟩ and |−⟩ → |+⟩, respectively.At weak disorder they can be evaluated from the Fermi golden rule class.