Floquet engineering with quantum optimal control theory

Floquet engineering consists in the modification of physical systems by the application of periodic time-dependent perturbations. The search for the shape of the periodic perturbation that best modifies the properties of a system in order to achieve some predefined metastable target behavior can be formulated as an optimal control problem. We discuss several ways to formulate and solve this problem. We present, as examples, some applications in the context of material science, although the methods discussed here are valid for any quantum system (from molecules and nanostructures to extended periodic and non periodic quantum materials). In particular, we show how one can achieve the manipulation of the Floquet pseudo-bandstructure of a transition metal dichalcogenide monolayer (MoS2).


Introduction
When a physical system is continuously driven by a time-dependent periodic perturbation, its properties change with respect to those of the unperturbed state. By properties of a system, one refers both to its structure, i.e. the arrangement of its constituents, and to the manner in which it reacts to other perturbations, i.e. its response functions. If these are changed, one can effectively design a new system (material, molecule, optical lattice, etc) with desirable properties by tuning the shape of the periodic perturbation. The name coined for this procedure is Floquet engineering, since one often uses the mathematical tools of Floquet theory to define then quantum structure of stationary states of driven systems, and to simplify the study of periodically driven systems-although this is not mandatory and one could deal with the problems using, e.g. higher order perturbation theory or the straightforward solution of Schrödinger's equation.
Floquet engineering of materials was first proposed in the context of Floquet topological materials [1][2][3]: it was shown how the topology of the light dressed electronic structure can be changed from that of the groundstate electronic structure. The concept has been successfully realized in optical lattices [4] (e.g. the Haldane model in an honeycomb lattice, a model with topologically distinct phases, was realized [5]; Zenesini et al [6] achieved the dynamic control of the quantum phase transition between a bosonic Mott insulator and a superfluid). Signatures of Floquet topological materials have been observed also in solid state materials, namely the opening of a topological gap [7] and the appearance of an anomalous Hall current [8]. There have been numerous further proposals for Floquet induced changes in materials [1,9], among them control of topology in Weyl-and Dirac-semimetals [10,11], induced topology [12] and higher order topologies [13].
In most occasions, the periodic perturbations that have been tried have simple monochromatic shapes. However, the number of possibilities obviously explode if one allows for complex multi-frequency shaped drivings. The difficulty, in that case, lies in navigating the huge parameter space needed to describe the many possible periodic functions that could lead to the desired outcome. Recently, we have demonstrated [14] one way to use quantum optimal control theory (QOCT) [15,16] to overcome this difficulty: QOCT can be used to find the shape that best achieves a given target in terms of the properties that one wishes for the driven system to have. In [14], in particular, the system of choice was a simplified tight-binding two-band model of graphene, and the target was the creation of a driven system with tailored electronic Floquet pseudobands.
In this work, we discuss several different procedures to exercise this kind of quantum control on Floquet systems. In addition to describing in more detail the method that was used in [14], we also propose some other possibilities, that may be more suitable for other situations, depending on the problem definition and on the computational requirements.
The key difference of this control formulation with respect to the usual optimal control is the consideration of periodic solutions to Schrödinger's equation. In the usual control problem, one is normally concerned with an initial value problem: the state of the system is given at some initial time, one establishes a target functional of the system trajectory, and seeks for the control parameters that maximizes that functional (e.g. the occupation of some target state at the final time of the propagation). In contrast, in the control problem that we are attempting here, we do not have a fixed initial state, but we work with periodic solutions to the equations of motion. One also establishes a target functional of the system trajectory (albeit using Floquet pseudoenergies and modes as they are convenient computational proxies), and seek for the control parameters that lead to the periodic solutions that maximizes that functional.
We note that we restrict the discussion to isolated systems-except for, of course, the presence of the periodic perturbation. The presence of an unavoidable environment, and the concomitant dissipation effects, should be considered by an extension of these methods that will be considered in forthcoming publications. Therefore, in this work we will consider the possibility of well-isolated states, whose prethermalized lifetimes are sufficiently long compared to the relevant frequencies and periods. This situation will be easier to achieve experimentally for cold atomic systems in optical lattices than for materials. Section 2 presents and discusses various methods for optimization of Floquet systems-some of them have already been employed in previous publications (by ourselves and by other groups), and some of them will be exemplified here for the first time. Those numerical examples are shown in section 3. Concluding remarks are given in section 4.

Methods
Let us start by formulating the problem that we want to solve. We consider a periodically driven system characterized by a time-periodic Hamiltonian where T is the time period, and v = v 1 , . . . , v M is a set of parameters that determine the shape of the time-dependent Hamiltonian. Typically, but not necessarily, H(v, t) = H 0 + g(v, t)V, where g(v, t) is a real time-dependent function, and V is the coupling operator to the external perturbation. The parameters v can be, for example, the amplitudes of the Fourier expansion of g. This Hamiltonian determines the evolution operator at all times; in particular at the periodic time T, U(T) (we arbitrarily set the origin of time at t = 0). And, in turn, this evolution operator defines the Floquet modes and energies (the quasienergies or pseudoenergies), as its eigenstates and eigenvalues (or, more precisely, the arguments of the eigenvalues): or, equivalently: Note that the pseudoenergies are only defined modulo Ω = 2π T , i.e. ε α + mΩ ≡ ε α for any integer m. In the following, as it is usually done, we assume ε α to belong to the first Floquet-Brilluin zone [−Ω/2, Ω/2).
The eigenstates u 0 α thus defined are time-independent objects. One however normally speaks of the Floquet modes or Floquet states as the periodic time-dependent objects defined by: The gist of Floquet theory is the fact that any solution |ψ(t)⟩ of Schrödinger's equation can be expanded as a linear combination of these modes: Thus, the Floquet modes constitute a valid basis set for the expansion of the time dependent evolution of the system.
It is clear that modes and pseudoenergies are functions of the evolution operator from which they are computed via equation (2): we could write ε α = ε α (U) and u α = u α (U)-hereafter, it the time argument is not specified, U ≡ U(T). Since, in turn, the evolution operator depends on v (we could write it as U = U(v)), they can also be considered functions of the parameters v; i.e. we could also write ε α = ε α (v) and Any property of the driven system can be expressed in terms of the pseudoenergies and of the Floquet modes just as any property of the static system can be expressed in terms of energies and eigenstates of the Hamiltonian. Let us therefore consider any (real) function of all of these: This includes any possible observable of the system, response property, etc. Without loss of generality, we now assume that the goal is to design a driven material whose f function is as high as possible: this is the merit function, in control theory terminology. Although we will present and solve specific examples of control objectives later, we discuss one possibility here for the sake of clarity. In many applications of QOCT, the goal is to find a control function that leads to the maximization of the expectation value of some operator O, often at the final time of the propagation, ⟨ψ(T)|O|ψ(T)⟩. In the context of periodically driven systems, it would probably be more relevant to consider the time average, 1 T´T 0 dt ⟨ψ(T)|O|ψ(t)⟩. If we consider equation (6), and assume that only the lowest energy level α = 0 is occupied (f 0 = 1), one arrives to the following definition for function f : In this case, the definition is independent of the pseudoenergies ε.
In the examples that we will show below, the situation is the opposite: since the goal is the manipulation of the band structure of a material, we will define a function f that is independent of the modes u and only depends on ε.
Whatever the definition, since pseudoenergies and modes are functions of the evolution operator, we can consider a function of evolution operators defined as: And since the evolution operator itself is determined by the parameters v, we can finally define a function The quantum optimal control problem can thus be posed as find the parameters v that lead to the maximization of function G.

Gradient-free optimization
The first and most direct and obvious approach to solving this problem is using a gradient-free approach to the maximization of function G. This procedure was for example used by Zhang and Gong in [17]. Furthermore, in that case a gradient-based procedure could not be used, because the function G was in fact piecewise-constant and discontinuous: the goal was to find topological invariants, and in consequence the key object is the Chern number.
In order to use a gradient-free approach, one only needs of a method to compute G(v). One way to do that is by solving Schrödinger's equation for the evolution operator: Once the evolution operator at time T has been computed, one can obtain the pseudoenergies and modes, and compute function G using equation (10). There is however an alternative to the computation of the evolution operator by explicit propagation: one can reformulate the theory in the Floquet-Hilbert space in the frequency domain (see section 2.4) in order to compute pseudoenergies and modes.
In any case, given any procedure to compute G(v), one can use one of the many optimization algorithms that do not require gradients. For example, the above-mentioned work [17] used a particle swarm optimization method.

Gradient-based optimization: target function F defined directly in terms of the evolution operator
If the gradient of G can be computed at a reasonable cost, the algorithms that use the gradient will be more effective than the gradient-free ones. QOCT provides one expression for this gradient [18,19]: where the inner product '·' in operator space is the Fröbenius product: This expression requires, in addition to the propagation given in equation (11), the propagation of the so-called costate B: Equation (13) is similar to the equation for the real propagator, equation (11)-note that B is also an operator-like object-but the boundary condition given in equation (14) is different. The computation of the gradient of G requires, therefore, two propagations: one for U, and one for the costate B. And, in addition, it requires of the computation of the boundary condition (14), the gradient of F. This computation, depending on the nature of F, may be trivial or very difficult.
Let us consider one example: one possible goal may be to find a driven system whose pseudoenergies match some predefined target values, {ε α }. In order to achieve that goal, the first idea that comes to mind is to set function (7), for example, as: Upon a successful maximization of function f (that achieves its absolute maximum at f = 0), we would exactly get the target pseudoenergies. Note that this definition disregards the modes u α .
The corresponding definition of the target function F(U) would be: One would therefore need to compute the gradient of this function, needed for equation (14), which in turn would require the computation of the gradients of ε α (U). We will discuss this possibility in the next section.
However, there is an alternative route that avoids this computation, and is the route that was followed in [14]. If the goal is to find a field that produces a given set of target pseudoenergies {ε α }, we can directly set a target evolution operator whose decomposition has the corresponding eigenvalues: One way to set this operator as target is to define the target function as: This function is maximized if U is equal (or equivalent, i.e. equal up to a phase factor) to the target evolution operatorŨ. In this way, the calculation of the boundary condition (equation (14)) is trivial: By working in this way, we do not need to worry about how to compute the derivatives of the pseudoenergies and modes with respect to the evolution operator components. In addition, this procedure has the advantage that it was already tested: the application of optimal control methods for the creation of unitary operators was pioneered by Palao and Kosloff [18], even though with a different application in mind: the design of gates in quantum information processing devices.
Some results obtained with this method were presented in [14], and therefore we will not provide new examples here. The caveat of this method is that it may be too restrictive: it essentially asks to find a driven system whose full set pseudoenergies and Floquet modes coincide with a prescribed set. Notice that, in the example given above, where the goal is setting the pseudoenergies, regardless of the modes, one has nevertheless to prescribe a set of arbitrary modes in equation (17) -at least the 'static part' , u 0 α . For many target definitions, it may be better not to force a particular set of pseudoenergies and modes, but rather to directly optimize the target function f. Two possible routes for this purpose will be described in the following sections.

Gradient-based optimization: general targets
Let us now go back to the general target definition in the form: that may also be written as: defining F(U) = f(ε(U), u(U)), and recall that the QOCT formula for the gradient of this function is: Let us now assume that one cannot, as in the previous section, write F as an explicit function of U, which would lead to a simple expression for the boundary condition (24) such as equation (19). We face the problem of computing the gradient elements ∂F There are a couple of technical issues that must be considered. First of all, note that function F is a real function of a set of complex parameters (the components of the unitary operator U). We wish to compute the derivatives of the function with respect to its parameters. But we must remember that these derivatives cannot be normal complex derivatives as defined in complex analysis theory: any real function of complex arguments is not differentiable anywhere (unless it is a trivial constant function).
This difficulty is often solved by rewriting this kind of functions in terms of both the arguments and of the complex conjugate of the arguments, as if they were independent of each other. Then, one may take the derivatives with respect to the variables, or with respect to their conjugates-in fact, that is the meaning of the notation ∂F ∂U * ij . This procedure is not obvious if, as it is the case in hand, one cannot explicitly write the function as a function of both the variables and the conjugate of the variables.
In this case, one must go back to the formal theoretical background behind the apparently illogical procedure of using the variables and their complex conjugates as if they were independent. This is the concept of Wirtinger derivative [20][21][22]. The basic idea is that, even though the function F is not differentiable in the context of complex analysis, if we consider it a real function of the real and imaginary parts of its arguments (and therefore, it is a real function of real arguments), it may be perfectly differentiable. One may then define the following derivative: It is easy to check that this definition leads to the usual procedure of taking the derivative with respect to the complex conjugate of the variables, keeping the variables themselves constant. However the definition works even if it is not easy to identify the decomposition into variables and conjugate variables, and is in fact the real rigorous mathematical definition of ∂ ∂z * [20]. A second technical subtlety is the following: let us consider the computation of the previous derivatives using the finite-differences expression: where E (ij) is the matrix that has all elements equal to zero, except for the ij element which is one. The meaning of ε α (U) is the αth Floquet pseudoenergy associated to the unitary operator U. It only makes sense if U is a unitary operator. In that case, it can be diagonalized, and the eigenvalues are complex unit numbers e iθ . Likewise for the modes u 0 α , which are the corresponding eigenvectors. However, in general, U + ∆δ ij need not be a unitary operator. Therefore, the functions ε α (U + ∆δ ij ) and u α (U + ∆δ ij ) are not well defined. One way to solve the problem is to extend the definition of functions ε α and u α , so that they are defined for any operator U, and not only unitary ones: they can be defined in terms of the elements of the diagonal of the Schür decomposition (to define ε α ) and of the corresponding Schür vectors (to define the u α ). Those objects exist for any matrix U, and when the matrix is actually diagonalizable, they coincide with the eigenvalues and eigenvectors.
Once we have clarified those two technical issues, it remains to actually compute the derivatives in equations (26) and (27). While there are ways to work out analytical expressions for ∂ϵα ∂ReU ij , etc and proceed by using the chain rule, for relatively small systems it is feasible, as we will show in the example shown below, to simply use a finite difference formula. Note, however, that this can become computationally intensive if the dimension of the quantum space is large.

Gradient-based optimization based on perturbation theory
We finally describe a method for the computation of the gradient of function G(v) that somehow bypasses the standard QOCT, and uses instead perturbation theory on the Floquet-Hilbert space. In fact, the following procedure may be the one that is best suited for larger systems, as it avoids the explicit construction and manipulation of the full evolution operator.
In the previous equations, the calculation of the Floquet modes and pseudoenergies was assumed to be done through the explicit construction of the evolution operator of the system, and its diagonalization. However, the cost of this procedure grows unfavorably as the dimension of the system grows. For larger systems, the usual procedure is to move to the 'Floquet-Hilbert' space of time-periodic states, and compute the pseudoenergies and Floquet modes from the diagonalization of the extended Hamiltonian.
Let us recall the key objects to fix the notation: let T be the space of time-dependent periodic square-integrable complex functions defined in the interval [0, T]; we may define an inner product in this space as: An orthonormal basis of this state is formed by the Fourier functions: that verify (n|m) = δ nm and n |n)(n| = I, thus effectively defining the Fourier expansion of any function: If H is the Hilbert space of the quantum system, we may define the 'Floquet-Hilbert' space as F = H ⊗ T . It can be understood as the set of all (periodic) quantum trajectories t → |ψ(t)⟩; the inner product in F is given by: Note the double angles to distinguish the inner product in F from the inner product in H. Often, states in F are also denoted with a double angle: |ψ ⟩⟩ ∈ F.
The Floquet modes u α , defined with equations (2) and (4), verify the following equation: If we regard the Floquet modes as elements of F and call them |u α ⟩⟩, this equation can be rewritten as: where now H and i d dt are operators in F and, moreover, we have made explicitly the dependence on the periodic driving parameters v. In this way, the time-dependent propagation equation is transformed into an eigenvalue problem.
Let us go back to our optimization problem. The goal is to optimize some function of a set of parameters v that determine the Hamiltonian, and this function is defined in terms of the pseudoenergies and modes: One may compute the gradient components ∂G ∂vm using the chain rule, if there is a way to compute the gradients ∂εα ∂vm and ∂uα ∂vm . But these can be obtained using the previous equation (33): one can now apply perturbation theory in the Floquet-Hilbert space. Thus, for example, for the derivatives of the pseudoenergies, we would use an extension of the Hellman-Feynman theorem: And, for the derivatives of the Floquet modes one can use an extended version of Sternheimer's equation: The solution of equations (36) and (37) provides a practical way to compute the derivatives of the pseudoenergies and Floquet modes with respect to the control parameters v, and in consequence one way to do quantum control for Floquet systems [23] In practice, one also needs to expand the previous equations into a complete basis. First, let us consider a generic basis {|i⟩} of the Hilbert space. We can then consider the basis {|im⟩ ≡ |i ⟩ ⊗ |m)} for F. Then, equation (33) takes the form: where: Equation (38) is an eigenvalue equation; if we truncate the quantum Hilbert space to a space of dimension d, and the Fourier basis expansion to M functions, the dimension of the problem is dM.

Examples
In this section, we will study in detail one example of application using the two new techniques described above (those in sections 2.3 and 2.4). The methods described above have been implemented in the qocttools code [24], used to perform the following calculations.
For the sake of clarity and generality, the discussion of the QOCT method presented in the previous sections was for an arbitrary quantum system. In the following, however, we will work with periodic materials, which benefit from the decomposition of the Hilbert space into k-points that derives from Bloch's theorem. This entails a huge computational simplification-at the cost of slightly complicating the notation. Thus, many of the objects introduced above should now carry a k-point index: ε kα , u kα , H k (v, t), U k (v), etc. The generalization of the expressions given above is fortunately rather trivial, and it just sometimes involves summation over the Brilluoin zone (or over whatever set of points in a region of the Brillouin zone that is deemed to be relevant). For example, equation (12) becomes: To exemplify this technique, we have chosen to work with a monolayer of a transition metal dichalcogenide (TMD). We have used a three-band tight-binding model of a MoS 2 monolayer (depicted in figure 1), fitted to reproduce the DFT first principles electronic structure [25]. In particular, the model that includes up to the third nearest neighbor, and has been parametrized using the generalized gradient approximation. The field-free Hamiltonian at each k point is given by a 3 × 3 matrix H k , whose elements are parametrized functions of k. It gives rise to a band structure depicted in figure 2 (red solid lines). Notice the direct gap at the K point, of almost 2 eV.
We must now consider the monolayer to be subject to an external time-dependent electric field. It is given, in the velocity gauge, by a spatially homogeneous vector potential A(v, t), where v are the set of parameters that determine the particular temporal shape of the field, and whose variation we may use in order to modify the material properties. The time-dependent Hamiltonian of the driven system can then be modeled using a Peierls substitution: A(u,t) .
The control parameters u 1 , . . . , u P will be the Fourier coefficients that determine the shape of the vector potential A(t): A y (u, t) = M n=1 u 2M+2n cos(Ω n t) + u 2M+2n−1 sin(Ω n t) .
The frequencies Ω n = 2π T n (n = 1, . . . , M) are multiples of the fundamental floquet frequency Ω = 2π T . In our case, we have set this frequency to Ω = 2 eV, in order to make it slightly larger than the band gap. We allow for two electric field components along the x and y components. The application of such periodic external fields may then be studied in terms of the Floquet states and pseudobands. Let us consider, for example a monochromatic circularly polarized field: where λ is set such that the electric field has an amplitude of 4 V nm −1 . The corresponding Floquet pseudobands are also depicted in figure 2 (blue dots). Recall that the pseudoband energies ε kα are only defined module Ω, and one may consider the bands to be folded within the Floquet-Brillouin energy range [−Ω/2, Ω/2). Additionally, one may then define the sidebands, by vertically shifting them by multiples of Ω: ε m kα = ε kα + mΩ, for any integer m are also valid pseudoenergies. These are replicas of the first pseudobands, but contained in the displaced energy intervals [−Ω/2 + mΩ, Ω/2 + mΩ), and are also shown in figure 2. These sidebands are actually visible in angular-resolved photoelectron spectroscopy experiments [26]. In They are therefore most visible if the mth frequency component of the corresponding Floquet mode is larger.
It can be seen how the Floquet pseudobands, for this reference strong circularly polarized driving field, are significantly different to the field-free bands. Several sidebands are visible (⟨u m kα |u m kα ⟩ > 0) for various values of m, which means that the associated Floquet modes oscillate with several multiples of the fundamental 2 eV frequency.
Let us now concentrate on the region around the K point; the two lowest field-free bands (one is the valence band, the other one is a conducting band) approach there, although a gap of nearly 2 eV remains. Suppose that we wish to modify the material in such a way that the Floquet pseudobands are much closer, or even close the gap. The first step is the definition of function f (equation (7)), which in this case can be defined solely in terms of the pseudoenergy at the K point, regardless of the Floquet modes; for example we may try: Therefore, the function only depends on the pseudoenergies at point K; the function grows with increasing values of the lowest pseudoenergy ε K0 , and also grows with decreasing values of the highest pseudoenergy ε K2 . The goal is therefore for the shaped optimized fields to induce a larger lowest pseudoenergy ε K0 , and a smaller highest psuedoenergy ε K2 , thus approaching them and leading to a narrower gap.
The optimization seeks for the parameters u ≡ u 1 , . . . , u 4M that maximize function f ; we add to this definition of the problem a bound constraint: In this way, the amplitude of each frequency component is limited to a certain value. We have performed these optimizations using the two procedures described in sections 2.3 (the QOCT formula for general targets) and 2.4) (the formula based on perturbation theory in Floquet-Brillouin space). The function maximization algorithm that we have used is the sequential least-squares quadratic programming algorithm [27] as implemented in the NLOPT library [28]. The results obtained with the two procedureds are, of course, identical, since they provide the same gradient. Figures 3, 4 and 5 display the results obtained when constraining the electric amplitudes to 0.5 V nm −1 , 1 V nm −1 and 2 V nm −1 , respectively. In the first case, the allowed amplitudes are too low, and the Floquet pseudobands, even for the optimized field, barely differ from the original field-free bands. When allowing the electric amplitudes to be of up to 1 V nm −1 per component ( figure 4), the optimized bands start to show a visible closure of the gap. Finally, if the limit is set to 2 V nm −1 (figure 5), the gap is almost fully closed.
Note that, in the definition of function f(ε) = ε K0 − ε K2 , the meaning of ε K0 is the lowest pseudoenergy at the K point, whereas ε K2 is the highest pseudoenergy at the K point (given that there are only three levels). As a consequence, if the highest and lowest band approach each other, the middle band will, too. That is why in the plots we observe that all bands approach. Note also that the motif is repeated for all values of m, i.e. for all sidebands, although only the ones with m = 0 and m = 1 have significant weight in the Fourier expansion of the modes.

Conclusions
In this work, we discussed several practical schemes for the quantum optimal control in the context of the Floquet engineering, by maximizing a target function defined in terms of a set of control parameters that determine the temporal shape of the driving fields. First, we briefly discussed the gradient-free optimization, where one needs to maximize the target function without the knowledge of its gradient. Second, we discussed in section 2.2 a gradient-based optimization scheme, that can be used for a special case of target functions, that are directly defined with the corresponding evolution operator. This gradient-based scheme has been used in our previous work [14], where it was demonstrated that it is capable of properly optimizing the Floquet states and their properties.
In addition to those methods, we further developed two novel optimization scheme for the Floquet states, capable of dealing with general targets, in sections 2.3 and 2.4. In particular, the last one is based on perturbation theory, and it does not require the explicit construction and manipulation of the evolution operator. Since the computational cost of the explicit construction of the evolution operators grows badly with the system size, this optimization scheme is expected to be suitable for applications to large system models, such as when using full ab-initio simulations on solid-state systems (for example, using time-dependent density-functional theory [10]).
To demonstrate these optimal control schemes, we have applied the gradient-based optimization schemes discussed in sections 2.3 and 2.4 to the bandstructure engineering of a monolayer of TMD, MoS 2 . Here, we intended to reduce the bandgap of the MoS 2 at the K point by optimally tuning the parameters of the driving field. As a result of the optimization, we demonstrated that the optimized driving fields can indeed reduce the bandgap of MoSe 2 at the K point, and in the strong field regime, can even close the gap.
The optimal control methods for Floquet systems described in this work are practical prescriptions for tuning the parameters of the driving fields in such a way that the systems may exhibit desired properties and functionalities, beyond the equilibrium phase. This opens a path towards an efficient engineering of material properties on demand [29][30][31]. We finish by noting that there may be another interesting application of this version of control theory: since there is a deep analogy between Floquet states of periodically-driven matter and quantum-electrodynamics (QED) states of matter in an optical cavity, one may extend this Floquet optimal control procedure to the QED materials engineering [32].

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).