Dynamically generated quadrupole polarization using Floquet adiabatic evolution

We investigate the nonequilibrium dynamics of the $S=1$ quantum spin chain subjected to a time-dependent external drive, where the driving frequency is adiabatically decreased as a function of time (``Floquet adiabatic evolution''). We show that when driving the rhombic anisotropy term (known as the ``two-axis countertwisting'' in the context of squeezed spin states) of a N\'eel antiferromagnet, we are able to induce an overall enhancement in the quadrupole polarization, while at the same time suppressing the staggered magnetization order. The system evolves into a new state with a net quadrupole moment and antiferroquadrupolar correlations. This state remains stable at long times once the driving frequency is kept constant. On the other hand, we find that we cannot achieve a quadrupole polarization for the symmetry-protected Haldane phase, which remains robust against such driving.

We investigate the nonequilibrium dynamics of the S = 1 quantum spin chain subjected to a time-dependent external drive, where the driving frequency is adiabatically decreased as a function of time ("Floquet adiabatic evolution"). We show that when driving the rhombic anisotropy term (known as the "two-axis countertwisting" in the context of squeezed spin states) of a Néel antiferromagnet, we are able to induce an overall enhancement in the quadrupole polarization, while at the same time suppressing the staggered magnetization order. The system evolves into a new state with a net quadrupole moment and antiferroquadrupolar correlations. This state remains stable at long times once the driving frequency is kept constant. On the other hand, we find that we cannot achieve a quadrupole polarization for the symmetry-protected Haldane phase, which remains robust against such driving.

I. INTRODUCTION
In order to find quantum states with desired properties, we can look in various spaces: In the chemical space we can investigate the multitude of natural compounds. This space can be further extended by synthesis, metamaterials or the replacement of chemical bonds by magneto-optical traps in ultracold quantum gases. Another possibility is to exploit the additional dimension of time and engineer new states in nonequilibrium conditions.
A particularly simple way to create a nonequilibrium state is a quantum quench, where the system is suddenly evolved with a new Hamiltonian. This can be used to observe the melting of equilibrium order parameters, such as string order [1,2], and can lead to quasi-steady prethermalized states [3] before heating sets in, but does not offer optimal control and is not easily implementable beyond ultracold-atom systems.
Alternatively, one can drive the system out of equilibrium by a periodic external force, e.g., a continuous laser beam with frequency Ω = 2π/T and period T . In practical terms, this setup is well-controlled in the high-frequency limit, where one can find the effective Floquet Hamiltonian [4,5] by means of a Magnus expansion of the original Schrödinger equation. In leading order, this results in renormalized system parameters, so that the problem can be analyzed using equilibrium techniques (Floquet engineering). Heating to an infinite-temperature state should take place eventually, but is shown to happen on exponentially long time scales for large frequencies, leading to a stable prethermalized state similar to the case of quenches [6,7].
In equilibrium physics, the concept of adiabaticity is fundamental. In practical terms, it can be used to define and traverse phase diagrams or prepare complex ground states by adiabatically changing the couplings of a Hamiltonian, e.g. using quantum annealing. Extending this concept to Floquet engineering, one can attempt to adiabatically change the drive parameters to further improve the degree of dynamic control of the system [27,28].
In this work, we adopt the specific protocol of initiating the system by driving a term with Ω = ∞, followed by an adiabatic decrease of Ω [29,30]. This adiabatically propagated state is called the "Floquet ground state" [30], and Ω is freed up as an additional control parameter in the procedure. This has been first studied for the integrable transverse-field Ising model, in which case the state was seen to undergo topological phase transitions and Kibble-Zurek scaling was observed [29,30].
Fortunately, in the case of one-dimensional chains, this "adiabatic Floquet" protocol lends itself to an efficient simulation even for non-integrable systems using matrix-product states (MPS). The initial state is guaranteed to have low entanglement for the class of gapped chains in accordance with the area law. This is a key property that is exploited by the MPS formalism [31]. Furthermore, as long as the change of frequency is slow enough, the entanglement is expected to grow only slowly and long propagation times may be reached. This stands in contrast to quench dynamics, where the entanglement entropy increases linearly with time [32], while the MPS bond dimension (i.e., the number of variational parameters to represent the state) has to increase exponentially.
In this paper, we show that the adiabatic Floquet arXiv:2210.16088v2 [cond-mat.str-el] 14 Apr 2023 protocol can be used to convert a conventional Néel antiferromagnetic state into an unconventional antiferroquadrupolar state. More specifically, we apply the protocol to the non-integrable S = 1 spin chain, a system which is mainly interesting for its symmetry-protected "topological" Haldane phase and the potential for spin-nematic (quadrupolar) order. The latter is a state with nonvanishing anisotropic second-order expectation values of the type S α j S β j = 0, while having a vanishing first-order expectation S α j = 0 (where S α=x,y,z j is a spin operator). Thus, it is an interesting quantum state that carries no magnetic moment, but still breaks the spin-rotational symmetry via a more complicated order parameter. A spin nematic can be regarded as something between a ferromagnet and a spin liquid: While it lacks magnetic order like the latter, it still breaks the rotational symmetry like the former and has a preferred axis. The name derives from the physics of nematic liquid crystals, which in a similar sense constitute a phase between a liquid and a solid [33]. While quadrupolar order is in principle possible in the S = 1 chain, in the following section we discuss that it is not easily achievable in equilibrium, motivating an extension to driven systems.
Starting from the ground state of an initial Hamiltonian with Ω = ∞, we slowly drive the system from the high-frequency to the mid-frequency region, representing the wavefunction as a MPS. We find that if the initial state is in the symmetry-protected Haldane phase, it still remains remarkably robust against the drive and no new phase transition is found. On the other hand, if the initial state is in a trivial ordered phase, then we are able to induce an overall quadrupolar moment and enhanced correlations, eventually reaching a stable phase with long-range antiferroquadrupolar order, where the staggered magnetization is suppressed.
The Hamiltonian is an extended variant of the Heisenberg spin chain: where J is the exchange interaction parameter.
The compound Ni(C2H8N2)2Ni(CN)4 (NENC) under the action of a periodic driving in the x 2 − y 2 orbital of the Ni(II) atoms. The figure is inspired by Fig. 1 in Ref. [34], but we have added a possible experimental driving setup. The structure of the chain repeats in the horizontal direction. Molecular orbitals from the Ni and N atoms have been represented in the x, y, z geometry by the lobes (see legend). Each Ni(II) is attached to the next one by a NC-Ni(CN)2-CN configuration. The magnetic properties of such Ni compounds have been studied experimentally. A theoretical description of these compounds is proposed by effective Hamiltonians in the form of Eq. (1) representing the S = 1 chain of the Ni(II). The experimental realization sketched in this figure corresponds to the driving protocol given by Eq. (12) of this work, where the drive in the xy-plane (represented as a sinusoidal wave) continuously changes the probability distribution of electronic x 2 − y 2 orbital in the Ni(II) atoms.
(S x j , S y j , S z j ) represents the spin-1 operator at the j-th site, with S α=x,y,z j representing the different spin projections. As the local basis |σ , we take the eigenbasis of S z j and denote the eigenvectors as |σ = + , 0 , − with the eigenvalues of +1, 0 and −1, respectively. Finally, h is the staggered magnetic field, while D is the anisotropy in the z-direction (easy-axis for D < 0 and easy-plane for D > 0). We set = 1 and J = 1, thereby measuring all energies in units of J and times in units of /J. For D = h = 0, there exists a gapped phase (the "Haldane phase"), which (for periodic boundary conditions) has a unique symmetry-protected ground state with exponentially decreasing spin-spin correlations. The robustness of the Haldane phase has been a focal point of previous studies in equilibrium [58][59][60][61][62][63][64], where it was found that it can be characterized by a non-local string order parameter [65][66][67] O α=x,y,z and by a global twofold degeneracy in the entanglement spectrum. It is protected by a combination of inversion symmetry, time-reversal (in the sense S x,y,z → −S x,y,z ) and combined rotations of π about a pair of axes [63].
Since a finite h breaks all these symmetries at once, even a small value destroys the Haldane phase [36]. However, it remains robust against the anisotropy term, which does not break any of the above symmetries [62][63][64].
In this case, the Haldane phase is a thermodynamic phase, which is stable in an extended region of the phase diagram, eventually losing in competition to strong-D phases (see below) once D exceeds a critical value. An interesting question is thus how this robustness extends into non-equilibrium. The other limiting cases of Eq. (1) are as follows: For D → +∞, the ground state is given by a product state of local 0 projections. In the D → −∞ limit, the ground state is two-fold degenerate, given by the Néel state . . . + − + − . . . and the Néel state shifted by one lattice site: . . . − + − + . . .. For h → ±∞, the ground state is given by a unique Néel state.
The S = 1 chain is also arguably the simplest system that allows for quadrupolar exchange. The quadrupole operator is defined as the traceless tensor It has five linearly independent components that can be grouped into a vector: so that αβ Q αβ j Q αβ j = 2 Q j · Q j . Quadrupolar exchange thus requires a product of four spin operators and is usually discussed within the bilinear-biquadratic model [68,69], given by: Because of the identity the Hamiltonian Eq. (5) boils down to a competition of ordinary exchange interaction and quadrupolar exchange. In 1D, quantum fluctuations generally prevent spontaneous ordering unless the order parameter is conserved. In this case, a quadrupolar state may rather be defined via quadrupolar correlations that dominate over spin-spin correlations. Such correlations are found with a 3-site period for J, J q > 0 and J q /J > 1 [70]. For J < 0, a ferroquadrupolar order was initially predicted close to the ferromagnetic phase [71], but highly accurate MPS calculations demonstrate that it either does not exist (with a dimerized phase found instead), or only exists in a very narrow parameter regime [70,72]. In 2D, quadrupolar phases are better defined and a finite quadrupole moment generally arises for J q / J > 1 [73,74]. We note that the choice yields the famous Affleck-Kennedy-Lieb-Tasaki (AKLT) ground state [60,61], which belongs to the Haldane phase, but is exactly representable by an MPS with very low entanglement. From a technical point of view, it is thus a more convenient representative member of the Haldane phase than Eq. (1) with h = D = 0. A possibility to generate finite quadrupolar moments in one dimension is the explicit breaking of the spin-SU(2) symmetry via the so-called "rhombic single-ion anisotropy" [75][76][77]: where we have introduced the standard spin-flip operators S ± j = S x j ± iS y j . A more general coupling to the square of the spin operator is also possible [77,78]. In the context of squeezed spin states, the term in Eq. (8) is known as the "two-axis countertwisting" (TACT) [79] and there are several proposals of how to implement it [80]. In equilibrium, a strong value of E will induce a finite quadrupole moment in the direction of Q x 2 −y 2 , similar to how an external field induces a finite magnetization [75]. However, this requires finding a material with large E and small D at the same time.
Recently, attention has shifted to a different regime, where bond-nematic rather than local order might be found. For an S = 1/2 system in a strong magnetic field close to saturation, deviations in magnetization are given by magnons. If the effective interaction between these magnons is attractive, they may form pairs and condense, with non-zero quadrupolar correlations S + i S + j playing the same role as anomalous expectation values of fermion-pair creation operators c † i c † j in a superconductor [81]. The best-documented experimental example where this may occur LiCuVO 4 [82,83]. The corresponding experiments are quite challenging and need to be performed in magnetic fields of 45-50T. Summarizing, a spin-nematic state is an exotic nonmagnetic state with a higher-level order parameter. The minimal spin value to observe it locally is S = 1. A stabilization of this state requires strong biquadratic exchange or a strong anisotropy. However, both are expected to be weak in real materials or require finding a fine-tuned point [84][85][86], so that the part of the phase diagram where spin-nematic phases are predicted could not be explored in practice. Attention has therefore shifted to the different physical regime of magnon pairing in high magnetic fields for S = 1/2 materials, which has its own challenges. Here, we pursue an alternative idea, namely the purposeful enhancement of quadrupolar interactions using nonequilibrium driving, starting from an ordinary S = 1 system (Haldane chain or Néel antiferromagnet).

B. Floquet adiabatic protocol
We evolve a system in the presence of a time-dependent term Af (t, Ω(t))V : where H 0 is the unperturbed Hamiltonian, f (t, Ω(t)) is the periodic envelope function of the drive, Ω(t) is the drive frequency, A is its amplitude and V is the operator of the drive. It is convenient to choose an envelope function that averages to zero over one cycle, Specifically, we set This means that for Ω(t = 0) = ∞, the state at t = 0 can be obtained as the ground state of H 0 and is a legitimate Floquet state.
In order to observe quadrupolar order in our setup, we choose to drive the rhombic anisotropy Eq. (8), which is quadratic in the spin operators (cf. Fig. 1): For t > 0, we start with a frequency of Ω 0 that is large enough to be connected with the infinite-frequency state at t = 0. We then adiabatically decrease the frequency in the range Ω(t) ∈ [Ω f , Ω 0 ] (Ω 0 > Ω f ) on a time window of length t f (see App. A for more details). After that, we assess whether a steady state has been reached by evolving the state with constant Ω f for another multiple of t f , i.e., until the end time t end = αt f , where α > 1 is a scaling factor. In other words, the steady state is observed on a timescale of (α − 1) t f .
The meaning of these control parameters is as follows: (1) The target time t f controls adiabaticity, whereby larger values make the process more adiabatic, so that we would like to choose t f as large as possible. (2) For the final frequency Ω f , the most interesting regime is on the energy scale of the system or below it, i.e. O(10 0 ) − O(10 −1 ) in units of J. However, if Ω becomes so small that Ω −1 ∼ t f , the system cannot be considered locally periodic at any time, and the Floquet theorem loses its validity. Therefore, one should choose a range of Ω for which all values satisfy Ωt f 2π ∀Ω ∈ [Ω f , Ω 0 ] [29], which means that smaller Ω f have to be offset by larger t f .
Taking all these trade-offs into account, we set Ω 0 = 100.0, Ω f ≥ 10.0, t f = 62.82, α = 1.5, 2.5.  (11), with Ω f = 0.3. We use an intermittent measurement process, where the red circles correspond to the measurement times of the observables. On the top axis, we have represented the corresponding values of the driving frequency Ω, which is kept constant at Ω = Ω f for t/t f > 1.

A. Time-Evolving Block Decimation (TEBD)
We use the infinite time-evolving block decimation (iTEBD) with a fourth-order Suzuki-Trotter decomposition [31,87] to compute time evolution of a given MPS directly in the thermodynamic limit. The main control parameter is the bond dimension χ, which encodes the amount of entanglement.
The MPS |ψ has a unit cell of L cell = 2 for a finite staggered field in Eq. (1). We represent it in the standard Γ − Λ notation [88], with A and B denoting the sublattices, and σ i = {+, 0, −}. The AKLT state (Sec. IV A) can be expressed analytically as a MPS with a bond dimension of χ = 2. The Néel ground state (Sec. IV B) is computed from an imaginary time evolution with a fixed bond dimension of χ = 200 and a decreasing Trotter step size until convergence has been reached; we have checked that using χ = 120 or L cell = 4, 6 does not change the results, and also cross-checked with the Variational Uniform MPS (VUMPS) algorithm [89] as an independent benchmark. During the real time evolution, it is convenient to choose a variable Trotter step size ∆t that depends on the current period. We do not let it exceed the maximum value of ∆t = 0.01. Furthermore, we set an initially fixed bond dimension of χ = 200 (for the quench starting from the AKLT state, the initial bond dimension χ = 2 is increased such that that the state is encoded exactly until χ = 200 is reached). This is sufficient for an upper bound of the discarded weight (defined as the sum of discarded weights on each Trotter substep) of = 10 −5 . In other words, fixing = 10 −5 does not lead to an increase of the bond dimension beyond χ = 200, which seems reasonable due to the adiabatic nature of the quench. Further comments on this issue as well as a comparison with exact diagonalization data can be found in App. B.
B. Measurement events and the steady-state limit When we keep track of all observables in a continuous fashion, the micromotion of the Hamiltonian dynamics leads to the appearance of oscillations as a consequence of the drive given by Eq. (11). We ignore this micromotion by an intermittent measuring at some of the zeros of the drive, as indicated in Fig. 2. This subset of N measuring points is defined by: Furthermore, we define the steady-state limit (SSL) value of an observable O in the Heisenberg picture by averaging over the last N SSL measurement times: Due to computational limitations, the above definition is strictly speaking a quasi -steady state, in the sense that the reachable times are finite, but long on the intrinsic time scale: αt f /J = 1 (α = 1.5, 2.5).

A. Driving from the Haldane phase
We apply the protocol given by Eqs. (9) and (12) to the AKLT state, which is the ground state of the Hamiltonian Eq. (7) (H 0 = H AKLT ) and an ideal representative of the Haldane phase.
Since the drive does not break the protective symmetries of the Haldane phase, we expect the phase to be robust at least against small amplitudes. Numerically, we are able to go up to A = 2.5 and find that the drive is still unable to destroy the Haldane phase because it is unable to destroy the string order. For the initial AKLT state, the string order is isotropic with a value of O α=x,y,z string = 4/9. Figure 3 shows its evolution and we find that it remains finite in all directions. (Note that this stands in contrast with a sudden quench, which is able to destroy the order [1,2].) Furthermore, we find no noticeable increase of a quadrupole polarization or correlations (not shown).
Our conclusion is that the robustness of the Haldane phase is difficult to overcome in our setup and one needs to start from an initial ground state which should already incorporate the breaking of its protective symmetries. To this end, we explore the driving protocol of Eq. (12) and apply it to the initial Hamiltonian given by Eq. (1) with non-zero h and D.
B. Driving from a trivial phase

Choice of parameters
Relatively small anisotropies and interchain couplings are enough to induce an ordered phase for the S = 1 chain. Still, in a lot of Ni-based compounds both are weak enough to keep them in the Haldane phase or at the edge of the phase transition line [37]. An exception is CsNiCl 3 , which shows Néel order below a critical temperature, and for which h = 0.051, D = −0.038 have been deduced [35,40] (though not without uncertainty [39]). A very large anisotropy of D ∼ −1.5 was recently proposed for BaMo(PO 4 ) 2 [90].
To demonstrate the general principle, we set H 0 = H Heis [Eq. (1)], D = −0.2 and h = +0.1, which puts the initial system firmly into the ordered Néel phase: The ground state has a finite staggered spin polarization of (−1) j S z j ≈ 0.74 in the z-direction. Spin-spin correlations of the z-projection show antiferromagnetic (AFM) long-range order, while quadrupolar correlations are short-ranged (see Fig. 5(b) for t = 0).

Suppression of AFM Neél order and emergence of quadrupole polarization
The time evolution of the staggered magnetization (−1) j S z j and the quadrupole component (−1) j Q xy j are presented in Fig. 4. The initial staggered magnetization of the Néel state is always reduced by the drive and is completely suppressed for A ≈ 8. This coincides with the emergence of a staggered net quadrupole moment in the Q xy -direction, which becomes finite and survives in the steady-state limit, where Ω f = 10.0 is kept constant. There is an optimal value of A where the quadrupole polarization is highest, i.e., it tends to be smaller for both very small and very large A.
We note that the polarization is obtained in the third component Q xy of the quadrupolar operator Eq. (4), while the driving term only couples to the first component Q x 2 −y 2 . This axial enhancement is somewhat analogous to the behaviour observed in Ref. 25, where switching on a rotating magnetic field in the xy-plane induces a finite polarization of the z-component of the spin, i.e., perpendicular to the plane of the driving.

The antiferroquadrupolar correlations
The above result suggests the dynamical emergence of a spin-nematic state due to application of the drive, where the quadrupole polarization dominates over the magnetization. To further investigate this state, we fix A = 6.0 and show the quadrupole correlations Q xy j Q xy j+r in Fig. 5. The appearance of an antiferroquadrupolar (AFQ) state is characterized by negative (positive) correlations at odd (even) sites r that remain stable for times t/t f > 1 in the steady-state limit. Note that there is a sign switch as a function of time for r = 1. In Fig. 5(b), we show the profile of the quadrupole correlations characteristic of the AFQ state in the steady-state limit, compared to their initial values at t = 0. We point out that our resulting state is distinct from the large-E phases reported in Ref. 75, where the polarization is in the direction of Q x 2 −y 2 , i.e. along the applied field E. If the ground state can be chosen to be real-valued, then Q xy j = 0 must necessarily hold in equilibrium. A polarization in the direction of Q xy therefore results from the time propagation during the drive, which gives an imaginary part to the wavefunction. Thus, our resulting state generally has no equilibrium analogue for systems described by Eqs. (1) and (5), even if the rhombic anisotropy term Eq. (8) is also present, unless spontaneous breaking of time-reversal symmetry results in a complex-valued ground state. We discuss this in some more detail in Sec. IV B 5.

Phase diagram in the steady-state limit
To better characterize the dynamically obtained AFQ state and explore its stability at longer times, we study the phase diagram in the steady-state limit (SSL) for the two order parameters of the staggered magnetization (−1) j S z j and the quadrupole polarization (−1) j Q xy j , as a function of the target frequencies and driving amplitudes (Ω f , A). The spin-nematic regime is reached when (−1) j S z j SSL ≈ 0 and (−1) j Q xy j SSL = 0 for an average over the steady-state regime Eq. (15). We define the crossover parameter between the two regimes as:   6. (a) Phase diagram in the steady-state limit (SSL) following a quench from the trivial Néel phase for the crossover parameter rSSL defined in Eq. (16). The three regimes are (1) magnetized: rSSL < π/4, (2) hybrid: rSSL π/4, (3) spin-nematic: rSSL ≈ π/2. (b) The SSL values of the emerging quadrupole moment for different Ω f as a function of the driving amplitude A. Reaching smaller Ω f enhances the net quadrupole moment requiring smaller amplitudes A. Strong values of A tend to suppress the overall enhancement of the quadrupole moment, cf. Fig. 4. (c) The SSL value for the staggered magnetization at different Ω f as a function of the driving amplitude A. A smaller Ω f requires smaller A to suppress the Neél order. Straight lines between data points in (b) and (c) are a guide to the eye. equal strength, which we call the hybrid region. Finally, the straight lines delineate the spin-nematic region The effect of reducing the target frequency Ω f on both the staggered magnetization and the quadrupolar moment is shown in Figs. 6 (b), (c). In the small-A region, the systematics of the AFM order and AFQ order are opposite, so that decreasing the frequency leads to weaker AFM polarization and stronger AFQ polarization. In the large-A region, the systematics is the same and decreasing the frequency leads to a weakening in both.

The AFQ wavefunction
Despite the apparent simplicity of the final xy-AFQ state, we find that it is entangled and cannot be written down as a simple analytical wavefunction. However, we can heuristically attempt to write it down as a simplified MPS by restricting ourselves to the three largest eigenvalues of the entanglement spectrum Λ AB and Λ BA (for the sublattices A and B) in Eq. (13).
Analyzing the numerical result for the ground state, we find that it approximately has the following MPS structure: where a σ ij = Γ A,σ ij denotes the non-zero matrix entries.
After switching on the drive, we keep monitoring the pattern of non-zero entries and find that the xy-AFQ state shows the following structure: Moreover, we also observe: As a technical detail, note that the reduced MPS also needs to be renormalized. Figure 7 compares Q xy s=A,B obtained with the full and with the heuristically reduced wavefunction, with overall good agreement. To derive an analytic expression for Q xy s=A,B , we note that Q xy s can be written using Q ± s = (S ± s ) 2 , S ± s = S x s ± iS y s as follows: We find: Equation (21) shows that for the initial MPS in Eq. (17), Q ± s = 0, as the matrices Γ s,± do not share common matrix elements at equal row and column indices. Conversely, the driven state given by Eq. (18) does contain non-zero products between the matrix elements and yields a net quadrupole moment Q xy j . Furthermore, we see that any contribution to Q xy s must come from the imaginary part of Q ± s . This implies that the matrix elements of Γ s,± must contain imaginary parts, which are acquired in the course of the unitary time evolution (beyond a trivial global phase). This shows that the obtained xy-AFQ state cannot be engineered from an equilibrium configuration by means of an arbitrary variation of the parameters of the model, provided that time-reversal symmetry holds (i.e., the ground-state wavefunction can be chosen real-valued).

V. SUMMARY
We have demonstrated that it is possible to convert a conventional Néel antiferromagnetic state on the S = 1 quantum spin chain into an unconventional antiferroquadrupolar state by applying the adiabatic Floquet protocol (starting with an Ω = ∞ initial state and adiabatically decreasing Ω). We chose to drive the rhombic anisotropy Eq. (12), as it contains squares of the spin operators. Up to the time scales considered, the engineered state remains stable when the final frequency is kept fixed. The range of amplitudes A and final frequencies Ω f where such an AFQ state is observed is shown in the phase diagram of Fig. 6. While we are limited to relatively large Ω f in terms of numerics, we surmise that the spin-nematic region extends to smaller Ω f (and therefore smaller A according to Fig. 6).
In contrast to this, we find that an initial state belonging to the symmetry-protected Haldane phase is robust against such driving, as evidenced by the preserved string order parameter. Thus, while a lot of the previous interest in S = 1 chains was motivated by an experimental realization of Haldane physics, we find that its robustness limits the possibilities for nonequilibrium engineering. We therefore propose that conventional Néel-ordered compounds are more useful in this regard. In particular, we see that the easy-axis anisotropy is crucial in stabilizing a net quadrupole moment.
Overall, the adiabatic Floquet protocol presents a controlled approach to driven systems in order to manipulate and engineer valuable quantum states. Due to the adiabatic changes, we can achieve much longer numerical propagation times than in the case of sudden quenches.
From a technical perspective, we have employed matrix product states in combination with the iTEBD algorithm; we have also used exact diagonalization (see App. B) to back up our claims. Due to the adiabatic nature of the quench, a bond dimension of χ = 200 is sufficient to bound the discarded weight from above. In fact, we provided analytical arguments that the drive from the Néel state can essentially be encoded with χ = 3. For small discarded weights shown in App. B, however, we observe a fast increase of the bond dimension and thus in principle an uncontrolled error on this scale. This might be attributed to non-generic features of the micromotion, though we lack a clear understanding of this issue, which is beyond the scope of this paper.
An experimental realization for the driven S = 1 spin chain considered here can be either using lasers with real materials [25,57,78] (cf. Fig. 1), or in cold-atom systems [91]. would like to thank the Publication Fund of the TU Braunschweig for financially supporting the freely accessible publication of this article.
Appendix A: Details of the time discretization To perform an adiabatic variation of the driving frequency Ω we proceed as follows: We logarithmically discretize [92] the interval [Ω f , Ω 0 ] into N cycles + 1 points {Ω j } j=0,...,N cycles , where N cycles determines the total number of cycles completed by the drive in the time window [0, t f ]; thus, t f is completely determined by the N cycles parameter.
The adiabatic limit is realized with N cycles → ∞ (equivalently t f → ∞). Numerically, we fix N cycles = 300. After a single cycle is performed with period T j from times [t j , t j+1 ] ∈ [0, t f ], we decrease the frequency by δΩ j from Ω j → Ω j − δΩ j , with the corresponding increment in the period T j → T j + δT j . For times t > t f , we keep the frequency constant and equal to Ω f , evolving the state to a final time t end = αt f .

Appendix B: Numerical tests
In Fig. 8, we show the quadrupole polarization (−1) j Q xy j following a quench from the Néel phase for different values of the discarded weight and maximum Trotter steps ∆t (the setup is analogous to Fig. 4(b) which was obtained using = 10 −5 and ∆t = 0.01). The corresponding bond dimension is shown in the inset. Once a maximum value of χ = 900 has been reached, the bond dimension is no longer increased due to computational limitations, and the time evolution is carried out using a fixed χ. We note that the error is no longer strictly controlled in this regime.
At a maximum ∆t = 0.01, the bond dimension does not increase for = 10 −5 , but it increases mildly for = 10 −6 and rapidly for = 10 −7 , where its maximum value of χ = 900 is reached quickly and the error is no longer controlled. This issue persists even if the quench is made more adiabatic (N cycles = 500). Physical quantities, such as the quadrupole polarization, however, seem converged. It is reasonable to assume that our adiabatic protocol can be modelled accurately using a small χ (in fact, χ = 3 seems sufficient, see the analytic discussion in Sec. IV B 5). A possible explanation is that the blow-up of the bond dimension at small is related to features of the micromotion.
We add supportive evidence by comparing with exact diagonalization (ED) data obtained system sizes of L ≤ 12 with periodic boundary conditions. The results are shown in Fig. 9. One can see that both approaches agree qualitatively and even quantitatively unless Ω f is small or A is large. In particular, the ED data indicates a finite steady-state value of the quadrupole polarization (or a trend in this direction as L is increased).
Finally, since time propagation using ED poses no entanglement problem, we check the robustness of the quasi-steady state for very long propagation times, using N cycles = 1000, so that t f ≈ 190, t end = 2.5t f ≈ 475 and finite systems of up to L = 10. The result is shown in Fig. 10. We observe that the quasi-steady state remains robust in this parameter range.

Appendix C:
In this section, we present general expressions for the Floquet adiabatic evolution that can be used as a starting point for methods beyond the direct time propagation used in this paper.
Given a purely periodic Hamiltonian H(t + T ) = H(t) with period T , the Floquet Hamiltonian H F is formally defined by the unitary operator over a single cycle in the time interval [t 0 , t 0 + T ]: whereT is the time-ordering operator. Note that in general, the Floquet Hamiltonian is dependent both t 0 , and T ; i.e. H F ≡ H F (t 0 , T ). In most of cases determining H F in an exact way is not possible due to the complex time ordering appearing on the right hand side. A common approach to approximate H F for a fixed frequency Ω = 2π/T is to truncate the Magnus expansion to a given order or to employ an effective Hamiltonian; however, the validity of these approaches is normally restricted to the region of high Ω values.
In the adiabatic Floquet approach described in Sec. II B, we start with the initial frequency Ω 0 = ∞, where H F is known exactly. Due to the slow decrease in Ω, H F undergoes infinitesimal changes (and so does the unitary evolution operator in Eq. (C1)). To determine the adiabatic variation of H F , we look at two neighbouring cycles with the drive period differing by an infinitesimal amount δT , and employ in their respective intervals the definition given by Eq. (C1). This way, one does not need to solve any of the Floquet Hamiltonian problems individually, as only the difference between the series expansions of Eq. (C1) needs to be taken into account. We discretize the time interval [t 0 , t f ] as a set of points t k=0,...,N =f . The values of t k will be called the switching times, representing the times where an infinitesimal variation of the driving period δT takes place. In each interval [t k−1 , t k−1 + T k ], the value of the drive frequency is fixed and given by Ω k = 2π/T k , and the drive function f (t) satisfies: We consider a time dependent Hamiltonian of the form: where H 0 is the time-independent part and X j are generic local operators. To simplify calculations, we further impose the following restriction on the drive function within each finite interval: On each time interval [t k−1 , t k−1 + T k ], we associate a Hermitian operator H (k) F defining the Floquet Hamiltonian in that interval: The last identity follows from the connection property of the evolution operator. Within a given time interval of fixed period T k , the series expansion for the unitary operator is: where the carets indicate that the operators are in the interaction representation, i.e.V k (s) = e iH0s V k (s)e −iH0s , and we define V k (s) ≡ f (s, T k )X (see Eq. (C3)) for a fixed value of the period T k . The term in brackets contains a product of the operatorsV k (t j ) whose order is irrelevant due to the time ordering operator acting over the bracket.
We consider now two adjacent unitary operators: with δT ∼ 0 an infinitesimal variation of the period, and δV k representing an infinitesimal change of the Floquet operator when a decrease in the driving frequency from Ω k to Ω k+1 takes place. The first interval of the evolution is [t k−1 , t k−1 + T k ] and the second is [t k−1 + T k , t k−1 + 2T k + δT ]. The initial form of the Floquet Hamiltonian is known in the Ω → ∞ limit: The variation with respect to the period is given by: The two neighbouring terms are explicitly written in their series expansion: By changing the variables of integration in the U k+1 F tõ t n = t n − T k , we obtain two different contributions when subtracting both series ∆U (k) where we have defined: Note that all summed up terms start at n = 1, since the identities appearing on the right hand side have cancelled out. For the k + 1 term we make use of the following expansion: We stress that the order of products inside the brackets is irrelevant due to the action of the time-ordering operator. The second term in the right-hand side of Eq. (C12) is identically zero when Eq. (C4) is satisfied over the infinitesimal interval of integration. For the first series, the term in brackets at a given order n is given by: where we have abbreviated a = t k−1 and b = t k−1 + T k , and defined: We can arrange t 1 a total of n times, with the remaining n−1 operators still subjected to time-ordering. Proceeding to all orders, Eq. (C12) becomes: