Periodically and Quasi-periodically Driven Dynamics of Bose-Einstein Condensates

We study the quantum dynamics of Bose-Einstein condensates when the scattering length is modulated periodically or quasi-periodically in time within the Bogoliubov framework. For the periodically driven case, we consider two protocols where the modulation is a square-wave or a sine-wave. In both protocols for each fixed momentum, there are heating and non-heating phases, and a phase boundary between them. The two phases are distinguished by whether the number of excited particles grows exponentially or not. For the quasi-periodically driven case, we again consider two protocols: the square-wave quasi-periodicity, where the excitations are generated for almost all parameters as an analog of the Fibonacci-type quasi-crystal; and the sine-wave quasi-periodicity, where there is a finite measure parameter regime for the non-heating phase. We also plot the analogs of the Hofstadter butterfly for both protocols.


Introduction
The rapid development of experimental techniques in cold atom systems draws a lot of attention to interesting quantum dynamics. An important scenario is to periodically drive the system. When the lattice is modulated periodically in time, one could generate topologically non-trivial band insulators [1][2][3][4][5][6][7][8] or introduce gauge field with non-trivial dynamics [?,9,10]. However, the quantum dynamics are usually hard to solve. Special cases where large simplification occurs due to symmetry reason is especially valuable.
The driven dynamics of Bose-Einstein condensates (BECs), which has been realized in [11][12][13][14], is one of such examples. As pointed out in [15] and later studied in [16,17], the dynamics of the Bose-Einstein condensates under the Bogoliubov Hamiltonian exhibit an SU(1, 1) structure, which, more generally speaking, is a consequence of the bosonic Bogoliubov transformations. The dynamical evolution of the system can then be described by the multiplication of a sequence of SU(1, 1) matrices, 1 acting on two components pseudospin (a, a † ) for bosonic creation and annihilation operators. The observable and the phase diagram can then be obtained by analyzing the corresponding SU(1, 1) matrices. The techniques applied in this work is a close analogy of those in the driven conformal field theory (CFT) in 1 + 1-D [22][23][24][25][26][27]. In the CFT case, the system has the Virasoro symmetry and for specific driving Hamiltonians, one could focus on the SL(2, R) ∼ = SU(1, 1) subgroup. Another analogy is the dynamics of the scale-invariant gases in the harmonic trap [28,29], as would be discussed in more detail later.
In this work, we consider two types of protocols for both the periodic and the quasi-periodic drivings. The two types correspond to modulating the scattering length 2 , which is equivalent to the coupling strength in our BEC setting, by a square-wave or a sine-wave in time. For the square wave protocol, in each period the evolution is determined by multiplying a finite number of matrices while for the sine-wave protocol we in general need to solve a Mathieu's equation, however simplifications can be made by Magnus expansion in the near-resonant limit [15]. In both protocols, the long-time dynamical phase diagram is then determined by the matrix element of the SU(1, 1) evolution matrix/Bogoliubov transformation in a single period as in [25]. We then turn to the driving dynamics with quasi-periodicity. For the square-wave quasiperiodic driving, the excitations grow exponentially for almost all parameters, while special points with power-law heating also exists [25] and will be discussed. For the sine-wave type quasi-periodicity, we consider the scattering length to be a superposition of two sinusoidal functions with incommensurate frequencies, where the evolution equation is a quasi-periodic Mathieu's equation [31,32]. 3 The paper is organized as follows: In section 2, we explain the SU(1, 1) group structure of the Bogoliubov Hamiltonian for a driven Bose-Einstein condensate focusing on the Heisenberg evolution of creation/annihilation operators. In section 3, we study the two protocols when the 1 Similar strategy for the case with SU(2) spin dynamics has been discussed in Ref. [18][19][20][21]. 2 One can also consider a time-dependent Zeeman field as in e.g. Ref. [30]. 3 See e.g. [33] for a review from the perspective of applied mathematics. scattering length a s (t) is a periodic function in time and the quasi-periodic case is discussed in section 4. Finally, we summarize our work and comment on the Bogoliubov approximation in section 5.

Dynamics of Bose-Einstein Condensates
We consider dilute single-component bosonic atoms and keep the s-wave interactions only. Taking the time-dependence of the scattering length into account, the Hamiltonian reads: Here we have set the mass of atoms m = 1 for simplicity. At the bare level, the scattering parameter g(t) can be related to the scattering length a s as When study dynamics around a specific state with large number of particles in the condensate, we expand ψ(r) = √ n 0 +δψ(r), where n 0 is the density of condensate particles. To the quadratic order, this gives the Bogoliubov Hamiltonian: where the constant term is only a function of total particle number N = dr(n 0 + δψ † δψ) and therefore conserved. The validity of using a simple Bogoliubov Hamiltonian to study the dynamics of condensates has been verified by comparing to a careful numerical study which includes the evolution of condensate wave function [15,34,35]. The Hamiltonian can also be written in the momentum space via Fourier expansion δψ(r) = 1 √ V k =0 e ik·r a k with V the system size. Thus, we have There are two ways to recognize the SU(1, 1) structure in this system (for each k):

Using the notation
, we can rewrite the time-dependent Hamiltonian as a linear combination Further introducing A y k = 1 2i (a † k a † −k − a −k a k ), we have the closed commutation relation that reassembles the su(1, 1) algebra [15,16] 2. We can also directly study the evolution of excitation operators a k and a † k . For a general g(t) and the corresponding unitary U (t) = T exp(−i t 0 H(t )dt ), the time-evolution generates a Bogoliubov transformation between a † k and a −k as follows In the first step, the unitary U (t) and U † (t) act on the operators a k and a † −k , while in the last step, the 2 × 2 transform matrix U k acts on the basis vector (a k a † −k ) (where we have used the inversion symmetry α k = α −k , β k = β −k to rewrite the second column of U k ). Further note that the (bosonic) commutation relation a k , a † k = 1 demands |α k (t)| 2 − |β k (t)| 2 = 1, therefore the Bogoliubov transformation matrix U k is an SU(1, 1) matrix.
Moreover, when the evolution consists of several steps U (t) = U n U n−1 ...U 1 , we have: That is to say, the full evolution U k (t) = U k,1 ...U k,n−1 U k,n .
In this paper, we focus on the periodic and quasi-periodic driving of the Bose-Einstein condensate |ψ(t) = U (t)|ψ(0) , and measure the total number of excitations n k (t) in the final state. For concreteness, let us choose the initial state |ψ(0) to be the fully condensed state with a k |ψ(0) = 0, 4 then after the evolution, we have which is directly related to the matrix element of the SU(1, 1) matrix U k (t). Besides |β k (t)| 2 , one could also consider measuring other elements of U k , by adding additional evolution right before the final measurement.
With the driving, particles in the condensate can be pumped into excitations with finite momentum. We call the system in the heating phase at certain momentum k when the occupation number grows exponentially as n k (t) ∼ e λ k t in the long-time limit, and we call the correspond exponent λ k > 0 the heating rate. From (9), this implies |β k (t)| 2 ∼ e λ k t , so does |α k (t)| 2 ∼ e λ k t due to the constraint |α k (t)| 2 − |β k (t)| 2 = 1. In other cases, n k may oscillates or increases polynomially, which we will refer as non-heating or critical point respectively.
As a side note, we would like to point out that the present discussion can be directly applied to specific dynamics of scale-invariant atomic gases in a trap [28,29]. To see the connection, let us consider a single harmonic oscillator 5 with the trapping frequency ω = 1 and a = 1 Then we can define the following set of generators that satisfies the aforementioned su(1, 1) algebra: Consequently, for Hamiltonians given by a linear superposition of these operators, the dynamics can again be studied using the SU(1, 1) group transformation. In the original basis, we can regroup the above generators into p 2 , x 2 and xp + px. The generalization to many-body scaleinvariant atomic gases, e.g. the unitary Fermi gas, can then be accomplished by noticing that these operators are a subgroup of the non-relativistic conformal group [36], where p 2 , x 2 and xp + px are generalized to time translation, special conformal transformation and dilation.

Periodic driving
In this section, we consider periodic drivings where g(t) = g(t + T ). We will discuss the phase diagram that is experimentally relevant. Another purpose for this section is to setup a base for the discussion in the quasi-periodic case. Let us assume the total evolution time t = nT with integer n, we have Here H F is the standard Floquet Hamiltonian. From the perspective of the Heisenberg evolution of excitation operator (8), we have correspondingly U k (nT ) = U k (T ) n , and we are interested in the large n asymptotics of U k (nT ), which determines the dynamical phases of the condensate, namely the heating, non-heating phase and the phase boundary (critical phase) between them.
The determination of the phase diagram has two steps.
1. We compute U k (T ) for a given protocol specified by g(t): The Heisenberg equation for a k and a † −k is given as follows or equivalently, with anti-time ordering operator T . In addition to performing the exact evaluation of U k (T ), several standard approximations such as Magnus expansion [15] can also be applied, as discussed later in this section.
2. Analyze the large n behavior of U k (T ) n based on a single U k (T ).
Assuming we have worked out Then the matrix element α k (nT ) and β k (nT ) of U k (nT ) can be related to the α k (T ) and β k (T ) via following formulas Here we have defined Note From (15), we conclude that the β k (nT ) grows exponentially only when |Tr(U k (T ))| > 2, and the exponent reads When |Tr(U k (T ))| < 2, the η becomes a pure phase and n k (nT ) oscillates periodically with frequency ω k = 2|δ k |/T . When δ k /π = p/q is a rational number, we have n k (nqT ) = 0. At the phase boundary with |Tr(U k (T ))| = 2, n k (nT ) grow quadratically if we have β k (T ) = 0. In the following, we consider two concrete protocols, and determine the phase diagram in terms of the experimental parameters.

Protocol 1: square-wave
In this protocol, the scattering length is tuned to be a square wave as shown in Figure 1(a). We then have To avoid instability, we assume both g j 0. Using (12), we have where We then compute U k (T ) = U k,1 (T 1 )U k,0 (T 0 ) and find When |k| is large, we expect the effect of g being small and no particles can be excited. This is consistent with (20), where we have Tr(U k (T )) ≈ 2 cos(k 2 (T 0 + T 1 )/2) 2 and β k (T ) ≈ 0. For general |k|, we compute (20) numerically and plot |Tr(U k (T ))| in Figure 1 (b) which can be served as a phase diagram. To be concrete, we fix T 0 = T 1 and g 1 = 0. At small g 0 , the heating phase appears near the resonance k 2 T = 2nπ, where 2 particles can be excited by the periodic modulation.

Protocol 2: sine-wave
In this protocol, we let g(t) be a single or a combination of sine functions. We first consider the case where g(t) = 2g 0 cos(ω 0 t) and ω 0 = 2π/T . Using (12), we further write out the set of equations satisfied by α k (t) and β k (t): This leads to [16] Regarding t as the spatial dimension, this is the Schrödinger equation for a particle in 1D with mass m = 1/2 and energy E = k 4 4 , moving in a potential with V (t) = −2k 2 n 0 g 0 cos(ω 0 t). The system would be in the non-heating phase if the corresponding energy lies in some bands and otherwise it will be in the heating phase. The detailed phase diagram has been worked out in [16]. Here we instead focus on the case with large ω 0 g 0 n 0 and near resonance ω 0 −k 2 ω 0 , where we could perform a rotation-wave approximation to (12) and find Here ∆ k = k 2 /2 − ω 0 /2 and the average is performed over time, i.e. O(t) := 1 T T 0 O(t)dt. Systematic improvements can be made by taking into account the 1/ω 0 correction. In the rotating frame, U k (T ) now has the same form as (19), with A k = ∆ k , B = g 0 n 0 and E k = ∆ 2 k − g 2 0 n 2 0 , leading to: When ∆ k < g 0 n 0 , the system is in the heating phase, with a heating rate λ k = ∆ 2 k − g 2 0 n 2 0 . On the other hand, for ∆ k > g 0 n 0 , the system is in the non-heating phase, with ω k = g 2 0 n 2 0 − ∆ 2 k .
We then consider a combined protocol where we divide each period into two halves as in the square wave case, as shown in Figure 2 (a). In the first half period, we modulate the scattering length as g(t ∈ [0, T /2)) = 2g 0 cos(ω 0 t) and in the second period we add an additional phase shift g(t ∈ [T /2, T )) = 2g 0 cos(ω 0 t + φ). We assume T = 4πn/ω 0 , with n ∈ Z. Due to the additional phase, we now have Therefore the trace of U k (T ) = exp(−i M k,1 T /2) exp(−i M k,2 T /2) is expressed as follows, We plot the numerical result of |Tr(U k (T ))| in Figure 2 (b). For ∆ k = 0, the system transits into the non-heating phase at φ = π for large ω 0 , as observed in the experiment [13]. Note that at this point, the second half of the evolution is canceled exactly with the first half, but this exact cancellation is not the general picture for a generic phase boundary (Tr(U k (T )) = 2). When including 1/ω 0 corrections, the phase boundary would shift, as studied in [15].

Quasi-periodic driving
Now we turn to the quasi-periodically driven condensate where g(t) is deterministic but has no periodicity. We consider both square-wave and sine-wave protocols as in the periodic driving case. The square-wave protocol is a time domain analog of the Fibonacci model for the one dimensional quasi-crystal (but generalizing to an arbitrary irrational number following [37,38]), while the sine-wave protocol corresponds to the quasi-periodic Mathieu's equation, whose lattice version is also known as the Audry-André model [39].

The set-up and the trace map
For this square-wave quasi-periodicity protocol, we define a sequence of 0/1 bits with an irrational number α ∈ (0, 1) Here y is the floor function which gives the greatest integer less than or equal to y. 6 For example, for the inverse golden ratio α = which is also known as Fibonacci word. 6 This definition is equivalent to the other definition that is commonly used in the quasi-crystal literature where χ [1−α,1) (t) = χ [1−α,1) (t + 1) is a periodic characteristic function with period 1, namely χ [1−α,1) (t) = 1 if 1 − α t < 1, and χ [1−α,1) (t) = 0 if 0 t < 1 − α. See e.g. [40,41] for the Fibonacci model when α is the inverse golden ratio.
Next, we evolve the system according to the sequence: if the bit is "1/0", we apply the unitary U 1/0 respectively, with U 1 = T exp(−i T 0 dtH 1 (t)) and U 0 = T exp(−i T 0 dtH 0 (t)). Again, for the inverse golden ratio example, we have ) or in terms of U k,1 and U k,0 : Note the order is reversed according to the definition (8). The evolution can be described more efficiently if we only probe the system at specific time. This simplification relies on the continued fraction representation for α: where {a n } n=1,2,3... are a sequence of positive integers, and a 0 = α could be negative or 0, for our case α ∈ (0, 1), we have a 0 = 0. Any real number has an unique continued fraction representation, for rational numbers, the continued fractions terminate at finite n, while the irrational numbers have infinite sequences [a 0 , a 1 , a 2 . . .].
Continued fractions are useful in finding the best Diophantine approximations. Operationally, we can truncate the continued fraction at order n, which produces a rational number (known as n-th principal convergent) where integers (p n , q n ) satisfies the recursion relation q n = a n q n−1 +q n−2 and p n = a n p n−1 +p n−2 . For example, for the inverse golden ratio α = √ 5−1 2 , we have a 0 = 0, a n>0 = 1 and q n are given by the Fibonacci numbers: q 1 = 1, q 2 = 2 and q n = q n−1 + q n−2 .
These rational numbers p n /q n provide best rational approximations of the irrational α in the sense that the product q n α is closer to an integer than any smaller q < q n (see e.g. [42] for more explanations). Consequently, one can show that 7 for V j defined in (27). This relation says the first q n elements of the sequence can be determined by copying the first q n−1 elements. Thus, we have the following recursion relation An important analytic tool that allows us to efficiently determine the heating rate is the trace map [40,37], which is a recursion relation among the following traces Tr (U k (q n T )U k (q n−1 T )) , y n k = 1 2 Tr (U k (q n T )) , z n k = 1 2 Tr (U k (q n−1 T )) . The the heating rate λ k T with fixed g 0 n 0 T = 2 determined by taking a period of q 50 ∼ 2 × 10 10 steps. The results show self-similar structures as in [25].

sin[arccos(y)]
is the Chebyshev polynomial of second kind (see appendix for a derivation). Using the trace map (37), one can straightforwardly show that is independent of n, hence an integral of motion (Fricke-Vogt invariant). When I > 0, for almost all initial conditions, (x n k , y n k , z n k ) flows to infinity as n increases.
Recall that we have k = k 2 /2, A k = k + g 0 n 0 , B = g 0 n 0 and E k = A 2 k − B 2 . In this protocol, the heating rate λ k can be approximated by taking the large n limit of periodic driving with a period q n T : we could approximate α ≈ p n /q n and use (27) to define the driving sequence within a period. Following (17), this gives with sufficiently large n. Generally q n grows exponentially with n and we can obtain accurate results for n being a few tenth.

Patterns in phase diagram
To proceed, we also need to choose the irrational number α, namely choose the modulation pattern that is determined by V j sequence. We will exam two specific examples in this subsection (1) α = √ 5−1 2 , which is the most popular choice in the study of 1D quasi-crystal and has the name "Fibonacci model" [40,41]; (2) α = π − 3.
(1) Fibonacci Driving. Let us start with the inverse golden ratio α = √ 5−1 2 , which is called the Fibonacci driving in [25]. In this case, since a n>0 = 1, we have One can extend the recursion relation (37) to n = 2 by defining U k (q 0 T ) ≡ U k,0 , which gives and consequently The I = 0 case where either sin(E k T ) = 0 or sin( k T ) = 0 corresponds to a single quench (shown as black dashed line in Fig. 6) and we should focus on the I > 0 case. It is known that when I > 0, the phase with zero heating rate λ k forms a Cantor set of measure zero in the phase diagram [38] and the system would generally be heated up exponentially in time. This is consistent with the numerical results shown in Figure 3. As shown in (b), the heating rate is non-zero almost everywhere, although the magnitude can be arbitrarily small with a self-similar structure as shown in (c-d). Later we will comment, with a comparison to the π-driving, the self-similarity in the Fibonacci driving is a manifestation of the self-similarity in its continued fraction representation Although the measure of phases with λ k = 0 is zero, special points with λ k = 0 are known and interesting. When two of the initial conditions (x n k , y n k , z n k ) become zero, y n k is then a periodic function of n with a period of 6. As an example, for (x 1 k , y 1 k , z 1 k ) = (a, 0, 0) we have It is straightforward to show that the number of excited particles oscillates periodically with respect to n at time q n T and the period is 3. On the other hand, if one probe the system at non-Fibonacci number times, there are still excitations and the number of excitations grows polynomially. An explicit example is shown in Fig. 4 (a) (with k T = 1.5π and E k T = 2.5π) compared to the exponentially heating (b) if we slightly perturb away from the special points (with k T = 1.5π + 0.01 and E k T = 2.5π).
(2) π Driving. Now we consider the choice α = π − 3 ∈ (0, 1), the point is that such irrational number has a rather irregular continued fraction representation, unlike the case for the inverse gold ratio, whose continued fraction representation has self similar pattern. We will show numerical that the patterns in the phase diagram reassemble the patterns in the continued fraction of π − 3: a 1 = 7, a 2 = 15, a 3 = 1, a 4 = 292 . . .
that is to say π − 3 = 1 7 +  (e-f). The heating rate for rational approximation x = p 3 /q 3 = 16/113. By comparing (b-f), we find p 1 /q 1 and p 3 /q 3 governs the heating rate at different magnitude scales. Figure 6: The non-heating parameters for g 0 n 0 T = 2. We consider all rational α = p/q with q 40. Black dashed lines correspond to I = 0.
In Fig 5 (a), we plot the density of log(λ k T ) as a function of two parameters (g 0 n 0 T, k 2 T ). In Fig 5 (b), we show a specific cut at g 0 n 0 T = 2 and demonstrate the expectation that the regime with zero heating rate λ k = 0 forms a Cantor set of zero measure. 8 We further zoom in a small regime in Fig 5 (c) and show that the self-similarity feature that was observed in Fibonacci driving is absent in the π-driving. On the contrary, what we see here in the π-driving is that the heating rate shows different pattern at different scales.
To further understand the formation of the structure at different scales, we show in Fig 5  (d) the heating rate distribution with p 1 /q 1 = 1/7, namely the first principal convergent that approximates the irrational α = π − 3. As excepted, the large scale structure in Fig 5 (b) is captured while details are lost. In Fig 5 (e) and its zoom (f), we use the third principal convergent p 3 /q 3 = 16/113, and find the patterns are well reproduced, even at the scale with magnitude of λ k ∼ 0.01.
A simple understanding of the above phenomenons can be achieved via analogy with the band-theory of lattice Hamiltonian. Starting from some p n−1 /q n−1 , we have a band structure where the in-gap states corresponds to the heating phase. When we increase n, the band folds and new gap is opened. Intuitive, after we repeat this for infinite many times, the λ k = 0 regime (corresponds to band) breaks into a Cantor set. Now, from the definition (32), we know that when some a n is very large, q n is much larger than q n−1 . Consequently, p n−1 /q n−1 becomes a particularly good approximation to α at certain magnitude scale: the correction when taking a n into account comes from folding the Brillouin zone a n times, which is expected to be small when a n 1. 2 . (c) The non-heating parameters for g 0 n 0 /ω 0 = 0.5 and g 1 n 0 /ω 0 = 0.1. We consider all rational ω 1 /ω 0 = p/q with q 40

Hofstadter Butterfly
Knowing that the system is almost always in the heating phase for an irrational α, we now give an overall picture for the phase diagram with varying α ∈ (0, 1). In Fig 6, we fixed g 0 n 0 T = 2 and draw all the non-heating region with α = p/q and q 40. The result resembles the Hofstadter butterfly [43]. Note that in our set-up, there is no reflection symmetry α → 1 − α, due to the asymmetry between U 0 and U 1 .

Protocol 2: sine-wave
Now we turn to the second type of quasi-periodic driving. We consider g(t) = 2g 0 cos(ω 0 t) + 2g 1 cos(ω 1 t). The driving is quasi-periodic when ω 1 /ω 0 is a irrational number. Again, the evolution is governed by (22), which we copy here: with the initial condition The model now corresponds to a 1D quantum mechanics in quasi-periodic potential [31,32], also known as quasi-periodic Mathieu's equation in applied mathematics, e.g. see [33] for a review. Now we fix an irrational ω 1 /ω 0 = √ 5−1 2 and g 0 n 0 /ω 0 = 0.5. We numerically solve (49) and fit the long time behavior to get the heating rate λ k . The result is shown in Figure 7 (b). For g 1 = 0 without quasi-random potential, we have non-heating phases around k 2 /ω 0 = 0, 1.6, 3, which corresponds to bands in a sine potential. When g 1 becomes finite, we expect a finite k 2 g 1 is needed for states to heat up. This is consistent with the existence of extended blue regions in Figure 7 (b) around k 2 = 0 and g 1 = 0.
Interestingly, we also see strip structures around the non-heating regions. To further explore structure, we focus on g 1 n 0 /ω 0 = 0.1 and change different ω 1 /ω 0 . We consider different rational numbers ω 1 /ω 0 = p/q, and plot all the non-heating regions on the ω 1 − k 2 plane. The result is shown in 7(c) which resembles some feature of the Hofstadter butterfly [43] for small but finite k 2 .

Summary
In this paper, we study two protocols of the periodically and quasi-periodically driven dynamics of Bose-Einstein condensates. We determine the phase diagram in terms of whether the system is heated or not, whether the heating phase has an exponential growing number of excitations. For the quasi-periodically case with a square-wave modulation, phase with λ k = 0 form a Cantor set of measure zero in the parameter space. We also exam a special non-heating point and find the number of excitations grows algebraically instead of exponentially. On the contrary, for the sine-wave quasi-periodic driving case, there is a finite measure regime where the system is in the non-heating phase. We also find analogs of the Hofstadter butterfly for both protocols. We expect that these experiments can be carried out with minor modification using the experimental system in [13].
Finally, we would like to remark on the validation of the Bogoliubov approximation. By keeping the condensate wave function as a constant, we are assuming the number of particles in the initial condensate is much larger the number of excitations. When the system is in the heating phase, this would ultimately break down. However, since this only happens in the late time, one expect the phase boundary receives little correction. This is again similar to the CFT case [25], where the validation is checked by comparing to a lattice simulation.
The third equation here is merely a definition, we just need to prove the first two equations. We have U k (q n T ) = U k (q n−1 T ) an U k (q n−2 T ) = S an−1 (y n−1 k )U k (q n−1 T ) − S an−2 (y n−1 k ) U k (q n−2 T ). (57) Taking the trace leads to y n k = S an−1 (y n−1 k )x n−1 k − S an−2 (y n−1 k )z n−1 k . Similarly, we consider U k (q n−1 T )U k (q n T ) = U k (q n−1 T ) an+1 U k (q n−2 T ) = S an (y n−1 k )U k (q n−1 T ) − S an−1 (y n−1 k ) U k (q n−2 T ), whose trace leads to x n k = S an (y n−1 k )x n−1 k − S an−1 (y n−1 k )z n−1 k .