Summary comments:

As a physicist my comments focus on the physics and the data interpretation, rather than the math used.

The approach presented is interesting and may be useful, but the more information on the context for the present work is needed. For example, the relationship of the errors and uncertainties in the present and previous approaches is not sufficiently clear. While some of the figures suggest much smaller errors (apparently 10K)  than in the current GOLD data products (30-40K near solar max) the temperatures retrieved for the May 2024 storm deviate dramatically from the released data which show unambiguous consistency between the temperature and composition fields, an indication that the temperature structure is real but is not apparent when the approach described in the paper is applied. The geophysically interesting structure observed during the storm seems to be obscured by the approach presented, possibly indicating biases inherent in the proposed approach when analyzing rare, atypical observations.

As noted by the authors, the approach presented may be useful for data with low signal to noise. The assumption of a normal distribution in GOLD’s current temperature retrieval, while appropriate for most of the dayside where the difference between normal and Poisson statistics are negligible, is not appropriate for all solar zenith angles. However, the authors mischaracterize the current solar zenith angle restrictions in the GOLD data products. Retrievals are possible with lower signal to noise data, but the current retrieval introduces a cold bias (as discussed by Evans et al. 2024). While the figures presented by the authors show that the proposed approach also becomes subject to biases at the larger solar zenith angles where signal to noise is lower (Figure 3), it may be productive to explore the use of Poisson statistics for data with low signal to noise. The use of data at large solar zenith angles is further complicated by the increases in peak emission altitude that occur at large SZA (as discussed in Evans et al. 2024).

The paper has multiple assumptions in the work presented and in comments on previous work that need to be clarified or corrected. An example is the comment in line 98 that a fully calibrated instrument model. Such a model is unnecessary for any of the approached that researchers have used, but a relative calibration is necessary for all the approaches, especially when using multiple emission bands. Other cases are noted in the annotated .pdf file.

A related issue is the possible geophysical implications of assumptions made, e.g., cap harmonics, by the authors. There is an underlying asymmetry in the morning versus evening or the northern latitudes versus the southern latitudes distorted by the symmetry within the harmonics. Another example is the significance of being dependent on the posterior distribution, e.g., in the May 2024 storm results.

The limitations of the Cantrall are noted very late in the paper (line 312?), after the extensive discussions that are likely to be unfamiliar to most readers. Seems to deflate the significance of the work. If introduced earlier, could you compare the present results with the previous Cantrall Matsuo results and include some discussion of the role of Poisson statistics in the differences? That could clarify the significance and possible limitations of each new fundamental assumption in your retrieval, Poisson statistics and Bayesian probability.

This manuscript presents a hierarchical Bayesian inversion framework for retrieving thermospheric neutral temperature from ratios of photon counts in two FUV spectral sub-bands (LBH). The method explicitly accounts for Poisson counting statistics and introduces spatial structure through a permanental (squared-Gaussian) process prior. The approach is demonstrated using both simulated GOLD disk emission data (November 2018) and real GOLD L1C observations, including the May 2024 geomagnetic storm. The results indicate improved robustness at low count levels and extended usability toward higher solar zenith angles compared to existing ratio-based techniques.

Overall, the paper is technically sound and addresses a relevant problem in FUV remote sensing. The methodology is novel in this application and well motivated. However, several aspects of the presentation, validation, and physical interpretation would benefit from clarification and strengthening. Below you may find my suggestions to improve this manuscript.

Major comments: 

Line 55: The paper treats the ratio of Poisson intensities as the primary geophysical variable. While this is mathematically well justified, the physical meaning of the inferred ratio field should be discussed more explicitly. In particular, it would be helpful to clarify under what conditions the linear relationship Z~mTeff + z0 remains valid (e.g., dependence on viewing geometry, vibrational population assumptions, or background contamination). Although it may be further explained in Cantrall and Matsuo, 2021; it would be helpful to provide a brief description in this manuscript to provide sufficient context to the reader.



Line 180: The spherical cap harmonic (SCHA) kernel is well motivated, but the practical procedure for selecting the smoothness parameter ν is not fully specified. The manuscript concludes that ν=1+10⁻⁸ performs best, but this appears empirical. A clearer description of how this choice would be made in an operational context (e.g., cross-validation, heuristic rules) would improve reproducibility.

Section 3.3: The reliability diagrams are a strong point of the paper. However, the text would benefit from a more explicit interpretation of what “overdispersed” means in this context and how this impacts scientific use (e.g., assimilation or model validation). It may also be useful to report typical uncertainty magnitudes in Kelvin for representative conditions.



Section 4.2: The comparison with TDISK is mainly qualitative. Quantitative metrics (e.g., mean bias, RMS difference) between the new retrieval and TDISK, in addition to CRPS relative to WAM truth, would make the evaluation more convincing. This is particularly important for the May 2024 storm case.



Minor comments

- Several symbols are introduced quickly (e.g., alpha, beta, q, p in the beta-prime distribution). A short table summarizing distribution parameters and their meaning would help readability.

- Figure 1 would benefit from explicitly labeling the local time and storm phase in the caption. Additionally, I suggest to plot the difference (in K or %) between the reconstructed temperature and the ground truth. It would provide the reader additional information regarding the spatial variability of the error.



- Line 10: There is an equation (E(N(R)) with no number. Also, could you add a general description of variable s (bold s) here (e.g., does it represent a 3-D dimension?). If I understood correctly, it is linked to the bold s used in the kernel, and there is a definition for bold s here. 



- In Eq. (12), the condition α_a>1 for the MAP estimate should be highlighted earlier, since it affects interpretability when counts are very small.



-The Laplace approximation introduced in Section 2.2 plays a central role in enabling the computational efficiency of the method, as it leads to a closed-form posterior and avoids the need for sampling-based Bayesian inference. While this connection is implied later in Sections 3 and 5.1, it may be helpful to state this link more explicitly when the approximation is first introduced.

The manuscript reports a method to derive thermospheric column integrated temperature or effective temperature by considering the impact of Poisson distribution in the FUV photon counts, such as those in the GOLD observations. The temperature is verified with simulated FUV disk emissions for a weak storm and GOLD observations for the same storm and May 11, 2024 super-storm. While this approach is worth to try, there are still some details that need to be included for a better understanding of the method.

• The 1^st half of the abstract below to introduction. The term “temperature” is not clearly defined in the abstract. Does it refer to a temperature at a fixed altitude, pressure level or height integrated (effective) temperature? There is no definition of specific bands for the ratio calculation. The last statement is on a broader solar zenith angle (SZA). What is the SZA range that this manuscript deals with? Most of the statements in the abstract do not provide specific information.

• Equation (2). The column-integrated or effective thermospheric temperature depends on wavelength of FUV emissions. Will the method give the same Teff, regardless the wavelength used in the temperature retrieval?
• How do readers use the Teff to refer the thermospheric condition? Since T increases with altitude continuously in the thermosphere, the Teff is expected to be equal to the T at a specific reference altitude. Can this reference altitude be estimated? If there is no such information of the reference altitude, how Teff is useful for understanding the thermospheric dynamics and space weather application?
• Equation (2). The volume emission rate and O2 observation and Tn.
• The manuscript lists many equations with different approximations (e.g. Eqs 3-12). However, there are no data-based plots to support the approximations.
• Lines 165-166. Could a plot of the simulated Poisson-distributed photon counts included in the manuscript?
• Lines 149-150. Due to the construction of the intensities as spatial processes,

the parameters of the marginal posteriors are calculated using spatial information and thus endow spatial structure on the resulting temperature field.

             Does this mean that the potential improvement in Teff estimation is associated with a lower spatial resolution? If this is true, the method introduces a bias where localized small scale temperature fluctuations or enhancements (e.g. TAD) is added to the background temperature. Smoothing of TDISK in Figure 5 could result in the retrieved temperature in Figure 5.  

• Figure 1. The Teff estimation depends on values of γ. How is the γ selected when there is no ground truth data to validate it?

• Figure 1. Is the center of the Teff maps for the nadir looing direction? If this is true, the estimated Teff apparently depends on the look angle and/or solar zenith angle. This makes the Teff map difficult to use for science investigation.
• Figure 5. The new T map shows a smooth plot while there many localized fluctuation in the Tdisk map. This is hard to judge that the smooth map is better than the Tdisk as many local temperature fluctuations. Furthermore, it is also better to compare Teff and Tdisk over an entire storm (including quiet time before and after the storm).
• Figure 6. The TDISK in the 2^nd column shows a nice vortex in enhanced temperature in the southern hemisphere which nicely reveals the Coriolis effect. However, such a vortex structure is not seen in the Teff.
• The manuscript extends the work by Clayton and Matsuo (2021) with a different approach to derive column integrated temperature. However, there is no comparison between the Teff.
• The algorithm is partially validated by using simulated FUV spectra based on WAM neutral profiles. How are the Teff validated using observed GOLD spectra?