Controlling quantum many-body systems using reduced-order modelling

Quantum many-body control is among most challenging problems in quantum science, due to computational complexity of related underlying problems. We propose an efficient approach for solving a class of control problems for many-body quantum systems, where time-dependent controls are applied to a sufficiently small subsystem. The approach is based on a tensor-networks-based scheme to build a low-dimensional reduced-order model of the subsystem's non-Markovian dynamics. Simulating dynamics of such a reduced-order model, viewed as a ``digital twin"of the original subsystem, is significantly more efficient, which enables the use of gradient-based optimization toolbox in the control parameter space. We validate the proposed method by solving control problems for quantum spin chains. In particular, the approach automatically identifies sequences for exciting the quasiparticles and guiding their dynamics to recover and transmit information. Additionally, when disorder is induced and the system is in the many body localized phase, we find generalized spin-echo sequences for dynamics inversion, which show improved performance compared to standard ones. Our approach by design takes advantage of non-Markovian dynamics of a subsystem to make control protocols more efficient, and, under certain conditions can store information in the rest of the many-body system and subsequently retrieve it at a desired moment of time. We expect that our results will find direct applications in the study of many-body systems, in probing non-trivial quasiparticle properties, as well as in development control tools for quantum computing devices.


I. INTRODUCTION
The remarkable experimental capabilities have led to the advent of quantum technologies and inspired intense efforts to develop optimal control methods for quantum systems across several fields (for a review, see Ref. [1]).Finding optimal control sequences for quantum many-body systems is a particularly important, but challenging task.Generally, full simulation of a many-body system dynamics for a given choice of control parameters requires resources exponential in the number of degrees of freedom; in gradient-based optimization methods over a large parameter space, such simulation has to be repeated many times, which makes such methods prohibitively demanding.
To overcome this challenge, gradient-free optimization methods combined with tensor-networks-based dynamics simulation were applied to control problems including many-body ground state preparation [2][3][4][5].Recently, methods using reinforcement learning techniques have been proposed [6,7]; such approaches can also be viewed as gradient-free optimization since they do not use the gradient of the reward function.The gradient-free optimization based control methods, however, are generally expected to be less efficient than gradient based methods [8].Recently, Ref. [9] demonstrated the advantage of gradient-based methods for the problem of ground state preparation by crossing the superfluid to Mott-insulator phase transition in the Bose-Hubbard model.
Here, we propose a different approach to a class of problems where a sufficiently small subsystem S of a many-body system is subject to time-dependent controls.Focusing on time evolution of degrees of freedom in S, we represent the rest of the many-body system E by its lower-dimensional "twin", or a reduced-order [10][11][12][13] model.Such a reduction effectively keeps track of the relevant degrees of freedom in E, discarding the ones which have little or no influence on dynamics of S. The reduced-order model may involve effective Hilbert space dimension that is orders of magnitude smaller than the one in the original problem.This allows us to use the powerful toolbox of gradient-based optimization methods, using automatic differentiation techniques [14] to calculate the gradient of the loss function.
Practically, to build a reduced-order model, we employ tensor-network techniques, that are widely used for dimensionality reduction in quantum many-body physics [15,16] and applied mathematics [17,18].Unifying the developed reduced-order modeling scheme with gradient based optimization yields an efficient method for quantum many-body control.
We use our approach to automatically design protocols for manipulating information propagation in strongly interacting systems.First, we consider an one-dimensional (1D) XYZ quantum Heisenberg chain with extra fields that break integrability.For simplicity, we choose a single spin as the subsystem where time-dependent controls are applied.Our algorithm is able to find sequences that restore quantum information locally, or transmit to another end of the chain.Physically, the identified pulse sequences of local operations inject and re-absorb long-lived quasiparticles in an optimized way.
Further, we apply the optimization method to the many-body localized phase.We are able to find control protocols for local dynamics inversion that outperform existing spin-echo-type protocols for many-body localized systems [19].Thus, our approach enables automated discovery of optimal generalized spin-echo sequences in interacting systems.
The method describe here can be readily applied in experiments with the current generation of noisy, intermediatescale quantum (NISQ) devices [20][21][22].Various quantum computing platforms, including programmable Rydberg simulators, trapped ions, isolated spin impurities in solids, and superconducting circuits arrays realize 1D spin chains [23][24][25][26][27][28][29][30] with the possibility to control qubits individually by means of optical or microwave pulses.Our approach shows that the non-Markovianity of many-body environment can be employed for generating excitations and information spreading across the system.We expect that a modification of our approach may also be used for many-body state preparation.
Here we are focused on the the realization of the control method in the coherent phase of quantum many-body systems assuming specific techniques to avoid fast thermalization.In the thermalized phase, our approach is not applicable and, at the same time, there are no reasons to expect here a possibility to maintain controllable coherent dynamics required in applications.
Our work is organized as follows.In Sec.II, we describe our general approach to building a reduced-order model.We illustrate this methods for a quantum spin chain, our primary quantum many-body model of interest, in Sec.III.In Subsec.IV A, we discuss the way to designing control protocols on the basis of reduced-order models.We illustrate several control protocols using the proposed method: in Subsec.IV B, we demonstrate the inversion of dynamics of the system in time; in Subsec.IV C, we discuss controllable information propagation across the system, in particular we show a possibility to realize information transmission from Alice to Bob locating at different ends of the chain.Finally, we conclude and discuss potential next steps in developing the proposed method in Sec.V.

II. BUILDING A REDUCED-ORDER MODEL
Consider a many-body quantum system that consists of two parts.Assume, that the dynamics of the first part is of our interest, while the dynamics of the second part is not.This induces a natural separation of a many-body system into the target system (the first part) and the environment (the second part) in spirit of the theory of open quantum systems.Our first goal is to build a low-dimensional effective model of the target system dynamics, whose numerical simulation is much faster in comparison with the original model.This is the key step allowing one to run thousands of optimization iterations at the control signal adjustment stage.Since dynamics of the environment is not of out interest, one can reduce its dimension in such a way that dynamics of the target system remains almost the same.We utilize tensor network techniques for this purpose.We split the entire system dynamics into a sequence of unitary transformations U and represent it as a tensor network shown in Fig 1a .Next, we use a singular-value decomposition (SVD) to decompose each U that is seen as a tensor with four entries into two subtensors ending up with a decomposed tensor network shown in Fig 1b .As one can note, the entire tensor network breakes up into two subnetworks, namely an environment network that is highlighted by the blue color and a system network that is highlighted by the red color.Note, that both networks have the form almost identical to a matrix product state (MPS) with only difference in the last dangling edge.The environment network by construction describes all the environmental effects in the system's dynamics.By reducing its dimension we end up with an effective low-dimensional environment network that leads to almost the same dynamics of the target system.We utilize the standard MPS truncation technique to build the effective low-dimensional environment network, this procedure schematically represented in Fig 1c .We fix a desirable accuracy of the truncation and the truncation algorithm provides an effective environment network with bond dimension r(k) depending on discrete time k.By taking the convolution between the effective environment network and the system network, as it is sketched in It is worth to notice, that the environment network is closely connected to an influence functional [31,32] that has been recently applied to a numerical simulation of quantum dynamics in variety of contexts.The most widely spread use case of the influence functional is the numerical simulation of a dynamics of an open quantum system coupled with a harmonic bath.Early approach to this problem [33,34] cuts off long-time memory effects by removing multipliers from the exact analytical form of the influence functional of a harmonic bath making it tractable for FIG.1: a Representation of a coupled system and environment dynamics as a tensor network.b The same tensor network after decomposition of each unitary block into two blocks.The entire network brakes up into two sub-networks namely a system network and an environment network shaded by red and blue regions correspondingly.c Schematic representation of the truncation procedure of the environment network.d The truncated tensor network that approximates the exact system and environment coupled dynamics.e Connection between the environment network and the influence functional.As one can see, the influence functional is a convolution of the environment network with complex conjugated version of itself.numerical treatment.More recent approaches such as TEMPO [35] and its variations and improvements [36,37] use low-rank tensor-network representation of the influence functional to improve accuracy and include long-time memory effects in the consideration.Another recent approach [38] aimed on replacement a complex harmonic bath by a simple one allowing numerical simulation of the system's dynamics.The core idea of the approach is to find such a surrogate bath, which has the same two-time correlation functions.In case of a Gaussian bath this is equivalent to the equality of influence functionals.The combination of tensor network techniques with the theory of the influence functional were recently applied to develop new analytical and numerical methods for correlated spin systems dynamics analysis and simulation.In particular, self-consistency equation for the influence functional has been introduced in [39][40][41] which allows one to study long-time thermalized dynamics of spin systems both analytically and numerically.It has been shown, that influence functional admits exact disorder averaging [42] making it possible to study many-body localization (MBL) phases rigorously.
Approaches based on the influence functional are especially successful in case of irreversible processes with weak memory effects.In these cases temporal entanglement is weak and a low-rank MPS can efficiently approximate an influence functional.Whereas such systems are fundamentally very interesting, they are not well controllable due to the high information loss rate.Systems with long memory effects and weak information loss are better controllable and more interesting from the optimal control perspective.For describing environmental effects in such kind of problems the environment network suits better.To justify this claim let us consider the connection of the environment network with the influence functional that is represented in Fig. 1e.As one can see, the environment network with bond dimension r corresponds to the influence functional with bond dimension r 2 .This means that the environment network can describe higher temporal entanglement and more complex memory effects.But on the other hand, the environment network can not efficiently describe irreversible loss of information.Any loss comes at the cost of increasing bond dimension.Therefore, environment network suits well for describing processes that are well controllable, i.e. with small information loss and high temporal entanglement and complements the influence functional based approaches.This motivates our choice of the environment network for the optimal control purposes.

III. REDUCED-ORDER MODELING OF A QUANTUM SPIN CHAIN
We begin validation of the proposed reduced-order modeling technique from building a reduced-order model of a discretized in time XYZ quantum Heisenberg chain with external magnetic field and with open boundary conditions.The dynamics of a spin chain within this model is described by two-spin unitary operators that read where # ∈ {left, mid, right} and τ is a time step size.The corresponding two-spin Hamiltonians take the following form where σ k i denotes a k-th Pauli matrix (k ∈ {x, y, z}) acting on a spin number i, J k is a coupling constant, h k is a component of the external magnetic field.In what follows we use tensor network diagrams in order to represent the final state of the spin chain after unitary evolution.One can think of U # as a four-way tensor that is represented graphically as a block with four edges.Combining blocks U # and initial spin chain's state into a tensor network, one represents the state of the spin chain at time T = N τ as it is shown in Fig. 2a, where n is a number of spins, N is a discrete time or a number of unitary layers, |ψ(0) = n i=1 |ψ i is an initial state of the spin chain.In the limit τ → 0, τ N → T one restores the standard continuous in time XYZ quantum Heisenberg model in an external magnetic field.
As a target system we choose either the first or the middle spin of the spin chain and the rest of the spin chain we treat as the environment.Throughout the paper we use l to denote the number of the target spin.For large spin chains the reduced-order modeling technique presented in Sec.II is not directly applicable, since it requires explicit manipulations with the exact environment, whose dimension grows exponentially with number of subsystems.However, one can take advantage of the environment structure, i.e. in the given case it is either a single chain connected to the first spin or two chains connected to the middle spin.If environment has a chain-like structure, one can build the effective environment by iterative adding more subsystems to it and truncating it when necessary.Such an iterative approach does not require explicit manipulations with the complete environment and therefore scalable and applicable to large chain-like environments.
We perform such an iterative model order reduction procedure and end up with an effective low-dimensional model describing dynamics of the spin of interest.The transition to the effective model describing dynamics of the first spin is represented in terms of tensor diagrams in Fig. 2b, where by the red color we highlight the first spin whose dynamics is of our interest.Similar illustration can be build for the case of the middle spin.The dynamics of the target spin withing the reduced-order model reads where l is the number of the target spin (first/middle), |ψ eff (k) is the joint target spin and effective environment state at discrete time k and l (k) is the state of the target spin.Note, that in general case the dimension d eff k = 2r(k) of |ψ eff (k) increases with time.An explanation of this effect is that the target spin gets entangled with more spins of environment with time, and therefore one needs to include more degrees of freedom of the environment into consideration.
We compare the exact dynamics of the target spin with the dynamics simulated by use of the reduced-order model.We choose the following parameters of the model J x = 0.9, J y = 1, J z = 1.1, h x = 0.2, h y = 0.2, h z = 0.2, τ = 0.15, and the following initial states of the target spin and the environment where the number of spins is n = 27, and the chosen model parameters correspond to a weakly non-integrable dynamics in the continuous in time case.We built the effective environment network by iterative adding spins to it, and truncating it each time when its dimension exceeds r max = 512.We set the truncation accuracy (see Appendix A and B for more details) to be = 0.01.The comparison of the exact target spin dynamics simulation with the dynamics simulation based on the reduced-order model is given in Fig. 3a, b.One can note that the dynamics of the reduced-order model matches perfectly the exact one.
In order to demonstrate that the reduced-order model is capable of the external control signal response prediction, we apply a random control signal (a sequence of random one-qubit unitary gates) to the target system and compare its exact dynamics with the reduced-order model based one.The comparison is given in Fig. 3c, d.As before, one can see that the exact and the reduced-order model based dynamics match each other.
In order to study how the dimension of the reduced-order model, i.e. 2r(k), scales with time for various numbers of spins and demonstrate the effect of dimensionality reduction we plot it in Fig. 3e against time for different n and the target spin fixed to be the first spin.
We also use the exact simulation of the whole spin chain dynamics in order to estimate the number of spins that are covered by the light cone propagating from the first spin.This number shows how many spins are involved in the dynamics of the first spin and thus the dimension of the Hilbert space of those spins is an upper bound of the reduced-order model's dimension.We plot how this upper bound evolves with time in Fig. 3e.One can see that the dimension of the reduced-order model grows much slower compared to the light cone based estimation of the effective dimension.While the reduced-order model's dimension for 27 spins reaches ≈ 10 4 at the final time step, the light cone covers the entire spin chains, which means that all 26 spins of the environment are involved in the first spin dynamics.This is an evidence of the proposed reduced-order modeling technique efficiency because the naive light cone based estimation of the effective dimension results in d eff N ≈ 2 27 ≈ 1.3 * 10 8 while the reduced-order modeling technique results in d eff N ≈ 10 4 .However, the model reduction soonly becomes intractable with increasing simulation time, since r(k) grows exponentially, and one can not simulate thermalization of the target spin properly.Nevertheless, our goal is to build the reduced-order model suitable for further simplification of different control problems and for this purpose the proposed technique suits well.
In the following section we move forward and apply the developed reduced-order modeling scheme to various optimal control problems.

IV. MANY-BODY OPTIMAL CONTROL: METHODOLOGY AND NUMERICAL EXPERIMENTS A. Reduced-order modeling based optimal control
The developed reduced-order modeling technique gives rise to a new class of optimal control methods in quantum many-body physics.Indeed, the main difficulty towards efficient quantum many-body optimal control is the necessity of running dynamics simulation thousands of times.This difficulty is substantially mitigated via the reduced-order modeling.The overall optimal control scheme breaks up into two steps: 1.One builds a reduced-order model of a quantum many-body system.Now dynamics simulation can be run thousands of times within a reasonable time; 2. One formulates an optimal control problem as an optimization problem written in terms of the reduced-order model and resolves this optimization problem using some optimization method.
It remains unclear what kind of optimization method to use.A typical control problem written in terms of the optimization problem takes the following form: where L is the loss function that is written in terms of the reduced-order model and it measures how good a control signal is (the smaller value of L is, the better control signal is), {u i } ∆N i=1 is a ∆N -steps sequence of unitary control gates applied to the system, i.e. it is a control signal that needs to be optimized, 1 is the identity operator.Note, that Eq. ( 5) is the constrained optimization problem.Since control gates are unitary, an optimization technique of our choice must preserve u † i u i = 1 for all gates.To solve the given optimization problem we found a Riemannian optimization algorithm [43,44] namely Riemannian ADAM optimizer [45,46] to be efficient.It performs a gradient-based search of the optimal point on a manifold defined by the constraints, in our case on the manifold of unitary matrices (a special case of the complex Stiefel manifold [47][48][49][50]).We calculate gradient of L w.r.t {u i } ∆N i=1 utilizing the automatic differentiation technique [14].Riemannian ADAM optimizer performs descent procedure towards the optimal point on the manifold of unitary matrices until convergence evaluating the gradient of L typically ten thousands times.Note, that without the use of the reduced-order model, even a single calculation of the gradient becomes extremely memory demanding since automatic differentiation requires to keep all intermediate data in memory.
Let us consider a simple example.Suppose we are allowed to control the first spin of a spin chain.The transition to the reduced-order model in this case is schematically represented in Fig. 2d.The exact model and the reduced-order one are interchangeable in terms of the control response prediction.Suppose the optimal control goal is to have the initial and the final state of the first spin the same for any initial state.In this case the loss function can be written as F , where Id is the identity channel, Φ(N |u 1 , . . ., u N ) is the channel that maps the initial first spin state to the final state.Here Φ(N |u 1 , . . ., u N ) can be written in terms of the reduced-order model as it is shown in a tensor diagram Fig. 2c and thus the loss function is cheap to evaluate.Below we provide three concrete examples of protocols for quantum control based on our approach.

B. Dynamics inversion via optimal control
Here we consider the first quantum many-body optimal control problem.Suppose that one has access to a disordered spin system in the many-body localized (MBL) phase [51][52][53][54] and it is allowed to apply a control signal to a dedicated single spin.One needs to design such a control protocol that runs dynamics of this spin "backward" in time.A typical example of such a protocol is a spin echo protocol [19,55] that runs dynamics "backward" in time in a sense that the controlled spin recovers information lost in the rest of the spin system after the spin flip operation.Our goal is to design an alternative control protocol that leads to better information recovery.
We start from a brief introduction into the origins of the spin echo protocol in MBL systems.Following the works [55,56] this effect can be explained by use of the phenomenological model of the MBL phase.The MBL phase in the thermodynamic limit can be characterized by an infinite number of local integrals of motion, which can be thought of as effective spin-half operators τ z i .In this terms the MBL Hamiltonian takes the following form where couplings J ij , J ijk , . . .fall off exponentially with separation with a characteristic localization length ξ.All terms of this Hamiltonian commute with each other.Therefor the total evolution operator U MBL (t) = exp (−itH MBL ) factorizes into the product of commuting exponents of individual terms.Consider one of those exponents that includes the participation of the first spin, it takes the following form: Taking into account the following relation: which follows from the Pauli algebra, one finally ends up with where U MBL/1 denotes the MBL evolution operator that describes the MBL dynamics of all spins but first spin and acts trivially (as the identity operator) on the first spin.It means that if the MBL system evolves for some time t then one applies the spin-flip control gate τ x to the first spin, then the system evolves for the same time t again, and finally one applies the spin-flip control gate τ x to the first spin again, one ends up with the completely the same state of the first spin as its initial state, i.e. the recovery of the information about the initial state of the first spin takes place.This is the essence of the spin echo protocol.The same consideration is valid for an arbitrary spin from the system.However, the phenomenological model Eq. ( 6) works well for the deeply localized phase.For the weakly localized phase, spin echo may barely be observed.Nevertheless, one can use the reduced-order modeling based optimal control technique from Subsection IV A to design an alternative multistep spin echo protocol suitable for a weakly localized phase.The multistep spin echo protocol consists in application of a sequence of unitary gates {u 1 , . . ., u ∆N } instead of a single σ x gate, where ∆N is the duration of the protocol (number of control gates), to the target spin at the middle of the dynamics observation.We designed this protocol for one of the models experiencing MBL dynamics.The dynamics of this model is driven by the following Floquet operator [57] where per-spin magnetic fields are random and sampled from the uniform distribution Uniform(0, 2π).The state of the whole system at discrete time k reads |ψ(k) = F k |ψ(0) .It is known that this model is in the localized phase for J < J * ≈ 0.4 [57].In our numerical experiments we consider system consisting of n = 21 spins, with coupling J = 0.3 corresponding to the localized phase, the target spin being the middle/first spin and compare the spin echo based dynamics inversion with a multistep spin echo based dynamics inversion designed by the proposed technique.
We set a particular quenched disorder, i.e we picked a particular configuration of external magnetic fields from the distribution Uniform(0, 2π).We slightly generalized the spin echo protocol in order to make it better suitable for a particular quenched disorder.Instead of the instant swap of the spin by σ x at the middle of observation we apply an instant unitary gate u that is optimized to achieve the best performance by using the proposed method.In other words, the generalized version of the spin echo protocol is the multistep spin echo protocol of duration ∆N = 1.As the initial state of the environment (all spins but the target spin) we take |ψ E (0) = n−1 i=1 |↓ .For the total number of discrete time steps N = 151 we built a reduced-order model describing the dynamics of the target spin.For the multistep spin echo protocol we turn control "ON" in the time interval from k start = 50 to k stop = 101, i.e. the total protocol duration is ∆N = 51 discrete time steps.For the spin echo protocol we turn control ON only for the single discrete time moment k = 76.
To adjust the control signal for getting the best echo effect at the end of the dynamics, one needs to formulate the control problem as the optimization problem.The initial and the final states of the target spin are connected via the quantum channel Φ(N |u 1 , . . ., u ∆N ) that can be defined by means of the reduced-order model.The closer Φ(N |u 1 , . . ., u ∆N ) to some unitary channel is, the better echo effect one has.The mutual information I(N |u 1 , . . ., u ∆N ) between subsystem of the corresponding Choi matrix Ω(N |u 1 , . . ., u ∆N ) reaches its maximum when Φ(N |u 1 , . . ., u ∆N ) is a unitary channel (see Appendix C for more details).Thus, maximizing I(N |u 1 , . . ., u ∆N ) one maximizes the echo effect.Therefore, the solution of the following optimization problem provides the optimal control signal maximize {u1,...,u ∆N } This optimization problem is solved by using the technique from Subsection IV A.
After getting the optimal control sequence, for both protocols we also run exact simulation of the entire spin chain in order to study the information flow under control and compare protocols with each other and with the case of control absence.Using the results of the exact simulation, we visualized information flow showing how the information about the initial state of the target spin spreads across the spin chain.For this purpose we utilize mutual information I l→m (k) introduced in Appendix C that shows how much information about the initial state of l-th (target) spin is kept in m-th spin at discrete time moment k.We also separately plotted I l→l (k) in order to demonstrate information revivals of the target spin.The results are given in Fig. 4. The three main conclusions could be made out of the Fig. 4. First of all, we observe the information revival at the end of the evolution for both spin echo and multi step spin echo protocols.This means, that information about the initial target spin state is being reconstructed at the end of the evolution.Second, by looking on the density plots we note, that at the second half of the evolution information about the initial target spin state starts to propagate backward towards the target spin for both control protocols.This means that the dynamics inversion takes place.Finally, one can see that the multistep spin echo protocol outperforms the spin echo protocol in terms of revival amplitude.Therefore, the multistep spin echo protocol works better than the standard spin echo protocol.
To check that the conclusions above are still valid after averaging over disorder, we performed averaging over ten different disorder realizations for J = 0.2 and all else parameters being the same.The same plots but for averaged quantities are show in Fig. 5.One sees that all the features we observed for a particular disorder realization are also valid in average.

C. Controllable information propagation in a quantum spin chain
In this subsection, we apply the proposed control technique to the control of information propagation in the discretized in time XYZ model discussed in Sec.III.We pick all the same parameters of the model as in Sec.III with the number of spins ranging from 13 to 27 and consider two control tasks aimed on controllable propagation of information through the spin chain.Within the first task we chose a bit artificial but complicated control problem causing non-trivial information flow under optimal control.The problem is formulated as follows: one needs to find such a control sequence {u 1 , . . ., u N } applied to the target spin that l (0) = l (N ) and l (N/2) = I 2 , where l is the density matrix of the target spin.In other words, we want the information about the initial state of the target spin to be completely absorbed by the environment at time T /2 and completely reconstructed back at the end of the dynamics.This control problem has the following formulation in terms of optimization: where Φ (k|u k , . . ., u 1 ) is a quantum channel that maps the initial state of a target spin to the state of the first spin at discrete time k, Id is the identity quantum channel, ∆ is a quantum channel that maps any state to the completely mixed state 1  2 I.We resolved this optimization problem using the technique from Subsection IV A. As before, we also did the exact dynamics simulation under the optimal control and without control in order to study how the information about the initial state of the target spin propagates in the spin chain.The information flow in all cases is visualized in Fig. 6.One can see that the optimal control sequence achieves the desired information flow, i.e., at FIG. 6: Density plots representing dynamics of I l→m (k) with turned OFF/ON control (x-axis is time, y-axis is the number of a spin) for different values of n and l being equal to 1 (the first spin) and n−1 2 + 1 (the middle spin).
the intermediate time, information about the initial state of the target spin dissolves in the rest of the spin chain; however, at the end of the dynamics, it is concentrated back in the target spin.Interestingly, the optimal control sequence uses reflection of the information flow from borders of the spin chain as a resource when it is possible, i.e. when the information flow has enough time to reflect from a border and get back.Note, that this is a non-trivial effect of many body echo, that is recognized and by the optimization algorithm with only use of the reduced-order model.
Within the second task, we apply the proposed method to design a control protocol allowing one to transfer quantum information through a spin chain.Let us assume that Bob prepares one of the spins in some state.The goal of Alice, who has access to one of another spins, is to apply such a control sequence {u N , . . ., u 1 } to her spin, that after time T Alice has her spin in the state as close as it is possible to the initial state of the Bob's spin.In other words, Alice has to "catch" information propagating from the Bob's spin and reconstruct the state of Bob's spin from this information.To formulate this task as an optimization problem, let us fix four linearly independent initial quantum states of the Bob's spin {|φ i φ i |} 3 i=0 whose corresponding Bloch vectors lie at the vertices of the tetrahedron, i.e.
k are components of vectors s (i) 3 i=0 that read Being able to pass these four states through the spin chain from Bob to Alice is enough to pass an arbitrary single spin state.For the fixed initial state of the Alice's spin (in our case |↓ ) one can formulate the problem of transferring states {|φ i φ i |} 3 i=0 through the spin chain as the following optimization problem minimize where (u N , . . ., u 1 , |φ i ) is the final state of the Alice's spin given the initial state of the Bob's spin and the control sequence.For each initial state |φ i of the Bob's spin we build a separate reduced-order model describing dynamics of the Alice's spin and utilize it to compute (u N , . . ., u 1 , |φ i ).The optimization problem Eq. ( 14) as previous ones is solved by using the technique from Subsection IV A. In order to address the performance of the obtained optimal control sequence {u N , . . ., u 1 } we compare initial states of Bob's spin with final states of the Alice's spins and study how the information about the initial state of the Bob's spin propagates through the spin chain.The results are given in Fig. 7.One can see, that the optimal control sequence applied to the Alice's spin is able to partly reconstruct the initial state of the Bob's spin.One can also observe how Alice "catches" the light cone that propagates from the Bob's spin and preserves the information about Bob's spin up to the end of dynamics by using the optimal control sequence.
V. DISCUSSION AND OUTLOOK In the present paper, we have proposed a new method for many-body quantum control that is based on the reduced-order modeling scheme accelerating a numerical simulation of many-body quantum systems in many orders of magnitude.This acceleration makes it possible to run tens of thousands iterations of a gradient based control signal search in reasonable total time.We have validated the proposed method on number of control problems including controllable information spreading across a spin chain and dynamics inversion in the MBL phase.
The proposed method gives rise to a new class of many-body control methods that have not been investigated before.Their field of applications varies from the development of new methods of error mitigation and noise suppression in quantum technologies to automatic search for new quantum materials, phases of matter and collective effects in many-body physics.
The proposed method can be generalized in various ways.For instance, instead of the iterative scheme for building the effective environment proposed in the paper, one can use tensor networks renormalization techniques such as introduced in Refs.[58][59][60][61][62].They are not restricted by chain like environments and one can try to build reduced order models for 2D or even 3D many-body quantum systems and systems with irregular topology that are common in the field of quantum chemistry.Another possible generalization lies in the transition from the control of local observables and partial density matrices to macroscopic observables, e.g. total energy, total polarization, etc.Indeed, together with the environment dimensionality reduction one can "renormalize" macroscopic observables leading to reduced-order models of macroscopic observables dynamics.This opens new possibilities for steering quantum manybody systems between different phases of matter via external control.The transition from "local" to "macroscopic" is possible not only for observables but also for control signals.For instance, instead of applying a control signal to a single spin one may want to apply the same time-dependent magnetic field to all spins.In this case, design of the reduced-order model is definitely more involved, but with the great development of the tensor networks toolbox it may be possible.The next interesting generalization consists in extraction of a reduced-order model from observed experimental data.It is often the case that one has access to an experimental setup with possibility to measure the response of a quantum system of interest on an external control signal.The question is whether is it possible to build the reduced-order model of a system of interest in this case based purely on observed data?With use of algorithms such as tensor-train cross approximation [17] one can try to do that efficiently and adaptively.Finally, the presented approach can be improved by unifying it with the influence matrix approach [39,40,63] allowing one to simulate long-time subsystems dynamics.The great development of tensor networks and dimensionality reduction techniques makes it possible to unify all the further generalizations of the proposed method into a universal framework opening great possibilities for automatic discovery of new quantum devices, phases of matter and quantum collective phenomena.
VI. ACKNOWLEDGMENTS I.A.L., M.A.G., and A.K.F.acknowledge the support by the RSF Grant No. 19-71-10092 (studies of the many-body control approach) and the Priority 2030 program at the National University of Science and Technology "MISIS" under the project K1-2022-027 (applications to spin chains).

VII. CODE AVAILABILITY
The code for all the numerical experiments is available via the link https://github.com/LuchnikovI/Quantum-manybody-dynamics-reduced-order-modeling.This tensor network can be splitted into two parts, a system network |S i1...i N and an environment network |E i1...i N , they read The diagramatic representations of both networks are given in Fig. 9a, b.The final joint system and environment state in terms of networks |S i1...i N and |E i1...i N can be written as follows Note, that both networks have a form very similar to the matrix product state (MPS) tensor network [15,64].The only difference with MPS is the additional dangling edge which is responsible for the final system (environment) state.
At this point we are ready to perform environment dimensionality reduction.The object whose dimensionality is being reduced is the environment network.Due to the relation Eq. (A5) the environment network is automatically in a so called left-canonical form, that is the starting point of the standard MPS truncation algorithm [17,18,65].This gives rise to an efficient environment network truncation technique, that is equivalent to the standard MPS truncation algorithm.It is easy to formulate this technique purely in terms of the environment dynamics induced by the quantum channel Ψ introduced in Eq. (A7).Consider discrete in time dynamics of the environment under the action of the quantum channel Ψ, i.e.
The central object we care about is the spectrum (λ 1 (m), λ 2 (m), . . ., λ dE (m)) of E (m), where eigenvalues arranged in the non-ascending order.If the spectrum of E (m) is mostly concentrated in r(m) largest eigenvalues, then one can project E (m) on a principle subspace that is the span of r(m) dominant eigenvectors, i.e. eigenvectors with largest eigenvalues.This leads to the truncated version of the density matrix that reads where π(m) is the orthogonal projector on the principal subspace.To gain more intuition about how the principal subspace is determined, we schematically plotted a typical spectrum in Fig. 10a.The relative error of the projection (truncation) reads where F stands for the Frobenius norm.In words it means that the error is equal to the square root of "mass" of eigenvalues in the spectrum tail that is cut and colored by red in Fig. 10a.By establishing a desirable error threshold dE j=r(m)+1 λ j (m) ≤ m one determines the principal subspace dimension r(m) and the principle subspace itself as the column space of the matrix ω(m) whose columns are r(m) dominant eigenvectors.The principle subspace is seen as a low-dimensional effective environment Hilbert space at time m.Typically E (m) gets nosier in time, i.e. its spectrum gets wider.In order to preserve the truncated spectrum tail "mass" the same, one needs to increase r(m) with time.Therefore, r(m) typically grows in time.This is schematically illustrated in Fig. 10b where one can see how r(m) grows with m due to the widening of the spectrum.To obtain the truncated environment network it is enough to insert projection operators π(m) = ω(m)ω † (m) in between of neighboring blocks B im+1 and B im for all m as it is shown in Fig. 10c.This results in truncated blocks that read  Finally, having the truncated environment network, one can build the reduced-order model of the system dynamics.The reduced-order model is defined by "effective" gates U eff m of size d S r(m) × d S r(m − 1) driving joint dynamics of channel Φ l→m (k) that maps the initial state of l-th spin to the state of m-th spin at time k.Its diagrammatic representation is given in Fig. 11a.This quantum channel fully characterizes correlations between the initial state of l-th spin and the state of m-th spin at k-th discrete time moment.To quantify correlations by a single value one can turn to the corresponding Choi matrix Ω l→m (k) that is represented in terms of tensor diagrams in Fig. 11b.This Choi matrix is seen as the density matrix of two-component quantum system and thus the mutual information between those components I l→m (k) is well defined and reads     l→m (k) are represented in terms of tensor diagrams in Fig. 11c.I l→m (k) suits well for our visualization purposes, it shows how information about l-th spin propagates in discrete time k and space m.Indeed, there are other quantities that may suit better for this role, e.g.quantum capacity [66][67][68], but we chose mutual information since it is easy to calculate.
Fig 1d, we get the effective low-dimensional model of the target system dynamics of time-depended dimension d eff k = dr(k), where d is the target system dimension.The natural formulation of the truncation technique in terms of discrete in time open quantum dynamics is available in Appendix A. Formal algorithm with its justification is given in Appendix B.

FIG. 2 :
FIG. 2: a Tensor diagram representing the spin chain state after unitary evolution.Total number of spins is n, N layers corresponds to the total evolution time T = τ N .b Tensor diagram illustrating the transition to the reduced-order model describing dynamics of the first spin.Varying edges thickness represents an increase of the effective environment dimension with time.The dimension of the joint system and effective environment state is depicted as d eff k .c Tensor diagram representing the quantum channel driving the dynamics of the target spin in terms of the reduced-order model.Asterisk symbol depicts complex conjugate.d Tensor diagram showing how one applies control signal {u1, . . ., uN } to the target spin (in the given case target spin is the first spin).The reduced-order model and the exact model can be used interchangeable, because they lead to the same response of the target spin to a control sequence.

FIG. 3 :
FIG. 3: a Comparison of the exact and the reduced-order model based dynamics of the first spin in the discrete in time XYZ model with external field consisting of 27 spins.The total number of discrete time steps is N = 50 that corresponds to total simulation time T = 7.5.σi = Tr(σi ), where i ∈ {x, y, z}.b All else being equal for the 14-th (middle) spin as the target system.c All else being equal for a random control signal.d All else being equal for a random control signal and 14-th (middle) spin as the target system.e Comparison of the reduced-order model dimension with the light cone based estimation (black line) of the effective dimension for different numbers of spins n. = 0.01 is an accuracy threshold used to build a reduced-order model (see Appendix B for the formal definition of ).

FIG. 4 :
FIG. 4: Information dynamics visualization for J = 0.3 and a first spin being the target spin and the multistep spin echo protocol, b first spin being the target spin and the spin echo protocol, c middle spin being the target spin and the multistep spin echo protocol, d middle spin being the target spin and the spin echo protocol.Density plots represent how the rescaled mutual information log(I l→m (k) + 10 −2 ) propagates across the spin chain, i.e. they show the flow of information about the initial state of the target spin.The blue band represents a region in time domain when the optimal control is applied.Line plots show how I l→l evolves in time for the turned ON/OFF control protocol.The rise of I l→l corresponds to the revival of information at the end of the evolution.

FIG. 5 :
FIG. 5: Averaged information dynamics visualization for J = 0.2 and a first spin being the target spin and the multistep spin echo protocol, b first spin being the target spin and the spin echo protocol, c middle spin being the target spin and the multistep spin echo protocol, d middle spin being the target spin and the spin echo protocol.Density plots represent how the rescaled averaged mutual information log( I l→m (k) + 10 −2 ) propagates across the spin chain, i.e. they show the flow of information about the initial state of the target spin.The blue band represents a region in time domain when the optimal control is applied.Line plots show how I l→l evolves in time for the turned ON/OFF control protocol.The rise of I l→l corresponds to the revival of information at the end of the evolution.

FIG. 7 :
FIG. 7: Density plots represent how information about Bob's spin propagates through the spin chain.Points in the Bloch ball show input states (Bob's spin initial states) and output states (Alice's spin final state).Plots correspond to the following parameters of the problem: a) n = 21, Bob's spin being the third spin Alice's spin being the spin number 19; b) n = 21, Bob's spin being the first spin Alice's spin being the last c) n = 15, Bob's spin being the first spin Alice's spin being the last spin; spin;

FIG. 9 :
FIG. 9: a Tensor network representing the joint unitary dynamics of the system and the environment.b The same tensor network with decomposed U tensors according to Fig. 8. c System network.d Environment network.
B i |ψ E , otherwise, (A13)and in the truncated environment network |E trunc i1...i N error introduced by the whole procedure is bounded above as follows (see Appendix B and Ref.[18])|E trunc i1...i N − |E i1...i N F |E i1...i N F≤ one require the error to be less or equal to some upper bound it is enough to set m = √ N which leads to|E trunc i1...i N − |E i1...i N F |E i1...i N F ≤ .(A15)Varying the value of one achieves a trade off between accuracy of approximation and effective environment dimension r(m).

FIG. 10 :
FIG.10: a A typical spectrum of the environment density matrix.The tail of this spectrum contains only a tiny fraction of the total eigenvalues "mass" and thus can be truncated.b A typical time evolution of the environment density matrix spectrum.It gets wider in time and therefore the principle subspace dimension r(k) grows in time.c The environment truncation scheme: (i) One inserts projectors ω(k)ω † (k) on the principle subspace in between of blocks forming the environment network.(ii) One contracts tensors in dashed red boxes and ends up with new yellow blocks that form truncated environment network.

( 1 )
l→m (k) is the first component density matrix and

( 2 )
l→m (k) is the second component density matrix.Both

FIG. 11 :
FIG. 11: a Diagrammatic representation of the quantum channel Φ l→m (k).b Diagrammatic representation of the corresponding Choi matrix Ω l→m (k).Note, that the only essential difference between Φ l→m (k) and Ω l→m (k) is in the multiplier 1 d S .c Diagrammatic representation of

( 2 )
l→m (k).Note, that due to the TP property of Φ l→m (k),

( 2 )
l→m (k) is proportional to the identity operator.