Multimode Fock states with large photon number: effective descriptions and applications in quantum metrology

We develop general tools to characterise and efficiently compute relevant observables of multimode $N$-photon states generated in non-linear decays in one-dimensional waveguides. We then consider optical interferometry in a Mach-Zender interferometer where a $d$-mode photonic state enters in each arm of the interferometer. We derive a simple expression for the Quantum Fisher Information in terms of the average photon number in each mode, and show that it can be saturated by number-resolved photon measurements that do not distinguish between the different $d$ modes.


I. INTRODUCTION
Photonic states with a large and fixed number N of photons play a crucial role in quantum technologies but are extremely challenging to prepare experimentally. The paradigmatic example are single-mode Fock states, |N ∝ (a † ) N |0 , where all the photons share the same spatio-temporal mode (a † ), and which are the basis of many quantum metrology protocols [1][2][3][4]. Nowadays, the most widely used method to generate them is based on combining with post-selection heralded single photons emitted in spontaneous parametric down-conversion processes [5][6][7][8][9]. This method, however, suffers from an exponential decrease of efficiency with N , hindering its application for large photon numbers. Single-mode Fock states can also be emitted naturally from entangled atomic states in ensembles with many more atoms than N [10]. However, exciting such atomic states is highly non-trivial because of the linear energy spectrum of such systems [11][12][13][14].
A way of circumventing these limitations is the use nonlinear systems for the generation of such photonic states. These type of systems appear in many different contexts, such as in cavities with Kerr-type non-linearities [15], in multilevel quantum dots due to biexciton binding-energies [16,17], or even in atomic ensembles simply because an atom can not be doubly excited or by exploiting Rydberg blockade [18,19]. These mechanisms ultimately translate into either non-harmonic energy splittings (Kerr cavity QED) or non-harmonic decay rates (saturation), which can be harnessed for multiphoton emission. An illustration of that is the proposal put forward by us in Ref. [20] to generate multiphoton states with quantum emitters coupled to photonic waveguides [21][22][23][24][25][26][27][28][29][30]. There, N excited emitters interact with the waveguide in the so-called mirror configuration [28,31], such that its dynamics is described by the well-known Dicke model [32]. In that situation, the emitters experience a nonlinear decay process, known as superradiant decay, which enhances the probability of emitting the photons into the waveguide as compared to other decay channels. Beyond the collective enhancement, the non-linearity has another effect: the * marti.perarnau@mpq.mpg.de photons released into the waveguide have an inherent multimode structure [13]: where a † ki is the creation operator of a waveguide photon of momentum k i . The coefficient A {k} = A k1,k2,···kn characterizes the multimodal structure of the wavepacket, and will be non-factorizable (A {k} = √ N !A k1 · · · A kn ) for photons emitted from any type of non-linear system (e.g., nonharmonic energies or non-linear decay rates). This non-trivial multimode nature of the emitted wavepackets forces one to revisit the results derived for single-mode Fock states as they are not necessarily valid anymore. For example, the multimodal structure poses limits on the scalability as Fock state sources from spontaneous parametric down conversion processes [33][34][35], is required to accurately predict the scattering of quantum pulses [36,37], or, as we showed in our recent manuscript [20], renormalizes the results of single-mode quantum metrology protocols.
Motivated by these observations, the goal of this article is to develop general tools to deal with multimode states generated in non-linear decays, both from the point of view of its characterisation as well as applications in quantum metrology. We start by considering a general N -dimensional emitter decaying in a non-linear fashion with a waveguide. The wavefunction of the emitted photonic state (given by A {k} in (1)) can be obtained through the techniques of [38], and it involves N ! terms due to bosonic symmetrisation. Given this highly non-trivial state, the main contributions of this article are: tential for quantum metrology [1][2][3][4]. Building on our previous work [20], our goal is to extend well known results in quantum optical interferometry [39,40] to the presence of a nontrivial multimode structure within the input photonic states, as in Eq. (1). For that, we consider phase estimation in a Mach-Zender interferometer, where the input state of each arm of the interferometer is a generic d-mode state with a fixed total photon number. Then, our main contributions are: 1. We show that the quantum Fisher information [41] (QFI) Q takes the particularly simple form where n i , m i are the average photon number in the ith mode of the two incoming wavepackets.
2. We show that the QFI (2) can be saturated by numberresolved measurements which cannot distinguish between the different d modes.
3. Finally, we also discuss the effect of photon loss in the interferometer and in the measurement devices given the proposal of [20] for quantum-enhanced metrology with twin Dicke superradiant states.
The paper is structured as follows: We start presenting the general multimode structure of the photons emitted from non-linear systems into waveguides in Section III, whereas the tools to characterize such photonic states are developed in Section III. These tools are applied to superradiant photonic states in Section IV. In Section V we consider quantum metrology with multimode states, and finally we summarize our findings in Section VI.

II. NON-LINEAR SYSTEMS DECAYING IN 1D WAVEGUIDES
In the interest of generality, we consider the emission process coming from a N -level system (|0 , |1 , ..., |N ) with energies ω j . Its free Hamiltonian is then given by (taking ≡ 1): with σ ik = |i k|. The system is coupled to a 1D waveguide, described by a one-dimensional and chiral photon bath with a linear dispersion (both the chirality and linearity assumptions can be relaxed obtaining similar results) ω q = cq. Taking c ≡ 1, its Hamiltonian read: Finally, the system-bath interaction Hamiltonian is assumed to be given by: FIG. 1. A the non-linear system (in blue) with N levels with energy ωn, couples to a 1D waveguide. The coupling with the waveguide induce single-photon transitions n → n − 1 (in red arrows) at a rate γn. Around the non-linear system, we depict the two example of non-linear system that we consider along the manuscript, that are, equally spaced atomic ensembles and anharmonic cavities.
where γ j denotes the decay rate of the transition from the jth to the (j − 1)-th level. The global Hamiltonian describing the emission process is then given by the sum of the three terms: H = H S + H B + H SB . The whole physical set-up is illustrated in Figure 1. We consider that initially the waveguide B is in the ground state, whereas the system S initially contains N excitations. When the excitations decay into the waveguide, the photonic state is described by the wavefunction (naturally extending the considerations of [13]): where T stands for time-ordering, and O t satisfies Furthermore, the bosonic creation and annihilation operators a s , a † s satisfy the standard commutation relation The use of this general light-matter Hamiltonian H allows to capture the physics of very different models, such as: • Saturated atomic ensembles in the atomic mirror configuration. As explained in Refs. [13,20], within the Markov approximation the coupling of the ensemble with the waveguide can be described by a single collective dipole operator. This can be effectively described as N -level system with equally spaced energy levels, ω n = nω 0 , but non-linear decay rates: γ n = Γ 1d n(N − n + 1). In the text, we shall call these emitted states superradiant states, or superradiant photonic states.
• Anharmonic cavities. In the case of a non-linear resonator, with a Kerr non-linearity of the form H B = ω a a † a + U a † a(a † a − 1), but coupled to the waveguide in a linear fashion such that its reduced dynamics is given by the standard Lindblad form Γ 1d (2aρa † − a † aρ − ρa † a)/2, the energy levels are now non-linear: ω n = nω a + n(n − 1)U , while the decay rates are harmonic γ n = nΓ 1d .
Along this manuscript, we will focus on the characterization of the first ones, as they will be the most relevant for quantum metrology. However, all the formalism developed is valid for any combination of {ω n , γ n }.

III. TOWARDS EFFECTIVE CHARACTERISATIONS OF MULTIMODE STATES
We start this section by developing techniques that enable us to compute relevant observables of |φ (N ) efficiently for large N . This is motivated by noting that the analytical form (6) contains N ! terms due to the time ordering in (7), making such an expression challenging to handle beyond low N .
We first consider a single emitted photon (with γ ≡ γ 1 − γ 0 and ω ≡ ω 1 − ω 0 ): where we have defined B Using (8), the commutation relation follows. Given such single-mode operators b γ,ω 's, we have developed techniques to compute efficiently observables of the form The computational techniques for dealing with (11) are rather involved and are developed in detail in the Appendix I. Here, we instead explain the main ideas and implications: 1. The computation of (11) is developed by expressing it as a recurrence relation, which is then transformed into a matrix multiplication.
2. The solution is exact as the integrals in (6) are carried out analytically. Yet, in practice it is convenient to perform the matrix multiplication numerically in order to access large N .
3. The size of the involved matrices is at most (4 n (N − n + 1)) 2 , although usually this can be reduced if there are some symmetries in the calculation (i.e. if some of the x i 's and y i 's in (11) are the same). In practice, this means that for n = 1 (corresponding to average photon number), one can easily reach up to N ≈ 1000, whereas for n = 2 (corresponding to the variance) one can reach N ≈ 100. This should be contrasted to the naive calculation of (11) from (6) which involves (N !) 2 integrals.
4. Since the b † γ,ω 's form an (overcomplete) basis of the Fock space spanned by a † t ∀t, one can in principle compute arbitrary observables through this approach.

IV. CHARACTERISATION OF SUPERRADIANT PHOTONIC STATES
We now apply the machinery developed in the previous section to the characterisation of superradiant photonic states emitted when N excited atoms are placed next to the waveguide in the atomic mirror configuration, as previously proposed by us in [20]. This corresponds to taking |φ (N ) in (11) with γ n = Γ 1d n(N − n + 1) and ω n = nω 0 .
Before proceeding to its characterisation, let us mention that there are two main features that make this proposal particularly appealing [20]: • In the absence of photon loss (i.e., when the atomwaveguide set-up is perfectly isolated) and assuming perfect control on the system's Hamiltonian, the protocol is deterministic and scalable. That is, simply by placing more atoms next to the waveguide, one can generate a larger N -photon state.
• In the presence of photon loss in free space, the probability of success scales as ≈ 1−ln(N )(Γ * /Γ 1d ), where Γ 1d is the decay rate into the waveguide and Γ * into free space, for Γ * /Γ 1d 1. This slow decrease with ln N arises due to the enhanced collective decay and should be contrasted with the standard success probability ≈ 1 − N (Γ * /Γ 1d ) obtained for N independent decay processes (see Ref. [20] for more details).

A. Average photon number
We first consider the photon number We start by fixing the decay rate γ and varying ω around ω 0 , corresponding to the frequency of the two-level emitters. Figure 2 shows how the average photon number is centred at ω 0 , as one would have expected physically. Next we consider n(γ, ω) for a fixed ω = ω 0 and varying γ. The results are shown in Fig. 3, where we compute n(γ, ω 0 ) for N = 100, 200, 300. Interestingly, note that the maximum of n(γ, ω) appears at x = N/ ln N . To understand  this, note that the time scale τ of decay of the photons is proportional to γ −1 . In particular, for the superradiant decay, the jth collective exitation decays with a time scale γ −1 j , with γ j = Γj(N − j + 1), and the average decay time is given by Hence, we realise from (6) that the choice γ = N/ ln N corresponds to single-mode photons decaying with the average decay time of the superradiant photons. This provides an heuristic explanation for the optimal choice of γ that maximises n(γ, ω 0 ), i.e., the number of photons in the mode b γ,ω .

B. Most relevant modes
The results of Figures 3 and 2 suggest that a rather large proportions of photons (between 50% and 60% for N = 100, 200, 300) can be contained in a small set of modes centred around the frequency ω and with inverse decay rate γ = N/ ln N . This motivates us to consider a set of D modes with frequency ω 0 and varying Γ = jN/ ln N , with j = 1, ..., D; that is, {b jx,ω0 } D j=1 with x = N/ ln N . Note that these modes are not orthogonal due to (10). A set of orthogonal modes can be constructed by solving the generalised eigenvalue equation, where T and R are matrices of size D 2 whose elements are given by This leads to a set of D bosonic modes The usual commutation relations are guaranteed by (14) (plus appropriate normalisation).
Let us now define the number operators (for a given D), and the corresponding ratio of photons, In Figure 4, we show C d as a function of N for D = 10 (i.e. the modes c k 's are a linear combination of 10 modes b kx,ω0 's). Note that C d quickly saturates with N . The values of C d for different D's are shown in the following Table, which is evaluated at N = 100: These numbers can be slightly increased by considering larger d's or D's. It is remarkable that with only 2 (3) modes we can cover 98.4% (99.6%) of the photons, and that more than 90% photons live in a single mode. Hence, although the superradiant state |φ (N ) may naively appear as a highly multimode state, it can be described by means of only a few modes. At the same time, it is worth noticing that although we reached this result by a rather heuristic method, these descriptions are almost optimal: since the small set of considered modes already contains around 99.9% of the photons, it is not possible that by considering a much larger (possibly infinite) set of modes the effective descriptions can considerably change. In other words, while a few modes seem to suffice to describe |φ (N ) with high accuracy, it is not possible to describe it by a single one.

C. Variance
Let us further characterise the fact that most photons are contained in a few modes by studying the fluctuations of the number operators n k 's given in (18). In particular consider the variance with ... ≡ φ (N ) |...|φ (N ) . To compute such expressions, first we expand them using the commutation relations (17) as and similarly for higher n 2 j . Each term can be evaluated by using (16) and then computing .., D and save the results. In the numerical simulations of this section, we take D = 8, for which we already cover 99% of the photons as shown in Table IV B.
The results are shown in Figure 5 for d = 1, 2, 3 and D = 8.
As expected, σ n d decreases as d increases, and for d = 3 and D = 8 the variance is rather small: less than one photon for N = 40. This confirms the previous results that a few modes can cover most of the photons contained in the superradiant state. Furthermore, we also observe that σ n d grows linearly with N , which is compatible with the proportion of photons staying constant with N as shown in Figure 2.

D. Effective descriptions
The main insight of the previous sections is that 2 or 3 modes suffice to describe most of the photons of |φ (N ) . We can use this insight to build effective descriptions of |φ (N ) within these subspaces. This can be achieved by computing the overlaps for all {c j } d j=1 such that d k d = N , and where d is the number of considered modes. For that, we first express (22) as a combination of 0|b k1 x (16). Then, in Appendix I, we develop a recurrence relation method to compute 0|b k1 x In order to compute (22), it is convenient to first compute all possible 0|b k1 x such that j k j = N and save such coefficients. There are O(N D ) such coefficients, which provides the necessary space to carry out the computation. Overall, the whole computation is challenging, both in terms of the number of operations and the complexity of each one, but becomes feasible for small N and D ≤ 10, which is enough to capture most of the photons (see Table IV B).
We consider the following unnormalised two and three mode state: or N = j,k |α j,k | 2 for three, indicates the overlap between two or three mode subspace is. In Figure 6, we show how N is close to 1 for low N ≤ 10, especially for d = 3 where it stays above 0.97. This implies that |φ (N ) can indeed be well described by two or three mode states as in (23). For d = 2, the coefficients are given in Figure 7 and the coefficients for d = 3 are provided in Appendix I. These approximate states of the form (23) become handy in calculations for quantum information tasks, as we will later illustrate for quantum metrology, and importantly our techniques enable us to quantify how close are our approximate states to the real |φ (N ) . As a final remark, we note that the linear decrease of N with N is also compatible with our previous results where we observed that the proportion of photons in a given mode stays constant.

V. APPLICATIONS IN QUANTUM METROLOGY
We now apply the different insights and techniques developed in the last section to quantum optical interferometry, in particular by considering phase estimation in a standard Mach-Zender (MZ) interferometer. We consider that each arm of the MZ inteferometer is described by a set of modes {a k } d k=1 and {b k } d k=1 , respectively (see Figure 8). When d = 1, we recover the standard two-mode optical interferometry [4], and our goal is precisely to extend well-known results to the multimode regime where d > 1. We start this section by introducing some basic concepts of quantum metrology, as well as describing the set-up we consider here in detail.
We consider the estimation of a parameter ϕ by measuring N photons which encode information about it. The (possibly entangled) N -photon state is described by |ψ ϕ . Let us assume that we apply a measurement M on |ψ ϕ , the statistics being described by a probability distribution P (s j |ϕ) where {s j } are the possible outcomes of the measurement given the value ϕ of the unknown parameter to be estimated. If this process is repeated ν times, in the limit ν 1 the Cramer-Rao bound guarantees that the mean-squared error ∆ 2φ of any unbiased and consistent estimatorφ of ϕ is lower bounded by [42,43]: where C is the classical Fisher Information (CFI), The CFI quantifies the best resolution of a particular estimation scheme as defined through |ψ ϕ and a measurement M . In their seminal work, Braunstein and Caves optimised C over all possible quantum measurements M [41]. The resulting quantity is the Quantum Fisher Information Q (QFI) , which satisfies, and, for pure states, is given by withψ ϕ = ∂ ϕ ψ ϕ . The QFI quantifies the potential of a particular state |ψ ϕ for quantum metrology [41].
Let us now discuss how we encode ϕ in the standard Mach-Zender (MZ) interferometer, but extended to the multimode regime. We consider as an initial state |ψ = |φ A , φ B , where with and where {α i } ∈ C(n) iff i α i = n, and with n + m = N . That is, we consider that two independent (but generic) states of m and n photons, each of them described by d modes, enter the two arms of the interferometer. For d = 1, the states are single-mode states and hence Fock states. We also assume bosonic commutation relations, which can always be satisfied by a proper choice of the initial modes through the generalised eigenvalue equation, as in (14). The initial state |ψ = |φ A , φ B goes through a balanced beam splitter U BS , gains a relative phase ϕ when traveling through the two arms of the interferometer, and finally enters another beam splitter; the final state then reads with whereã j ,b j are the output modes. Using this transformation we easily obtain, On the other hand, |ψ ϕ can be explicitly written as: Finally, it will be useful to note thaṫ A. Quantum Fisher Information for pure multimode states When |ψ = e −iϕH in (28), then we have the convenient expression for the QFI: Q = 4 ψ ϕ |H 2 |ψ ϕ − ( ψ ϕ |H|ψ ϕ ) 2 . Using (31) and (34), we obtain where we defined the average photon numbers n j ≡ φ A |a † j a j |φ A and m j ≡ φ B |b † j b j |φ B . In the particular case n i = m i , e.g. for twin states, we finally obtain From this expression one can immediately recover the QFI of twin Fock states Q = N (1 + N/2) [39], which corresponds to d = 1 and n 1 = N/2. For twin superradiant states, from our considerations of Sec. IV B we have that 2n 1 ≈ 0.90N , 2n 2 ≈ 0.08N and 2n 3 ≈ 0.02N , from which we obtain Q ≈ 0.41N 2 + N , hence recovering our previous results [20]. More generally, (38) provides a simple and clear expression for the potential of a particular multimode state for optical interferometry, and from it we learn that 1. To obtain the QFI of a twin multimode state with N/2 photons in each arm, it is enough to compute the average photon number of the internal modes of each arm.
2. When the multimode state is generated through a nonlinear decay as described in Section II, then our techniques enable us to compute the QFI for large photon number (N up to N ≈ 1000).
3. Heinseberg scaling (i.e. Q ∝ N 2 ) is possible when the number of relevant modes is independent of N , and quantum-enhanced scaling (i.e. Q ∝ N 1+α with α > 0) when the number of relevant modes grows sublinearly with N .

B. Number resolved measurements and QFI
Although the QFI provides the maximal sensitivity of |ψ ϕ to ϕ, it is equally important to understand how to achieve it in practice. For that, in this section we consider number-resolved measurements in both outputs of the interferometer. In particular, we consider two types of measurements • Number-Resolved (NR) measurements. We also consider standard photon counting measurements that are not able to distinguish between the different modes. In this case the corresponding probability distribution to obtain n and m photons in each output of the interferometer reads In both cases, we assume that the detector frequencybandwidth is larger than the photonic wavepackets linewidth, which in the superradiant case scale as γ ∝ N/ ln N (see Sec. IV A and Fig. 3).
3. n ± m photons in one output and m ∓ m photons in the other one, with m ≥ 2. Then we have, P αβ (ϕ) = O(ϕ m+1 ), which again does not contribute to (26).
Putting together these considerations we obtain: where in the third line we used that |ψ ϕ=0 has only support in the subspace of n ± 1 photons in one arm and m ∓ 1 in the FIG. 9. Classical Fisher information for MNR measurements (in orange) and for NR measurements (in blue) when two copies of |ψ1 in (44) (case (a)) or two copies of |ψ2 in (45)  other one, and in the fourth line we used ψ ϕ=0 |ψ ϕ=0 = 0. Hence, we conclude that C = Q around ϕ = 0 for MNR measurements. Crucially, the derivation (43) follows analogously for NR measurements, i.e., for the coarse-grained distribution (41). This has the very important consequence for practical implementations that experimentally one does not need to distinguish between the different modes to saturate the QFI. These conclusions hold around ϕ = 0, but there is in principle no reason why they should also hold for other ϕ. To address this point, we consider two illustrative states: and compute the CFI for the corresponding twin states (i.e. |ψ 1 , ψ 1 or |ψ 2 , ψ 2 are the input states of the MZ interferometer). In Figure 9 we plot the CFI for (44) and (45) with n = 5 and given MNR and NR photon measurements. We observe that C = Q for MNR measurements whereas for the NR measurement the equality is only saturated around ϕ = 0 (or multiples of π/2). It is worth stressing that one can always add phase shifters to compensate for ϕ = 0 during the estimation process in order to guarantee that C = Q for NR measurements. These are crucial considerations to take into account in implementations of quantum metrology with multimode states. It is also worth noticing that the multimode structure of the state leads to a non-trivial change of the CFI for NR measurements, as illustrated in Figure 9. Finally, we also compare NR and MNR measurements for twin-superradiant states as shown in Figure 10, which also illustrates the optimality of choosing ϕ ≈ 0. To obtain the results of Figure 10, we have used the effective descriptions obtained in Section IV D which enable us to describe superradiant states through effective descriptions using a low number of modes (here we use two-mode descriptions, d = 2, for simplicity).
While our techniques enable us to compute the metrological properties up to hundreds of photons, we note that current experimental implementations of photon-number detectors are limited to few photons resolution (N ≤ 10), see e.g. [45,46]. It is also worth mentioning that recent theoretical proposals show that atoms coupled to the waveguide could also be used for photon detection up to considerably larger photon numbers [47][48][49], which hints to the exciting possibility that atoms suitably coupled to waveguides can be both used to generate and detect photons. In the next section we discuss a proposal for quantum metrology with superradiant states in the presence of photon loss for realistic photon numbers, namely N ≤ 10.

C. Photon loss in the interferometer and the measurement device
Besides finding specific measurement schemes to saturate Q, in practice it is also crucial to consider imperfections in the interferometer and in the measurement devices. In fact, quantum-enhancements in metrology are largely affected by imperfections (either in the interferometer or due to imperfect measurements), as photon loss prevents Heisenberg scaling for sufficiently large N [50][51][52][53]. To intuitively understand why Heisenberg scaling is lost, note that for obtaining C = Q in (43) requires photon-measurement detectors that are capable of distinguishing a single photon (i.e. between N photons and N ± 1 photons in the outcomes of the MZ interfereometer); yet, in the presence of any finite photon loss (see Fig.  8), this is no longer possible when ηN 1, where η is the probability of losing a photon. Still, quantum enhancements in the presence of photon loss can appear as a better prefactor in Q ∝ N [51][52][53] as classical schemes are limited by the shot-noise limit Q ≤ N . This advantage however highly depends on the state into consideration: for example, GHZ states, which are optimal in ideal conditions, are known to quickly lose any sensitivity to ϕ in the presence photon loss. Twin Fock States (TFS), on the other hand, are known to be a good candidate for quantum metrology even in the presence of photon loss η (quantifying the probability of losing a photon in each arm of the interferometer), as in this case the QFI becomes which is half of the optimal one in the limit of large N [51,54]. It is also important to stress that although these considerations (e.g. (46)) are derived in the asymptotic limit (ηN 1), they provide valuable insights already for moderate N [51,54]. To extend (46) and the general considerations of [50][51][52][53] to multimode states is certainly an interesting but also challenging endeavour, as the multimode structure of the state makes it difficult to diagonalise it to be able to compute Q. In this Section, we instead pursue a more humble goal: We compare (twin) superradiant and Fock states for some specific case-studies of MZ interferometery and find that they perform similarly even in the presence of photon loss (which is expected as superradiant states contain 0.9N photons in a single mode).
We consider photon loss by adding a beam splitter with transmitivity η before each measurement apparatus (this is equivalent to placing the beam splitters before the second beam splitter of the MZ interferometer as the losses are symmetric). This is implemented by adding orthogonal modes {e j } D j=1 , {f j } D j=1 and implementing the transforma-tionsŨ BS : and We again characterise the input superradiant states through the effective descriptions obtained in Section IV D with d = 2. These effective descriptions enable us to easily account for photon loss, which would be highly challenging through the continuous descriptions (6). Still, our results are limited to low N , both because we can only obtain the coefficients of the state for N < 10 (see IV D), and also because the possible outcomes of the experiment (i.e. the size of the probability distribution P (MNR) α,β (ϕ)) grows as O(N 4 ) when dealing with two-mode states in each arm of the interferometer, hence making it difficult to compute the analytical expression (26) for high N .
In Fig. 11 we show C in the presence of photon loss in the measurement devices for different configurations: twin Fock states of N = 8 (orange) and N = 10 (blue) photons and NR measurements, and twin superradiant states with N = 10 and NR (green) and NMR (red) measurements. We observe how twin superradiant states perform close to Fock states, and how NR measurements for superradiant (and hence multimode) states become optimal around ϕ = 0, as expected from our previous considerations. It is worth pointing out that as η increases, the optimal value moves away from ϕ = 0, which happens for NR and NMR measurements, and for single-mode and multimode states. These observations are confirmed by further numerical results for other values of N and η, which are shown in Appendix II. From these results, we conclude that twin-superradiant states behave fairly similar to twin-Fock states in terms of their metrological performance also in the presence of photon loss; this conclusion is indeed expected given that superradiant states contain ≈ 0.9N photons in a single mode.

VI. CONCLUSIONS
To sum up, we have developed the theoretical tools to characterize (i.e., compute observables) of a wide class of multimode photonic states coming from the emission of a general non-linear level structure. Besides, we provide a construc-tive way of capturing the properties of these multimodal states with few-mode descriptions. To illustrate the potential of these tools, we have applied them to the case of superradiant photonic states, showing, for example, their observables can be captured efficiently already with two or three modes up to a large number of photons. Finally, we applied these ideas to a phase estimation proposal based on twin superradiant states and number resolved measurements. Our results suggest that twin superradiant states of N photons are a promising candidate for quantum metrology, as they perform approximately as a twin Fock state with ≈ 0.9N photons, which are in fact the number of photons contained in a single mode. The crucial difference is that twin superradiant states can be generated in a deterministic and scalable manner (assuming no photon loss and perfect control on the system's Hamiltonian to ensure the collective decay), in contrast to standard probabilistic methods to create Fock states, whose success probability decays exponentially with N even in ideal conditions. We hope these ideas motivate its experimental implementation in nanophotonic waveguides coupled to atoms [21][22][23][28][29][30], artificial emitters [24,26,27], or molecules [25].
Our results in quantum optical interferometry are in fact general and can be applied to arbitrary multimode states with a fixed photon number. Indeed, we have considered phase estimation in a Mach-Zender interferometer and derived a simple expression for the quantum Fisher information obtained when the input state of each arm of the interferometer is a generic n-photon d-mode state, and shown that it can be saturated by number-resolved measurements that cannot distinguish between the different modes. This is a crucial observation for experiments, as it shows that single-mode proposals for quantum metrology, e.g. [39,40,44,51,[54][55][56], can be naturally extended to the multimode regime without requiring extra resources in terms of the measurement devices.

APPENDIX I. RECURRENCE RELATIONS
Here we build recurrence relations to compute As the derivation is rather non-trivial, for clarity we will start by computing simple but relevant cases where n = 0, 1, 2 before deriving a general relation for (1). Along the derivation we will use that from (8) it follows that and similarly

A. Normalisation
It is instructive to first check that |φ (N ) is normalised. We want to compute, Using the symmetry of A t1...t N under permutations over {t j } and the commutation relation (8), we arrive at In order to solve this integral, which includes a time-ordering operation T , we split the integral as a sum of integrals using There are N ! such integrals, and they are equivalent. Hence we have that This integral can be easily worked out using ∞ s dse −at = e −as /a, which leads to the desired result B. Average photon number We first consider f (x 1 , y 1 ,x 1 ,ỹ 1 ). Before proceeding to its calculation, first note Using (8) and (2) and the symmetry of A t1...t N over permutations, we first have with This integral is challenging to compute because of the time-ordering T . We will compute the integral through a recurrence relation. The rough idea is as follows: first one splits the integral as a sum of integrals using (6). Since the number of integrals increases exponentially with N , it is crucial to use that when any of the t j 's is integrated, the resulting integral is the same due to the symmetry of A under permutations. This allows us to keep the computation efficient, as the number of integrals grows linearly with N . Let us implement this idea by developing a recurrence relation where each step corresponds to an integration through one of the variables of integration (u, v, t 1 , t 2 ...). Let us start by defining the integrals with j = 1, ..., N . Note that We can then compute F (N −1) 1,1 by noting the following recurrence relation (with j = 1, ..., N ): which provides the desired result. When computing this numerically, it is convenient to compute instead which directly provides I (N ) 1 (x 1 , y 1 ,x 1 ,ỹ 1 ) in (14).

D. Higher order terms
Given these previous considerations, it is in principle not difficult (but quite tedious) to extend these techniques to higher-order correlators of the form φ (N ) |b x1,y1 ...b xn,yn b † x1,ỹ1 ...b † xn,ỹn |φ (N ) . Essentially, following the previous considerations we need to compute integrals of the form In complete analogy with the previous considerations, we can define and the following recurrence relation can be derived (a natural extension of (25)), and c (N −k−n+1) This recurrence relation (33) can be computed through a matrix multiplication of M O(N ) in analogy with the previous sections, where M is now a matrix of size (2 2n (N − n + 1)) 2 . Hence we notice that the complexity of the calculation grows exponentially with the order of the correlator. where we used (16). As in the previous section, we can compute each integral by solving the following recurrence relation: where θ[x] is the step function (θ[−|x|] = 0 and θ[|x|] = 1), together with the initial condition and the coefficients This provides the desired solution as: Following the same logic as in the following sections, this integral can be computed by a product of matrices in a space of dimension which is approximately bounded as In Figure 12 we illustrate these ideas by computing all 0|c k1 1 c k2 2 c k3 3 |φ (N ) with N = 8, D = 8 and d = 3.

II. NUMERICAL RESULTS ON QUANTUM METROLOGY WITH PHOTON LOSS
This section shows more numerical results on the CFI with photon loss and considering as input states twin Fock states and twin superradiant states, as shown in Figures 13 ,14 and 15.