{'fraction_non_alphanumeric': 0.03314638250140213, 'fraction_numerical': 0.007646289025986166, 'mean_word_length': 4.71546105353136, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Defining a well-ordered Floquet basis by the average energy

At the moment, the most efficient method to compute the state of a periodically driven quantum system is using Floquet theory and the Floquet eigenbasis. The wide application of this basis set method is limited by: a lack of unique ordering of the Floquet eigenfunctions, an ambiguity in their definition at resonance, and an instability against infinitesimal perturbation at resonance. We address these problems by redefining the eigenbasis using a revised definition of the average energy as a quantum number. As a result of this redefinition, we also obtain a Floquet-Ritz variational principle, and justify the truncation of the Hilbert space.


I. INTRODUCTION
The works of Shirley and Sambe [1,2] introducing Floquet theory to quantum mechanics, have enabled the efficient calculation of the time-periodic Schrödinger equation and kicked off renewed interest in periodically driven quantum systems (Floquet systems). The Floquet method has since become common place when studying model systems with periodic driving [3][4][5]. Applications beyond model systems, and towards realistic systems derived from first-principles, are still limited by a lack of a proper definition of the ground state and a variational method which can derive it efficiently.
We condense the common problems of the Floquet method in the left side of fig. 1, along with an often disregarded, yet intuitive way to solve them in the right side of the figure. Just like the energy levels of static systems, Floquet systems are characterized by the quasi-energies n and their corresponding Floquet eigenfunctions u n (t).
where H(t) = H(t + T ) denotes the time-periodic Hamiltonian of the system we wish to study. Complications arise from the appearance of countably infinite equivalent solutions to eq. (1), with quasi-energies offset by a multiple of the driving frequency ω = 2π T , but representing the same physical state, i.e. the same solution of the original Schrödinger equation. Thus it is impossible to label the eigenstates and define a ground state by the quasi-energy alone, as it is made apparent in the left side of fig. 1. On the other hand we find hints of ordering when we consider the eigenstate's energy spectrum, as defined by: We show this in the right side of fig. 1 1) and (2)) of a typical Floquet system. The energy spectrum amplitude of each state is depicted by varying thickness and opacity.
are chosen to depict three common features present in a typical Floquet system. Firstly we can see that a ground state can be uniquely defined (state Ψ 0 ), given a robust enough criterion, e.g. ordering by the average of the spectrum. Secondly the energy spectrum is plagued with near resonant states which creates an ambiguity in how to define the eigenstates. This is similar to the degeneracy problem of static systems, e.g. quantization axis of the Hydrogen's 2p 1 states. In the Floquet system here, resonant states Ψ 1 and Ψ 2 are not uniquely defined from the solutions of eq. (1). This ambiguity can be resolved by using additional quantum numbers, and the ones based on the energy spectrum are the most intuitive. Lastly there are numerous states Ψ n with relatively similar quasi-energies, but well separated energy spectra compared to the low-lying/active states. Such states should not affect the physical system significantly and we can truncate the Hilbert space and significantly simplify the calculations needed.
In this paper we propose the average energy to be the central parameter to label and define the Floquet eigenstates. We show that this definition has a variational derivation, is stable against infinitesimal perturbations and allows for the truncation of the Hilbert space. Before diving into this investigation, we present a short history of the problem of defining a Floquet ground state in section I A, where the reader can find various alternatives that are currently used. We summarize the potentially ambiguous notations in section I B for the reader to refer back to if the equations become unclear. At the beginning of section II we present the basic derivation of the Floquet method with special emphasis on the problems concerning the eigenstates definition. The infinite time average energy and its corresponding eigenstates are defined in section II A. For practical applications we derive the observed average energy in section III A, with the reasoning explained at the start of section III. Using that definition we derive an equivalent Ritz variation method in section III B. We use the two-level system with a circularly polarized driving as the minimal example system and exemplify the prior discussions in sections II B and III C.

A. History of the problem
The Floquet method has been thoroughly studied for decades and the problems of labeling the eigenstates and defining a ground state have been known from the beginning [1,2,6]. So far there have been three ways of solving this problem, by labeling the eigenstates using: a Floquet-Brillouin zone in quasi-energy space [2,4], an adiabatic continuation [7][8][9], or a perturbation method [2,[10][11][12]. These methods are commonly applied to finite model systems where in principle the full Hilbert space is accessible, and there is minimal need to label the eigenstates. On the other hand, these methods fail when applied to realistic systems, where the Hilbert space has infinite degrees of freedom. The incompatibilities in the realistic system boil down to the resonance condition of infinitely many eigenstates and their avoided crossing problem [8,13]. While there have been various adiabatic methods proposed to overcome this problem [8,9,14], there is still a need for eigenstate labeling methods which are efficient, do not rely on an adiabatic continuation and can be applied to realistic systems.
Another common approach is to circumvent the need to identify any particular eigenstate and extend the system to an open quantum system [6,[15][16][17]. This is analogous to changing the problem from calculating the ground state to calculating the thermal equilibrium, which in this case is commonly referred to as the Floquet steady-state. This is often the most achievable and physically relevant state. Since the steady-state does not generally have a Gibbs distribution [6,15,[18][19][20][21], at the low temperature limit the system does not generally reduce to a pure ground state, which can diminish the significance of labeling the eigenstates and defining a ground state. However, even here, the eigenstate labeling problem is not fully circumvented, since there is still no criterion for truncating the Hilbert space. Thus the computation cost become unfeasible in realistic systems.
To give an example of the problems coming from the lack of unique labeling, we look at the attempts of formulating Floquet ab-initio methods [22,23]. These methods were constructed to calculate a poorly defined ground state, for which the Ritz variational principle is assumed. One of these methods, the Floquet Density Functional Theory [22], was shown to be fundamentally flawed due to the ambiguous definition of the ground state [24][25][26], an argument that can be extended to the other variational based ab-initio derivations as well. One method that remains valid is the Floquet Density Matrix Renormalization Group [27] and similar iterative methods. But even there, the solution and its convergence is dependent on the selection criterion used at each iteration. A starting point for reformulating these ab-initio methods is to uniquely redefine the ground state of a general Floquet system, particularly one that is applicable to realistic Hamiltonians.
The potential usage of an average energy to label the eigenstates has been considered before, either indirectly [13] for selecting relevant adiabatic states, or directly [28] for deriving an effective Gibbs distribution. It should be pointed that the average energy definition there differs from the one we propose in this paper. This usage is often disregarded since other observables offer a better agreement, e.g. regular energy [28]. However this does not rule out the potential usage of the average energy to approximate and select the significant eigenstates for describing the steady state solution or truncating the Hilbert space. Using the average energy definitions from [28,29] can still be problematic as we will explain in more detail later on.
To resolve these problems we derive a more robust definition of the average energy.

B. Notations
In this paper we use the following notations in order to keep the equations more compact. The details are clarified throughout the paper as it becomes relevant. Readers can skip this section and come back if any notations appear unclear. These notations are not in any particular order, so we just present them below as is: Fourier decomposition: Average expectation: Energy difference The time-dependence will often be dropped, unless explicitly needed for clarity. In general we reserve the superscripts to indicate different variations of the functionals, eigenstates, etc.. These can often be mixed with each other. If the superscript is missing any definition is applicable and/or the perturbed/untruncated definition is used, depending on the context.
E T Observed average energy, E i Truncated average energy at step i.
Subscripts are generally reserved to label the eigenstates, eigenvalues, etc.. We reserve the following subscript (or superscript in the case of truncation) notations: U is reserved for the transformation matrix from an eigenstate representation to an arbitrary basis, and is constructed from the vector columns of the Floquet eigenfunction Transformation matrix Hilbert spaces are denoted like H, and the equivalent eigenspace specific to each Hamiltonian H is E H . Eigenspace: (15) The Floquet interaction picture with respect to interaction v (H = H 0 + v) is denoted by the subscript I .
Non-interacting propagator: Interaction picture propogator Interaction picture operator Finally we use the usual δ notation to indicate functional variation in arbitrary direction, derivative, etc., and we add a constrained variation notation. The notation is generally omitted if we consider the full Hilbert space.

II. EXACT FLOQUET SYSTEMS
We start off with a thorough examination of the Floquet theory in the closed quantum systems. Consider the Schrödinger equation of a quantum system described by the time-periodic Hamiltonian H(t + T ) = H(t), with a corresponding driving frequency ω = 2π T Using Floquet theory we can solve this equation as an eigenproblem resembling the solution of a static Hamiltonian From these eigenfunctions we can derive the propagator and the wavefunction in eq. (20) at arbitrary times. This is analogous to the time-dependent wavefuntion solution of a static Hamiltonian.
Here C n is the usual overlap at a known time point t = 0, and will be used throughout this paper. In order to preserve the norm of the wavefunction, the summations in eqs. (21) and (22) and further on, are limited to different Floquet eigenstates, which will be clarified shortly. It is important to note that the quantum system and all of its observables are strictly determined by the wavefunction Ψ(t), and not by the solutions of the Floquet Hamiltonian H F (eq. (1)) Unfortunately we cannot directly calculate this wavefunction Ψ(t) or propagator U (t) efficiently. Instead we calculate them indirectly through the eigenpair of quasi-energy and Floquet eigenfunction ( n ; u n (t)) as in eqs. (21) and (22). Although these eigensolutions do not have a direct physical meaning, they are computationally accessible, and sufficient to describe arbitrary states. Deriving the Floquet eigenstates is the only reliable and efficient way of describing the system and it is thus the main objective of this paper Early on, Sambe has formalized the Floquet system's extended Hilbert space HT ∋ u(t), upon which the Floquet Hamiltonian operates and defines the eigenspace E H [2]. This extended Hilbert space is the tensor product of the original Hilbert space H on which H(t) operates at any given time, and the countably infinite Fourier space T which guarantee the periodicity of the Floquet functions, in accord with Floquet theory [30]. In this extended Hilbert space, the Floquet Schrödinger equation has a systematic block matrix form: The Floquet eigenfunctions diagonalize this Hamiltonian, and form an orthonormal complete basis set spanning HT.
Sambe also showed that a variational principle on the quasi-energy is possible, with the stationary points corresponding to the Floquet eigenfunctions.
This variational principle offers an alternative method of deriving the Floquet eigenfunctions, which is more computationally efficient and necessary for many theoretical formulations. From here on, we will assume that the exact eigensolutions { n , u n (t)} are known. The infinite dimension of the Floquet Hamiltonian in Fourier space T suggest that there are a countable infinite number of eigenpairs { ′ n ; u ′ n (t)}. These infinite solutions can be grouped into subsets, with the solutions related to each other by harmonic shifts (eq. (32)), and describing the same wavefunction Ψ n (t) (eq. (33)).
We refer to the wavefunction Ψ n (t) as the physical Floquet eigenstate, to distinguish it from the ambiguous Floquet eigenfunction u n (t). The physical Floquet eigenstates form a complete basis set and fully describe the propagator in eq. (22). These basis sets span the Hamiltonian's Hilbert space H instead of the extended one HT.
Because of this ambiguity, only a subset of the eigenfunctions u n (t) are needed in summations like eqs. (21) and (22) in order to span the Hilbert space H and describe the propagation of arbitrary wavefunctions. The choice of eigenfunctions is arbitrary, as long as the Hilbert space H is fully spanned, or equivalently there are no pairs related to each other by a harmonic shift (eq. (32)). We refer to these choices in eigenbasis as a choice of quasi-energy shift. Any observable, interaction, etc. are independent of this choice, except for the quasi-energy which is not a true observable according to eq. (25). From here on we assume that an arbitrary choice of quasienergy shift is performed and labels m, n describe distinct physical Floquet eigenstates.
Another important consequence of the expanded Hilbert space, is that the Floquet eigenfunctions are not uniquely defined within the subspace of resonant eigenfunctions. That is to say, we can rotate the basis {u a (t)} within the resonant subspace a = b = , and obtain an- In this paper we reserve the labels a, b to describe resonant Floquet eigenfunctions, and we assume that we have selected the subset of eigenfunctions so that the resonant eigenfunctions have the same quasi-energy a = b = .
For each resonant set with quasi-energy , HT denotes the resonant Hilbert subspace to which it belongs. This ambiguity problem is analogous to the degeneracy problem of static systems, where for example in the Hydrogen 2p 1 degenerate space, we can arbitrarily choose the axis of quantization, and any choice gives valid Hamiltonian eigenbasis. However there are two caveats in the Floquet system compared to the static case. Firstly, the number or density of resonant Floquet eigenfunctions is practically infinite in realistic systems, e.g. the Volkov states [36,40], while in general the static system is finite, with the exception of flat band systems which are still resolvable. Secondly, the energy spectrum as defined in eq. (2) differs between different resonant Floquet eigenfunctions ( fig. 2b), and is dependent on our choice of eigenbasis.
A consequence of this energy spectrum difference is that, if we include thermodynamic effects, resonant eigenstates interact differently and are not equally distributed at equilibrium [41]. In the extreme case of highly separated energy spectrum, we can find a resonant eigenbasis where only a few states are occupied in the steady-state, so that we can truncate the Hilbert space without affecting the description of the steady state. If resonant eigenfunctions share the same energy spectrum, these are indeed degenerate.
In the case of a static system, we can use an additional quantum number or an adiabatic continuation to label the eigenstates at and around the degeneracy point. In limited cases [8], we can also use the adiabatic method to define the Floquet eigenfunctions around the resonance (e.g. fig. 2a). On the other hand, defining the eigenstates using an additional quantum number would be more efficient, so we aim to find an appropriate parameter that would result in reasonable and/or truncatable eigenbasis (e.g. fig. 1). The most natural choice is one based on the energy spectrum, and the first possible choice is the average energy.

A. Infinite time average energy
We redefine the average energy as simply the timeaveraged expectation value of the Hamiltonian of an arbitrary state/initial wavefunction up to some time T , to be clarified later in the paper. This definition gives a proper observable (eq. (25)), and as such, it is independent of our choice of the quasi-energy shift.
This definition is closely related to the energy spectrum average (eq. (40)), and can be adapted for various theoretical formulations.
Before defining the averaging time T , we will address the previous definition of the average energy [28,29]. There the average energy has been defined as the average energy expectation value over a period T , only of a Floquet eigenfunction u n (t) (derived from eq. (1)).
This value is also independent of our choice of quasienergy shift (eq. (32)) and thus could be used for labeling the eigenstates. An obvious flaw of this definition is that it is not defined at resonance, since the Floquet eigenstate itself is not uniquely defined there. Another problem is that, the generalization of this definition to arbitrary states (eq. (43)) is incompatible with the variational principle, and is not defined by the Floquet eigenstate's average energies¯ n (regardless of their definition). Nevertheless, we refer to this generalization as the effec-tive average energy, to be used in later derivations.
Effective average energy: This suggest that we cannot use this average energy definition to derive a Floquet eigenstate directly, and we need a different average energy definition.
Coming back to the definition in eq. (39), we can simplify the equation by using the exact propagator (eq. (21)) expanded in an arbitrary choice of eigenfunctions: In this section we choose the averaging time T to be the limit at infinity, and refer to the resulting average energy as the infinite time average energy, and equivalently the eigenstates derived from it.
For now we assume there are no resonant eigenfunctions (mod( n − m , ω) ≠ 0). The infinite time average energy functional trivially simplifies to the weighted sum of the eigenstate average energies Where due to non-resonance the eigenstate average energies are unambiguously defined as either the effective or infinite time value¯ T n =¯ ∞ n =¯ n . Immediately we see the benefit of this definition that we can derive the eigenstates from the variation of the functional form. The variational principle in this case is over the initial wavefunction Ψ(0), as opposed to the previous variations over the Floquet functions u(t).
In theory this variation is sufficient to define and find all eigenpairs The superscript in E 0 H indicates that only the solutions at t = 0 are calculated and the Hilbert space is limited to H, not to be confused with the notation E H 0 . From these we can also derive the time-dependent Floquet functions u n (t) and quasi-energies n from the Schrödinger equation, and subsequently the propagator and any observable of arbitrary states at arbitrary times.
In practice however, the exact propagator U (t) in eq. (39) is inaccessible, and the variation with respect to the initial wavefucntion Ψ(0) is computationally unfeasible. Fortunately as long as there are no resonance conditions, we can use the quasi-energy variation principle (eq. (30)) to derive the Floquet eigenstates and the effective average energy (eq. (43)) which coincide with the infinite time average energy for the non-resonant Floquet eigenstates.
Thus we can define a well ordered eigentriplet ., up to a trivial quasienergy shift. The correspondent eigenspace E ∞ H spans the extended Hilbert space HT (having included equivalent solutions eq. (32)), and the propagator is directly obtained from eq. (21).
From these eigentriplets we can uniquely define the ground state ( 0 ,¯ ∞ 0 , u 0 (t)), up to a trivial quasi-energy shift. This ground state can be derived variationaly through the initial wavefunction Ψ(0), Floquet function u(t), or using a Lagrange minimization method.
The Lagrange multiplier λ and functional derivative δ [u] δu are vectors spanning the Hilbert space dimensions, and guarantee the minimization in¯ T [u] is taken over the Floquet eigenfunctions. The Lagrange multiplier¯ guarantees the normalization constraint is satisfied. Although this form is more complicated with the addition of λ, it can facilitate the derivation of other approximations, e.g. a Hartree-Fock variant, and we include it here for future reference.
As for the resonant region, we cannot use the Floquet eigenfunctions defined by eq. (1) to simplify eq. (46) to eq. (48), since the eigenfunctions themselves are not uniquely defined in this region. This also means that the eigenstate average energies are not well defined from eq. (42), and we have to redefine them here. We limit ourselves now to the resonant Hilbert subspace HT .
Choosing an arbitrary orthonormal basis set in this subspace ⟨⟨u a u b ⟩⟩ T = δ ab , the infinite time average energy has the following form: We point out that the effective average energy func-tional¯ T (eq. (43)) is equivalent to the infinite time average energy functionalĒ ∞ (eq. (56)) within the resonant Hilbert space, which is why we refer to it as effective average energy.Ē HT ∀u a (t) ∈ HT . (58) Equation (56) can be represented in a matrix form, whose eigenvalues¯ ∞ n are independent of the basis set. Diagonalizing this matrix we find a unique eigenbasis, from which we uniquely define the eigentriplet ( ,¯ ∞ n , u n (t)).Ē The coefficients c an are the corresponding eigenvector projection onto the arbitrary basis {u a (t)}. We label the eigenstates here with m, n as to indicate that the ambiguity coming from the resonance is resolved and the eigenstate average energies are defined from the diagonalization of eq. (59). Considering the relation to the energy spectrum, this diagonalization gives us the basis set with the least overlap between their energy spectra and maximally separated average energies ( fig. 2b). The average energies themselves can have real and accidental degeneracies, which can be further resolved by expanding this procedure to the energy variance and other higher order moments. Combining all of the resonant and non-resonant subspaces, we can define the well-ordered eigenspace E H of an arbitrary Hamiltonian H(t). The theoretically sufficient eigenspace defined from eq. (50) remains unchanged upon the combination of the Hilbert subspaces, where now the eigenstate average energy¯ ∞ n is uniquely defined by the infinite time average energy variation/diagonalization: As for a computationally accessible form, we have to rely on the quasi-energy variation. Since the infinite time average energy and effective average energy are equivalent in the resonant space, we can derive the variational principle in two steps over quasi-energy and effective av-erage energy functionals.
We have to be careful of the order of the variations so as to properly limit the variation space of the effective average energy. The variation in eq. (65b) is not valid in the full extended Hilbert space. Equivalently we can absorb these complications into a Lagrange multiplier, and extend the method in eq. (53c).
where crit indicates the search of a critical point, not only a minimum. So far we have a unique description of the Floquet eigenbasis up to a trivial quasi-energy shift E ∞ H = { n ,¯ ∞ n , u n (t)}. This basis set describes the propagator (eq. (21)) uniquely, is compatible with the variational principle, is well ordered through the average energy, and is unambiguous at resonance.

B. Exact two-level system example
The simplest toy model that shows the properties discussed above is the two-level system driven by a circularly polarized interaction. The exact solutions can be found in various textbooks [36], derived in a number of ways. We choose the adiabatically connected solution with varying driving strength.
Here the parameters ω 0 , ω, δ, V, Ω correspond to the natural oscillation frequency, driving frequency, detuning, driving strength and Rabi frequency respectively. We choose a blue detuning ω = 1.5ω 0 , and we normalize all the parameters to ω 0 = 1. We plot the energy spectra and infinite time average energies of the Floquet eigenstates in fig. 2. The corresponding analytic formulas are as follows: In this example the eigentriplets { ± ,¯ ∞ ± , u ± (t)} are defined and labeled according to the adiabatic continuation with respect to the driving strength V . Generally this would not be available due to the ambiguous definition of adiabaticity [8,9]. For the sake of argument we assume that the eigentriplets are still undetermined at the resonance points and we calculate them using eqs. (64) and (66). At the first resonance V = V 1 (such that Ω = ω), we expand the infinite time average energy matrix in the basis of eqs. (69a) and (69b).
The infinite time average energy matrix is already diagonalized, so the Floquet eigenbasis defined from eq. (66) coincide with the adiabtically connected one (eqs. (69a) and (69b)). We can better see the effect of this diagonalization in the eigenstate energy spectrum ( fig. 2b), where we vary the definition of the Floquet eigenfunction by mixing the Floquet eigenfunctions with quasi-energy = 0, or equivalently varying the initial wavefunction Ψ(0) = u(0) in eq. (74).
The average energy minimum and maximum occur at the eigenfunctions of the matrix eq. (75), which are equivalent to the adiabatically connected ones. We note that at these points the spectra are maximally separated from each other.

III. INFINITESIMAL PERTURBATION PROBLEM
So far we have a method of defining and calculating the Floquet eigenstates of an ideal Hamiltonian using the infinite time average energy, including the resonant region. However, within this framework, the eigenstates, propagators, etc. are sensitive to perturbations, particularly at and near resonance (similar to the perturbation problem of static degenerate systems). As a result, we cannot confidently model real systems because of the inevitable computational or modeling errors, even though the infinitesimally small perturbations should not have any measurable effects. Thus we have to refine the eigenstate definition to reliably model the real systems regardless of any infinitesimal perturbation.
In real experiments, one of the fundamental limitations of any time-dependent quantum system is the observable timescale, which we will use to quantitatively define and assess the accuracy of the Floquet eigenstates. Associated with this observable timescale is a minimum energy resolution ξ, which from here on is the desired degree of accuracy and for simplicity we imply its appropriate dimension transformations in upcoming equations.
To better understand the effect of the infinitesimal perturbation, we assume that the real/unperturbed Hamiltonian H 0 (t) are known, along with the eigenstates we wish to approximate 0 n ,¯ 0 n , u 0 n (t) .
A model/calculated Hamiltonian H(t) is riddled with infinitesimal undesired perturbations/errors v, and in principle, it is hard to distinguish between the perturbation and the actual components of the real Hamiltonian H 0 . So the only computationally accessible solutions would be these perturbed Floquet eigenstates { n ,¯ n , u n (t)}, corresponding to the following Floquet Schrödinger equation.
For simplicity, we will only consider the weak static perturbation v ≪ ξ ≪ H 0 , ω, although these methods can be generalized for other perturbation methods, e.g. high frequency expansions [10,37].
At the near resonance regime (mod( a − b , ω) < ξ), we can use the degenerate perturbation theory method to derive the perturbed eigenstates [38,42]. In this case the quasi-energy difference is bounded ( − 0 = O(v)), while the eigenfunctions u(t), their energy spectrum P ∞ (E) and infinite time average energy¯ ∞ can differ drastically. Even when the perturbation becomes infinitesimal v ≪ ξ, these definition difference can be beyond the acceptable resolution ξ. Thus the infinite time average energy method becomes unstable and impractical for realistic applications, e.g. defining a ground state.
On the other hand the perturbation has minimal observable effects if we limit ourselves to small enough timescales t < T U max . Using the interaction picture we can quantify the differences between the model/perturbed and real/unperturbed systems. The propagators U (t) are sufficient to characterize any observable difference, and the interaction picture propagator U I (t) extracts this difference by definition.
We can ignore the effects of non-resonant elements (mod( n − m , ω) > ξ), as the interaction picture propagator does not diverge from unity at any time-scales.
As for the resonant parts we can express it using the transformation matrix U(t): In this form we can see that U I (t) only starts to diverge from unity at timescales t ≳ (O(1 v ), since the transformation matrix is unitary at any give time UU † = 1.
We can define this as the timescale T U max , up to which any observable of the model and real system are indistinguishable.
Thus the model can approximate the real system for weak enough perturbations v ≪ ξ, and we want to find an eigenstate definition which can describe both systems within the timescale limitations t < T U max . Since the infinite time method goes beyond this limit, we have to choose a different quantum number to describe the observable effects.

A. Observed average energy
We define the observed average energy as the average energy expectation value up to a finite time T (eq. (39)), within reasonable experimental timescales T < T U max . Our goal is to find a definition of the average energy functionalĒ T which closely approximates the real/unperturbed oneĒ 0T within a predefined acceptable accuracy range ξ. With this we can trivially approximate the unperturbed eigenstates using a variational principle with tolerance ξ.
For simplicity we assume the real system does not have near resonance pairs, only exact resonances and far from resonance states. In this case the real system is well described by the unperturbed infinite time eigenstates E ∞ H 0 , and the unperturbed observed average energy quickly converges to the infinite time one.
Including near resonances in the real/unperturbed system, the observed average energy and its eigenstates can still be defined for the real system by reapplying the discussions in this section, in which case the observed unperturbed eigenstates E T H 0 would be the target we wish to derive instead of the infinite time solutions E ∞ H 0 . We explain more in appendix B, and here we only concentrate on deriving the infinite time solutions of the real/unperturbed system through the observed average energy.
Decomposing the definition in eq. (84) using the unperturbed eigenstates 0 n ,¯ 0∞ n , u 0 n (t) , we find two timescale boundaries TĒ min and TĒ max within which the observed average energy satisfies: The derivation of these boundaries and related discussions are presented in appendix A. Roughly these boundaries are related to the timescale where the unperturbed eigenstates become relevant T > TĒ min , and the timescale before which the infinitesimal perturbations v can be ignored T < T U max ≲ TĒ max . These boundaries only become relevant in further theoretical derivations based on the Floquet average energy, and in practice we are only concerned if a reasonable timescale can be defined within these boundaries. Since for any averaging time within this region, we get an equally good approximation to the exact average energy functional, we can choose an arbitrary timescale to define the appropriate observed average energy functional, and in most cases T ≈ 1 ξ would be a good choice for this.
With this definition of observed average energy alone we can derive a close approximation of the unperturbed eigenbasis E 0 H 0 from its variational principle around the initial wavefunctionĒ T As for the computationally accessible forms, we rely on the observation that the effective average energy¯ T [u] in the near resonant susbspace HT ξ retain the same properties of eqs. (57) and (58), up to the acceptable accuracy ξ.
[u] − < ξ δ [u] The variational procedure is thus analogous to the one presented in section II A, but in a more extended search space.
where eq. (95a) implies the quasi-energy resonance is not lifted if the difference is within the acceptable error ξ.
This accounts for the infinitesimal resonance lifting previously presented, so that we get consistent solutions near the real/unperturbed eigenstates 0 n ,¯ 0∞ n , u 0 n (t) for arbitrary weak interaction v ≪ ξ.
The equivalent Lagrange minimization method remains roughly the same, but with appropriate change of Lagrangian multiplier λ ξ to account for the finite resolution of δ : The same variational procedure are used when the real system has infinitesimal near resonance (appendix B). We also point out that, although in these derivations we have assumed we know the exact/unperturbed average energy and eigenstates we wish to approximate, in practice these are not necessary. We can discuss whether the observed average energies and their eigenstates of the model/perturbed system closely apporximate the real system by varying the timescale T or equivalently the acceptable accuracy ξ, and observe the stability of these solutions.
Thus, we have a more robust description of the Floquet eigenbasis E T H = n ,¯ T n , u n (t) , which has the same benefits of the definition in II A, but also is robust against infinitesimal perturbation. In this form unaccounted perturbations or numerical errors will not break the definition of the basis set.

B. Floquet-Ritz variation principle
The power of the Ritz variational principle is in the approximation of the ground state on a much smaller Hilbert subspace. For this we have to investigate how the ground state approximation changes as we increase the Hilbert space.
First it should be noted that in the full extended Hilbert space HT, both the infinite time and observed average energies are lower bounded by those of the ground states,¯ ∞ 0 or¯ T 0 , respectively. Thus at the limit of the Hilbert space expansion, we have a Ritz-like variation principle.Ē For the latter variation to be applicable, we assume that the two lowest observed average energies are sufficiently well separated¯ T 1 −¯ T 0 ≫ ξ, so that the global minimum gives a good approximation of the ground state . In order to have a proper Ritz variational principle, we assume we have a well behaved expansion series of the Hilbert space HT i , e.g. one derived from a Davidson algorithm [43,44], and we can find a lower-bound in the expansion i that closely approximates the exact solution. In other words, for a given accuracy ξ, the ground state { 0 ,¯ 0 , u 0 } HTi , or equivalently the average energy func-tionalĒ[u] HTi , evaluated in the Hilbert subspace HT i are within the acceptable accuracy ξ of their exact counterparts, and subsequent expansions will not improve the accuracy more than that.

=Ē[u]
HT In the simplest case we can model such a Hilbert space expansion with a monotonic decrease in the coupling v i between two Hilbert subspaces, HT i and its complement HT i⊥ , corresponding to the Hamiltonians H i and H i⊥ respectively. The Hilbert subspace HT i is where the current ground state approximation resides, and increases after each iteration. The full Floquet Hamiltonian can thus be decomposed as follows at any given step i in the expansion.
Whether a given algorithm does yield such a property is still up for debate. For now we are only concerned if such an algorithm can yield a Ritz-like variational principle and what ground state definition it should follow.
First let's consider the behavior at a weak coupling threshold v ic ≲ ξ. In this case we can use the previous arguments in section III A to find that the observed average energy ground state of the full Hilbert space T 0 ,¯ T 0 , u T 0 HT ∈ E ξ H is closely approximated by the decoupled solution of eq. (102): which in this case is the observed ground state HTi c of the Hamiltonian H ic (t) in the limited Hilbert space HT ic , assuming the correct subspace is chosen. Similarly all of the expansion points above this threshold i > i c have the same property, as long as the expansion series satisfies eq. (103). We can simplify this statement using the observed average energy functional evaluated in the truncated Hilbert space HT i , compared to the exact one evaluated in the full Hilbert space HT: The same is not true for the infinite time ground state definition, where we can find cases where the asymptotic limit does not correspond to the full Hilbert space solution, provided that such a limit can even be found.
(106) Therefore a Floquet Ritz variational principle can be formulated based on the observed ground state definition. The condition for the convergence and the Ritz inequality are as follows: The Floquet Ritz method has similar properties to the static method. If the initial guess is not good enough we risk convergence to an excited state rather than the desired ground state (eq. (108)).
There is however a caveat to this method. The functional used in evaluating the average energyĒ T HT i ′ at intermediate steps, namely the effective average energy functional¯ T [u] HT i ′ , or the propagator U (t) HT i ′ , do not necessarily give a good approximation to the average energy functional in the full Hilbert spaceĒ T HT . It is thus possible to find "average energies" well below the observed ground state average energy in the full Hilbert space.
This does not contradict the Ritz variational principle in eq. (108), since the solution at these points do not satisfy eq. (107) and are not convergent yet. However this can still create problems depending on the algorithm used to approach the convergence. E.g. if we select the ground state at each step to determine subsequent Hilbert spaces HT i+1 , it is possible to flip-flop between approximations of various different eigenstates ( T n ,¯ T n , u T n ), and even to converge to an excited state despite starting from a good initial guess.
Nevertheless, in principle a Floquet-Ritz variational principle is possible on limited Hilbert spaces HT i , as long as we use the observed average energy labeling, either explicitly or implicitly. This is the main point we want to state regarding a possible Floquet-Ritz theory.

C. Perturbed two-level system
A minimal example including the perturbation effect presented above is the two-level system in eq. (68) perturbed by a weak static coupling. Physically this can be regarded as a stray static electric field contamination or computationally it can be a numerical error.
For consistency we consider the unperturbed Hamiltonian (eq. (68)) to be the real system and we project the perturbed Floquet Hamiltonian on the unperturbed Floquet eigenfunctions (eqs. (69a) and (69b)): We will focus on the first resonance point Ω = ω. Here the driving frequency ω is much lager than the static components ( 0 , so we can approximate the Hamiltonian by ignoring the time-dependent components in the projected Floquet Hamiltonian (eq. (111)). The exact perturbed Floquet eigenstates can then be approximated to:Ĥ These eigenstates correspond to the infinite time eigenstates with the average energies almost coinciding with each other.¯ ∞ ± ≈ ±v Here we see the limitations of the infinite time average energy method. As we approach the limit v → 0, the infinite time average energies difference of the model/perturbed systems vanishes¯ ∞ + −¯ ∞ − → 0, while the real/unperturbed solution at v = 0 does not (¯ 0 + −¯ 0 − = ω 0 ≠ 0). Such mismatch in energies does not occur in the static case and is specific to the resonance of Floquet systems. We also see that the eigenstates (eq. (115)) do not change as we approach the limit v → 0, even as they are defined by the infinite time average energy.
We plot the infinite time Floquet eigenstates and average energies in fig. 3 with the perturbation exaggerated to v = 10 −3 ω 0 for clarity. This system is simple enough that we can compute the exact propagators.  As shown in fig. 3a, the infinite time eigenstate solutions of the perturbed system are closely approximating the unperturbed eigensolutions ( fig. 2a) at all the non-resonant points, with the exception of a trivial label swap caused by the difference in the adiabatic continuations. However at the resonance point these solutions differ drastically, which is best seen in the resonance energy spectrum ( fig. 3b) plotted against the initial wavefunction with the definition given in the previous example (eq. (76)). Here we notice that the range of the average energy variation is much smaller than that of the unperturbed case ( fig.2b). Also, the eigenstates of the unperturbed and perturbed systems correspond to different values of θ (dashed and continuous gray line). Thus we confirm that the labeling of the Floquet eigenstates by the infinite time average energy can be inconsistent with the inclusion of infinitesimal perturbations.
Having shown the limitations of the infinite time approach, now we compare it to the observed average energy method. First we want to confirm that our original premise for deriving the observed eigenstates is valid, i.e. any observable of the perturbed and unperturbed system can be within an acceptable error at some timescale. Since the exact propagators are known at the first resonance point (eqs. (120) and (121)), we can directly find the timescales T where propagators in the real and model systems are equivalent, for arbitrarily chosen accuracy ξ. For our purposes it is sufficient to expand the model propagator U (t) up to first order in v and derive the relevant timescale T from there.
As long as the observation timescale and the timescale for determining the average energy are within this range, we can conclude that the model/perturbed system does indeed well approximate the real/unperturbed system. Here we have used the dimension transformation of ξ with the typical energy scale as ω 0 . For the two-level system this rough approximation is sufficient, especially since the most interesting physical system is where all of the interactions are of the same order O(ω) = O(ω 0 ) = O(V ) and the time-periodic interactions are most prominent.
Within these timescales T and acceptable accuracies ξ, we can derive the observed average energy functional E T and its resulting eigenstates directly from eqs. (39) and (89). But first we look at the observed average energies boundaries (TĒ min and TĒ max ), which we have proposed to determine the region where the model eigenstates approximate the real ones. Since we know the exact Hamiltonians, we can directly derive them from eqs. (A3) and (A6): We see here a good agreement between the two higher bounds (TĒ max and T U max ) considering that O(ω) = O(ω 0 ) = O(V ), in accord with our derivation. The lower boundary in this system, assures that the observed average energy functional approximates the infinite time average energy of the unperturbed system, and not another functional with different critical points.
In the limit of v → 0 we find the higher bounds diverge to infinity, suggesting that at arbitrary timescales, we would not be able to find an observable difference between the model and real system, perfectly in accord with the physical intuition. Comparing the lower and higher bounds (TĒ min and TĒ max ), we find that for small perturbations v we can always find an acceptable error ξ, up to which the model Hamiltonian approximates the real one, although this is often too large for practical applications.
For this example system we choose the acceptable accuracy to be ξ = 10 −2 ω 0 so that we can define an observable timescale between the boundaries. The effect of these boundaries are best seen in fig. 4, where the energy functionalĒ T [Ψ(0)] is evaluated for varying timescales resulting into three distinct regions. These regions represent the timescale where the energy functionalĒ T [Ψ(0)] approximates: (I) the instantaneous energy; (II) observed/unperturbed infinite time average energy; (III) perturbed infinite time average energy. We can see the accuracy of these descriptions in fig. 4b, as well as the deviation from the acceptable regime outside the boundaries (with the exception of region (III), which has not yet converged due to numerical limitations).
Next we look at the derivation of the observed average energies and the eigenstates. Using the unperturbed Floquet basis, we can compute the exact quasi-energy, observed average energy and effective average energy at the near resonance for different Floquet functions u θ (t), other than the eigenstates. For sake of simplifying the equation, the Floquet functions are constructed by mixing the unperturbed Floquet eigenfunctions, which we wish to derive. u θ (t)⟩ = cos θ u 0 We see that the effective average energy¯ T is indeed a close approximation to unperturbed average energyĒ 0 up to the accuracy of ξ. Since the unperturbed average energies are well separated¯ 0 + −¯ 0 0 ≫ ξ, and the quasienergies are relatively stationary (eq. (130)), we can get a close approximation to the unperturbed basis set from the variation of eq. (131). We plot this observed average energy in fig. 5, where we use a numerical approximation of the energy spectrum (eq. (41)) with T = 1 ξ to simulate the limited observable accuracy ξ and differentiate from the exact energy spectrum in fig. 3b.
Including the finite accuracy effect on the energy spectrum, we find it to closely resemble the unperturbed spectrum fig. 2b. The stationary point of the effective average energy (eq. (131)) as well as that of the spectrum average taken from fig. 5, closely agree with the unperturbed solution. We depict the acceptable accuracy range calculated from eq. (134), as the vertical gray areas, with the inner most area corresponding to an error of v instead of ξ. The error in the eigenstate wavefunctions is normalized by the average energy difference as follows: (134) Finally we show the possibility of the Ritz variation principle by adapting Davidson algorithm to the Floquet system, similar to Sahoo et. al. [44]. This example system is simple enough that we can use this formulation directly, while more complicated systems would require a more refined method to take into account the average energy variation as well. Starting from the initial guess of the undriven ground state u i0 , we calculate the ground states u in as defined in this paper and the residue vector r in at each step i n , which represent the current approximation of the ground state and the convergence condition, respectively. We have executed the procedure equivalent to [44] and we will not go into detail about this procedure. We only note that we've changed the algorithm of selecting the ground state and we focus on how the algorithm converges and the corresponding solution.
The results converge to an error of ξ within two steps, without having to create the infinite matrix H F . The convergent solution is also well within the accepted accuracy to the observed ground state. If we disregard the acceptable accuracy ξ and instead choose an accuracy ≪ v, we would eventually converge to the exact infinite time ground state.

IV. CONCLUSION
In this paper we have proposed a robust method of defining the Floquet eigenstate in an ordered fashion using the average energy. The main differences from the previous methods are: we are able to uniquely order the eigenstates independent of the quasi-energy shift; we can uniquely define the eigenstates near resonance where it would otherwise be ambiguous; the method is robust against infinitesimal perturbations; we can systematically cut off the Hilbert space and retain the accuracy within a reasonable timescale. Based on this definition we derived variational methods of approximating the eigenstates, which in principle are more computationally efficient. For this method additional consideration has to be made to the physical timescale that we wish to investigate.
The physical significance of the average energy was not presented in this current work. Some intuition can be found by evaluating the steady state of the open-quantum system. We conjecture that for reasonable systems, this choice of basis set gives a good approximation of the exact steady state, i.e. the density matrix is close to diagonal and occupied by a few low lying states in this eigenstate representation, even as the Hilbert space is truncated.
Having a variational principle for deriving the ground state, various method such as Floquet Hartree-Fock could be adapted to approximate the ground state. Since the ground state is generally not sufficient to describe a physical steady-state, additional excited states are required, which can be derived in a similar variational manner, with similar Excited Floquet Hartree-Fock. With sufficient calculations (estimated from the resulting energy spectrum), the computation of the physical Floquet steady-state is relatively trivial, at which point we have a good approximation of a real periodically driven system at a long enough time-scale to have equilibrated with the environment.
Ĥ 0 (t),v I (τ ) dτ u n (t)⟩ dt + O v 2 . (A1c) The first term (eq. (A1a)) quickly converges to the unperturbed average energyĒ 0 [Ψ] (regardless of its definition) within a few cycles of the driving T > O(T ). This timescale is much lower than the timescales we will be discussing and will be ignored. Thus the goal is to find the timescale T where the remaining terms (eqs. (A1b) and (A1c)) vanish for arbitrary initial wavefunction Ψ(0), or in this experssion, arbitrary coefficients C 0 n .
From the second term (eq. (A1b)), we obtain the lower boundary TĒ min which is independent of the perturbation v of the model H(t). Depending on which definition of the average energy and eigenbasis we wish to approximate in the unperturbed system (appendix B), the lower boundary TĒ min can change. For simplicity we assume the unperturbed system has no finite near resonances, so that we approximate the infinite time solutions of the unper- Even when the unperturbed Hamiltonian has exact resonances (∃ ω 0 mnl = 0), the lower boundary remains finite due to the diagonalized definition of the infinite time average energy (eq. (59)) and the infinite time eigenstate definition. Similarly in the case of infinitesimally small, but finite near resonance (∃ ω 0 mnl ≪ ξ), the lower boundary can be within the acceptable timescale TĒ min < T ∼ 1 ξ, depending on whether we treat these states as resonant states or not (appendix B). We can further lower this boundary if we limit the average energy functional approximation to a smaller Hilbert space, e.g. only approximating the functional near the ground state. A rough physical understanding of the boundary TĒ min is the timescale from where the weakest significant interaction or avoided crossing can be resolved up to the accuracy ξ.
The higher boundary TĒ max is derived from the remaining terms in eq. (A1c), and are specific to each perturbation v. We can safely ignore the perturbation effects on non-resonant states and only consider the effects of resonant and near-resonant states (mod( a − b , ω) < ξ). Calculating the condition for eq. (85) to be satisfied we get a form of the higher boundary as follows: This higher boundary TĒ max closely approximates T U max which defines up to which timescale the propagator and any observable, are closely approximated in the model/perturbed system and the real/unperturbed one. So a rough physical intuition of this higher boundary is the timescale up to which the effects of the infinitesimal resonance lifting can be ignored. Depending on what we define to be the perturbation v, and what the unperturbed eigenstates are, the lower and higher boundaries could cross (TĒ max < TĒ min ), in which case we would not be able to find an average energy definition that satisfies Ē T −Ē 0T < ξ, and the Floquet eigenstates of the two systems could differ significantly for any observed average energy definition. Different choices of perturbation v or unperturbed eigenstate defintion E T H 0 could recover this condition (e.g. appendix B). Otherwise it could simply be that the perturbation v is not weak enough to be ignored and the model does not closely approximate the real system.