Straightforward derivation yields the expression of OTOCs C t = C 1 t + C 2 t − 2 Re C 3 t , (5) where the three terms in right side are defined by C 1 t = ⟨ ψ R t 0 | p 2 m | ψ R t 0 ⟩ , (6) C 2 t = ⟨ φ R t 0 | φ R t 0 ⟩ , (7) and C 3 t = ⟨ ψ R t 0 | p m | φ R t 0 ⟩ , (8) with |ψ _R (t ₀)⟩ = U ^†(t _n , t ₀)θU (t _n , t ₀)|ψ(t ₀)⟩ and |φ _R (t ₀)⟩ = U ^†(t _n , t ₀)θU (t _n , t ₀)p ^m |ψ(t ₀)⟩.

To get the state |ψ _R (t ₀)⟩, one needs three steps: 1) forward evolution t ₀ → t _n , i.e., |ψ(t _n )⟩ = U (t _n , t ₀)|ψ(t ₀); 2) action of the operator θ on |ψ(t _n )⟩, i.e., | ψ ̃ ( t n ) 〉 = θ | ψ ( t n ) 〉 ; and 3) backward evolution t _n → t ₀, i.e., | ψ R ( t 0 ) 〉 = U † ( t n , t 0 ) | ψ ̃ ( t n ) 〉 . C ₁ (t _n ) (see Eq. 6) is just the expectation value of the p ^2m taken over the state |ψ _R (t ₀)⟩. For the numerical calculation of the state |φ _R (t ₀)⟩, one should first construct the operation of p ^m on the initial state |ψ(t ₀)⟩, i.e., |φ(t ₀)⟩ = p ^m |ψ(t ₀)⟩. Then, forward evolution from t ₀ to t _n yields the state |φ(t _n )⟩ = U (t _n, t ₀ )|φ(t ₀)⟩. At time t = t _n , the action of θ on the state |φ(t _n )⟩ results in a new state | φ ̃ ( t n ) 〉 = θ | φ ( t n ) 〉 , starting from which the time-reversal process t _n → t ₀ yields the state | φ R ( t 0 ) 〉 = U † ( t n , t 0 ) | φ ̃ ( t n ) 〉 . The norm of |φ _R (t ₀)⟩ is just the C ₂ (t _n ) (see Eq. 7). As the two states |ψ _R (t ₀)⟩ and |φ _R (t ₀)⟩ are available at the end of time reversal, one can calculate the C ₃ (t _n ) according to Eq. 8.

It is known that in the PT -symmetry breaking phase, the norm of quantum state N ψ ( t n ) = ⟨ ψ ( t n ) | ψ ( t n ) ⟩ exponentially increases for both the forward and backward time evolutions. To eliminate the contribution of norm to OTOCs, it is necessary to take the normalization for the time-evolved state. Specifically, for the forward evolution t ₀ → t _n , we set the norm of the quantum state equals to that of the initial state, i.e., N ψ ( t j ) = ⟨ ψ ( t 0 ) | ψ ( t 0 ) ⟩ with 0 ≤ j ≤ n. The backward evolution starts from the time t = t _n with the state | ψ ̃ ( t n ) 〉 , whose norm N ψ ̃ ( t n ) = ⟨ ψ ( t n ) | θ 2 | ψ ( t n ) ⟩ is expectation value of θ ² with the state |ψ(t _n )⟩. It is evident that the value of N ψ ̃ ( t n ) is important information encoded by the operation of θ on the state |ψ(t _n )⟩. Based on this, we take the normalization of the quantum state in the backward evolution t _n → t ₀ in such a way that its norm equals to N ψ ̃ ( t n ) , i.e., N ψ R ( t j ) = N ψ ̃ ( t n ) . One can find that for both the forward and backward evolutions, the norm of a time-evolved state always equals that of the state which the time evolution starts from. The same procedure of normalization is applied in calculating C ₂ (t _n ). Therefore, we have the equivalence N φ ( t j ) = ⟨ φ ( t 0 ) | φ ( t 0 ) ⟩ and N φ R ( t j ) = ⟨ φ ̃ ( t n ) | φ ̃ ( t n ) ⟩ (0 ≤ j ≤ n) for the forward evolution and time reversal, respectively.

We rewrite the C ₁ as C 1 t = 〈 ψ R t 0 | p 2 m | ψ R t 0 〉 = 〈 p 2 m t 0 〉 R N ψ R t 0 , (9) where N ψ R ( t 0 ) = ⟨ ψ R ( t 0 ) | ψ R ( t 0 ) ⟩ is the norm of the quantum state |ψ _R (t ₀)⟩ and ⟨ p 2 m ( t 0 ) ⟩ R = ⟨ ψ R ( t 0 ) | p 2 m | ψ R ( t 0 ) ⟩ / N ψ R ( t 0 ) indicates the exception value of p ^2m of the state |ψ _R (t ₀)⟩ with the division of its norm. We numerically investigate both the forward and backward evolutions of the norm N , and the mean values ⟨θ⟩ and ⟨p⟩ for a specific time, e.g., t = t ₁₀. It should be noted that we define the expectation value of observable Q as ⟨ Q ⟩ = ⟨ ψ ( t ) | Q | ψ ( t ) ⟩ / N ( t ) with N ( t ) = ⟨ ψ ( t ) | ψ ( t ) ⟩ . It is evident that such kind of definition eliminates the contribution of norm to mean value. Figure 3A shows that the norm is in unity during the forward time evolution (i.e., t ₀ → t ₁₀) and remains at a fixed value, i.e., N ψ R ( t 0 ) ≈ θ c 2 during the backward evolution (i.e., t ₁₀ → t ₀). For t ₀ → t ₁₀, the value of norm equals to that of the normalized initial state, so N ( t j ) = 1 . For the time reversal t ₁₀ → t ₀, our normalization procedure results in the equivalence N ψ R ( t j ) = ⟨ ψ ( t n ) | θ 2 | ψ ( t n ) ⟩ . Interestingly, our numerical investigations in Figures 4A, C, E show that for the forward evolution, the initial Gaussian wavepacket rapidly moves to the position θ _c = π/2; it should be noted that the initial Gaussian wavepacket has not moved to the position θ _c before t ₄. This is the reason why C(t) decays sharply before t ₄. During the time reversal, it remains localized at θ _c with the width of distribution being much smaller than that of the state of forward evolution. Correspondingly, the mean value ⟨θ⟩ has very slight differences with θ _c (see Figure 3A). Since the quantum state is extremely localized at θ _c , one can get the approximation N ψ R t 0 = ⟨ ψ t n | θ 2 | ψ t n ⟩ ≈ θ c 2 . (10)

Figures 4B, D, F show the momentum distribution of the state during both forward and backward evolutions. For the forward evolution, the quantum state behaves like a soliton which moves to a positive direction in momentum space, resulting in the linear increase of the mean momentum, i.e., ⟨p⟩ = Kt (see Figure 3B). The mechanism of the directed acceleration has been unveiled in our previous investigations [29,56]. Intriguingly, at time t = t ₁₀, the action of θ yields a state with a power-decayed shape, i.e., | ψ R ( p , t 0 ) | 2 ∝ ( p − p c ) − 2 (see Figure 4F). Most importantly, during the backward evolution, the quantum state still retains the power-law decayed shape, for which the center p _c decreases with time and almost overlaps with that of the state of the forward evolution. This clearly demonstrates a kind of time reversal of transport behavior in momentum space.

In the aspect of the mean momentum ⟨p⟩, we find that the value of ⟨p⟩ linearly decreases during the backward evolution and is in perfect symmetry with that of the forward evolution, which is a solid evidence of time reversal. In the end of the backward evolution, the quantum state |ψ _R (t ₀)⟩ is localized at the point p = 0 (see Figure 4B). By using the power-law distribution |ψ _R (p, t ₀)|² ∼ p ⁻², it is straightforward to get the estimation of the expectation value of p ^2m , i.e., ⟨ p 2 m ⟩ ψ R = ∫ p − N / 2 p N / 2 p 2 m | ψ R ( p , t 0 ) | 2 d p ∝ N 2 m − 1 . Taking both ⟨ p 2 m ⟩ ψ R and N ψ R ( t 0 ) in Eq. 10 into Eq. 9 yields the relation C 1 t ∝ N 2 m − 1 θ c 2 , (11) which is verified by our numerical results in Figure 5. As an illustration, we consider the cases with m = 1, 2, and 3. Our numerical results of the late-time saturation values of C ₁ are in good agreement with Eq. 11. It is now clear that the scaling of C(t) with N originates from the power-law decay of the state |ψ _R (t ₀)⟩. The reason for the formation of power-law decayed wavefunction has been uncovered in Ref. 65.

We proceed to evaluate the time dependence of C ₂(t) in Eq. 7, which is just the norm of the state |φ _R (t ₀)⟩ at the end of backward evolution. According to our normalization procedure, the value of C ₂ equals to the norm of the state | φ ̃ ( t n ) 〉 = θ | φ ( t n ) 〉 , hence C 2 = 〈 φ t n | θ 2 | φ t n 〉 = 〈 θ 2 〉 N φ t n , (12) with N φ ( t n ) = ⟨ φ ( t n ) | φ ( t n ) ⟩ and ⟨ θ 2 ⟩ = ⟨ φ ( t n ) | θ 2 | φ ( t n ) ⟩ / N ( t n ) . We numerically find that the state |φ(t _n )⟩ is extremely localized at the position θ _c during the forward evolution (see Figure 6). Then, a rough estimation yields ⟨ θ 2 ⟩ ∼ θ c 2 . The norm N φ ( t n ) equals that of the initial state |φ(t ₀)⟩ = p ^m |ψ(t ₀)⟩. By using the initial Gaussian wavepacket ψ ( p , t 0 ) = ( 1 / σ ℏ eff 2 π ) 1 / 4 ⁡ exp ( − p 2 / 2 σ ℏ eff 2 ) , one can straightforwardly obtain N φ t n = ∫ − ∞ ∞ p 2 m | ψ p , t 0 | 2 d p = 2 m − 1 ‼ 2 m α m , where α = 1/(σℏ ²) and (…)!! denote a double factorial. Taking both the ⟨θ⟩ and N φ ( t n ) into Eq. 12 yields the late-time saturation value C 2 t ∼ θ c 2 2 m − 1 ‼ 2 m α m , (13) which is in good agreement with our numerical results in Figure 5.

The value of C ₃(t) depends on both the states |ψ _R (t ₀)⟩ and |φ _R (t ₀)⟩ (see Eq. 8). Figure 7 shows the probability density distributions of the two states in both the real space and momentum space. For comparison, the two states are normalized to unity. One can find the perfect consistence between |ψ _R (t ₀)⟩ and |φ _R (t ₀)⟩. Then, we roughly regard C ₃ as the expectation value of the p ^m taking over the state ψ _R (t ₀) or φ _R (t ₀), i.e., C 3 ( t ) ≈ ⟨ p m ( t 0 ) ⟩ ψ R N ψ R ( t 0 ) N φ R ( t 0 ) , where according to aforementioned derivations N ψ R ( t 0 ) = θ c 2 and N φ R ( t 0 ) = C 2 ( t ) . By using the power-law decayed wavepacket |ψ _R (t ₀)|² ∝ p ⁻², one can obtain the estimation ⟨ p m t 0 ⟩ ψ R ≈ ∫ p − N / 2 p N / 2 p m | ψ R p , t 0 | 2 d p ∼ 0 for odd m , N m − 1 for even m . . (14) Accordingly, the C ₃ is approximated as C 3 t ∼ 0 for odd m , η N m − 1 for even m . . (15) with the prefactor η ∝ θ c 2 ( 2 m − 1 ) ! ! / 2 m α m 1 2 .

We numerically calculate the absolute value of the real part of C ₃. Interestingly, our numerical results of |Re [C ₃]| is in good agreement with the analytical prediction in Eq. 15 (see Figure 5), which proves the validity of our theoretical analysis. We further numerically investigate the |Re [C ₃(t)]| at a specific time for different N. Figure 8 shows that for B = p, the value of |Re [C ₃(t)]| is nearly zero with varying N, which is consistent with our theoretical prediction in Eq. 15. For B = p ³, the value of |Re [C ₃(t)]| has slight difference with zero for large values of N, signaling the derivations with Eq. 15. This is due to the fact the quantum state |ψ _R (t ₀)|² is not exactly symmetric around p. In order to quantify such asymmetry, we numerically investigate the difference of the sum of the probability between the positive and negative momentums Δ ρ = ∑ n = 0 N 2 − 1 | ψ R ( t 0 , p n ) | 2 − ∑ n = − N 2 − 1 | ψ R ( t 0 , p n ) | 2 and find that it is non-zero Δ _ρ = 0.13. Interestingly, for B = p ², the value of |Re [C ₃(t)]| increases linearly with increasing N, which is clear evidence of the validity of our theoretical prediction.