Polynomial filter diagonalization of large Floquet unitary operators

- The manuscript proposes a new diagonalization scheme for Floquet operators

- The proposed method is convincingly demonstrated to have an advantage over existing exact diagonalization schemes, both in terms of time and memory

- No new physics, no results regarding the physical properties of the studied system

- The complexity of the algorithm is not discussed

The manuscript "Polynomial filter diagonalization of large Floquet unitary operators" introduces a method of exact diagonalization of unitary Floquet operators. The method relies on an observation that one can formulate an efficient multiplication of vectors by Floquet operators for physically interesting systems. A combination of a polynomial spectral transformation and Arnoldi algorithm that utilizes the efficient matrix-vector product is then used to calculate eigenvalues of the Floquet operator.
The manuscript shows that the proposed method indeed provides eigenvectors of the considered Floquet operator and convincingly argues that the proposed method significantly outperforms both full exact diagonalization and shift and invert method. The Floquet systems play a major role in our understanding of non-equilibrium dynamics. Hence, the paper fulfils the two expectations of submissions to SciPost Physics: i) it details an important computational discovery, ii) opens a new pathway in an existing research direction. As such, it should be published in SciPost Physics. However, I recommend a revision of the manuscript according to the remarks specified below.

1. The paper presents no results on the considered physical system. In fact it is even not clear what precisely is the considered system: are the two-site unitaries drawn with the Haar measure? While the aim of the paper is to introduce the method, including results showing some basic physical properties of the calculated eigenstates (for instance the system size dependence of entanglement entropy) would widen the scope of the paper.

2. From the present manuscript it is not clear what is the complexity of the proposed method. A more thorough analysis of the scaling of the various steps of the algorithm with matrix size should be included. Those scalings could then be compared with results shown in Fig. 3.

3. Could the efficiency of the method be improved by varying the coefficients of the polynomial in the spectral transformation (6)?

4. Can investigations of Floquet systems of size L=20 bring new physics as compared to systems of size L=15 accessible to full exact diagonalization? I believe that referencing the recent debate about many-body localization transition [Phys. Rev. E 102, 062144 (2020), EPL 128, 67003 (2020), Phys. Rev. Lett. 124, 186601 (2020), Annals of Physics 427, 168415 (2021), etc.] would illustrate well the importance of finite size effects in non-equilibrium phenomena.

1-Useful extension of standard ED methods for (non-)Hermitian matrices to unitary matrices, which can (as argued) be immediately relevant for Floquet systems.

2-The presented methods are clearly correct and their advantages over other methods are supported by strong numerical evidence.

3-The main ideas are clearly presented.

1-The text may be a bit technical and to concise for the general reader.

2-No data is provided. No code is provided – not even pseudo-code.

The paper introduces an adaptation of general ED methods to unitary matrices. In particular, a polynomial function is introduced that maps a complex number on the unit circle to a general complex number, whose norm depends on how close the initial number is to a specified target (on the unit circle). The polynomial in question is a finite geometric sum, which could also be understood as a sort of truncated delta function around the target, hence the close values are strongly enhanced. This function allows to compute the eigenvalues and eigenvectors of a unitary matrix close to a set target using only matrix-vector multiplications (due to the polynomial form). Since relevant unitary matrices will be dense in general, this can be far more efficient than shift-invert methods as argued and demonstrated in the paper. In between, it is also described how an efficient matrix-vector multiplication of a unitary based on a “brickwork” circuit can be performed.

After introducing these concepts, the paper mostly presents numerical evidence of the superiority of the polynomial method as opposed to the shift-invert method and full diagonalization. For this, 50 eigenpairs close to 1 are computed for different sizes and with different parameters. The evidence shows clearly that there is an advantage in terms of runtime and memory requirements. Given that polynomial methods are well established also for general systems, this is really not too surprising, hence there is no reason to doubt the reported findings.

Therefore, the question of publication boils down to the question of whether the results offer a significant numerical advantage and whether they are interesting to the general reader. One could argue that it is more or less a straightforward application of known polynomial methods to unitary matrices – this is also the feeling I got when first reading the text. However, given that unitary matrices are very important, not only in Floquet but also in Quantum Physics generally, and given that work on similar methods for Hermitian systems has been published and is widely known, I believe that the paper should be published in SciPost Physics.

However, I feel that some changes (apart from using the SciPost layout) should be made as listed below in the requested changes section. A change that I do not request, but would suggest the author to consider, would be to add some discussion of how computing eigenpairs for larger sizes could be useful. I think that already a short discussion of concepts like level statistics and their meaning in the study of phase transitions or the study of heating in Floquet systems could increase the readability and usefulness for the general reader significantly.

Acceptance criteria (https://scipost.org/SciPostPhys/about#criteria):
In my opinion (see the discussion earlier) the paper meets the SciPost expectation

“3. Open a new pathway in an existing or a new research direction, with clear potential for multipronged follow-up work;”.
.

It also meets the general acceptance criteria except for

“5. Provide (directly in appendices, or via links to external repositories) all reproducibility-enabling resources: explicit details of experimental protocols, datasets and processing methods, processed data and code snippets used to produce figures, etc.;”

While I am confident that the results can be reproduced in principle given the paper, no code for the computations or figures and no data is provided.

0-(formal) Provide the necessary code and data to fulfill the reproducibility requirement. Use the SciPost layout.

1-The expressions in eq. 6 and 7 should be discussed in a bit more detail – at least in an Appendix. For instance the functional form $g_{k}(e^{i \theta}) $should be provided, such that one sees the equation of the “apple”. Also the derivation of the mapping of an arc to the “apple” should be be given.

2-Concerning Fig. 2: I find the color scheme here rather confusing. Not only is the scale different for every plot, so that one can not compare the sizes based on color (to be fair this is done in Table 1 and Fig. 3), the relative ranges are also completely different. For instance in the top figure the numbers range from (presumably) 1 to around 11 (?), so about a factor of 10, while in the bottom they range from 600 (?) to 1250 (?), so about a factor of 2. As a result there seems to be a lot of variation in the color of of the bottom plot, while the top plot seems almost fully blueish. Thus the plots look visually different but there is no real meaning behind it, or it is not explained well – in any case there should be some changes in the figure or discussion thereof.

3-A clarification as to the exact definition of the random circuit might be useful: is every unitary drawn from a circular ensemble or is there a finite set of unitaries and they are chosen at random? Is there an averaging involved in the benchmarks?

4-Is the “fast” matrix-vector multiplication used in all benchmarks (if applicable)?

5-Curiosity: Do you know why the residuals for full diagonalization in Tab. 1 are so bad (compared to the other methods)?

6-Is there an insight as to why one gains around four sites, or is this a purely empirical observation?

7-In the last sentence of the conclusion possible optimizations are hinted at, perhaps at least a sketch of how these could be achieved could be given here.

8-Finally, the references should be extended. While the bibliography is complete in the sense that all statements are referenced, it seems that there are extensive references for some topics like discrete time crystals (13 out of 40 – by the way references 21-23 are arguably not related to time crystals but rather focus on prethermalization and energy absorption in generic systems),
while for other topics only a single reference is used with more well known works being available. As an example, [1] is used as a reference for heating to an infinite temperature state in generic Floquet systems, but 10.1103/PhysRevX.4.041048 is a well known source for the same statement. Other works I can see fitting (off the top of my head) are: 10.1103/PhysRevLett.116.120401, 10.1080/00018732.2016.1198134, 10.1103/RevModPhys.89.011004, 10.1080/00018732.2015.1055918, 10.1103/RevModPhys.91.021001.