We thank the referee for detailed comments that we address below.
We have made significant changes to improve the manuscript in response
to the referee and marked all changes in blue. We have now made the item
numbers consecutive, so all items now go from 1 to 20.

> 1. Definition of “bipolar” region: How do the authors define the bipolar
> region ? What does “weak”
> bipole mean in Table 1 and 2 ? For instance, for models with theta = 0 and
> 30 in Fig.7, i see there are a lot of magnetic concentrations which are not
> much larger than the size of turbulent eddies, albeit the region with the
> vertical field is confined in a relatively narrow band. When filtering at
> the certain wavenumber, does the large-scale bipolar region remain even in
> these models ? The authors should clarify the definition of the “bipole” in
> this study.

Bipoles are opposite polarities of spots whose size is that of several
correlation lengths. A definition is now given in the introduction on page 1
on the left column. Here we also now specify the meanings of "weak", although
the more qualitative attributes (yes, no, and weak) are now discussed in the
penultimate paragraph of Section 3.2 and they are also now indicated in Fig.9.
The Fourier-filtering is indeed done for the purpose of averaging over
the turbulent scales and retaining the large-scale magnetic structure
of the bipoles, as discussed in Section 2.1 and 3.2

> 2. Differences between DNS and MFS: The authors try to distill the essence
> of the DNS with the MFS. However, it is hard to say that the MFS could
> successfully reproduce the characteristics of the DNS because there are a
> lot of differences between them; the amplitude of the growth rate, its
> latitudinal dependence, the pattern of the vertical magnetic field organized
> at the surface, and its dependence on the Coriolis number and latitude. The
> authors should discuss much more what causes these differences, how the
> coronal layer affects the result at least qualitatively, and why the
> narrower band-liked structure is developed in the DNS but not in the MFS.
> Why is the mode with the higher wavenumber more preferably excited in the
> MFS ?

We agree and are particularly concerned about the lack of a banded
appearance of the BRs in the MFS in comparison with DNS. We have now
mentioned this shortcoming at the end of the new Section 4.5. There is
also the mismatch in growth rates, which can partly be explained by us
having neglected the vertical anisotropy, which was called q_g in our
earlier papers. This is also now addressed in the paper.

> 3. Differences between DNS and linear theory: The WKB approximation might
> be applicable to the system because k_f H_rho ~ 30 is assumed in the
> authors’ DNS, at least in the linear stage. Then, according to the linear
> theory given in Losada+12, the max growth rate of the NEMPI is expected to
> decrease with the colatitude, i.e., higher in the pole and lower in the
> equator because the frequency of the inertia wave becomes higher in the
> larger colatitude. Nevertheless, it seems that the growth rate of the NEMPI
> measured in the DNS does not show such a monotonic behavior (e.g., Figs.6
> and 13). The authors should clarify why it is so different from the
> prediction of the linear theory. I guess the fastest growing mode of the
> NEMPI may not be included in the range k < kf/6 (cut-off wavenumber) for
> some models , OR, maybe not fully-resolved by the DNS because the
> horizontal wavenumber of the fastest growing mode should be varied with the
> latitude. It might be better to do the convergence check at the lower
> latitudinal region.

Our section 4.2 has now been significantly expanded. In particular, we
have now addressed the concerns in the last paragraph of Sect. 4.2 we
discuss the difference in growth rates between vertical and horizontal
magnetic fields. Both cases play a role in the simulations, which is
one of the reasons why a direct comparison with analytic theories
is not possible.

> 4. Interpretation on the rotational dependency: Based on the linear theory
> of Losada+12, the increase of the rotation rate enhances the frequency of
> the inertia wave, resulting in greater suppression of the NEMPI. The result
> of the DNS is however not quite so simple. It shows non-monotonic
> dependence of the NEMPI on the Coriolis number. It may be one of
> interesting findings in the authors’ results that the relatively-moderate
> Is this a consequence of the coronal layer or of the rotational shrinking
> of the turbulent eddies ? Is the rotational dependence a natural result of
> the competition between the enhancement of the NEMPI due to the rotational
> shrinking of eddies and the rotational suppression ? The authors should
> give some physical interpretation on the non-monotonic rotational
> dependency at least qualitatively.

At Co ~ 0.002-0.003, the growth rate just decreases. Or is the referee
talking about the increase of the growth rate at Co > 0.05 seen in the
green line? This is caused by dynamo action and is now also discussed
in the expanded Sect. 4.2; see the penultimate paragraph.

> 5. S3.3: The authors discuss in the first paragraph of S3.3 the dependence
> of the BR formation on the initial field strength despite the latitudinal
> dependence should be presented here according to the title of the
> subsection. With including this, the authors should re-edit the content
> ``drastically” to maintain its integrity.

We have now reworked this entire section and put things in the right
place.  The subsection numbers are now shifted by one, because we have
now also added here what used to be in the appendix.

> 6. Fiducial model & convergence study: The authors explain ``main”
> properties of BR formation with ``Run A4 of Appendix”. It makes the logical
> flow of the manuscript difficult to follow because the diagnostics
> quantities of Run A4 cannot be found in Table 1 of the main body. It may be
> better to combine the convergence study of Appendix into the main body if
> the authors use Run A4 to explain the basic characteristics. Otherwise, the
> results of Run A2 should be used to explain them.

We have now moved the appendix to the main part of the paper and inserted
it as Section 3.1. We also now discuss mostly Run A2 instead of A4 and
use A2 for most of the relevant figures.

> +++++++++++++++++++++++++
> Other relatively minor comments:
> +++++++++++++++++++++++++

> 7. S3.1:
> - what is the definition of E_M^z in Fig.3 ? Is the energy contained in the
> z-component of B-field ?

This is now Figure 5. E_M^z is the magnetic energy spectrum based on
the z component of the magnetic field. This is now explained in the
third-last paragraph of Section 2.1.

> 8. - what is the physical meaning of (2k_*E_M^{z,max})**1/2 in body text &
> Tables. is it corresponding to the something like compactness parameter of
> BRs ? what should we read from the value ?

This is also now addressed this in the third-last paragraph of
Section 2.1. It is the energy per logarithmic wavenumber interval
and is evaluated at k* = 2*k1.

> 9. S3.2:
> - caption of Fig.4 : black —> blue

Corrected now. This is now Figure 7.

> 10. - Fig.5: for the fair comparison, the authors should exhibit the surface
> distribution of the Bz-field with the different Co at the ``same
> colatitude”.

This is now Figure 6, except that we show there now a run with Co=0.0076.
For Co=0, we should not have written theta=0, because all values of
theta are equivalent. So in this sense it is a fair comparison.
We have now added a comment in the caption.

> 11. S3.3:
> - this subsection is poorly organized. it needs a considerable revision.
> - In the body text, the authors mentions the models with Co = 0.0015 and
> 0.0075. But, in Fig.6, the latitudinal dependences of the models with
> 0.0012, 0.0023, and 0.0061 are shown. In addition, the caption of Fig.6
> seems to disagree with the figure.

We have now moved the first part of this subsection to section 3.6 as
the second paragraph. The correct numbers are those in the table,
i.e., 0.0012, 0.0023, and 0.0076. We have made them now consistent.

> 12. - Why is the error bar required for Fig.6 ? If so, it should be added to
> Figs 4, 8, and 13 as well.

Figure 6 is now Figure 8. As the referee points out elsewhere, there are
significant non-monotonicities, and one wonders to what extent they can
be explained by the finite length of the time series. If the time series
is too short, the time average will not be statistically meaningful.
To assess this, we have added error bars. In some cases, the error bars
themselves were unreliably short, which has now been corrected. We have
now added a comment about error bars in the penultimate paragraph of
Section 2.1. We have also now added error bars in all the other figures.

> 13. - In Fig.7, why is the bipolar structure pronounced at
> theta=60 deg? why is the inclination of 45 degree is preferred?

This is now Figure 9. We have now added a remark in Section 3.5 that the
inclination angle will be discussed again in Section 4.3 when we study
MFS. We are there able to conclude that it is not an artifact of the
finite domain size, and that even very close to the equator (10 degrees
away from it), one observes this 45 degree angle, while at the equator
it is 90 degrees. We thus state that we remain puzzled about this finding.

> 14. S3.5:
> - is there any clear relationship or scaling law between the critical
> Coriolis number which gives the rotational quenching and the initial field
> strength ?

No, we have not yet obtained any clear scaling law between the
critical Coriolis number which gives the rotational quenching and
the initial field strength.

> 15. S3.6:
> - Fig.10: Since the suppression of the turbulent pressure by the large-sale
> B-field is the primary cause of the NEMPI, it would be natural that the
> amplitude of the negative effective magnetic pressure, i.e., |P_eff|,
> increases with B_0. However, for me, it is hard to understand the reason
> why P_eff is almost constant when varying the Coriolis number because the
> kinetic energy of the turbulence should be varied when changing the
> rotation rate. The authors should give some physical interpretation on this
> at least qualitatively.

A detailed theoretical study of the effect of rotation on the function
P_eff(B_0) is the subject of a separate ongoing study. So far we have not
obtained any clear scalings. The current results have now more clearly
been illustrated in the new Table 4 on page 10. We also note that the
kinetic energy is not strongly influenced by rotation in this setup with
forced turbulence, unlike for convective turbulence.

> 16. - Fig.11: Where are the red and blue curves coming from ? The solution of
> eqs (13),(14) & (17) ? If so, the authors kindly mention that in the body
> text explicitly.

Yes, it is a fit to Eq.(15) and we have now provided the details of how
the fit parameters are obtained; see the new paragraph around Eq.(21).

> 17. S4:
> - How many mesh points are used for the MFS ? same as DNS ? The authors
> should show it explicitly because the properties of the NEMPI depends on
> the resolution as the authors studied.

We have now added text in the new section 4.1 where we discuss the number
of meshpoints in the second paragraph.

> 18. - As already mentioned in the major comment 2, I again recommend that the
> authors discuss much more about the difference between the DNS and MFS and
> its cause.

We now discuss two major issues. One concerns the q_g parameter that
is found to be negative and reduces the efficiency of NEMPI. We have
now added this in the MFS; see the revised Eq.(12) and the newly added
parameterization in Eq.(13). Results are now included in Fig.14. Regarding
the lack of banded structures, we have now added a new subsection 4.4.

> 19. - In S3.6, the authors show that the parameter sets q_{p0} = 13, beta_* =
> 0.65, and beta_p = 0.18 are better for describing the DNS in the viewpoint
> from the negative effective magnetic pressure. Then, it should be better to
> adopt this parameter set mainly, not the parameter set found by Losada+13,
> to compare the MFS with the DNS. Is there any physical reason to adopt the
> set of Losada+13 ?

We have now addressed this in the third paragraph of the new subsection
4.1 where we mention various uncertainties in the MFS modeling as the main
reason for adopting a broader view when it comes to specific parameter
choices in the mean-field models.

> 20. - In the DNS, I can see the region with the strong vertical field is very
> localized and isolated in some models, e.g., Fig.5 (rightmost panel) &
> Fig.7 (left bottom). I mean that there is a clear separation between the
> magnetic concentration and the other region in some DNS models. In
> contrast, in the MFS, the vertical field prevails in the surface region.
> Even in Losada+12, I cannot find such a compact isolated magnetic
> structure. Is it possible to reproduce the localized and isolated magnetic
> structure by the NEMPI even in the MFS when adopting arbitrary parameter
> set ?

This is a concern and has been addressed in the response to the previous
point and also in connection with our response to point 2 of the referee.
We have also now added a new figure 21 to show more clearly the appearance
of isolated magnetic structures in the nonlinear phase.