Photon counting statistics in Gaussian bosonic networks

1) The paper is overall well written and is also pedagogical.
2) The overview of the field, i.e. list of references, is satisfactory.
3) The steps of the theoretical analysis and the flow of the calculations can be followed very straightforwardly. An useful compact notation (matrix form) is used.

Examples of the developed method are provided exclusively for the simplest cases, specifically for two or three coupled bosonic modes. Some of the quantities analyzed, such as correlation functions and entanglement for the case of two coupled modes, could, in principle, be derived without employing the proposed method.

The authors investigate the full counting statistics of N coupled bosonic modes, considered a many-body open quantum system.
Two of the authors have expanded upon their earlier work (Ref. [35]), extending the analysis from a single resonator to the case of N coupled resonators.
The work is intriguing and, in my opinion, warrants publication, provided that the authors satisfactorily address my questions and properly resolve the concerns outlined below.

1)
The model Hamiltonian in Eq. (1) is presented in the rotating frame of the driven resonators. The authors state that the resonators "are all driven at the frequency $\omega_d$" which is not accurate. The squeezed term appearing in Eq. 1 ($\sim {\hat{a}^{\dagger}}^2_i$) originates from the parametric pumping at $2 \omega_d$. If this were not the case, after transforming to the rotating frame, the parametric pump would manifest as fast-rotating and, by applying the rotating wave approximation (RWA), it would need to be disregarded. Therefore, the authors should revise the text accordingly.

2)
Since the system is purely linear, the parametric pumping at the drive frequency $2\omega_d$ with $\omega_d \sim \omega_i$ can yield divergent parametric oscillations (self-induced oscillations) unless the system is sufficiently detuned from the region of such instability (Arnold's tongue) and/or the strength of the parametric drive $r_j$ is significantly smaller than the damping $\gamma_i$ (numerical coefficients may vary based on the model's notation). Could the authors provide a comment on this point?

3)
There is another implicit assumption in their model that has not been explicitly mentioned: The frequencies of the resonators must be nearly equal, within a specified range determined by the resonator losses, i.e., $\gamma_i$, namel $\omega_i \sim \omega_d$.
If this condition is not met, one resonator will be driven out of resonance, which effectively means it is not driven at all, particularly if the detuning is much larger than the linewidth of its response function. Technically, this implies that, by applying the rotating wave approximation (RWA), the coherent drive $f_i$ for the resonator with $|\omega_d - \omega_i| \gg \gamma_i$ will appear as fast oscillating and would need to be disregarded in the Hamiltonian.

4)
The authors should emphasize that the Lindblad equation in Eq. (3) describes a drive-dissipative open quantum system, specifically a system out of thermal equilibrium. This is linked to the observation that the operators appearing in the dissipator in Eq. (3) are not associated with the "true" eigenstates of the quadratic Hamiltonian that describes the system, therefore the system cannot reach the thermal equilibrium. In this regard, in addition to citing Ref. [42], which pertains only to equilibrium conditions, the authors may consider citing, for example: (i) [https://arxiv.org/abs/2409.10300] and/or (ii) [https://arxiv.org/abs/2407.16855].

5)
Related to the previous issue, the authors discuss the validity of the Lindblad equation that includes the dissipator operators identically to a system of non-coupled/non-interacting resonators, which is also known as the "additive approximation".
They assert that this approximation is valid for 'weak coupling in the network compared to the dissipation rates $\gamma_i$', namely $(g_{jk}, \lambda_{jk})\ll\gamma_i$. However, the authors present results with $g = \gamma$, which contradicts their assertion.
However, the regime $(g, \lambda) \ll \gamma$ is fully dissipative and is opposite to the strong coupling regime $(g, \lambda) \gg \gamma$ , where the coherent interaction between the different modes is larger than the dissipative rate, potentially leading to interesting effects (e.g., entanglement, etc.).
I am not sure whether the condition provided by the authors is the right one. For instance, in circuit QED (one qubit coupled to a resonator), the 'additive approximation' is frequently employed even in the "strong coupling regime", and it functions effectively. The system proposed here seems analogous: if we consider a transmon as a qubit (i.e., a nonlinear resonator) coupled to a cavity (a second resonator), I see no fundamental difference between the two systems. The "additive approximation" surely breaks down in the Ultra-Strong Coupling (USC) regime in which we have $(g,\lambda\sim \omega_i)$, see a recent paper Phys. Rev. Lett. vol. 132, p.106902.
Although it has been well-established that decay operators derived from an uncoupled system can result in unphysical effects when applied to a coupled system [J. Phys. A 6, 1552 (1973)], their use is frequently considered a valid approximation.
Could the authors analyze better this point?

6)
The authors have developed a powerful theoretical framework, and it would be valuable to see how this approach can be applied. I understand that the first two examples analyzed (two coupled modes and entanglement) are pedagogical in nature, intended to illustrate the implementation of the method. The last example, involving the three-mode circulator, is particularly intriguing due to its broken time-reversal symmetry.
In the final part of this section, the authors consider the case without the coherent drive. In their calculations, they assume (or fix) the net current flowing between cavities 1 and 2, and then they present the probability distribution of the current. However, why is it not possible to 'derive' the net current flowing through the system and demonstrate the broken time-reversal symmetry instead of assuming it from the outset?
The proposed method should be capable of directly computing the current. Could the authors provide a clearer explanation of this issue?

Finally, I would like to comment on the authors' last statement at the end of the paper: 'It may be possible to extend our framework to non-Gaussian systems with Kerr or other nonlinearities.'
This seems improbable a priori. The ansatz of a Gaussian state is heavily utilized, particularly for transitioning from the Lindblad equation to time-dependent equations for the correlators, averages, and shifts.
In principle, non-Gaussian states could lead to a hierarchy of equations for higher-order cumulants. Could the authors provide better insight into how it might be possible to extend the method to accommodate non-Gaussian systems?

- After Eq. (17), there is an abrupt discontinuity in the theoretical derivation. Specifically, the quantity $\rho(\vec{n}, \vec{m}, t)$ is neither defined by a formula nor through an equation. Although the authors refer to Refs. [54, 55], for the sake of completeness, it would be helpful to include a formula or equation here to define $\rho(\vec{n}, \vec{m}, t)$.

- In Eq.(14) it appears a sum over the modes $\sum_{j=1}^N$ which I don't understand. The quantity $P(t;\vec{n},\vec{m})$ already includes all the modes so why do we need to sum? or maybe I missed something?

Ask for major revision

1- the reader can obtain an in-depth understanding of photon counting and the method
2- clear potential for follow-up work in a variety of directions
3- well-written and illustrated

1- missing connection to an alternative approach

The authors present a method to obtain the photon counting statistics of Gaussian bosonic networks where every bosonic subsystem is embedded in an individual thermal environment. The time evolution of the full network is described by a Lindblad master equation that is at most quadratic in the bosonic ladder operators. The main method the authors employ are Gaussian states which are naturally applicable (the Liouvillian is quadratic). The authors introduce necessary concepts to understand the connection between the Lindblad equation, probabilties for photon emission and absorption, and the (cumulant) generating function of the photon current. The cumulant generating function can be obtained by solving three coupled differential equations. Additionally, the two-time correlators ‘waiting time distribution’ and the ‘second-order coherence’ are discussed; both quantities describe the correlation of the radiation. After explaining how to obtain the photon counting statistics at long times, the authors finish by applying their method to three examples. There, they also discuss a potential entanglement witness that can be obtained from the photon counting.

The paper is well-written and given the usefulness of the method, the length of the manuscript is definitely justified. All key ingredients to understand the subject are stated in the main text. Further, additional information and details can be found in the appendices such that every reader can obtain an in-depth understanding of the topic.

In the introduction, the authors state that there does not yet exist a method to obtain the counting statistics of Gaussian bosonic networks. Regarding this, the authors need to be aware of the concept of ‘third quantization’ which also works on quadratic Liouvillians. Relevant references would be, e.g., [arXiv:1007.2921] (Lyapunov equation for covariance matrix), [arXiv:2304.02367] (generating functions), [arXiv:2302.14047] (connection to Wigner function). What is the connection between the two methods? What are respective advantages and disadvantages? I suggest to add an appendix to discuss this (see points 1 and 2 in the requested changes).

This (and other points) have to be appropriately addressed in the revised manuscript.

1- What is the direct connection between Gaussian states and third quantization? I believe that it is related to the fact that the algebraic Riccati equation can be solved by an eigenvalue decomposition (Wikipedia).

2- What are advantages and disadvantages of both methods? Is there something that can be calculated in one framework easier than in the other? Include an appendix to discuss point 1 and 2.

3- For the TMS-system, or nondegenerate parametric oscillator, the full counting statistics have already been obtained in the past, e.g., in [Phys. Rev. A 46, 395](not cited). Please include relevant references and compare the results.

4- After the submission of the paper, advancements in the detection of single microwave photons have been made, e.g, [Phys. Rev. Lett. 133, 076302]. Maybe update the list of references regarding this point.

5. How is the waiting time distribution impacted by a finite detection efficiency? Maybe include a numerical simulation for such a case.

6- While it is mentioned and cited, the authors might consider to explicitly include the connection between the second-order coherence and the Fano factor to provide a more complete picture and draw a nice connection between the short- and long-time statistics.

7- I did not find a definition of $\Theta_{N,0} $. Is $\Theta_{N,0} =\Theta_0 – I_{2N}/2$ (page 11) where $\Theta_0$ is obtained from Eq.10?

8- While the Hamiltonian of the system can be deduced by the matrix $\mathcal{A}$ it might be useful for the reader when the Hamiltonians are written down in section 4.

Ask for major revision