Requested changes:

1. “The inline nonzero moment defined just after Eq. (3) is confusing since the romanized d seems out of place. If the hbars from Eq. (3) were written with the dt terms, then the reader would quickly understand the notation associated with the inline equation. You could also write the dt/hbar after the delta.” Response: As suggested, we have moved the dt term to the end of the in-line expression. We agree that this improves the clarity.

2. “The reader could be helped with the inclusion of a sentence in the introduction just following the two sentences about PGPE and SPGPE. At this point in my read, I now know your tools and your focus. However, if I am not familiar with the tools, then I won't really see the contrast between them and the classical field approach. A summarizing sentence would set the reader up for the discussion of your methodology.” Response: We have added a brief description of the c-field methodology, and the PGPE and SPGPE in the introduction as requested.

3. “I understand why 2.2 and 2.3 are separated out from 2.1 but as single paragraphs, they dangle. Perhaps it would be better to just incorporate them?” Response: Thank you for the suggestion, and we have considered this. However, following the report of the second referee, we have extended 2.3, and so have decided to retain the current structure.

4. “In 2.2, is there a reason for working with the region defined by 95% of the radius? I didn't find intuition with a quick skim backward. I have several ideas about why one might/should do this but it would be good to know the author's perspective.” Response: We only consider vortices in the region defined by 95% of the radius due to numerical considerations. In the low-density regions near the boundary of the system, the phase can become filled with numerical singularities (‘ghost’ vortices) that are picked up by our vortex detection algorithm. These ghost vortices do not enter the system or affect the dynamics in a meaningful way, so we therefore remove them from consideration by ignoring any singularities detected in the region r > 0.95*R. We have included this explanation in the text.

5. “The start of 3.2 and its footnote was dense/technical. Consider reformatting into a trichotomized list. The technical matters could be unpacked more but in the end, the reader needs to walk away with the three regimes. This gets a little lost in the succinct but technical discussion which left me feeling split between re-reading the technical statements or pushing forward with the trichotomy in mind.” Response: We appreciate the suggestion to reformat the start of Sec. 3.2 as a trichotomised list, and to expand on the explanation therein. We have updated the text, and hope that the readability of this section has been improved as a result.

6. “I found the first sentence of 3.4 awkward and felt that it could be clearer with a re-organization of content to get rid of several parenthetical statements.” Response: We have rewritten the opening sentences of Sec. 3.4, and hope that they are now clearer than they were previously.

7. “In the second to last paragraph of section 3, the authors state that the system may eventually return to the Gaussian fixed point. Why does this argument exist? It seems to me that clusters might have a similar sort of topological protection and so I could imagine that they may not totally annihilate each other, at least at zero temperature as the authors note. This idea seems to split over this paragraph and two others in section 4. As it is an important and interesting question, perhaps a more focused discussion would be better. Alternatively, the authors could add on a sentence that hints at the forces at play that would prevent long-time protection of the clusters and the openness of this question.” Response: To clarify the prediction that the system should eventually return to the Gaussian fixed point, we have added the following text to the end of Sec. 3.4.2: "It is predicted that the system would eventually return to the Gaussian fixed point, because vortex–phonon interactions give rise to gradual vortex diffusion [61], which serves to break up vortex clusters and encourage vortex–antivortex annihilation. However, it is possible that the time required for the system to cross over to the Gaussian fixed point may approach infinity as T->0." Naively, one would expect that a nonlinear equation such as the Gross-Pitaevskii model should exhibit ergodicity, and therefore eventually 'forget' the details of its initial state and thermalise by finding the maximum entropy state available under the constraint of energy and norm conservation. In this case, the thermal state would correspond to a system without vortices, as this would maximise the entropy. However, it is possible that the topological protection of the vortex clusters would prevent the system from ever reaching this state at the lowest temperatures. We have also added the following sentence in the conclusion to make this point: "Although the vortices should gradually diffuse toward the fluid boundary due to interactions with phonons, perhaps this effect is not strong enough to overcome the topological protection of the vortices."

8. “I like the rapid-fire question list at the end of the conclusion. It really helps punctuate their program/results. However, since I was enjoying the read, I wanted to hear about the author's opinions about whether the methods of the paper are enough to address these questions. If not, then what changes/hurdle do they expect? Completing in this way would keep me interested in their next work or becoming more acquainted with this/their research program.” Response: Thank you for the suggestions. We have expanded the conclusion section to elaborate on some of our questions and comments, and to incorporate some recent work. We hope the referee finds these additions valuable.

Requested changes:

1. “The SPGP is valid if the temperature is higher than half-Tc. Here, Tc is the BEC transition temperature. For example, the review paper [Advances in Physics 57 (5), 363-455 (2008)] describes this point. However, their numerical simulations consider a low temperature region 0<T/Tc<0.3, so that their use of the SPGP cannot be justified. Also, both SPGP and PGP methods introduce a cutoff parameter, but they have not given any discussions on how to determine it. The authors should address the issues and describe them in the revised version.” Response: While the referee states that c-field methods are only valid at finite temperature for T>~ 0.5 T_c, this is in fact a rule-of-thumb estimate based on the typical numbers of atoms in a 3D harmonic trap in alkali gas experiments. Here we consider a two-dimensional flat-bottomed trap, which has a significantly different density of states to a 3D harmonic trap, and so the temperature region of validity is different. In fact, due to the enhancement of fluctuations in lower dimensions, the validity regime is actually extended. The discussion of validity of c-field methods in the Blakie et al. review paper the referee quotes is quite detailed, and describes several conditions that are not always consistent with one another. One often-quoted validity condition for c-field methods, also discussed in the Blakie et al. review, is that the number of particles per mode should be significantly larger than one, in order that the neglection of quantum fluctuations is valid. If we normalise our system to have O(10^5) atoms in the c-field region (an experimentally feasible number), on plotting the momentum distributions for n_0 = 1 and n_0 = 0.75 [see insets of new Fig. 4(c,d)], then we can see that all the modes in all simulations in fact satisfy this condition reasonably well. For completeness, we have attached to the response an additional figure (Fig. 1 in the attachment) showing the data for n_0 = 0.9 and n_0 = 0.85 (details as for Fig 4, with n_0 = 0.9 in the left column, and n_0 = 0.85 in the right). A paragraph providing a discussion of the choice of cutoff and the validity of the c-field methodology has been added to the end of Section 2.

2. “Discussion on the relation between turbulent decay and NTFPs. Usually, universal relaxation of NTFPs can be characterized by dynamical scaling of spatial two-pint correlation functions [Phys. Rev. Lett. 101, 041603 (2008)]. The nice review is given in [J. Berges, arXiv:1503.02907 (2015)]. Despite the importance of the scaling, the authors do not show time evolution of the correlation functions at all and there are no discussions on the dynamical scaling. In my opinion, the present results are insufficient for judging whether the turbulent decay is close to the Gaussian or anomalous nonthermal fixed points. I recommend the authors to confirm the dynamical scaling and estimate the power exponents.” Response: This was an excellent suggestion. Based on this, we have performed additional analysis of the dynamical scaling of both the two-point correlation function and the mode occupation spectrum, and we have included this in the revised manuscript. We have also extracted the correlation length from the correlation function, which shows evidence of power-law scaling in time. We have added a new figure (Fig. 4) to the manuscript, which shows the collapse of correlation functions and occupation number distributions, when the axes are rescaled with the correlation length. We have also added a second subfigure to Fig. 6 (previously Fig. 5) in which the correlation length is plotted as a function of time, for comparison with the mean nearest-neighbour vortex distance. We have accordingly updated the text in Sec. 3.4 to describe the new analysis, and provide further discussions of the results. We hope our additions provide a more complete picture of the statistical evolution of the system.

Minor questions:

1. “Relation between Onsager vortex formation and inverse cascades. In the introductory part, Onsager vortex formation and inverse cascades are mentioned. However, I cannot understand a relation between two concepts in the 2D GP model in the present introduction. It would be nice to clarify it.” Response: We have added to the sentence in the introduction which draws the association between the inverse cascade and the clustering of same-sign vortices in 2D quantum turbulence.

2. “Two rate equations for vortex number. Is it possible to fit the numerical data using Eq.(5)? If possible, we cannot judge which rate equations are correct. How do you answer it?” Response: As stated in the text, these equations should be regarded as entirely phenomenological. As such, we would not claim that either model is 'correct' or 'incorrect' based on our fitting procedure, unless one provided a significantly better fit than the other. It is indeed possible to fit our data to Eq. (5) and obtain values of Gamma_1, Gamma_2, Gamma_3, Gamma_4 as a function of temperature. We state in the text that: "When fitting the two rate equations (5) and (6), we find that both describe the data equally well. However, since the latter involves fewer free parameters, we opt to use it over our earlier model." We have attached a figure to our response in order to clarify this (Fig. 2 in attachment). This figure is identical to the main frame of Fig. 3 of the manuscript, except that we have additionally overlaid the fits from Eq. (5) as red dotted lines. The values of Gamma extracted from these fits are also given in the table in the inset. As can be seen, the difference between the quality of the fits to Eq. (5) and (6) is minimal. However, the Gamma values extracted from Eq. (5) do not provide a particularly clear interpretation due to the large number of free parameters in the model. There is a general trend in the data as the temperature is changed (i.e. the four-body term dominates over the one- and two-body terms at T=0, while the opposite is true at higher temperatures); however, the picture is less clear than that provided by Eq. (6). We also emphasise that the quantitative values of Gamma should not be regarded with too much significance; rather, we are interested in the qualitative trend in the Gamma values as a function of temperature. We argue that fitting to Eq. (6) gives a clearer interpretation of how the N_v(t) curves change as a function of temperature (as seen in the insets of Fig. 3 in the text), and hence we have chosen to use that model in the manuscript.

3. “One-body particle loss Realistic ultra-cold gases suffer from one-body particle loss. Does it disturb the evaporative vortex cooling?” Response: One-body loss is expected to have a relatively unimportant effect on vortex dynamics in 2D. At the first level of approximation in a flat-trap, it will simply reduce the density, which will cause a small increase in the size of the vortex cores, and speed up the time-scale for vortex dynamics, but not have any other effect.