We sincerely thank the editor and referees for considering our manuscript. We found that the detailed and in-depth referee reports were very helpful in improving our manuscript. We have addressed all the comments and revised our paper accordingly (see authors' reply to each referee report).

Reply to Report 1

We deeply thank the referee for taking the time to read our manuscript and to give the in-depth review.

We have addressed all the given remarks and requested changes in the revised version. Below we provide the point-by-point responses to the referee comments.

Response to Referee's remarks

a) The high-expansion analysis of Sec. 5.2 is performed on the hamiltonian evolving $W_k$ (Majorana fermions), and it is not explained how to relate it to a high-frequency expansion for the spin density matrix.

We regret that the previous version did not explain the relationship between the spins and Majorana fermions. Actually, it is difficult to give an explicit connection between the density matrices in the spin and Majorana-fermion representations due to the nonlinear nature of the transformation. However, in the Majorana-fermion representation, we can address the energy dispersion relation for the elementary spinon excitations.

In the revised version, we have expanded the model description in Sec. 5.2, consolidating the connection between the spin and Majorana-fermion representations.

b) The phase diagram in Fig. 1 changes drastically from 0th to 2nd order expansion. Is an expansion to 2nd order sufficient to get the full picture of the physics involved?

This is an important question, but we could not perform the 4th order perturbation theory. What we can say at present is that, in the vicinity of $T=0$, the phase diagram resembles that obtained by a different approach in Ref. [83]. Once the 4th order HF expansion formulas are available, we would like to investigate this point.

In the revised version, we have added, for clarity, a sentence "We leave it an open question to examine how the phase diagram changes for higher-order calculations, which require considerable effort."

c) It is not immediately clear to me why Eq. (50) implies property (ii).

We meant to point that $\mathcal{D}$ in Eq. (50) is exactly the same as the original Eq. (47). This is a special property of this model; The dissipator in $\mathcal{L}_\mathrm{eff}$, in general, proliferates over the entire chain. Given that the dissipation is confined on the edges even under drives, the dissipation is expected not to influence the bulk property much.

We would not state that Eq.(50) is proof of the property (ii). Instead, we think it provides some supporting information for why the property was seen in Ref. [83].

We have remarked the above points in the revised version.

d) In Eqs. (51)-(55), the dissipator $\mathcal{D}$ commutes with $\mathcal{S}_z$ independently of the sites on which the jump operators act. Therefore, while it is true that the drive prevents an edge dissipator from propagating into the bulk, Eq. (50) would still hold even for a dissipator acting on the entire chain. Would property (ii) still be true in that case or would the phase diagram change?

The referee is correct. The property (ii) holds even if the dissipator acts also on the entire chain. Namely, the dissipator does not propagate in the superoperator commutators.

Meanwhile, in such a case, it becomes difficult to analyze the NESS quantitatively using the Majorana representation. This is because $\sigma_j^\pm$ entail string operators nontrivially acting on $i(<j)$ except $j=1$ and $N$. We think it is still an open question how the phase diagram changes in such cases.

In the revised version, we have added a paragraph to remark this.

e) Would $\gamma_{\alpha\beta}(\epsilon_m-\epsilon_n-k\omega)$ in Eq. (82) have a simple representation in terms of Feynmann diagrams of the dissipative process?

Yes. We have added a new figure (Figure 6 in the revised version) for illustration.

Response to Requested changes

1 - Please see our response to Remark d).

2 - Please see our response to Remark c).

3 - Please see our response to Remark a).

4 - We corrected the character code in the figure files.

5 - We reorganized and enlarged the figures.

6 - Corrected. We appreciate the careful reading.

7 - Specified. We totally agree this specification improves accessibility.

8 - Please see our responses to Remarks b) and e).

Reply to Report 2

We gratefully thank the referee for taking the time to carefully read our long manuscript and to give valuable comments, which we found very helpful for revision. Below we respond to each of the given comments.

(1) I think in the abstract and the introduction it does not become sufficiently clear that the proposed high-frequency (HF) expansion does actually not directly provide the non-equilibrium steady state (NESS) of the system, but that it is rather employed for approximating the micromotion operator and the effective Liouvillian, from which the NESS still has to be computed numerically.

The referee is correct. We have remarked that numerical computations are necessary in the abstract and introduction in the revised version.

(2) In the paragraph following the one containing equation (25), it is stated that in case the one-cycle evolution operator has a negative real eigenvalue, both the Floquet Liouvillian L_F as well as the effective Liouvillian L_eff are not of Lindblad form. However, one should be aware that there are also situations, where L_F is not of Lindblad form, despite the fact that the one-cycle evolution operator has no negative real eigenvalue (see Refs. 64 and 65). In this case, it can happen that L_eff is of Lindblad form, while L_F is not. An example for such a situation is described in a recent preprint [arXiv:2107.1005] that appeared parallel with the present manuscript and that equally uses van Vleck's high-frequency expansions.

We thank the referee for pointing out the recent preprint (we assume it is arXiv:2107.10054). In the revised version, we have remarked that such an example exists with citing the preprint.

(3) The proof of Lemma 1 is rather brief. I think it would be worth of giving a more detailed explanation either in the main text or an appendix.

We appreciate this suggestion. We have extended the proof in the main text.

(4) Regarding the second example of a three level system: I am not sure whether the bath model used is meaningful. Namely, the dissipator is derived microscopically using the energy eigenstates of the undriven system. However, at the same time the Hamiltonian of the system has also a time-periodic term that should be taken into account when deriving the dissipator (as it is done in a later section). The authors provide the high-frequency expansion of both the Hamiltonian and the micromotion operator. I think it would be interesting to see also the high-frequency corrections to the dissipator.

Physically, the referee's claim is correct. The bath model in Sec. 5.3 is an approximation and a more precise treatment is given in Sec. 6.2, where we incorporate the high-frequency corrections in the dissipator.

We note however that the approximate bath model in Sec. 5.3 is useful in understanding the influence of the dissipator correction. In fact, such time-independent dissipation is often used in a broad class of physical systems including diamond NV centers, quantum optics systems, magnetic resonances, semiconductors, etc. The time-independent dissipator is known to work well at least at a semi-quantitative level in many systems. Therefore, it should be important to validate this phenomenological approximation anyhow. Our comparison between Secs. 5.3 and 6.2 provide theoretical support for the approximation (see e.g. the similarity between Figs. 3 and 6).

In the revised version, we have added a paragraph emphasizing that the bath model in Sec. 5.2 is widely used in the NV-center community and turns out to be a good approximation comparing Secs 5.2 and 6.3. We have also added two sentences in Sec. 6.3 to remark the similarity between Figs. 3 and 6.

Minor point: One might also want to write the driving term as a numerated equation like all the other terms appearing in the Liouvillian.

We agree and have displayed the equation.

(5) At the end of section 5 it is stated that the approach does not require numerical time integration, but uses only linear algebra. This is true. But one at the same time one still needs to compute an eigenstate of the effective Liouvillian, which is not necessarily a simple task. I think it would be good to include a discussion, in how far computing the steady state from L_eff can be easier than computing the steady state from time evolution. A place for such a discussion might also be the introduction.

We regret we did not give the important discussion on the computational merits. While both approaches are feasible in principle, the approach with $\mathcal{L}_\mathrm{eff}$ is easier in two ways.

The first point is the number of matrix-vector multiplications. The numerical time integration approach requires numerous matrix-vector multiplications before $\rho(t)$ reaches the NESS. Here the matrix (vector) means the Lindbladian $\mathcal{L}t$ as a matrix and $\rho(t)$ as a vector. The number of multiplication is roughly estimated by $\sim1/(\gamma_0 \Delta t)$, where $\gamma_0$ is the dissipation strength and $\Delta t$ is the time stepping. To obtain high-accuracy results, we have to decrease $\Delta t$ and need numerous multiplications. On the other hand, in the HF-expansion approach, we are to find the eigenstate $\eta$ of $\mathcal{L}_\mathrm{eff}$ with the largest real part. For this purpose, we can use the famous Lanczos algorithm, where the required number of matrix-vector multiplications are greatly suppressed. So we can reach the NESS more efficiently with the HF-expansion approach within the HF approximation.

The second point is the efficiency in obtaining expectation values of observables in the NESS. To obtain $A(t)=\mathrm{tr}[\rho_\mathrm{ness}(t)A]$ for an $A$ in the time evolution approach, one needs to evolve the NESS density matrix over one driving period. Here again, one needs many matrix-vector multiplications. On the other hand, in the HF-expansion approach, we need a smaller number of them. This is because, thanks to $\rho_\mathrm{ness}(t) =e^{\mathcal{G}t}\eta$ and the expressions of $\mathcal{G}t^{(k)}$, we know the exact $t$-dependence in the form of $e^{-im\omega t}$. Thus, once we perform a few matrix-vector multiplications and obtain $\rho_m$ in $\rho\mathrm{ness}$ as $\sum_m \rho_m e^{-im\omega t}$, we can immediately calculate $A(t)$ at an arbitrary $t$ in the driving period.

We have added these discussions at the end of Sec. 5 and in Sec. 4.

(6) I have a few comments regarding section (6), where the authors discuss the Floquet-Lindblad equation, as it is obtained for a driven system coupled to a thermal bath by employing the Floquet-Born-Markov approximation in combination with a rotating-wave approximation. The section is titled “time-dependent dissipator for weak thermal contact”. However, as the authors write, this problem can always be mapped to a time-independent problem in the interaction picture. And actually, the authors always treat this time-independent problem. Thus, emphasizing the time-dependence of the dissipation in the section title might be slightly misleading.

We agree with the referee. In the revised version, we have changed the section title to "Microscopically-derived dissipators by weak thermal contact." Correspondingly, we have also changed the title of Sec.5 by adding "Phenomenological."

The approach pursued in this section is slightly different from the one described in the previous sections. Namely, the authors employ the high-frequency expansion to the system-Hamiltonian without dissipation, in order to compute the Floquet states, which are then used in a second step to derive the Floquet-Lindblad equation, which is then mapped to a time-independent problem. This is an interesting and valid approach, however, as far as I see, it is different from what has been described in chapter 3. I think the authors should either point out that now a different approach is derived or they should explain in more detail how this approach is related to the previous one.

We thank the referee for this clarification. We regret that our previous version was not clear enough.

As the referee pointed out, we here derive a different approach making use of the time-independent frame. This approach only applies to Floquet-Lindblad equations (FLEs) derived from RWA but is particularly useful for this class of problems.

However, in more general FLEs, this approach does not work, and we need to consider using the general framework described in chapter 3. While we did not, we could have applied the framework to this section as well without using the time-independent frame.

In the revised version, we have extensively rewritten the beginning of Sec. 6 to stress the aim of this section.

Appendix D and section 6.1.1 almost entirely correspond to the results presented already in Ref. 41 and other papers. While this is pointed out at the beginning, during section 6.1.1. the authors always refer to appendix D, when presenting results, rather than referring to the original papers. Since the manuscript is already very long, the authors might consider to shift the long section 6.1.1.. into appendix D, since it mainly reviews known results.

We agree with the referee and have shortened Sec. 6.1.1 by moving derivations to Appendix D while keeping definitions. We have also cited Ref. 41 appropriately in the main text.

We would like to stress that our derivation is a generalization to allow degeneracies in quasienergy differences and is completely given in modern notations (earlier studies omit derivations). So we believe that Appendix D should deserve publication.

(7) The authors might consider removing section 6.3, containing the example of the inverse Faraday effect in an open driven Heisenberg chain, from the present manuscript and transforming it into an independent second publication. Namely, while most of the paper is focusing on a method, the high-frequency approximation for open quantum Floquet systems, this section puts a strong focuses on the physics of this particular effect. This is interesting, but after examples 1, 2 and 3 not necessarily required to illustrate the method. Moreover, in this way the paper becomes extremely long. However, this point (7) is only a suggestion.

We thank the referee for this suggestion. Although we have considered this possibility, we come to the conclusion that we keep section 6.3 as is because our emphasis is also on the test of the Floquet-Gibbs state as well as the specific phenomenon. We have learned that too many examples lengthen the paper, which we shall keep in our minds for future publications. We thank the referee for the advice.

All changes made in the revised version are highlighted in blue.