First of all, I apologize to the authors and the editor(s) for the delay in producing this referee report.

The manuscript describes a new implementation of the Monte Carlo (MC) method to solve the vector radiative transfer equation. The manuscript confirms recent evidence that MC methods can combine both flexibility and accuracy. The first had long been clear, and indeed MC methods had been the usual choice for complex geometries. The second attribute, however, had been less obvious to the radiative transfer community. The current paper gives further evidence that MC methods are "exact" and, when properly implemented, can achieve arbitrary accuracies if one is willing to pay the cost of running sufficiently long simulations. In my opinion, this is the main point of the paper and I believe it's an important one.

The manuscript provides enough context to understand the implementation of the MC method, and demonstrates sufficiently its performance with comparisons against standard benchmark cases. This said, my main criticism of the manuscript, however, is that  the implementation is based on "intuition", i.e. on tracking what happens to the simulated photons at each collision event. That makes it difficult to understand, for example, why neglecting the polarization state of the photon may lead to solutions that fail to convergence. In the same vein, it's difficult to understand why the implementation neglects variance reduction techniques. The integral formulation described in Garcia Munoz & Mills, 2015 (based on O'Brien, JQSRT, 1992) gives a more straightforward framework to tackle these issues. 

The manuscript is very well written. References in the manuscript seem complete.

I didn't check, but the MC implementation appears to be publicly available, which will help further disseminate the ideas developed here. 

I recommend the manuscript for publication, after addressing the following very minor comments.

Methods. It took me a while to figure out that the MC implementation is valid only for plane-parallel atmospheres, perhaps because the introduction mentions the application of polarimetry to disk-integrated observations. I suggest making it clear in the Introduction (#1) or Methods (#2) sections. 

* line ~16. w.r.t. is not standard, I think.

* line ~20. "at the top of at the atmosphere". Language issue.

* Fig. 1. In my printout, angles theta' and phi' become corrupt. It may just be a problem with my pdf reader though. Further, it would help if you could show the actual geometry of the problem, with the cells in both the x and y directions. This would make it clearer that the current geometry is plane parallel. 

* In Methods section, please state how V>0 is defined. 

* line ~150. "w will become zero upon TOTAL absorption".

* line ~166. Having the code decide between scattering by a gas molecule or an aerosol particle (or different types of them, if different aerosol types coexist) is inefficient in a MC method. It would be notably faster if the optical (scattering, absorption) properties are properly averaged and assigned to "average" particles. Could the authors please comment on this?

* pg. 8. Footnote. Setting w=1e-16 to terminate a simulation seems unnecessarily stringent. My guess is that the calculation converges for w~1e-6. Could the authors please confirm this statement?

* pg. 9. Just wondering, how easy would it be to extend the surface treatment to non-Lambertian reflection?

* line ~261. "In the backward mode, every photon that is not completely absorbed contributes to the result, assuming that the light source is at an infinite distance." I was a bit puzzled to read this. My understanding is that the backward model works best when the observer has a narrow entry cone or, equivalently, when the observer is very far. For the backward implementation, it is no problem if the star subtends a finite solid angle.

* line ~324. Did the authors "switch off" the direction sampling based on the polarization state of the photon and, in that case, did they run into problems to converge the model? I feel curious to know the answer.

General comments: 

This paper presents a new 3D atmospheric radiative transfer model based on Monte-Carlo Method.

It presents the methods commonly used to track the photons through the atmosphere and to compute the reflecances. Several validation cases are then presented showing that the MONKI code can correctly reproduce previous results. Finally, a simulation of the Earth and Venus atmosphere were computed and analysed, demonstrating that the code is able of handling different optical properties.

The paper is well written and clear. However, it does not really contribute to anything new. 

In the conclusion of the article, three contributions of the MONKI code are listed, but are not convincing  as they stand : 

1- line 481 “common limitation of Monte Carlo codes is their failure to converge in optically thick and strongly polarising atmospheres”. Can the authors  be more specific about the values beyond which the convergence does not occur and give an exmaple of such atmosphere. 

2- “The current version of MONKI does not use variance reduction techniques” . I suppose  other codes may choose wether or not to use  variace reduction technique, which are often necessary to achieve convergence in complex 3D atmospheres. So, it does not appear  to be a real advantage.

3- “In order to speed up MONKI’s simulations, we implemented and presented a precomputation approach for optical paths to the detector.”. This may be a new method for obtaining converging results more quickly but quantitative values showing the gain in computational speed have to be added to the article. Note that adaptative grid was already introduced in Villefranque et al. (2019)

4 - “MONKI can import different particle scattering matrices and can, in principle, also simulate light scattered by non-spherical aerosol particles and ice clouds”. The same applies for other code as long as optical properties are known.

To increase the novelty of the paper, I would recommend to go deeply in the analysis of the methods used to sample the scattering direction at each interaction presented in appendices E and F. I may be wrong, but it seems to me that the rejection method is not commonly used in Monte-Carlo Atmospheric radiative transfer. A deeper analysis of this method comparing to the conventionnal one (tabulated cumulative distribution function) in terms of accuracy, convergence and computation time, could be valuable to increase the interest of this publication for the community.  

In summary, I  understand that the authors wish to publish this article to have a reference to cite it for further study or when it will be used by other scientific teams as  the code is freely available, which is important to mention. Nevertheless, I would recommend to try to improve the paper by taking the above comments into account.

Specific comments : 

    • Equation of line 137 : This is correct only in the principal plane, ie when the scattered radiation is in the same plane of the incident radiation.

    • To avoid confusion, in Equation (10) and (13), authors should be more precise by indicating the direction of the incident radiation and of the view (for eq 10) or scattered (eq 13) direction.

    • Section 2.6, line 274 : does it means that the  χ  matrix need to be stored and each scattering in the reverse mode to be apply at the end of the incident Stokes vector ?

    • Section 3.1, Table 1: The results are in forward or bacward mode ? Can you add the absolute and relative difference in the Table

    • Section 3.3, Figure 6 : Why do you use the standard deviation of 3DMCPOL and not the one of MYSTIC that it is used to calculate the difference 

    • Section 3.3, Figure 6: Why divided the results by 1000/E0 and not using the same quantity than in the IPRT comparison (Emde et al. 2018) ?

    • Section 4 : I do not know Venus atmosphere, is it common to have clouds between 50 and 70 km ? Do you have some references to support the configuration ?

    • Appendices E and F describe how the rejection method is used for the Rayleigh and Mie scattering instead of the inverse CDF method using tabulated values as is usually done. Have the two methods been compared in terms of accuracy and computation time. These can be an interesting result to add in the paper. See my general comments.  

    • Line 505. MONKY is designed to support the preparatory studies for the aforementioned spectropolarimetric missions, which means calculating radiances at different wavelenghts. Howerver, no information is given on how the multispectral properties of the atmosphere including gaz absorption and particles scattering will be accounted for. 

Typos : 

line 58 : delete “and”

Reference : 

Villefranque, N., Fournier, R., Couvreux, F., Blanco, S., Céline, C., Eymet, V., et al. (2019). A path-tracing Monte Carlo library for 3-D radiative transfer in highly resolved cloudy atmospheres. Journal of Advances in Modeling Earth Systems, 11. https:// doi.org/10.1029/2018MS001602