Thank you for supporting publication of the present manuscript.
As detailed below, we have taken into account the comments and suggestions of the reviewers that we believe have overall improved the quality of the present paper.}
{In my understanding methods that are solved through non-linear equations display multiple solutions due to the dependence on the starting point and presence of multiple local minima.
While the presence of these multiple solutions may be connected to the disappearance of denominators, the problem is much more general.
The problem of multiple solutions exists even if fully self-consistent GW is solved in a finite temperature formalism (on imaginary axis - which effectively removes the problem with vanishing denominators).
Different solutions can be then illustrated by considering their physical properties for example $S^2$.
In such cases, the potential energy curves are smooth within a region where one solution is dominant and one can show properties of these solutions as a function of $S^2$, see J. Chem. Phys. 155, 024101, and J. Chem. Phys. 155, 024119.}
\\
\alert{
The reviewer is correct; these issues are more general and do not originate solely from the vanishing denominator.
For example, we have recently studied multiple solutions in another non-linear method, coupled-cluster theory [J. Chem. Theory Comput. 17, 4756 (2021)], where their origin stems from the non-linear nature of the equations, only.
However, in the case of $GW$ at least, we have clearly observed that these multiple solution issues appear when denominators vanish.
More precisely, in most $GW$ methods, because one only considers the quasiparticle solution at each iteration while discarding the other solutions (known as satellites, discontinuities are observed.
This has been studied in details in a previous study [J. Chem. Theory Comput. 14, 5220 (2018)].
In a fully self-consistent scheme where one takes into account all these solutions at each iteration, these issues do not appear, as we have recently discussed in a recent paper [Front. Chem. (9, 751054 (2021)].
The references provided by the reviewer has been added in due place and this point has been clarified in the revised version of the manuscript.}
{Convergence accelerators such as DIIS, KAIN, LCIS, can be successfully used and smooth PES can be produced when one is relatively close to a local minima, one should not expect that DIIS will help with convergence where multiple close lying solutions exists and when two solutions are competing, see J. Chem. Phys. 156, 094101 (2022).}
Indeed, convergence accelerators such as DIIS can be used to ease convergence but they will not make these discontinuities disappear as their origin is more profond.
This particular case is discussed in a recent paper [J. Chem. Theory Comput. 14, 5220 (2018)] where we have provided our implementation of DIIS within $GW$ methods (see Appendix).
This point has been further stressed in the revised version of the manuscript.}
{It is my understanding that the calculation on the illustrative example H2 in 6-31G basis illustrates the existence of multiple such solutions.
The issue of introducing a regularizer to avoid vanishing denominators is more of an issue for GW on real axis.
Of course these issues of vanishing denominators and multiple solutions due to the non-linear character of the Dyson equation are connected but I think it is worth mentioning explicitly how different versions of GW are affected by these issues.
Frankly, I find it infuriating that in the GW community many times it is not clearly spelled out which flavor of the so-called GW is used and how the equations are solved.
\alert{We agree with the reviewer that, historically, the literature has been unclear on how to solve the $GW$ equations, but we would like to stress that our group has made a clear effort to provide all the working equations and details necessary to solve these equations in the different cases.
We have added additional details and references in order to guide the readers to the relevant papers when one can find all the necessary details regarding the implementation of $GW$ methods.}
{The authors should specify on which axis (omega) they perform their evaluation since this will not be clear to a regular chemistry audience.
The convergence characteristics on imaginary or real axis is very different.
Please, specify it. It may seem trivial to a regular GW practitioner, it is not trivial for a regular theoretical chemist not very familiar with Green's functions.}
{It would be reassuring if the authors mentioned that they use $O(K^9)$ scheme just to illustrate multiple solutions and this should not be done in regular calculations.
A reader not familiar with GW algorithms may conclude that this what they advocate that should be done.
On the same note, it would be good if the authors mentioned that the scheme that scales as $O(K^6)$ is again not what one does in most GW schemes that are in use for calculations of larger systems.}
\alert{In Fig.~1 we plot the quasiparticle solution for each orbital, i.e, the solution with the largest spectral weight which is obtained using the $G_0W_0$ scheme.
{In Fig.~3, it is unclear to me that the traditional shift and the shift introduced by Evangelista result in demonstrating that Evangelista shift is superior.
It seems to me that the difference between the values of qp energies before and post shifting are of the same order of magnitude for both regularizers.
Could authors elaborate what they see differently and what my untrained eyes could not see?}
\alert{We agree that this part was incomplete and, therefore, we have made extensive efforts in order to improve it.
In particular, we have changed the notations regarding the various regularizers that we have studied.
We now use $\eta$ for the traditional regularizer and $\kappa$ for Evangelista's regularizer.
We have also included additional graphs for different values of $\eta$ and $\kappa$ which shows how the quasiparticle energies are altered by the choice of the regularizing function and the values of $\eta$ and $\kappa$.
\alert{Indeed the HOMO and LUMO orbitals do not show discontinuities along the dissociation coordinate so no need for a correction. Thus, it is an important feature that the regularization introduces only a small correction for these orbitals. It is also true that the regularization introduces a correction of few eVs for the LUMO+1 ($p=3$) and LUMO+2 ($p=4$) orbitals but we have to note that the quasiparticle solutions of Eq.~2 for these orbitals appear at the poles of the self-energy. So the regularized self-energy has to do a large correction which leads to large error on the quasiparticle energies. Moreover, it is also essential to notice that we talk here about the $G_0W_0$ scheme but in case of a partially self-consistent scheme then the use of regularization seems critical.
{In Fig.~5, I find very surprising that authors find that the curve is so "bumpy" already at the equilibrium distance.
Usually in the quantum chemistry community, it is a standard practice to pass a solution from the previous geometry to the subsequent one as a starting point.
This removes "bumpiness" in regions where denominators do not vanish.
Are indeed the denominators already vanishing here at the equilibrium?
If they do, then I find it highly unusual.
Or are these competing solutions that could be removed by passing the starting point between neighboring geometries without a regularizer?}
\\
\alert{
}
\item
{In my opinion, the novelty of this work is in showing how the unfolding scheme introduced by Booth group can be used to access different quasiparticle solutions.
I am less convinced about the goodness of the regularizer since I suspect that the solutions depend quite a bit on the magnitude of $\eta$.
Can this paper be published as a rapid communication?
Yes, but a plot showing that the total energy is relatively insensitive as a function of regularizer $\eta$ would make the argument more convincing to broaden the applicability of the partially self-consistent Green's function on the real axis.}
As stated above, the regularizer section has been improved and we are happy if our manuscript is published as a Regular Article (as recommended by the editor).
\noindent\textbf{\large Authors' answer to Reviewer \#2}
{The article by Monino and Loos presents a way to avoid numerical difficulties arising from GW calculations in moderately correlated electronic systems, where a small perturbation in a particular parameter of the problem, such as the interatomic distance, leads to a sudden jump in the main quasiparticle peak, as computed within the GW approximation. The article is written clearly and includes new, original insights that could be useful for an audience of specialists interested in electronic-structure calculations. I support the article for publication, but ask the authors to address my points below: }
{In standard GW calculations in solids, one assigns a quasiparticle peak to the solution of Dyson's equation that is associated with the largest value of the spectral function.
This is the only general method that works in solids as, rigorously, there is no strict solution for Eq.~2 in a generic system - unless one lets $\omega$ be complex.
Making this connection, especially around or before Eq.~8, would clarify the method to a readership that is not only interested in molecular systems.}
{I am not very happy with the term "intruder state" and with calling the procedure a "regularized GW method".
These terms give the impression that there is a spurious solution that doesn't belong to a particular part of the Hilbert space and which leads to a sudden jump of the quasiparticle energy.
Instead, my understanding is that the physics behind the jumps has nothing to do with spurious states, but one's insistence on making a one-to-one correspondence of the mean-field energy with a peak in the interacting spectral function.
For strongly correlated systems, the spectral function may have pronounced satellite peaks, whose intensities may vary significantly with a small change in an external parameter.
This is a point that the authors explain well in the text, but the terminology of intruder state confuses this otherwise simple physical picture.
In addition, the term "regularized GW method" gives the impression that the GW approach has unphysical features.
Instead, this procedure is nothing but a small broadening parameter that smooths out sudden jumps between neighboring peaks in the spectral function. }
We understand the point of the reviewer but "intruder state" and "regularization" are well-defined terms in the electronic structure community which re not linked with the appearance of spurious solution.
The intruder state problem is well documented in multireference perturbation theory and comes usually from a poor choice of the active space.
By definition, an intruder state has a similar energy than the zeroth-order wave function and should be then moved in the model space; this is exactly what is happening in the case of $GW$.
{How does the result depend on the details of the "regularization" function?
The authors should be aware that traditional GW calculations performed in solids often use a small broadening, and hence such a "regularization" is naturally captured by many codes - whether or not on purpose.}
As detailed below (see the answers to Reviewer \#1), we have thoroughly modified and expanded this section to test the effect of the regularization function and its parameter.
Shouldn't the brown dashed curve with $\eta=0.01$ go to the solution with the main quasiparticle peak for the left and right limits of the plots?
My naive intuition is that a small $\eta$ means that we are doing less smoothing to the self-energy/spectral function (or, in the language of the authors, performing a weaker "renormalization").
So, why doesn't the curve with $\eta=0.01$ jump abruptly to the solutions computed without the smoothing/"regularization" procedure?}
