Dear Claudio Attaccalite,

thank you very much for managing the review process for our manuscript. We are happy that both reviewers are positive about our work.

Below, we address the reviewers' comments and remarks on a point-by-point basis and detail how we modified the manuscript accordingly.

Regarding the code, we followed all suggestions and now include explicit examples, auto-tests, improved the code documentation, and added a benchmark comparison.

The compilation error found by one of the reviewers occurred on a type of computer that we are not supporting. Indeed, the new ARM ("M1"/"M2") chips used by Apple do not, at this time, support quad-precision for Fortran. We will monitor the situation and add support once this becomes possible. The reason for the HDF5 error found by the same reviewer was the use of the outdated version 1.10.5 of the HDF5 library (we require >=1.12.1).

We hope that our code and manuscript are now to your satisfaction and that of the reviewers.

Thank you for your time and consideration M. Pickem, E. Maggio, J.M. Tomczak

REPORT 1:

"The work by M. Pickem, E. Maggio, and J. M. Tomczak presents a new code to solve the Kubo formula in a semi-analytical way. The manuscript is well written and the code needed in the community as it connects the low temperature quantum incoherence scattering to the high temperature resistive one. It certainly deserves publication."

We thank the referee for the positive evaluation of the merit of our code and methodology.

"However there are quite some points that I think the authors should address first: * The authors discussed and compare their approach to the BTE in the relaxation time approximation. However nowadays the de-facto calculations mode is the iterative BTE solution in software such as EPW, Perturbo, Elphbolt or Abinit (Refs. 13,14,17,19 - you might be missing PRB 102, 094308 (2020) for the Abinit). Some of the statements made by the authors do not hold with the iterative solution and need to be discussed."

We thank the referee to bring PRB 102, 094308 to our attention. We have added the publication to the list of references. Regarding the "iterative BTE solution", we would like to point out that our approach does not require solving self-consistent equations. Indeed, iterative BTE solvers compute effective scattering rates and resulting changes in distribution functions. Our methodology, instead, takes scattering rates as a given input. In other words, the scattering problem, e.g., electron-phonon scattering, must have already been solved elsewhere, before it is post-processed by our code. Therefore, a comparison of our methodology to BTE in the relaxation-time approximation is more natural.

"* The LinReTraCe software rely on a linearlization of the electron selfenergy. However such linerization (e.g. Eq. 22) is know to be numerically unstable for all but the band edges (e.g. interband transition with low lying energy minima with similar energy). A Dyson-Migdal solution might be more appropriate. How do the authors addressed this instability ?"

Contrary to Boltzmann solvers, we do not deal with differential equations. We are linearizing a given self-energy around at a chosen energy to extract and post-process energy shifts, quasi-particle weights and scattering rates. The linearization could be around the Fermi level for a Fermi liquid, or the bare conduction (valence) band minimum (maximum) for a semiconductor. While the linearization is an approximation, it does not lead to any instabilities in the algorithm.

"* At the end of Sec. 2, the authors mentioned that the generalization of the band-curvature to the Wannier basis has not yet been derived. What about Eq. 60 in Sec. 3 ?"

In Eq. 60 we indeed compute band-curvatures in the Wannier basis before transforming them into the Kohn-Sham basis. Like the velocities, curvatures are not gauge-covariant, so it matters in which basis to perform the momentum derivatives. What (to the best of our knowledge) has not been properly derived by anyone are expressions for, say, the Hall conductivity sigma_{xy}^B starting from the continuum. Indeed the standard derivation for sigma_{xy}^B starts by coupling the vector potential via the Peierls substitution. This leads to group velocities of the lattice (k-derivative of the dispersion). Curvatures (second derivative w.r.t. k) then emerge, technically, through a partial integration. The standard conductivity sigma_{xx}, instead, has been derived, both, from the tight-binding/lattice model (Peierls substitution approach --> group velocities) and the continuum (minimal coupling --> dipole matrix elements), see, e.g., PRB 80, 085117 (2009) for details. It will be interesting to derive sigma_{xy}^B starting by including the vector potential via the field-theoretical minimal coupling prescription.

We have reformulated the respective paragraph along the above lines.

"* In Fig. 1, there might be a typo in two blue boxes with "[,""

The typing was intentional to indicate that the third list entry is optional.

"* In the footnote 17, page 30, the authors indicate that Boltzmann codes typically use fixed chemical potential or would exhibit "massive" numerical instabilities in gapped system. I disagree. As far as I know, all the BTE codes mentioned above compute the chemical potential using, e.g., a bisection method. It is very stable for the relevant temperature range (e.g. above 50 K)."

While one of the merits of our code is indeed a reliable determination of the chemical potential, we might have oversold the comparison to some of the Boltzmann codes. We have rephrased the relevant sentences. We insist, however, that Boltzmann RTA codes (e.g., BoltzTraP) typically perform chemical-potential scans, not searches. Moreover, the advantages of the Kubo formalism over Boltzmann are most pronounced at low temperatures. To reliably find the chemical potential in a semi-conductor down to single-digit Kelvin, a bisection method directly applied to the Fermi function will fail, for the fundamental reasons explained in Sec. 4.2.

"* In section 5.1, the author mentioned that LinReTraCe is more stable that in Boltzmann codes. Could the authors elaborate on that ? In particular, I think it would be very useful to compare the electron/hole conductivity/mobility of a simple semiconductor (e.g. Si) as a function of temperature between LinReTraCe and any of the established BTE codes (they all roughly agree with each other) that computes the transition probabilities/scattering rates. Then accuracy, efficiency and stability can be discussed in details."

We specified that the better stability happens for the convergence with respect to the number of momenta. The reason is simple: The polygamma functions include the lifetime broadening and are, hence, smoother than the (derivative of the) Fermi function. Thus, this is a consequence of the kernel functions used, not of implementational details. However, in the case of BoltzTraP, instead of a direct evaluation of transport functions, a "transport distribution function"---a sort of transport density-of-states---is computed. Like for electronic-structure DOS (when not using the tetrahedron method) this requires a purely numerical (non-physical) broadening that smoothes the energy-dependence. LinReTraCe, instead, only uses broadenings that are physical.

We have reformulated the respective paragraph and now include, in the new Sec. 5.1. a chemical-potential scan for silicon, benchmarked against BoltzTraP2. The silicon calculation also serves as the new "make test" check of the installation.

"* The SciPost code acceptance criteria can be found at: https://scipost.org/SciPostPhysCodeb/about#criteria In particular they specify that benchmarking tests must be provided. I have cloned the code repository athttps://github.com/linretrace/linretrace.git As far as I can tell, there are no benchmarking tests, automatic test-suite with reference data etc. In addition, the userguide.pdf located in linretrace/documentation is very limited (5 pages). It is also indicated that "At least one example application must be presented in detail". Even if the authors presented some examples in the manuscript, no such example is available with the code and can be run."

In the revised version of the code and the manuscript, we now include a benchmark that also serves as a full example, as well as automated checks for the availability of prerequisites (make validate) and the compilation (make test). We understand that the SciPost article itself is the userguide. We have therefore renamed the additional pdf document with minimal instructions into "tutorial.pdf". Indeed, the latter simply provides a full summary of the installation procedure from the manuscript and another step-by-step guide for a simple (tight-binding) example.

"I have also looked at the Fortran source files located inlinretrace/src_linretrace The main.F90 is relatively well documented (comments in the code) but for other subroutines, it might be improved."

We thank the referee for this in-depth checking. We have further documented all subroutines in the code as mentioned, as well as explicitly documented all python classes and functions therein.

REPORT 2:

"Article 202206_00012v1 by Pickem Maggio and Tomczak is an interesting presentation of a new software package to evaluate transport quantities from electronic structure input data. The semi-analytical evaluation of the Kubo integrals comes from a recent PRB, and is implemented in the present software, going beyond standard Boltzmann approaches in allowing for more broadening effects and energy integration of the scattering processes between states. The code is a new and clean implementation, the presentation is quite thorough, and should be published. There are a number of useful general comments and physical observations about transport which will benefit many practitioners, though they are a bit buried in this very lengthy document."

We thank the referee for the positive assessment of our work.

"A few comments about the manuscript and the code: - My main complaint would be that the strongest advantages over normal Boltzmann approaches are not really demonstrated or showcased. The lifetime effects and scatterings possible (beyond the energy and T dependent models already in existing codes) do not jump out and shock me."

As the referee noted above, the new physics extracted from our methodology has already been published in two articles: Qualitatively different results (compared to Boltzmann codes) are particularly found in narrow-gap semiconductors with finite lifetimes. We now showcase this physics in Figure 10 by including, for comparison, results in the Boltzmann relaxation-time approximation (dashed line). Important new features beyond Boltzmann theory include, among others, the low-temperature saturation of the resistivity and the Hall coefficient---fully absent in semi-classical approaches. Our methodology is relatively new and its consequences have not been fully explored just yet. Nonetheless, in the spirit of open science, we intend to make our code already available at this early stage.

In the manuscript, we deliberately also included example systems for which Boltzmann codes give results comparable to Kubo linear response. We did so, because we think that practitioners of band-theory are one important target audience of our code. Indeed, also for band-like systems, our code may still hold technical advantages. For instance, LinReTraCe computes a few more physical observables than most other codes. Also, LinReTraCe can include temperature and state-dependent scattering rates more easily.

"The absence of a low T peak in TlPbTe is a disappointment, but further remains a mystery: it could be due to yet other extrinsic effects, and the proposed inverse engineering is a bit dangerous in this respect. Any rho(T) curve can be fit with many many scattering models if one allows arbitrary tau(n,k) variations..."

The absence of the low-T feature in the simulated resistivity is expected. Our code merely post-processes a given electronic structure and associated scattering rates. DFT only provides an approximation to the former. Having no ab initio results for the latter, we used a phenomenological scattering rate that yields good agreement at high temperature. The missing feature in the resistivity, paired with the good description of the thermopower (for which the influence of the scattering rate approximately cancels), then indicates that the scattering rate---its state- and T-dependence---is far more complex than our phenomenological ansatz. Hence, the disagreement with experiment is not a short-coming of our transport code, but instead calls for a better description of disorder-induced lifetimes using an electronic structure methodology beyond standard band-theory. We agree with the referee that fitting the scattering rate so as to reproduce experiment is dangerous. Fitting makes more sense, when the physics is established and only the coefficients in front of a known temperature-dependence (e.g., a Fermi-liquid like A times T^2) are determined through reverse-engineering.

"- The solutions of the quasiparticle equations are not always easy to linearize, as the authors require for their technique, especially if correlations are strong. At the end of 1.1 some discussion of the conditions of physical validity of this approach is needed. A "few kBT" can be huge in strongly correlated very narrow d or f bands."

The referee is correct. In strongly correlated systems, the self-energy may exhibit non-linearities already on frequency scales of a few k_BT. For instance, our method is indeed not applicable to Mott insulators, where there are poles in the self-energy in the vicinity of the Fermi level.

In correlated metals and band- or hybridization-gap insulators, self-energies are more benign. There, our methodology provides a conceptual improvement over RTA-Boltzmann technique, which, them, are only well-defined with static scattering rates and even neglect the linear term leading to the quasi-particle weight. Differences between RTA-Boltzmann and Kubo occur primarily at low temperatures, see Figures 10 in the manuscript. 50K corresponds to about 4meV. On that scale and temperature the self-energy can be linearly approximated to good accuracy for the target systems mentioned above. For strongly correlated materials at intermediate and high temperatures, a more costly full-Kubo evaluates should be preferred. Still, our methodology represents a significant step beyond semi-classical techniques.

In section 1.1 we now explicitly mention these limitations, as suggested by the reviewer.

"- p.13: "Generalizations of the band- curvatures to the Wannier basis have not yet been derived". I do not fully understand, though it sounds interesting. BoltzWann implements band curvatures to get the full transport quantities - what is missing or what more do the authors have in mind?"

For this point, see also the reply to referee 1. Band-curvatures can of course be computed from the Wannier Hamiltonian (we do so in Eq. 60). However, the occurrence of curvatures relies on the derivation of the transport function from the lattice/tight-binding/Wannier model, i.e. employing the Peierls substitution approach to the lattice Hamiltonian. A better way is to start from a field-theory perspective and to introduce the vector potential through the standard minimal coupling and to only afterwards express this continuum Hamiltonian in a local basis. For the conductivity sigma_{xx} this has been done and compared, e.g., PRB 80, 085117 (2019). For sigma_{xy}^{B}, instead, the continuum derivation has not yet been performed (as far as we know). Instead, the Peierls approach leads to the usual group velocities. The curvatures (second derivative with respect to momentum) emerge only through a technical step in the derivation, namely a partial integration, see, e.g., PhD thesis of Wenhu Xu (Rutgers University).

We have reformulated the respective paragraph along the above lines.

"- Figure 3: I have never had such massive problems in converging chemical potentials in insulators even with small effective temperatures. The proposed refinement is certainly interesting, but I suspect there is some instability which can be fixed with the usual algorithms."

As explained, there are fundamental limits to basic algorithms (e.g., bisection) for determining the chemical potential when using the Fermi function and a naive sum of occupations. Maybe there are more sophisticated algorithms that the referee considers to be usual. In any case, we have opted to use the refinement procedure discussed in the manuscript, which we demonstrate to work very well.

"- I have downloaded and installed the package. On my macbook pro the compilation fails with gcc 11 (cpu does not implement quadruple precision) - I did not debug for very long, but it would be good to have a more robust interface for the quad precision code. On our old cluster things compiled fine with gcc 10.2"

We are not providing support for the Mac ecosystem. Indeed, as the referee notes, the new ARM processors do not provide quad precision. There are numerous queries about this in Apple forums. Once/if a solution for these new type of processors is found, we will add instructions.

"- I tried a few systems from the github repo with examples - these should be mentioned more clearly in the manuscript."

We have followed the referee's advice and describe examples in more detail. In particular we have added a step-by-step example and benchmark for silicon. A further tight-binding example is discussed in the tutorial.pdf in the repository.

"The default config.temp files work but produce occasional errors (see below). A more thorough testing of the code and examples is essential. I agree with the other referee that a small test suite is important to ensure the code is well compiled."

In the revised code, we have included, both, standardized "make validate" and "make test" procedures. These checks will now tell the referee that the installed HDF5 package was an older version (1.10.5) than LinReTraCe requires (>=1.12.1), see section 3.1 in the manuscript. In previous HDF5 versions (1.8 and 1.10) the implementation of logical datatypes as enumerations (necessary as we interface HDF5 via the h5py library in Python3) is not working properly, hence the error message.

"The runtime for PbTe with 8 MPI threads was actually much longer than I expected (480s), comparing with BoltzTraP for example, though still bearable. The initial claim "as fast as Boltzmann" seems a bit overblown."

Per se, evaluating polygamma functions instead of derivatives of the Fermi function does not require significantly more effort. In the case of the PbTe (results in Figure 18), BoltzTraP will indeed be faster. The reasons are: (a) BoltzTraP evaluates, as alluded to above, a transport DOS once and for all. A change in temperature then only requires a new integral over this DOS. Our more sophisticated kernel functions instead have to be re-evaluated for each temperature point. (b) BolzTraP computes transport for a fixed chemical potential, while LinReTraCe searches for the chemical potential at each temperature point. When keeping the chemical potential fix, our calculations will speed up. Further, one could use the energy-truncation flag when creating the hdf5 input-file. The input data included in the repository used a rather generous window.

"Overall the manuscript English is excellent, but there are a few typos and a final read through would be beneficial. page 11 If the same quasi-particle renormalization (no s) page 13 the derivative in k should be outside a parentheses for ∂kU†(k)H(k)U(k)" "ensures" not "assures" absent in Eq 55 page 16 explain color code in figure 1 page 17 "the corresponding interface" page 21 The first column (no s) A generalization (no s) page 25 bottom: the phrase concerns a code snippet and not table 1, there is no continuity in the text. page 32 "identically"

We thank the referee for the careful reading and the spotting these mistakes, which we have corrected in the revised manuscript. Table 1 does show the program flow, but we have rephrased the respective sentence.

"M1Max chip compilation error cpsipg.f:11:43: 11 | COMPLEX(kind=selected_real_kind (32)) WPSIPG,W | 1 Error: Kind -1 not supported for type COMPLEX at (1) -> this means precision not supported by the cpu, but quadruple precision is certainly possible. An alternative coding might guarantee universal compilation"

M1/M2 chips are unfortunately ARM processors which are a very recent development: quad precision is required by our code, hence full support cannot be provided for these systems as of now.

"PbTe config.temp error: HDF5-DIAG: Error detected in HDF5 (1.10.5) thread 0:

000: H5T.c line 1754 in H5Tclose(): not a datatype

major: Invalid arguments to routine minor: Inappropriate type ??"

HDF5 version 1.10 is unfortunately outdated, we require v 1.12.1 or above, see implementation details. The new "make validate" will now alert the user if the required (versions of) libraries are not available.

Code:
- added "make validate" check of library/compiler prerequisites
- added "make test" suite to check compilation
- expanded in-code documentation
- renamed userguide.pdf to tutorial.pdf and updated it
- minor updates (v1.1.1 --> v1.1.7)

Manuscript:
- added reference PRB 102, 094308 (2020)
- added description of the new "make validate" and "make test"checks
- added Si-example (new Section 5.1) for benchmark comparison to BoltzTraP (Figs. 7, 8)
- added comparison of LinReTraCe and Boltzmann response in model of narrow-gap semiconductor (Fig. 10) to highlight merits of Kubo kernels.
- mention more clearly the limitations associated with the linear approximation of the self-energy in the introduction.
- reformulated comment regarding band-curvatures in Wannier basis (Section 2.4.2)
- corrected typos
- minor changes throughout