Ab initio quantum models for thin-film x-ray cavity QED

We develop two ab initio quantum approaches to thin-film x-ray cavity quantum electrodynamics with spectrally narrow x-ray resonances, such as those provided by M\"ossbauer nuclei. The first method is based on a few-mode description of the cavity, and promotes and extends existing phenomenological few-mode models to an ab initio theory. The second approach uses analytically-known Green's functions to model the system. The two approaches not only enable one to ab initio derive the effective few-level scheme representing the cavity and the nuclei in the low-excitation regime, but also provide a direct avenue for studies at higher excitation, involving non-linear or quantum phenomena. The ab initio character of our approaches further enables direct optimizations of the cavity structure and thus of the photonic environment of the nuclei, to tailor the effective quantum optical level scheme towards particular applications. To illustrate the power of the ab initio approaches, we extend the established quantum optical modeling to resonant cavity layers of arbitrary thickness, which is essential to achieve quantitative agreement for cavities used in recent experiments. Further, we consider multi-layer cavities featuring electromagnetically induced transparency, derive their quantum optical few-level systems ab initio, and identify the origin of discrepancies in the modeling found previously using phenomenological approaches as arising from cavity field gradients across the resonant layers.

While progress has largely been enabled by the "source driven revolution" [9] in the x-ray regime, it is likely that it will not be sufficient for the establishment of more general x-ray quantum optical methods including non-linear and correlated quantum phenomena. Instead, the development of novel control techniques, target optimizations as well as novel physical platforms will be paramount [12].
In particular, suitable photonic environments for the nuclei enable one to engineer enhanced couplings between x-rays and nuclei, to implement more versatile nuclear level schemes, and to simulate otherwise unavailable coherent control fields.
In this context, thin-film cavities with resonant nuclei embedded in the guiding layers are particularly promising [3,11]. Such cavities (see Fig. 1 for an example) have already facilitated the experimental observation of a multitude of quantum optical phenomena, including the collective Lamb shift [20], an electromagnetically-induced transparency (EIT) mechanism without an externally applied control field [21], and spontaneously generated nuclear coherences [22]. Certain schemes have been found to provide access to new nuclear and x-ray observables, such as interferometric phase detection via Fano resonances [23]. While the coupling between x-rays and single nuclei is weak compared to the rate of cavity leakage, the observation of collective strong coupling [25] and Rabi oscillations between nuclear ensembles [26] have recently been reported, and the experimental realization of slow light [24] further hints at the possibility of quantum systems with strong non-linearities [42] even at low light levels [43]. A key factor contributing to the success of the x-ray cavities is the fact that in the low-excitation regime, the total system of cavity and nuclear ensembles can equivalently be reinterpreted as an effective nuclear few-level scheme, with properties going beyond what is available in the bare nuclei [21,[44][45][46]. We note that while thin-film cavities have become particularly popular in the nuclear resonance scattering community, they also facilitate x-ray quantum optics with electronic resonances [47].
Together, these advances establish the field of nuclear cavity quantum electrodynamics (QED) with hard xrays, opening various avenues to connect to the impres- x-ray cavity QED. A side-on view of a thin-film layer cavity, typically consisting of cladding layers (red) and a guiding layer (green) doped with thin resonant layers (gray), is shown. In the example in the figure, the resonant layers contain atoms or nuclei (black) featuring ultra-narrow transitions in the hard x-ray range, such as those provided by Mössbauer nuclei (see zoom, for the example of 57 Fe, ω ≈ 14.4 keV and γ ≈ 4.7 neV). The system is probed spectroscopically at grazing incidence (wave vector k with parallel component k , perpendicular component k ⊥ and incidence angle θ) by x-radiation. Typical observables include off-resonant cavity reflection spectra as a function of incidence angle and resonant nuclear spectra (sketched plots).
sive quantum optical toolkit in the optical and microwave regime that is available throughout various platforms of resonator [48][49][50][51][52][53][54] and waveguide [55,56] QED. On the theoretical side, the nuclear resonant scattering community has largely employed semi-classical or mean-field methods based on perturbative scattering theory for the various transitions [57][58][59][60], including variants of linear dispersion theory [61,62] known as the layer [11,63,64] or Parratt's [65] formalism, Shvyd'ko's time and space picture [66] as well as Maxwell-Bloch equation treatments [28,66]. These approaches enable one to accurately model the related experiments, which so far essentially operate in the weak-excitation regime. However, they are not well suited to interpret the resulting spectra in terms of quantum optical phenomena. Motivated by this, a recently developed quantum model for the nuclei-cavity dynamics based on the input-output formalism [44,45] has crucially provided the possibility to interpret the observed features quantum optically. By expressing the dynamical scattering process in terms of effective quantum optical level few-level schemes, this theoretical approach has enabled the identification of various phenomena in the linear regime [22-24, 26, 45]. Beyond that, it may be able to serve as a predictive theory for non-linear and quantum phenomena at higher intensities [12].
Despite the practical success of the phenomenological input-output model for thin-film x-ray cavities with Mössbauer nuclei [44,45], there are various open questions connected to its applicability. In the multi-mode and multi-ensemble case [45], which is for example crucial to describe the EIT phenomenon observed in [21], two heuristic extensions of the model, including a dispersion phase and an envelope factor, are required to fit the cavity spectra. While the resulting model is able to reproduce the main EIT phenomenon, quantitative and qualitative deviations are found when the spectra are scanned against incidence angle [45]. In the linear regime, where the spectral observables can be modeled rigorously and with excellent empirical agreement by established semi-classical and scattering theory methods [11,60,63,64,67], these issues mainly limit the interpretation in terms of quantum optical phenomena. Beyond semi-classics and at higher excitation, however, these observations cast doubts on to what extent such phenomenological quantum models can be used as a predictive tool at larger driving field intensities, where no experimental results are available at the moment.
Here, we develop two ab initio quantum approaches to thin-film x-ray cavity QED. Throughout the analysis, we focus on cavities containing Mössbauer nuclei, but the general results equally well apply to cavities with electronic resonances. Our results not only resolve the previous discrepancies between the quantum optical model and the semi-classical descriptions, but also establish new theoretical tools that can be applied in the fully quantum sector and at higher intensities. The ab initio character of our approaches further enables direct optimizations of the cavity structure and thus of the photonic environment of the nuclei, by providing a predictive theory in which the involved approximations are well under control. As one application, this opens an avenue towards the design of effective quantum optical level structures including the control fields realized by the cavity.
Our first approach constitutes an ab initio version of the phenomenological few-mode model [44,45], based on our recent general framework for few-mode approaches to quantum scattering problems [74]. It removes the need for heuristic extensions and to fit the quantum optical parameters. Instead, the ab initio approach directly determines the parameters from the cavity geometry, and improves the prediction of scattering spectra to essentially perfect agreement. We further provide insight into the approximations involved in the model, providing a solid foundation for its application in the linear sector and opening the path for applications beyond this regime. At the same time, this progress promotes the quantum optical approach of constructing few-level models for xray cavity QED systems to an essentially exact and theoretically well-founded method.
Our second approach is based on a well-known Green's function technique [70][71][72][73], which in our case allows to derive a Markovian master equation for the nuclei coupling via the cavity environment directly from the cavity  Overview of theoretical approaches to thin-film x-ray cavity QED with Mössbauer nuclei, including the standard semi-classical layer formalism [11], which is based on a variant of linear dispersion theory [61,62,65], and a successful phenomenological few-mode model [44,45], which is based on the quantum input-output formalism of cavity QED [68]. In this paper, we present two additional ab initio quantum approaches to the problem, including an application of the ab initio fewmode theory [69] and of established Green's function techniques [70][71][72][73] to the layer geometry and the nuclear quantum optics setting. Besides completing the connections between the theories illustrated in this diagram and showing which approximations are involved in each case, the main results presented in this paper are indicated by the numbered magenta arrows. structure. This method provides an alternative to the ab initio few-mode approach, involving different approximations. Most importantly, the Green's function formalism provides a numerically efficient way to calculate effective nuclear level schemes.
Next to these general results, the examples we use to verify the two ab initio methods and to showcase their capabilities are chosen in such a way that they at the same time substantially advance the understanding of nuclear cavity QED, mainly in three respects. First, we demonstrate that the ab initio methods are capable of modeling the cavity-nuclei system not only in the vicinity of a chosen resonance, but across the entire range of probing photon frequencies and incidence angles. This opens perspectives for new applications beyond the current focus on individual cavity resonances. Second, we generalize the presently established quantum optical approaches to the practically relevant case of arbitrarily thick resonant layers. This is mandatory to achieve quantitative agreement for standard cavities used in recent experiments. Interestingly, we find that in the thick-layer limit the system is no longer equivalent to an effective nuclear few-level scheme, but rather to a set of self-coupled con-tinuum ensembles. Third, we resolve previously not understood discrepancies in the quantum optical modelling of multi-layer cavities featuring EIT spectra, and show that they originate from spatial field gradients across the resonant layers.
Finally, we provide a clarifying guide and overview of the existing and here presented theoretical approaches to thin-film x-ray cavity QED in Fig. 2.
The paper is structured as follows. Section II introduces and discusses the few-mode ab initio approach. As a starting point and for later reference, we summarize the main features of the phenomenological quantum optical model as presented in [44,45] in Sec. II A. In Sec. II B, we develop the corresponding ab initio theory. The approach is based on our recently developed method [69] connecting the input-output formalism to quantum potential scattering theory. In Sec. II C, we provide a general solution for linear scattering observables of the resulting ab initio model. Section II D disentangles the effects of the empty cavity and the contributions of the nuclei to the observables, which is exemplified in Sec. II E for well-studied cavities featuring Fano resonances. The effective nuclear level scheme represented by the cavity is derived in Sec. II F by eliminating the cavity degrees of freedom. In Sec. II G, we verify and illustrate the method for an analytically solvable example geometry. As an application, we generalize thin-film x-ray cavity QED to cavities containing thick resonant layers in Sec. II H.
The following Sec. III develops and discusses the Green's function approach as a second ab initio method. After introducing the relevant concepts of Green's functions and macroscopic QED in Sec. III A, we use the linear dispersion theory to derive an effective wave equation in Sec. III B that demonstrates the relation of this quantum theory to the established semi-classical approaches [11,63]. Section III C derives the nuclear Master equation using the Green's-function method. As a main result and advantage of this method, we show how to derive effective nuclear level schemes in Sec. III D, providing a numerically efficient method. Sections III E-III G discuss the general solution in the linear sector and further analyze the numerical efficiency. Finally, in Sec. III H, we apply the Green's-function approach to multi-layer cavity structures in the EIT configuration, as studied in [21,45], and identify field gradients across the layers of finite thickness as the origin of previously not understood inaccuracies of the phenomenological quantum optical models [45]. We further provide the first ab initio calculation of the effective nuclear level scheme in the cavity environment. Our method removes the need for an ambiguous fitting procedure and opens design opportunities for more complex cavity geometries.
Finally, Sec. IV discusses and summarizes our results. Details on derivations and parts of the calculations, as well as an in-depth practical comparison of the phenomenological and ab initio few-mode model are provided in the Appendices. In this Section, we briefly summarize the essential features of the phenomenological quantum optical model for thin-film cavities, which are doped with layers of Mössbauer nuclei and probed spectroscopically in grazing incidence by hard x-radiation [11,22], as illustrated in Fig. 1. Such setups have been studied in a multitude of experiments, featuring a variety of different cavities (see [20][21][22][23][24][25][26]47] for the most recent quantum optics experiments). For simplicity, we restrict ourselves to the unpolarized and unmagnetized case, which already features the essential ingredients. With polarization and magnetization included, one only obtains additional coupled equations [44,45].
The phenomenological model is given by the following equation for the density operator ρ characterizing theb The coupling constants in the model include the system-bath coupling ( √ 2κ R,λ ), the cavity mode decay constants (κ λ ) and the mode-nucleus coupling (g λln for nucleus n in ensemble/layer l). See also Table I for the functional dependencies of each model parameter. The input-output scattering matrix Sio describes the scattering process between bath modes via the cavity interacting with the nuclei, and is assumed to provide the full scattering information for x-rays being reflected of the system. For certain cavities, heuristic extensions of the model are necessary to provide more accurate descriptions of the cavity and nuclear response (see text).
cavity and the nuclei [45], The Hamiltonian contribution, consists of a cavity part H cav , a nuclear part H nuc , their interaction H int and an external driving term H drive . The two Lindblad terms describe the incoherent cavity decay and the spontaneous emission of the nuclei. The model is illustrated in Fig. 3. We work in natural units with = c = 1. The dynamics of the empty cavity in the interaction picture is given by [45] whereâ λ is the bosonic cavity mode operator with index λ,b in is an operator characterizing the driving field  Fig. 3, the modes feature additional Lamb shifts (blue arrows at κ11, κ22) and cross-mode couplings (κ12). Instead of a single drive, the bath is now an angle-dependent continuum with systembath coupling constant W λm (corresponding to √ 2κ R,λ in the phenomenological model). The input-output matrix Sio still contains the bath mode scattering information, however, to obtain the plane wave scattering matrix S, one has to employ an additional background scattering contribution (see text). Again, Table I lists the functional dependencies of each model parameter.
and 2κ R,λ is the drive coupling to mode λ. ∆ C is the cavity mode detuning relative to the external field with frequency ω. It depends on the incidence angle θ and is parametrized in terms of the resonance angles θ λ of the modes at the nuclear resonance frequency ω nuc , which is used as a reference. The detuning dependence on the external driving frequency can usually be neglected [22] since the nuclear line width is typically orders of magnitude smaller than the cavity resonance width (see also Sec. II D). Interactions with nuclei in the cavity are included via the Hamiltonian contributions [45] whereσ z = |e e| − |g g| is the Pauli z-operator of a single nucleus, and |e (|g ) is the corresponding excited (ground) state. Further, ω nuc,l is the nuclear transition frequency of ensemble l. The indices are chosen such that l corresponds to the different resonant layers or ensembles and n runs over all the N l nuclei within one such layer. g λln is the coupling constant of a nucleus to the cavity mode. The coupling is assumed to have the form g λln = g λl e iφn [45], including a position dependent phase φ n due to propagation along the waveguide. We note that since the prefactor is independent of the nucleus index n within one layer l, the above form includes the thin-layer approximation. That is it is assumed that the mode profile does not vary significantly across a layer that is considered as one ensemble, such that describing thick resonant layers requires the inclusion of multiple sub-ensembles dividing the thick layer into thinner ones (see Secs. II H, III H for an illustration and further details).
In addition to these Hamiltonian terms, the model features the incoherent contributions [45] where L cav [ρ] describes the cavity decay and L SE [ρ] the single nucleus incoherent decay. Here, ρ is the density matrix of the system and {·, ·} is the anti-commutator. In order to calculate output operators from the driving input and the cavity dynamics, the input-output relation [45] is invoked. Spectroscopic observables can then be calculated from the output operators, such as the reflection coefficient given by [45] Within the adiabatic approximation [44,45], the mode operator dynamics are given by [45] a λ = 2κ R,λbin − i ln g λlnσ − ln The nuclear dynamics can then be obtained from a nucleionly Master equation obtained by adiabatically eliminating the cavity modes in Eq. (1). The phenomenological few-mode model summarized so far has been employed for the analysis of various experiments in x-ray cavity QED [22-24, 26, 47], providing a basis for the quantum mechanical interpretation of the system.
However, for certain cavities, it was found that heuristic extensions to the model are necessary. These include a dispersion as well as an envelope factor resulting in the modified reflection coefficient [45] where r disp = |r disp |e iφ C is the dispersion factor that is introduced to account for a dispersion phase φ C as well as Pheno model Ab initio model Mode resonance ∆ C,λ (ω, θ)∆ C,λ (ω, θ) Cavity loss κ λ κ λλ (ω, θ) Bath/drive coupling 2κ R/T,λ W λm (ω, θ) Mode-nucleus coupling g λl g λl (θ) Background scattering -S bg (ω, θ) Heuristic extensions r disp , renv(θ) -TABLE I. Overview of the differences between the phenomenological model for thin-film x-ray cavities [44,45] and its ab initio counterpart developed in this paper. The quantum optical parameters (left) and their functional dependencies in each case (middle and right) are shown.
for far off-resonant modes and r env is the envelope factor that is introduced to account for total reflection behavior at grazing incidence. While with these heuristic extensions included, good agreement with the layer formalism and experiment is found in many cases (see e.g. [22-24, 26, 44, 47]), there are still quantitative and even qualitative differences for important setups, such as the EIT cavity [45] that is investigated experimentally in [21]. In addition, the heuristic extensions are neither justified nor derived from the quantum mechanical input-output theory [68,75], restricting the predictive power of the model. In our analysis below, we will remove these restrictions, clarify the relations between different theoretical approaches, and prove an accurate complete description of the system without the need for heuristic extensions (see Fig. 2 for an overview).
B. Ab initio few-mode theory for thin-film x-ray cavities In this section, we present an ab initio version of the phenomenological model for thin-film x-ray cavities, based on our recently developed general formalism for few-mode theories [69]. The ab initio few-mode theory for the given problem can be written in close analogy to the phenomenological version summarized in Sec. II A. The main differences are the functional dependencies of the coupling constants and the treatment of the inputoutput relation as illustrated in Figs. 3, 4 and Table I. Here, we outline these differences and define the structure of the ab initio theory. A detailed derivation is given in Appendix A. Unlike the phenomenological model presented above, we develop the most general form of the ab-initio model without invoking the Born-Markov approximation for the cavity decay or an adiabatic elimination of the cavity modes, in which case a standard Lindblad equation formulation is not sufficient.
Nevertheless, the Hamiltonian governing the dynamics is analogous to the phenomenological version in Eq. (2), The empty-cavity dynamics is governed by which generalizes Eqs. (3). Note that in the ab initio formulation, no explicit driving fieldâ in is specified. Rather, the coupling to all bath modesb m (ω) is included, with their free evolution term H field . For this reason, we also do not work in an interaction picture here, but with the bare frequencies such asω(ω, θ), which is an effective bath mode frequency at a given incidence angle or parallel wave vector (see Appendix A for details). For notational brevity, we omit the parametric angular dependence of the mode operators.
In comparison to the phenomenological model, we note that the angular and frequency dependence of the effective mode detuning∆ C,λ (ω, θ) may differ from the one assumed in the phenomenological model by ∆ C,λ (ω, θ) in Eq. (3c) (refer to Sec. II G and Fig. 6 for an example). Similarly, W λm (ω, θ) extends the in-coupling rate 2κ R/T,λ , where m is an index for the external channels, such as reflection and transmission denoted by R/T in the phenomenological model [22].
Interactions with the nuclei are included via the contributions which are essentially identical to their phenomenological counter-parts Eqs. (4), only now the mode-nucleus coupling is dependent on the incidence angle. The incoherent nuclear decay corresponding to Eq. (5) is still given by The cavity decay term in the phenomenological model L cav [ρ], however, has no direct counterpart in the ab initio few-mode theory, since the frequency dependence of the cavity-bath couplings induces non-Markovian dynamics of the cavity modes [76][77][78]. In particular for cavities with overlapping modes [69], which are realized in standard x-ray cavities as we will show in Sec. II G, this effect is not negligible and no Markovian Lindblad term for the cavity modes can be obtained. We note, however, that when tracing out the full cavity as a structured bath for the nuclei, the Markov approximation is applicable for many cavities due to the narrow line width of the nuclear response, which is the basis for the Green's function approach presented in Sec. III.
In the ab initio few-mode approach, instead of a Lindblad term for the cavity decay, one can write a frequency domain equation for the cavity mode operators [69] which is a direct generalization of Eq. (8) in the phenomenological approach going beyond the adiabatic and Markov approximations, with the effective mode detun-ing∆ C,λ (ω, θ) defined in Eq. (11c).
In comparison to the phenomenological model version in Eq. (8), the cavity loss rate is now replaced by a frequency and incidence-angle dependent cavity loss matrix κ λλ (ω, θ), corresponding to −iΓ λλ (ω, θ) in [69]. It has a matrix structure, as it includes cross-mode coupling terms, going beyond the single-mode cavity decay rates κ λ in the phenomenological model.
In order to calculate observables, the exact inputoutput theory from [69] allows one to obtain the frequency-domain input-output relation for the input [output] bath operatorsb generalizing Eq. (6) in the phenomenological model. In addition to the parametric dependencies of the input-output coupling, solving the input-output relation Eq. (16) does not yield the full scattering information, as shown in [69]. Instead, an additional background scattering contribution S bg (ω, θ) is necessary to translate the bath-mode scattering into plane-wave scattering (refer to Sec. II G and Fig. 7 for an example). This additional contribution is absent in the phenomenological model, but crucial to reproduce the response not only in the vicinity of the studied resonance, but also away from it. Most notably, it accounts for the fact that the empty-cavity response is treated exactly within the ab-initio few-mode approach [69] independent of the number of included cavity modes.
In summary, we see that the differences between the ab initio and the phenomenological version of the model are a modified frequency and angular dependence of the coupling constants, as well as the cavity decay rate becoming a matrix introducing cross-mode decay terms. We also note that the parameters may in general be complex quantities and a background scattering matrix should be included for accuracy over a large frequency range. The differences between the models are summarized in Table I and illustrated in Figs. 3, 4.

C. General solution in the linear regime
In this section, we provide a general solution to the multi-mode multi-ensemble ab initio few-mode theory for x-ray cavities with Mössbauer nuclei, which is outlined in the previous section, in the linear regime. By writing the equations in Heisenberg-Langevin form, we find a closedform solution for the linear scattering observables of systems containing an arbitrary number of modes and layers or ensembles. By employing the exact input-output formalism [69], multi-mode coupling, overlapping modes, non-Markovian, open system and scattering effects are included in the solution method.
For simplicity, we again consider a single polarization and unmagnetized samples. The method, however, can straightforwardly be applied to the full Hamiltonian of x-ray cavities with thin film Mössbauer nuclei, including polarization and magnetization [44].

Heisenberg-Langevin equations of motion
The Heisenberg-Langevin equations of motion for the nuclear operators defined by the ab initio model in Sec. II B read [45,68,79] We note that in this form, the equations are particularly simple, in the sense that the individual ensembles only interact with each other via the cavity modes. The manyparticle aspect apparent from the nuclear index n can be treated by introducing collective operators [22,45] or by working with single particle operators and employing the permutation symmetry of the model [79,80]. The latter can be employed since our Liouvillian is invariant under the relabeling (l, n) → (l, n ) for nuclei within each ensemble l, which is the case due to the coupling and decay constants being independent of n within the single parallel wave vector approximation. We adopt the permutation symmetry in the following, which then allows us to drop the index n on expectation values of single nucleus operators in Eq. (17) and to reduce the corresponding sums over nuclei within one ensemble. The resulting equation for the cavity mode expectation values is given by Eq. (14) and now simplifies to [69] where we have dropped the parametric angle dependences for brevity. We note that the number of nuclei N l in ensemble l appears explicitly in the equations due to the use of the permutation symmetry [79]. The dynamical equations (18) are completed by the input-output relation Eq. (16).

Solution of the multi-layer multi-mode equations in the linear regime
In [45], the linear equations of motion are solved for the phenomenological model within the adiabatic elimination approximation for the special case of a two-layer EIT configuration, which is notably investigated experimentally in [21]. Here, we solve the above ab initio equations in the linear regime without the need for adiabatic elimination, a Markov approximation or a specific mode-ensemble configuration. The adiabatic approximation will, however, be employed a posteriori to gain insight into the final solution and the resulting nuclear line shape.
The crucial approximation to obtain linear spectral observables is the weak excitation approximation, which can conveniently be performed in the Heisenberg-Langevin approach by setting σ z l ≈ −1 [81] for all ensembles l (see Sec. III B for further details). This results in the nuclear lowering operator equation of motion becoming linear and independent of the other nuclear operators. Physically, this way of solving the equations of motion corresponds to eliminating the nuclei instead of the modes, which can be done without further approximations in the linear regime (see also Sec. III B). An alternative way of performing this approximation is the Holstein-Primakoff method [82,83].
Within this low-excitation approximation, one obtains three linearly coupled matrix equations in the frequency domain, one for the nuclear lowering operators, one for the cavity modes and one given by the input-output relation. These equations can be solved straightforwardly by linear algebra. For the nuclear lowering operators we obtain Substituting Eq. (19) into Eq. (14), we find Thus, the nuclear ensembles effectively modify the cavity propagator matrix D. Finally, inserting Eq. (20) into the input-output relation Eq. (16), we obtain the scattering solution The scattering matrix, defined in matrix notation by [69] b thus reads As noted before, the input-output scattering is complemented by a background scattering contribution S bg (ω) [69], which translates the bath-mode scattering into plane wave scattering via Since the background scattering is independent of the nuclear ensembles, it can be computed for the empty cavity and multiplied with the input-output solution without additional effort when solving the interacting system [69]. These results provide a general solution for linear scattering observables in the ab initio version of the phenomenological few-mode model for thin-film x-ray cavities with Mössbauer nuclei.

D. Separating the contributions of cavity and nuclei
While the above Eq. (24) together with the background scattering contribution constitutes a complete solution of the problem, it does not allow for a straightforward separation of the individual contributions of the cavity scattering and the nuclei. To establish this separation, we note that the relevant frequency scale of the nuclei (e.g., ∼neV, quality factor of Q ∼ 10 12 for 57 Fe, see also Fig. 1) is typically orders of magnitude smaller than that of the cavity (∼keV).
In order to isolate the nuclear contribution, we approximate the frequency response of the cavity as constant on the scale of the nuclear response. If ω nuc is any one of the nuclear transition frequencies, we can then write That is in Eq. (24), the only relevant frequency dependence on the nuclear scale is the ω in the nuclear addition to D (int) in Eq. (21). This approximation essentially amounts to performing the adiabatic approximation [44,45] a posteriori on a spectral level and can also be considered as a type of Markov approximation [68,77,84], since it relies on the cavity's frequency response being flat on the scale of the nuclear response. We refer to Sec. II G and Fig. 9 for an explicit demonstration of the validity of this form of the approximation. In the experiments so far [20][21][22][23][24]26], the approximation was found to be possible due to the narrow nuclear response. We note, however, that when strong coupling effects are present [26] or for broader electronic resonances [47], the approximation may break down and it may not be possible to separate the cavity and the nuclear response.
We can now separate the nuclear contribution by employing the Woodbury matrix identity [85,86] to write the inverse of the modified cavity propagator as where the first addend is the (inverse) empty cavity propagator D −1 , such that the second term amounts to the nuclear contributions. The matrices appearing in Eq. (28) are defined by their components as whereg λl are the rescaled coupling constants of the collective states [45] (see also Appendix A 6 a). Inserting Eq. (28) into Eq. (24), the scattering matrix then reads This formula shows the desired separation into the empty-cavity background S [no nuclei] io (ω nuc ) and the nuclear contribution. The interpretation of the various contributions will become clearer in the following section discussing pertinent example cavities.

E. Example: Fano cavities
In this section, we illustrate Eq. (30) using standard cavity structures that are known to feature simple, yet relevant, nuclear line shapes in their frequency-dependent response. The basic element contributing to the response can be identified using the single-cavity-mode approximation, in which the inverse D −1 appearing in Eq. (30) can be identified via Eq. (27) as a Lorentzian. In the general multimode case, since g D −1 g † is a matrix coupling the different transitions, the nuclear contribution therefore is a superposition of coupled Lorentzian resonances interfering with themselves and the flat cavity background. As a result, the nuclear line shape can be rather complex for the general multi-mode multi-ensemble case. It is, however, well known that for certain cavity structures, the response reduces to a single Lorentzian interfering with the cavity background [23,87], giving rise to Fano line shapes [23,[88][89][90]. In the following, we use such cases to illustrate Eq. (30).

Single nuclear ensemble
For a single nuclear ensemble, the index l has only one value, such that Λ reduces to a scalar. Similarly, g and, for single channel scattering, W reduce to a vectors with elements corresponding to the different cavity modes λ. As a result, for the reflection channel, W † D −1 g † and g D −1 W become complex numbers independent of ω, thus characterizing the magnitude of the nuclear contribution as well as its relative phase to the empty-cavity contribution. For the spectral shape, it therefore remains to analyze where we have defined This is the Lorentzian line shape of a single nucleus, modified by the collective Lamb shift ∆ LS and the superradiance Γ S for the ensemble in the cavity. Due to the interference with the cavity background, this Lorentzian in general becomes a Fano line shape. Indeed, the above form of the nuclear line shape closely corresponds to what is obtained in the multimode phenomenological model (see in particular Eq. (46) in [45]). The main difference here is the form of the D-matrix, which is diagonal in the phenomenological case [45] and now includes the cross-mode coupling terms in the ab initio generalization (see Table I). In the special case of isolated cavity resonances [69], the D-matrix may well be approximated as diagonal, such that the phenomenological assumption is well justified.
We thus recover the corresponding predictions of the phenomenological model [23,45], at the same time establishing an improved ab initio representation of the complex level shift ∆ LS − iΓ S /2.

Single cavity mode, uniform nuclear ensembles
Another special case where the nuclear contribution is of single Lorentzian form arises when there is only a single mode and in addition the ω nuc,l of all nuclear ensembles are equal [23,87]. In this case, the index λ is absent, D reduces to a scalar, and W and g become column vectors with elements W m corresponding to the external channels and g l to the layers. As a result, we obtain from Eq. (24) Applying the scalar version of the Woodbury formula leads to Again, the result indeed corresponds to a Fano lineshape, recovering the phenomenological result [45] with A corresponding to G 2 (Eq. (48) in [45]). We note again, that the result is fully expressible in terms of the rescaled couplingsg l (see also Appendix A 6 a).

F. Effective few-level scheme realized by the nuclei in the cavity
As shown in [44,45], one can associate an effective nuclear level scheme to the nuclear cavity QED system, such that the spectroscopic response of the two agree in the low-excitation limit. This way, the spectra of the nuclear cavity QED system can be interpreted in terms of quantum optical phenomena. For instance, a spectral dip in certain two-ensemble x-ray cavities could be interpreted as an EIT phenomenon [21,45]. In [44,45], the effective level scheme is derived using an adiabatic elimination of the cavity modes, giving rise to equations for the nuclear system alone with effective couplings between the ensembles. The parameters of this level scheme are then obtained using fits of the model to numerical calculations of the QED system's response or to experimental data.
In this section, we derive an effective few-level scheme from the ab initio few-mode model given above and generalize it to arbitrary ensemble configurations. While the adiabatic elimination is not necessary to solve the linear scattering problem, it is still useful in order to obtain an effective level scheme. Due to the ab initio character of the theory, this approach further allows us to express the quantum optical level parameters and couplings directly in terms of the cavity geometry. Alleviating the need for parameter fitting not only removes possible ambiguities due to overfitting or degeneracies, but also amplifies the interpretational power in terms of the relevant cavity modes and nuclear ensembles responsible for the observed effects.
As a first step, we write Eq. (17) in terms of collective operators [87]Ĵ As a next step, we adiabatically eliminate the cavity modes by substituting Eq. (14) into Eq. (35) and employing the adiabatic approximation [44,45] on the spectral level Eqs. (26,27). The adiabatic solution for the mode operator Eq. (14) can then be transformed into the time domain to givê Substitution yields where the effective drive coupling matrix is given by and the effective level coupling matrix is These equations can be reformulated in terms of an effective Hamiltonian (Color online) Schematic illustration of effective level schemes equivalent to the total system of cavity and nuclei in the linear low-excitation regime. The structure of the level scheme is defined via the couplings given in Eqs. (40) and (41), Single nucleus decay terms are not depicted for clarity (see text). and the effective Lindblad operator where L SE [ρ] is the local spontaneous emission term defined in Eq. (A21). This effective level scheme again closely corresponds to the results in the phenomenological model (in particular Eqs. (23)(24)(25) in [45]), with the aforementioned differences in the ab initio coupling parameters summarized in Sec. II B. Eqs. (40) and (41) can be interpreted as a system of interacting collective spins. Together, they define an effective level scheme in their low-excitation subspace. Starting from the ground state |G defined as the state without photons in the cavity and with all nuclei in their ground state, the excited states of this level scheme are collective excitonic states of the many-body QED system, which can be defined as and the couplings within this effective level scheme are illustrated in Fig. 5. For example, for a single ensemble l the collective Lamb shift and the Purcell enhanced superradiance are given by respectively. With more than one ensemble, the coherent inter-ensemble coupling is given by ∆ ll = Re[G ll ] , and incoherent cross-ensemble decay terms with decay rate γ ll = −2Im[G ll ] are also present. The latter can be shown to yield spontaneously generated coherences in certain cases [22]. The various transitions are driven from an external channel m with the effective Rabi frequency Ω lm . The spontaneous-emission contribution in Eq. (41) requires additional care, because it a priori operates on the single-particle level, rather than on the level of the collective transition operators. But in the low-excitation-or more generally the linear regime, it simply adds to the collective decay rate [45], such that the system can fully be expressed in terms of collective operators as shown in Eqs. (37). However, Eqs. (40), (41) are also valid beyond the linear regime, where methods based on permutation symmetry [79,80] can be employed to tackle the spontaneous emission term to extend the collective operator treatment in [12]. We note, however, that the non-linear dynamics may break the permutation symmetry by coupling to other parallel wave vector modes.
The effective many-body couplings contain the bare coupling constants g lλ instead of the rescaled coupling constantsg lλ = √ N l g lλ , because the underlying equations contain collective operators defined as sums over single-particle operators. In the linear limit, those can be replaced by effective single-particle operators, which leads to the appearance of the rescaled coupling constants [12,45] (see also Sec. III D). In this sense, the linear limit is independent of the number of coherently interacting nuclei, see also Appendix A 5, as is also the case in the semi-classical layer formalism [11], where only the nuclear number density appears as a material parameter. However, beyond the linear excitation regime, the number of coherently interacting nuclei is crucial, since it determines the onset of non-linear effects [12].
We further note that by including the magnetic substructure of the nuclei, more versatile level schemes become accessible [22].

G. Analytically solvable example system
In this section, we illustrate the ab initio few-mode approach to nuclear many-body cavity QED using a specific example. We chose the simplest possible layer cavity geometry which features a resonance structure that resembles what is practically investigated in x-ray cavity QED [20][21][22][23][24]26], and also has found applications in various contexts, see e.g. [91][92][93][94][95]. It furthermore has the advantage that it can be analytically solved within our framework.
We present detailed comparisons to the semi-classical theory and to the phenomenological few-mode theory. In the process, we demonstrate the advantages and the potential of the method, showing differences to previous interpretations and completing the theoretical picture of the few-mode approach. As an application, we generalize the quantum optical model to cavities containing thick resonant layers.

Cavity layout
The cavity structure considered in the following is depicted in Fig. 6. We choose the resonant nucleus 57 Fe as the archetype Mössbauer nucleus and practically most used species for nuclear x-ray cavity QED experiments, and consider a single unmagnetized thin resonant layer inside the cavity, where collective Lamb shifts [20] and Fano line shapes [22] can be observed. As the off-resonant dielectric material for the guiding layer, we choose the corresponding isomer 56 Fe, which does not feature the nuclear resonance. Choosing two isomers has the advantage that the electronic contribution of the resonant layers is equal to the surrounding cavity material. Consequently, when the resonant layers are placed at varying positions, the electronic scattering properties of the cavity remain the same. This feature allows us to separate resonant effects, such as which modes the nuclear transitions couple to and which effective level scheme is realized, from cavity structure effects, such as which modal environment is realized in the cavity. For the substrate, we choose an idealized material with a high refractive index, which essentially mimics a perfect mirror. This choice is motivated by the fact that a single layer of dielectric alone usually features weak resonances. A reflecting mirror substrate is the simplest way to realize a better resonance structure.

Analytic solution
In order to employ the ab initio few-mode theory, one first has to choose a few-mode basis. While the optimal basis may be difficult to find, one can use a constructive approach by choosing a basis that becomes locally complete in the region of the interaction if infinitely many discrete modes are considered [69]. Note that different choices for the basis vary in the number of modes required to achieve convergence of the results.
Here, we choose Dirichlet modes, that is modes with hard wall boundary conditions at the top and bottom surface of the layer cavity, giving the transverse mode functions [69,86,97] where L is the cavity thickness. The few-mode basis is then constructed as a finite subset of these modes by choosing the relevant set of λ indices. By observing the convergence of the spectral observables, we select the significantly contributing modes. The corresponding one-dimensional problem for this set of few-mode bases has been solved in [69,86], giving the system-bath couplings and few-mode propagator compares results for the ab initio formalism, which equivalently can be obtained using the layer formalism [11] (solid blue), e.g., implemented in the software package pynuss [96], to fits using the phenomenological model. A five-mode model with heuristic extensions (dashed yellow) is necessary for good agreement in the spectral range up to 8 mrad, reproducing the reflectance over the range of four modes. Without heuristic extensions (solid green), the fit is poor already in the first mode and beyond 6 mrad, properly describing only the second and third modes. To fit a larger spectral range, the phenomenological model requires inclusion of more modes. This is a major drawback to the ab initio few-mode theory, where the off-resonant reflectance is always represented exactly for any mode number [69].
where L is the thickness of the cavity layer, ω 2 λ /2 = (π 2 λ 2 )/(2L 2 ) +Ṽ defines the system mode frequency, V (ω) = (1 − n 2 )ω 2 /2 is the energy dependent potential, n is the complex refractive index of the guiding layer material ( 56 Fe), β = ωL = kL is the scaled dimensionless momentum variable [86], k is the wave number in units of m −1 , α = (β 2 − 2Ṽ L 2 ) 1 /2 , s = λ∈Λ Q (2λ 2 π 2 )/(α 2 − λ 2 π 2 ), and λ is the system-mode index, which is summed over the chosen few-mode subset Λ Q of the locally complete Dirichlet basis. Here, we have written the hermitian conjugate of the system-bath coupling W explicitly, in order to demonstrate which quantities should be excluded from the complex conjugation operation (see Sec. A 5 for details). For example, α contains an imaginary part via the complex refractive index, but should not be conjugated within the non-hermitian Hamiltonian prescription.
Due to the perfect mirror substrate, the system reduces to a single-channel problem with reflection only, such that that the system-bath couplings W λ (ω) do not depend on a channel index m. The background scattering matrix then reduces to a scalar given by [86] At a given incidence angle or parallel wave vector, this solution can be applied to our x-ray system using as the index of refraction and scaled momentum variable, respectively, see Appendix A 1.
The interaction with the resonant nuclei is governed by the coupling constant where we have dropped the dependence on θ for brevity, d is the nuclear dipole moment and f LM is the Lamb-Mössbauer factor. Details on the light-matter coupling constants are summarized in Appendix A 6. These expressions for the quantum optical parameters can directly be substituted into the general solution in the linear regime derived in Sec. II C to obtain the observables which are investigated in the following.

Empty cavity
We start by analyzing the off-resonant reflectance of the cavity, i.e., in the absence of the nuclear resonances, which represents a property of the empty cavity. In practice, this quantity as function of the incidence angle is also known as the rocking curve and can be measured in a cavity containing resonant nuclei by tuning the incident x-ray beam slightly away from the nuclear resonance, which does not have any influence on the much larger scales of the cavity linewidths. Here, we study the reflectance as function of incidence angle and photon shows the full reflectance, which can be calculated via Parratt's formalism or equivalently via the ab initio few-mode theory. Note that the rocking curve in Fig. 6 is a section through this figure at fixed energy. The ab initio approach decomposes the total reflectance into an input-output and a background contribution, which are shown for a single mode theory (λ = 1) in panels (c) and (d), respectively. As expected, the input-output contribution extracts the first resonance of the spectrum. The corresponding phenomenological single mode theory (b) does not correctly reproduce the first resonance. This is most clearly visible from the position of the resonance in the phenomenological model, which is shown as a red line in (a) and only coincides at the energy 14.4 keV where the parameters of the model were fitted (white cross).
energy. We compare the results of the ab initio framework to corresponding calculations using the phenomenological input-output theory, and further use the semiclassical layer formalism [11] implemented in the software package pynuss [96] as a reference. The latter is a reimplementation of parts of the conuss software package [64] with substantial extensions focused on applications in nuclear quantum optics. The lower panel of Fig. 6 shows results for the rocking curve of the system, i.e. the incidence-angle dependent reflectance observed with incidence angle equal to the exit angle (θ − 2θ geometry). We find excellent agreement between the ab initio and the semi-classical approaches. Note that the different cavity mode resonances contributing to the rocking curve overlap strongly. Nevertheless, it is important to note that the ab initio few-mode theory is equivalent to the layer formalism calculation for any set of chosen system modes, a feature which is based on the exact result shown in [69] and which can be checked by comparing the analytic formulas given in Sec. II G 2 to the layer formalism result for the geometry [11]. This has the practical advantage that the few-mode approximation only affects the interaction with the nuclei, but not the empty-cavity properties. This way, different ap- On resonance, both the ab initio and the phenomenological approach yield excellent agreement, featuring deviations on the few-percent level and below. At higher incidence angles, however, the ab initio theory agrees significantly better. A more detailed comparison, including different approaches to fitting the phenomenological parameters, is presented in Appendix B.
proximations are disentangled, contributing to the understanding of the relevant processes.
For the phenomenological few-mode theory on the other hand, no such equivalence is guaranteed. Instead, the model parameters are fitted for a chosen mode number such that the rocking curves agree in a certain angular range. Thus, the few-mode approximation affects the interaction with the nuclei and the empty-cavity response. The model fits in the bottom panel of Fig. 6 shows that for a five-mode model with heuristic extensions [dashed yellow] (see Sec. II A for details), the agreement is good in the range of the four lowest resonance dips. Without heuristic extensions [solid green], however, the fit is poor already in the vicinity of the first cavity mode.
These results demonstrate a key advantage of the ab initio few-mode model, that the off-resonant cavity scattering and hence the rocking curve is treated exactly. The model further does not require heuristic extensions and therefore gives a clear picture of the quantum mechanical interpretation of the cavity resonances.
In order to investigate the origin of the difference between the phenomenological and the ab initio theory fur-ther, Fig. 7 shows a two-dimensional generalization of the rocking curve, where the incidence angle and energy are varied. We note that these spectra are calculated at a fixed off-resonant material refractive index for simplicity. These two-dimensional empty cavity spectra give theoretical insight into the resonance structure of the cavity, which forms the electromagnetic environment for the nuclei. The layer formalism [11] and the ab initio few-mode theory are again numerically equivalent, only in the latter, the scattering matrix is decomposed into a mode-number dependent input-output and a background scattering part [panel (a)]. For the single mode results (λ = 1) shown in the figure [panels (c,d)], we find that the input-output result captures the first resonance line across the whole two-dimensional space, including the total reflection cutoff at low incidence angles. The phenomenological single mode theory in (b), however, does not capture the resonance trajectory, except at 14.4 keV, where the model parameters are fitted to the one-dimensional rocking curve. The reason for this difference is the angular and energy dependence of the quantum optical parameters, in particular ∆ C (θ, ω), which differs between the ab initio and the phenomenological model (see Sec. II B). The phenomenological resonance trajectory can also be expressed analytically by solving ∆ C,λ (θ, ω) = 0 using the assumed form of the phenomenological cavity detuning Eq. (3c) giving The comparison of this phenomenological resonance trajectory (shown as a red line in Fig. 7) to the actual spectral minimum of the first resonance reveals the difference between the two descriptions.

Nuclear spectra and interacting systems
In Fig. 8, we turn to the interacting scattering problem. We consider the nuclear resonance spectrum of the cavity as a function of detuning and incidence angle (panel b), which is accessible experimentally [20][21][22][23]45]. Note that here, the scale of the detuning axis is on the order of neV (γ57 Fe ≈ 4.7 neV), in contrast to the keV scale of the driving radiation, which determines the empty cavity properties investigated in Sec. II G 3.
The figure compares the nuclear spectra for the phenomenological and ab initio few-mode approaches to the semi-classical layer formalism [11], which serves as a wellunderstood reference in the linear scattering regime. In the phenomenological case, we use a five mode model to fit the nuclear spectra. For comparability, the ab initio few-mode result is also calculated using five cavity modes of the Dirichlet basis (λ ∈ {1, 2, 3, 4, 5}).
In the fitting procedure for the phenomenological model parameters, we first fit the empty cavity parameters to the rocking curve, yielding the yellow dashed line in Fig. 6 as the best fit. This includes the heuristic dispersion phase, which is necessary to obtain a good result already in the empty cavity case (see Fig. 6). Due to the absence of a cladding, an envelope factor (see Sec. II A) is not necessary for this cavity. The mode-ensemble interaction parameters are then determined by fitting the nuclear spectrum, using the two-dimensional layer formalism result in Fig. 8(b) as the fit objective, without prioritizing a certain angular region. We note that other fitting procedures are possible, providing optimized representations for different spectral features. In Appendix B, three options, some of which were already suggested in [45], are investigated and compared. The above procedure is found to provide the best description of the spectrum across the whole angular range up to 8.5 mrad. In contrast, for the ab initio method, no fit is required, and a single model is adequate for the entire spectrum. Fig. 8(c) and (d) show the residual deviations of the ab initio few-mode result and the phenomenological result obtained using above fitting procedure, both calculated for five cavity modes. The residual deviation is defined as |R(θ, ∆)−R reference (θ, ∆)| with the layer formalism result as the reference. We find that while both approaches provide good fits across the whole range, at higher incidence angles the ab initio theory progressively performs better than the phenomenological approach. The excellent agreement of both theories at resonance with the first cavity mode is illustrated by Fig. 8(a), where the one-dimensional slice at θ = θ 0 is shown.
We conclude that both approaches are well suited to model the given cavity, with the ab initio approach yielding a better global description. In comparison to the phenomenological model, it provides the main advantage that due to the absence of a fitting procedure, the quantum optical interpretation is unambiguous. In addition, the model captures the off-resonant properties exactly for any mode number and can be systematically brought to convergence by including larger mode numbers. These features are important in particular when going to more complex cavity structures and nuclear level schemes.
Furthermore, increasing the number of modes provides a systematic way of improving the ab initio results. For 30 or more modes, the relative deviation is below 2% in the entire Fig. 8(c) and on this level is dominated by residual effects of the 0.5 nm layer thickness (see also Sec. II H). This systematic improvement is not guaranteed in the phenomenological model (see also Appendix B).
We note that these results do not invalidate the phenomenological approach, which may still provide numerical advantages for fitting unknown cavity structures, for example in the experiment, where the matrix elements required for the ab initio few-mode theory may be difficult to calculate (see, however, the Green's function approach in Sec. III for a numerically effective alternative). Instead, the applicability of the phenomenological approach is clarified and an upgraded version can be employed when necessary. This progress in the under- standing of the approach bears potential in particular for calculations beyond the linear regime [12], where the theoretical models are required to be predictive and hence well understood.
The investigation here is supplemented by a detailed comparison of the phenomenological and ab initio fewmode approaches in Appendix B, where different fitting routines and the convergence at higher mode numbers are considered.

Coupling constants and frequency dependence
In Fig. 9, we investigate the parametric dependencies of the coupling constants in the ab initio few-mode theory. In particular, we consider the nuclear complex level shift [86] given by F = ∆ LS − iΓ S /2 = g D −1 g † and the quantity F R = −2πiW † D −1 g † g D −1 W, which quantifies the nuclear resonance peak height in analogy to 2κ R for the cavity resonances [44]. Together, these two quantities contain the information to calculate the nuclear spectrum in Fig. 8 according to Eq. (30). Since we consider a single nuclear ensemble in a cavity with only the reflection channel being open, both of these quantities are scalars.
Results are shown in Fig. 9 as a function of photon frequency and incidence angle over the range of the cav-ity spectrum in Fig. 7. We note that all parameters vary significantly in the frequency and angular direction, featuring non-trivial functional dependencies that give rise to the practical differences to the phenomenological model (see Sec. II G 4). However, while there are significant variations in the energy direction on the keV scale, the parameters are essentially constant in the neV-µeV range around the nuclear resonant energy. This confirms that the adiabatic approximation performed in Sec. II D is valid for such cavities. The non-trivial angular dependence nevertheless persists. We note that since the cavity parameters are usually obtained by fitting to the rocking curve [45], they are particularly susceptible to the differences in the phenomenological description.

H. Generalization to thick layers of resonant nuclei
In this section, we develop the quantum optical description for x-ray cavities with thick resonant layers. The case of thick resonant layers is practically relevant as it has been studied extensively in nuclear forward scattering experiments [18,[98][99][100] and has been used in grazing incidence, for example, for novel spectroscopy techniques [94] or the study of magnetism in thin films, see, e.g. [92,95]. Placed inside cavities, thick resonant layers are particularly interesting from the theoretical perspective, since they are difficult to describe using a phenomenological approach. The reason is that in the latter, the mode functions are not known, such that one would have to fit a continuum of light-nucleus coupling parameters to describe the spatial variation throughout the thick layer.
Within the ab initio few-mode theory on the other hand, the description is rather straightforward and can be considered as the continuum limit of the general multilayer solution derived in Sec. II C. We further show that the thick layer system can be described as a self-coupled continuum ensemble, rather than the few-level schemes obtained in the thin-layer case.

Theory
With the general solution of multi-mode multi-layer systems available, the case of thick layers can be solved simply by partitioning the layer into finely spaced subensembles and taking the continuum limit of many such sub-ensembles.
The main structure-dependent quantity to compute in order to obtain the scattering matrix is which modifies the intra-cavity propagator according to Eq. (21). For a thick layer of resonant material, adjacent atomic layers of nuclei are usually separated by a distance much less than the variation in the mode profile. Dividing the thick layer into thin sub-ensembles, we can then take the continuum limit of the sub-ensemble sum, such that where ω nuc (z), ρ N (z) are the position dependent frequency and nuclear number density, respectively. For homogeneous resonant layers the z-dependence of the couplings is mainly determined by the mode profiles. z i [z f ] is the lower [upper] edge of the thick layer.
Evaluating the integral above gives the necessary quantity to compute scattering observables. Similarly, one can compute a continuum version of the effective level scheme given by the continuum limit of Eqs. (40,41) and the effective Lindblad term This yields a physical interpretation of the thick layer as a self-coupled continuum ensemble rather than an effective few-level scheme as in the thin layer case.

Thick resonant layer in the example cavity
As a practical example we again use the cavity from Fig. 6 and consider a variable thickness of the resonant layer, while keeping it at the center of the cavity.
Since the resonant layer material is homogeneous, the number density of resonant nuclei and the resonance frequency of the nuclear transition do not vary with z. In this case, denoting the thickness of the resonant layer by t res , we can evaluate Eq. (50) to obtain where the mode integral for the cavity under consideration is which can straightforwardly be evaluated analytically. Fig. 10 shows nuclear spectra of the example cavity at different layer thicknesses. The results from the continuum ab initio few-mode theory show excellent agreement with the semi-classical pynuss [96] calculation, confirming the validity of our quantum model for thick layer cavities.
Finally, to interpret the spectra, we note that the selfcoupled continuum ensemble features an EIT-like dip effect for thick layers. At layer thicknesses comparable to the cavity thickness, additional distortions are found, indicating the emergence of higher order spectral interferences.

III. GREEN'S-FUNCTION APPROACH
In this section, we apply a well-known Green's function technique [70][71][72][73] to the nuclear x-ray cavity QED system. This provides a second ab initio approach to the problem, which involves different approximations. The investigation completes the comparison of the three different ab initio approaches depicted in Fig. 2 and contributes as a numerically efficient method for calculating effective nuclear level schemes.
As a theoretical result, we show that the layer formalism is essentially analogous to an effective propagation equation that can be derived from the Green's function quantization by linear dispersion theory, up to applicable but unnecessary approximations that can be removed (such as the rotating wave approximation). Consequently, in the linear regime, all three ab initio approaches (layer formalism, ab initio few-mode theory, Green's function approach, see Fig. 2) are essentially equivalent in that they can be used as alternative methods to calculate linear scattering spectra. However, each of these methods provides different advantages with regards to the interpretation of the quantum system or the mode structure of the cavity due to the different underlying approximations. In addition, the methods presented in this paper go beyond the semi-classical linear dispersion theory in that they bear the potential for direct application in different sectors beyond the linear scattering and low excitation regimes.
The Green's function technique presented in this section most importantly provides a numerically efficient method to calculate effective nuclear level schemes for arbitrary complex layer stacks. It also removes the need for a fitting procedure that is necessary in the phenomenological model in [44,45], with the main advantage over the ab initio few-mode approach in Sec. II F being the straightforward numerical implementation even for complex systems. On the other hand, the Green's function approach facilitated here does not provide access to the dynamics of the cavity modes, which is naturally included in the ab initio few-mode approach. We demonstrate the validity and usefulness of the Green's function level scheme for practically relevant example systems, including the EIT and non-EIT cavities investigated in [21,45]. In the process, we resolve an open question with regards to the quantum optical interpretation of this system that was raised in its phenomenological model description [45], where the main EIT feature was reproduced, but qualitative and quantitative differences of the two-dimensional nuclear spectra were found.

A. Green's function and macroscopic QED
The quantization of absorbing dielectrics (see [72] for a review) has been studied extensively in the literature (see e.g. [70][71][72][101][102][103]). A macroscropic prescription based on the Green's function [71,72] gives the Hamiltonian in the dipole approximation [71,104] wheref (r, ω) are bosonic operators,Ê is the electric field operator and unlike in previous sections, polarization is included. The bosonic operators appearing in the Hamiltonian are related to the field operator by [71] where the Green's tensor is defined via and the dielectric permittivity is allowed to be frequency dependent. Approaches based on such Hamiltonians employing the classical electromagnetic Green's function of the system are known as macroscopic QED [72], which has been used extensively, for example, for the study of dispersion forces and related phenomena [103], as well as recently for the description of atom-light scattering in nanostructures [73,105] and atomic-waveguide QED [106]. The Green's function quantization provides an alternative to the normal modes approach in Appendix A that was used as the basis for the ab initio few-mode theory in Sec. II B. Besides its numerical efficiency that will be demonstrated later, it has the advantage that absorption is described rigorously without the need for the non-Hermitian Hamiltonian prescription in Sec. A 5. As a drawback, the approach does not offer an interpretation in terms of the modes or resonances of the cavity structure, which are often a driving force for the design of novel and exotic cavity structures, both in the x-ray regime [12,21,23,26] and at lower wavelengths [107]. Instead, all the information about the cavity environment is contained in a single function. In this context, we note that the Green's function can in principle be approximated in terms of mode parameters [73]. For overlapping modes structures, however, such an approximate treatment is non-trivial and results in similar problems as found in the phenomenological model. On the other hand, this provides a route towards direct numerical optimization of cavity structures via the Green's function approach (see e.g. [108]). The approach could also be complemented by various decompositions of the Green's function known from resonance theory [107,[109][110][111][112][113][114].
We further note that the macroscopic QED framework, that is treating the material as a refractive index, also has its limitations. Complementary approaches connecting to electronic structure theory [115] and quantum chemistry [116] are available for a variety of parameter regimes and, while being computationally demanding, allow to fully integrate associated effects.

B. Linear dispersion theory and relation to the layer formalism
The equations of motion for the Hamiltonian Eq. (55) are all linear, except for the usual non-linear term featuring a product ofσ z and the field operators (see Eq. (17) for the few-mode counterpart of this contribution). In [104], it is shown that analogously to the linear calculation in Sec. II C and the linear dispersion calculation in [62,69], these equations can be tackled by employing the low excitation approximation σ z (t) ≈ −1 and solving the resulting linearly coupled differential equations. The result can be expressed in frequency space as [104] Ê (r, ω) with [104] σ + ln (ω) = d l · Ê (r ln , ω) (ω + ω nuc,l ) , and Substitution into Eq. (58) and using Eq. (57) shows that the electric field obeys the effective wave equation where we have assumed real d for simplicity. The effect of the coupling to the nuclear transitions can thus be interpreted as a modification of the frequency dependent refractive index [69,117]. We note that in contrast to the formula resulting from the simple model in [69], the above approach includes polarization and a rigorous treatment of absorption. On the other hand, a factor of ω 2 nuc,l ω 2 is absent due to the E · r gauge that is adopted in [104], which is appropriate for our weak coupling scenario, but has to be adjusted at extreme coupling strengths [117][118][119].
We further note that by transforming the above effective wave equation into the time-domain, one can obtain a propagation equation that is of the form of Shvyd'ko's time and space wave equation [66], with the wave equation kernel expressed explicitly for the elastic scattering case considered here.
When one is interested in steady-state scattering properties, the above equation can be solved directly in frequency space using transfer matrices or Parratt's method [65], which has been generalized to the layer formalism in the nuclear resonance scattering community [11,64]. The above derivation thus unveils the connection to the semi-classical theories used in nuclear resonant scattering, clarifying its relation to the full quantum theory of the absorbing dielectric environment interacting with the nuclear transitions. The main insight is that the central approximation in these approaches is the low excitation approximation. In the linear excitation regime, where this approximation applies, the approaches are then analogous. We note that in practice, slight differences in the formulas can be found, since additional convenient approximations are made in the x-ray case, which are well applicable, but unnecessary from a formal perspective. Examples include the rotating wave approximation. We also refer to Appendix A 6 b for further context.

C. Nuclear Master equation
Starting from the macroscopic QED Hamiltonian in the rotating-wave approximation, the Born-Markov approximation can be used to derive an effective Liouvillian for the nuclei interacting with each other via the electromagnetic field [71,73]. Since nuclei and x-rays usually feature very weak coupling and the cavities in use are highly leaky, the Born-Markov approximation is applicable in most cases. Recently, systems featuring collective strong coupling have been reported [25], where the Born-Markov approximation may break down and the ab initio few-mode theory presented in Sec. II B may be advantageous.
Following the approach in [73,105] and again excluding polarization effects for simplicity, we can write the effective Hamiltonian in our ensemble notation aŝ and the Lindblad term as where J ll nn [Γ ll nn ] is the nucleus-nucleus coupling [decay] constant and E in (r ln ) = E in (r ln ) is the driving field [73] providing a driving strengthΩ ln = d * · E in (r ln ) for the nucleus n in ensemble l. Within the Born-Markov approximation, the nuclei are thus described as driven by the incident cavity field without resonant modification. The couplings and decay constants can be obtained from the Green's function of the system [73] by J ll nn = µ 0 ω 2 nuc,l d * l · Re[G(r ln , r l n , ω nuc,l )] · d l , (64) Γ ll nn = 2 µ 0 ω 2 nuc,l d * l · Im[G(r ln , r l n , ω nuc,l )] · d l .
In the above form, the effective nuclear Hamiltonian includes intra-layer couplings between individual nuclei, constituting an effective level scheme beyond the single parallel wave vector approximation employed in Sec. II F. The resulting many-body coupling scheme is illustrated in Fig. 11. In order to obtain an ensemble picture as in Fig. 5, one can introduce spin-wave operators [105,106,120] σ ± l (k ) = nσ ± ln e ±ik ·r ,n , which under the translational invariance assumption diagonalize the effective Hamiltonian [105]. Here, since we are mainly interested in the linear sector, where the parallel wave vector k of a given excitation is conserved, we can then again derive an effective Hamiltonian for the subspace at a single parallel wave vector (see also Appendix A 4), as we show in the following. Such an effective level scheme is of interest, since it corresponds directly to the effective level schemes derived in Sec. II F from the ab initio few-mode approach and to the corresponding interpretation of recent experiments [20-24, 26, 45, 47]. The Green's function approach provides an alternative ab initio method to calculate these level schemes, with different approximations involved bearing the potential for different generalizations. The main advantage of the Green's function approach over the ab initio few-mode theory is its numerically efficiency for calculating effective level schemes in the case of the layered cavity geometry, as we demonstrate in the following sections.

D. Effective nuclear level scheme in the low excitation subspace
In the low excitation regime, where σ z ln ≈ −1, the equation of motion for the lowering operator resulting from the Born-Markov Master equation reads where G(r ln , r l n ) =J ll nn + i Γ ll nn 2 and we have dropped the expectation value brackets for brevity. Using the approximate translational invariance of the x-ray cavity, we can write the Green's function as [121] G(r ln , r l n , ω) = d 2 k (2π) 2 G(z l , z l , k , ω)e ik ·(r ,n −r ,n ) , where G(z l , z l , k , ω) is the Fourier transformed Green's function and lattice offsets between the ensembles of less than a lattice constant are neglected. Similarly, we define G(z l , z l , k ) = µ0ω 2 nuc,l d * l · G(z l , z l , k , ω nuc,l ) · d l . Transferring to the spin-wave basis [106], the equations of motion simplify to where we have used n e i(k −k )·r ,n = (2π and used the approximation n E in (r ln )e −ik ·r ,ln ≈ N A d 2 r ,n E in (r ln )e −ik ·r ,n , which is valid for grazing incidence illumination, where the phase variation of the illuminating beam is small over the length scale of the lattice parameter, and which effectively neglects Bragg scattering. Importantly, the driving field above is also approximated as independent of n within one ensemble l. If all nuclei are located at the same z l , this is an exact representation of the ensemble in the considered geometry. Practically, however, one often wishes to interpret a resonant layer of finite thickness as a single ensemble [20,21,45]. The driving field is then approximated as constant over the layer thickness and taken equal to its central value at z l , which we denote as the thin-layer approximation. We refer to Sec. III H for a practical discussion of this approximation and its consequences, as well as to the previous discussion of thick layers in the context of the ab initio few-mode theory in Sec. II H.
We see that Eq. (69) provides a closed set of operator equations at a given parallel wave vector. The effective subspace Liouvillian corresponding to this linear equation of motion is given bŷ and the Lindblad term as where We see that this effective level scheme corresponds in close analogy to the one derived from the few-mode approach in Sec. II F as depicted in Fig. 5, with the driving field being expressed as a channel mode expansion in the In summary, in the linear regime, the nuclear dynamics for an incident field of a defined parallel wave vector are given by the effective level scheme Eqs. (71), (72). For an incident field containing multiple parallel wave vectors, the superposition principle applies in the linear regime.
We note that as also observed in the ab initio fewmode theory (see Appendix A 6), the effective linear level scheme parameters depend on the dipole moment scaled by N/A , such that only the nuclear number density is relevant for linear observables.

E. Linear solution in frequency space
Eq. (69) can be solved in frequency space to give [69,73] where F. Reconstructing spectral observables The Born-Markov Master equation or the effective nuclear level scheme given above define the dynamics of nuclear observables for a given driving field, which can be solved to obtain, for example, the linear regime solution Eq. (75). Scattering observables can be reconstructed using a generalized input-output equation [73], which is also valid beyond the linear sector. At a given parallel wave vector and in frequency space, it readŝ Substituting the linear frequency space solution for the lowering operator yieldŝ This formula describes the field at position z for a given parallel wave vector k and frequency ω, including the nuclear response. It is thus to be interpreted as a steady state solution for a driving field at frequency ω. E in (z, k , ω) is the corresponding field in the absence of the nuclei. It therefore is proportional to the mode profile and a driving amplitude factor α in where our notation is chosen close to [121]. The function E 0(n) q (z, k , ω) is the mode profile that is normalized to the surface of the cavity from where the radiation is incident [121], with 0 [n] indicating top [bottom] illumination. The polarization index q accounts for sand p-polarization. For details on the mode functions used to expand the Green's function [121] we refer to Appendix C.  12. (Color online) x-ray cavities with nuclei in EIT (a) and non-EIT (b) configuration, which were investigated experimentally in [21] and modeled theoretically in [45]. (a) and (b) show the cavity structure (for materials, see legend) and off-resonant field distribution for the EIT (a, cavity 1) and non-EIT (b, cavity 2) case, respectively. We note that over the width of the blue resonant layers, the field distribution shows visible variations. The white dashed lines indicate the boundaries of sub-ensembles (see main text) over which the field distribution does not vary significantly.
In grazing incidence experiments at synchrotron facilities [11], a common observable is the reflection spectrum as investigated in Sec. II G. In the Green's function approach, such spectra can be obtained by evaluating the total field at the surface of the cavity Similarly, the transmission can be obtained by wherez 0 [z n ] is the position of the surface boundary of the uppermost [substrate] layer (see Fig. 20). For the substrate having absorptive character via a non-zero imaginary part of the refractive index, the wave amplitude decays upon propagation through the medium, such that the transmission coefficient characterizes the amplitude ratio at the last layer surface [121].

G. Numerical efficiency for the layer geometry
Conveniently, the layered x-ray cavity geometry is one of the few cases [73,103] where the Green's function for the cavity is known analytically [121]. In particular, G(z, z , k , ω) can be expressed algebraically via an analytic recursion formula [121] similar to Parratt's formalism [65]. This feature makes the above approach  Fig. 12. The spectra are calculated using the Green's function approach in the third cavity mode, and the well-known layer formalism (blue) is compared to the Green's function approach presented in this paper. Without dividing the two nuclear ensembles into subensembles (red), the spectra fit well and reproduce the main spectral features, but fail to reproduce more subtle additional features (see insets). These are captured if each resonant layer is divided into three sub-ensembles (yellow), and therefore can be attributed to the layer thickness, resulting in a field gradient across the layer (see also Fig. 14).
highly numerically efficient. Compared to the ab initio few-mode theory approach to calculating effective nuclear level schemes presented in Sec. III D, this feature poses a major advantage when one is interested in the calculation of effective quantum optical parameters as depicted in Fig. 5, and, for example, opens optimization opportunities. The formula for the Green's function and its evaluation, as presented in [121], are summarized in Appendix C. In the following, we employ a numerical implementation thereof in order to benchmark the approach, and to demonstrate its usefulness. As a main result, we provide an ab initio nuclear level scheme for the EIT cavity configuration investigated experimentally in [21], resolving previous discrepancies in the quantum optical description of the system [45].

H. Application: Electromagnetically induced transparency cavities
The interesting cavity geometry mentioned above that has notably been studied experimentally in [21] features a sharp spectral dip in the spectrum, which has been shown to correspond to an EIT phenomenon, and a corresponding effective level scheme has been derived [21,45]. The system can be described by the phenomenological fewmode theory and the comparison to spectral observables from semi-classical calculations showed that the main EIT feature can be reproduced. However, there are unexplained quantitative and qualitative disagreements [45]. In addition, it is unclear to what extent the heuristic extensions of the model (see Sec. II A) that were found to be necessary [45] influence the interpretation in terms of an effective level scheme.
In this section, we apply the Green's function approach developed above to the two cavities studied in [11,45], resolving these open questions. In particular, we show that the qualitative disagreements in the spectra arise due the relatively thick resonant layers that were used in these cavities, causing the formation of cavity mode field gradients across the layers. The new approach can incorporate such gradients, and thus provides excellent quantitative agreement, improving the previous phenomenological description of the system significantly. In addition, we unambiguously calculate the effective nuclear level schemes using the ab initio method and investigate its quantum optical parameter trends.

Double resonant layer cavities and spectral benchmarks
The two cavity layer stacks under consideration are depicted in Fig. 12 and are identical to the geometries investigated in [45] to understand the EIT effect observed in [21]. They consist of platinum cladding layers enclosing a carbon guiding layer that is doped with 3 nm thick resonant 57 Fe layers at certain positions. In both cases (cavity 1 and cavity 2), one of the resonant layers is located at the cavity center, where the field distribution in the third cavity mode at incidence angle θ = θ 3 features an anti-node. As shown in Fig. 12, in the EIT cavity (cavity 1, panel a) the second resonant layer is placed in the adjacent field distribution node closer to the cavity surface, while in the non-EIT cavity (cavity 2, panel b), it is placed in the node closer to the bottom of the cavity.
We refer to the two structures as "EIT" ("non-EIT") cavities, because they feature (do not feature) a pronounced spectral dip in the respective nuclear spectra shown in Fig. 13 [21]. The two panels compare the spectra obtained using the Green's function approach to reference spectra calculated using the semi-classical layer formalism [11,64] implemented in the software package pynuss [96]. For the Green's function method, two different curves are shown in each panel. The red line corresponds to a model in which each of the two layers is treated as a single nuclear ensemble, neglecting possible field gradients across the layer. While the qualitative agreement is good and the major spectral features are reproduced, there are quantitative differences, most visible in additional small spectral features shown in the insets. To capture these features, we divide each layer into three sub-layers in our theoretical model (see white dashed lines in Fig. 12), to better account for the variation of the field intensity across the layers. The result is shown as the solid yellow line, which essentially agrees perfectly with the semi-classical calculation. The improvement can be understood since in Fig. 12 it is clear that the field distribution varies visibly over the thickness of each resonant layer, while it is constant to a good approximation over the thickness of each of the sub-ensembles.
We further note that the agreement between the spec-tra is much better than in the case of the phenomenological few-mode fits shown in [45], even when the subensemble partition is not considered (red line in Fig. 13). This general quantitative improvement is due to the absence of heuristic extensions and problems related to the fitting procedure [45], which are not needed in the ab initio theories reported here.

The EIT effect and thick layer sub-ensembles
Next, we extend the discussion to the two-dimensional spectra as a function of detuning and incidence angle, which were also investigated for the two cavities in [45] in the context of the phenomenological few-mode model. In this reference, it was found that in addition to the quantitative differences in the one-dimensional spectra (see Fig. 13 and the discussion above), there are qualitative features that are not captured by the phenomenologically fitted model even with the heuristic extensions included.
Results are shown in Figs. 14 and 15 for the EIT-and the non-EIT cavity, respectively. In each figure, panel (a) corresponds to the layer formalism calculation which again serves as a benchmark. Panels (b) and (c) show the Green's function results with and without sub-ensemble partition of the resonant layers, respectively. Finally, panel (d) shows the absolute value of the difference between the results in panels (a) and (b).
We see that while (d) demonstrates excellent agreement of the Green's function description with subensembles to the layer formalism, the corresponding results without sub-ensembles (panel c) are missing a spectral feature as indicated by the red dashed ellipse. Its absence shows that approximating thicker resonant layers as a single thin layer can lead to qualitative differences in the theoretical description.
For the EIT case in Fig. 14, we further find that the agreement of the model without sub-ensembles is still rather good close to the third mode minimum θ = θ 3 , justifying the thin-layer approximation at this incidence angle. However, the differences become sizable in the region of the red dashed ellipse, which can be understood by noting that the field distributions vary more rapidly as function of position in the cavity at higher incidence angles. Fig. 15 shows analogous results for the non-EIT cavity (cavity 2). In this case, the absence of the spectral feature in the two-ensemble model also causes the surrounding spectral structure to change, with an avoided mode crossing turning into a merging point (see region inside the red dashed ellipse).

Effective nuclear level schemes
Finally, we discuss the effective level scheme corresponding to the Green's function approach to x-ray cavities with nuclei. Fig. 16 shows the effective nuclear level scheme and the quantum optical coupling constants for   Fig. 16 as an example. The coupling constants can be calculated directly using the Green's function approach and are shown in the table in matrix form for the EIT (cavity 1) and non-EIT (cavity 2) case. As before, we consider excitation in the third cavity mode (θ = θ3). The drive vector has been normalized, since it is proportional to the applied x-ray field amplitude. The corresponding quantities as a function of incidence angle are shown in Fig. 17. We note that the layer ensembles are labeled by their position in the cavity from top to bottom, such that the central layer in the anti-node corresponds to index 1 in cavity 1 and to index 2 in cavity 2.
the EIT-and non-EIT cavities in Fig. 12. For simplicity, the case without sub-ensemble partitioning is shown, with the notation adopted from [45]. Within the Green's function approach, the couplings and other quantum optical parameters can be calculated directly and unambiguously from the cavity geometry. The resulting parameters at the incidence angle θ = θ 3 are tabulated in Fig. 16. In Fig. 17, we show the quantum optical parameters as a function of incidence angle, revealing the changes of the three level system over the different modes. In particular, the collective Lamb shifts, the superradiant decay rate enhancements and the level couplings between the two ensembles in the EIT cavity (cavity 1) are depicted. Around each mode minimum (θ 1 , θ 2 , θ 3 ), the collective Lamb shift and superradiance roughly show the typical behavior of real and imaginary parts of a complex Lorentzian [23,46]. For the coupling parameters δ 12 and γ 12 , similar structures are observed. Only in the third mode, where the EIT phenomenon is observed, a different functional dependence and a negative γ 12 is found.
These results show how complex cavity structures can be unambiguously interpreted in terms of quantum optical models, which paves the way for designing effective nuclear level schemes via tailored mode environments.
As a last consistency check, we show that the collective Lamb shift and superradiance calculated here indeed correspond to what these quantities have been associated with in the nuclear cavity QED literature so far [20,23]. To this end, we replace the second resonant 57 Fe-layer by its off-resonant 56 Fe counterpart and fit a generic Fano model [23,89,122,123] to the resulting line shapes at each incidence angle. The Fano profile fit function for the reflectance spectrum is given by [23] |r(∆) where (∆) = (∆ − δ 1 )/(γ 1 /2). The scale factor σ 0 , the complex Fano parameter q and δ 1 , γ 1 are the fit parameters. This fit provides an alternative way to obtain the collective Lamb shift δ 1 and the superradiance γ 1 from the spectra, which can be compared to the ab initio predictions. Indeed, the fit results (dashed lines in Fig. 17b) show excellent agreement with the Green's function calculation, confirming the interpretation of the effective nuclear level scheme derived from the Born-Markov Master equation in Sec. III D.

IV. DISCUSSION AND SUMMARY
In summary, we have presented two ab initio approaches to describe thin-film x-ray cavities doped with narrow resonances such as those provided by Mössbauer nuclei on a quantum mechanical level. Our results improve the previously introduced phenomenological inputoutput model [44,45] for the system and resolve multiple open questions in the theory, which up to now have hindered extensions towards new parameter regimes and cast doubts on the models' predictive capabilities at higher intensities. This progress to an ab initio and, in the linear regime, essentially exact theory provides qualitatively new value to quantum optical interpretations used in the field, in the same spirit as recent developments connecting other sectors of quantum optics and ab initio theory [111,115,119,124].
The two ab initio methods presented here are directly applicable to model linear scattering experiments in grazing incidence, which have been performed extensively in the platform of hard x-ray cavity QED with nuclei [20][21][22][23][24][25][26] and electronic resonances [47], and to interpret them quantum optically in terms of an ab initio effective nuclear level scheme. The ab initio perspective and the exact treatment of the linear sector put the quantum optical interpretation on firm grounds, establishing it as a tool in the quantum theory of nuclear resonance scattering. Beyond that, we provide clear theoretical connections between the existing approaches and outline the approximations involved in each case, paving the way for generalizations. As a consequence, one can now seamlessly switch between the different theories, depending on which formulation is most favorable for the respective purpose or regime. In this context, we emphasize that a main advantage of our method over the existing theories is their applicability in the fully quantum mechanical sector. That is, the approaches bear the potential to describe non-linear and correlated quantum dynamical effects in hard x-ray cavity QED, going beyond linear, low-excitation, semi-classical or mean-field phenomena.
Our first approach is based on a few-mode description of the system, and promotes the established phenomenological models [44,45] to an ab initio theory without the need for heuristic extensions, thereby resolving previously not understood discrepancies in the modeling. The ab intio theory features modified angular and frequency dependencies of the quantum optical parameters, which can now be calculated directly from the cavity geometry. We further presented a closed form solution to the resulting equations of motion for general multi-mode multiensemble systems in the low excitation regime, including the derivation of scattering observables and the effective nuclear level scheme.
For an analytically solvable example cavity that features strongly overlapping modes, we demonstrate the advantages of our theory over the previous phenomenological models, and show that it provides quantitative agreement with semi-classical approaches. We investigate the functional dependencies of the quantum optical parameters, showing where the improvements provided by our theory originate from and finding non-trivial behavior as a function of incidence angle. As an application, we for the first time extend the quantum optical modeling to the realistic case of resonant layers with arbitrary thickness. We find that the resulting effects can straightforwardly be captured in the ab initio few-mode theory, a task which is difficult in phenomenological approaches due to the large set of fitting parameters required. Our results show that thick layers lead to coupled continuum ensemble descriptions, a feature which may lead to the x-ray implementation of qualitatively different quantum optical model systems than the few-level systems realized so far [20][21][22][23][24][25][26]47].
As our second approach, we develop an alternative method based on well-known Green's function techniques [70-73, 103, 125]. The main motivation for this second approach is its numerical efficiency for calculating effective nuclear level schemes. We again demonstrate the theory's connection to the existing semi-classical nuclear resonant scattering literature [11,63], deriving an effective wave equation including the nuclear resonance dynamics in the linear regime. We then show how a quantum optical description analogous to the previously investigated effective nuclear level schemes can be obtained and demonstrate its numerical efficiency, which is ensured by an analytic solution for the Green's function of the relevant geometry [121].
To showcase the power of this approach, we apply the Green's function method to the case of the EIT cavities investigated in [21,45], where qualitative and quantitative discrepancies of linear spectra between the phenomenological few-mode model and semi-classical theory were found [45]. Our approach is able to calculate the quantum optical description of the system in terms of an effective nuclear level scheme without the need for heuristic model extensions or a fitting procedure which were necessary before [45], and further reveals the importance of sub-ensembles of the relatively thick resonant layers in the system which are responsible for the aforementioned qualitative deviations. We show that quantitatively, the approach yields essentially perfect agreement if the resonant layers are divided into a sufficient number of sub-ensembles. We further present the first ab initio calculation of the resulting effective nuclear level schemes and investigate trends of the quantum optical parameters. These results pave the way for tailoring mode environments in x-ray cavities to design nuclear quantum systems.
Beyond the layer geometry, the Green's function formalism is easily adaptable to alternative photonic environments for the nuclei. Such structures may become accessible with improving fabrication techniques and already reported examples include other dimensional waveguides [126], curved channels [127], nanowires [128] and periodically structured surfaces such as nanodots or nanodiscs [129]. While unlike in the layer geometry and some others [72,73,102], the Green's function is usually not known analytically in these cases, various numerical and modeling approaches exists to tackle such geometries [73,103,130].
In this context, we also note that in comparison to the ab initio few-mode approach, in the presented Green's function formulation the mode structure of the cavity is hidden and cannot be included in the quantum optical modeling. It remains to be seen if the Green's function can also be modeled in terms of the modes or resonances [73] to address the relevant case of overlapping modes characteristic of x-ray cavities, e.g., using momentum space representations of the Green's function.
As a whole, our results provide a comprehensive solution for the linear scattering regime of thin-film x-ray cavity QED with Mössbauer nuclei or ultra-narrow resonances, resolving previous discrepancies and revealing the connection between different existing theoretical approaches. As an outlook, we expect this progress in the understanding of the low-excitation sector to provide a solid theoretical foundation for describing phenomena in the non-linear and correlated quantum dynamics regime of this platform, which may be accessible at current and upcoming x-ray facilities [6,12] such as x-ray free electron lasers, where the first experiment on nuclear quantum optics has recently been performed [38].
where r is a displacement vector in the layer plane, k is the corresponding parallel wave vector, C is a normalization factor for the parallel mode component, and f m (z, k ⊥ , θ) is an effective one-dimensional mode function at incidence angle θ. The latter fulfills the effective one-dimensional mode equation The perpendicular wave component can be obtained from the mode energy via k ⊥ = ω sin(θ) and is to be understood as a scattering boundary condition, that is the incoming perpendicular wave vector in the far field (see Fig. 1 for an illustration).
In summary, the classical three dimensional problem for the layer geometry can be reduced to an effective onedimensional problem with an effective dielectric constant and wave energy These substitutions are employed in the analytical calculation in Sec. II G.

Cavity quantization and few-mode theory
The above wave equation can be quantized canonically [131], resulting in the Hamiltonian whereĉ m (k , k ⊥ ) are bosonic operators at a given parallel and perpendicular wave vector, with ω(k , k ⊥ ) their corresponding frequency. The operator equations of motion for this Hamiltonian are equivalent to the Maxwell wave equation [131]. The advantage of the quantized formalism is that we can treat interactions with resonant nuclei (see Sec. A 3) beyond mean-field Maxwell-Bloch [28,66] or semi-classical scattering treatments [11,63]. We note that the canonical normal modes quantization is valid for real refractive indices. To account for complex refractive indices, which are relevant in the hard x-ray regime, we employ a non-hermitian Hamiltonian prescription, which allows the above Hamiltonian to be used directly (see Appendix A 5 for details). For a complete treatment of absorptive processes within the framework of macroscopic QED see Sec. III, where an alternative approach is developed. The Hamiltonian Eq. (A7) features the typical normal mode continuum of open scattering systems [69,131]. In order to connect to the successful phenomenological model [44,45] for thin-film x-ray cavities with Mössbauer nuclei, the notion of a few resonant modes has to be introduced. To perform this step systematically without a model or fitting prescription, we employ the ab initio few-mode theory [69], a projection method that yields a few-mode basis transformation. As shown in [69], treating the problem in such a few-mode basis can provide a powerful approach to the light-matter quantum dynamics by extracting the relevant degrees of freedom of the light field.
Following this approach, we can partition the cavity Hamiltonian given above into a few-mode part and an external bath. If we choose the few mode basis to respect the translation symmetry, the full Hamiltonian can be written as where where we have employed the relabeling from Appendix A 1 andâ λ (k ) are bosonic few-mode operators at each parallel wave vector, which have frequency ω λ (k ). Similarly,b m (ω, k ) are the external bath operators, which couple to the few-mode operators with coupling strength W λm (ω, k ). We see that the Hamiltonian is a linear combination of one-dimensional few-mode terms at each parallel wave vector k . The parallel direction thus still features a continuum. Since there is no confinement or resonance structure in the parallel direction, this continuum can not easily be removed by another few-mode projection.

Nuclear resonant interaction
The interaction with nuclear transitions can be described within the dipole and rotating wave approximation by the Hamiltonian andσ +,−,z are the Pauli operators for the nuclear transitions. We have included the effect of multiple ensembles (indexed by l), each of which contains N l individual nuclei (indexed by n) and couples to the system modes (indexed by λ) with coupling constant g λln (k ) as well as to the bath modes with coupling constantg mln (ω, k ). These coupling constants can be calculated from nuclear transition and material parameters (see Appendix A 6 for details), which are also used in the layer formalism [11]. We note that if a sufficient number of system modes is chosen, the few-mode approximation can be performed [69], that is the direct interaction of the nuclei with the bath modes can be neglected.
We further note that in addition to these Hamiltonian terms which arise from the nucleus-light coupling, nuclear transitions can feature additional incoherent decay channels, such as internal conversion [11,60,132]. For example for 57 Fe where the α-factor is 8.56 [11], internal conversion makes up about 1 − (2α + 2) −1 ≈ 95% of the incoherent decay rate (see Appendix A 6 b for details). In our prescription, this effect can be accounted for by a Lindblad term [22] where γ IC is the internal conversion decay constant. We note that since the field continuum is still present in the theory, radiative incoherent losses are already included in the light-nucleus interaction Hamiltonian and should not be added as an additional incoherent Lindblad term. If the few-mode approximation is performed, a small residual decay contribution is added to the incoherent decay rate to account for the weak interaction with the removed bath modes, which tends to zero at large system mode numbers.
Within the dipole, rotating-wave and non-hermitian field Hamiltonian approximations (see also Appendices A 5, A 6 b), the Hamiltonian derived above is a fully general description in few-mode form, with the translational invariance of the system explicitly implemented.

Effective one-dimensional problem
We see that the derived few-mode Hamiltonian is already very similar to the phenomenological model [44,45], with the main difference being the continuum of parallel wave vectors. In order to complete the connection, we derive an effective one-dimensional description in this section, which is well applicable in the linear excitation regime.
We first note that in practice, the nuclear cavity QED system is often studied spectroscopically. That is the system is excited by a collimated and highly monochromatic * beam from a modern low-emittance x-ray facility such as a synchrotron or an x-ray free electron laser in conjunction with a high-resolution monochromator. At anticipated light sources such as x-ray free electron laser oscillators [6], similar setups are to be expected.
In such a setup, the exciting light defines a narrow range of incidence angles and consequently a narrow range of parallel wave vectors. Due to the assumed translational invariance of the system in the layer plane, the parallel wave vector is a conserved quantity in the lowexcitation regime (see Sec. III D for details, and [133] for a related discussion of higher excitations), such that each parallel wave vector forms an isolated subspace of the dynamics. Therein, we can then obtain the effective one-dimensional Hamiltonian with and H 1D int (k (in) ) = λ∈modes l∈ensembles n∈1,2...N l g λln (k (in) )â λ (k (in) )σ + ln + h.c. .
(A19) * On the spectral scale of the cavity, the beams are monochromatic. On the scale of the nuclei on the other hand, the exciting radiation has a broad spectrum.
In addition to the internal conversion Lindblad term, we obtain a radiative contribution to the incoherent decay (A20) where γ rad is the spontaneous emission rate of the nuclear transition in the cavity environment. The two incoherent terms can be combined to (A21) which is the single-transition version of the term used in [44]. In the case of 57 Fe, the full spontaneous emission rate γ SE is dominated by internal conversion, such that we can approximate the latter by the natural linewidth of the transition γ SE = γ IC + γ rad ≈ γ.
Noting that k (in) appears only as a parametric dependence now, we see that the above description provides the improved input-output model summarized in Sec. II B, where the parallel wave vector dependence is rewritten in terms of the incidence angle and the frequency of the transition energy.
We have thus provided an ab initio generalization of the successful phenomenological model [44,45] for thinfilm x-ray cavities with Mössbauer nuclei and clarified its origin as well as the involved approximations. The basic structure and its relation to the original phenomenological version [44,45] are summarized in Sec. II B.
We note that, as we show in the following sections, the improved model provides a number of qualitative advantages, allowing the quantum optical description to be applied to new systems and featuring essentially exact predictions in the linear regime. The ab initio character of the theory further provides a solid foundation, which, as a main motivation beyond the quantum interpretation of linear scattering experiments, will allow the method to be applied as a predictive tool in the non-linear and correlated quantum dynamics regime.

Complex refractive index in the ab initio few-mode theory
In the treatment presented in Appendix A 2 and also in the original development of the ab initio few-mode theory [69], a real refractive index is considered. x-ray cavities, however, feature significant material absorption and as a result a complex refractive index should be accounted for.
A rigorous treatment of material absorption and the resulting effective quantum theory can be obtained in various ways (see [72,103] for a review as well as [134] for a recent advance). In general, the resulting light-matter Hamiltonian is highly complex even if the resonant quantum interaction is not included.
Here, we resort to a simple approach, which is also employed in the standard nuclear resonant scattering literature including the semi-classical scattering theory [11,63]. In the approach, the Maxwell wave equation is directly coupled to the resonant nuclei or atoms, while neglecting quantization effects of the light field. For a real refractive index, the generalization to the quantum level is given by the canonical quantization scheme (see Appendix A 2). The absorptive character of the material is then included by using the complex refractive index to obtain a non-Hermitian Hamiltonian. Since the approximation is also included in the standard semi-classical xray scattering theory [11,63], previous experiments substantiate its validity at least for intensity observables at low driving fields. The approach can be complemented by the quantum jump formalism [135][136][137], similarly to recent work using Green's function techniques [106], where the absorptive bath is treated rigorously from the outset (see Sec. III for a detailed comparison of the approaches).
Practically, the non-Hermitian few-mode Hamiltonian approach is implemented by performing calculations as for a real refractive index, and then transferring to the non-Hermitian theory by substitution of the complex refractive index. We note that for numerical implementations, one has to ensure that complex conjugation operations, as they are used for example in the quantum scattering theory [86] underlying the projection scheme, are not applied to the refractive index.

Coupling constant in the ab initio few-mode theory
In this appendix, we derive the mode-nucleus coupling in terms of known nuclear resonance and material parameters.

a. Mode coupling and quantization area
The few-mode coupling in our approach can be written as [69] g λln (k ) = − idω nuc,l 2ω λ (k ) where χ λ (z, r , k ) is the three-dimensional system mode from Eq. A13 and z l the vertical position of the nuclei, which is by construction independent of n within the thin-layer approximation. If z l varies within one ensemble, the latter should be divided into multiple subensembles (see also the thick layer treatment in Sec. II H).
Since the system modes themselves also fulfill the translational invariance, we have whereχ λ (z l ) is the effective one-dimensional mode function evaluated at the nuclear layer position z l and C is a normalization factor. We therefore see that the mode normalization constant is important to determine the value of the coupling. In order to obtain the constant, we resort to the usual box quantization procedure [138] for the parallel direction. We require that where A is a parallel quantization area, such that we obtain Consequently, the coupling constant is given by The phase factor e ik ·r ,ln can be absorbed into an algebra preserving redefinition of theσ-operators in the effective one-dimensional Hamiltonian [87], amounting to a basis change that effects the nuclear dipole moment by a phase. The coupling is thus independent of the nuclear index n within one ensemble, such that we have We see that the coupling remains dependent on the quantization area A . While this feature may initially seem unphysical, it is to be expected, since the number of nuclei N l participating in the dynamics also depend on this area. As is also noted in the phenomenological model [44], in the linear regime the collective coupling constant is the relevant quantity, where ρ N is the number density of the resonant nuclei and f LM is the Lamb-Mössbauer factor encoding the fraction of nuclei which participate in the recoil free scattering or decay process [11,60]. t l is the thickness of the ensemble layer l and we assume the thin layer limit for a single ensemble by taking z l to be independent of the nuclear index n. We further drop the parametric parallel wave vector dependence for brevity. In the linear limit, the properties of the system thus only depend on the number density and not on how many nuclei participate in the dynamics, which is consistent with the semi-classical layer formalism [11,60,87].
Beyond the linear limit, however, the absolute number of nuclei and the quantization area become important and do not reduce to the number density in the equations. This feature can be understood by recognizing that in a fully translation invariant system, the number of nuclei is necessarily infinite. Physically, however, only a finite number participates in the dynamics due to a limited coherence volume and a limited size of the excitation beam. The quantization area A should therefore be chosen to capture these physical features. We further note that the translationally invariant description essentially neglects finite excitation size effects and excitation spreading at the edge of the excitation region. A phenomenological argument to include such contributions is given in [12].

b. Effective nuclear dipole moment
A remaining question is how to compute the effective single nucleus transition matrix element d from tabulated resonance parameters. The connection to quantities that are conventionally used for example in the layer formalism [11] can be obtained by a simple comparison of physical scattering observables.
The formula for the resonant contribution to the refractive index of an resonant nuclear medium is [87] n57 where we have taken 57 Fe as an example, and "magnetic splitting and polarization dependence" [87] have been neglected [87]. The quantities here are defined as in [87], including the Lamb-Mössbauer factor f LM , the resonant wave number k 0 , the spins for the ground (excited) state I g (I e ) and the internal conversion factor α.
A formula for the resonant contribution to the refractive index in the linear regime can alternatively be obtained from our effective transition theory, for example via a linear dispersion theory [69] calculation, which gives ε res (r) = − |d| 2 ω − ω nuc + i γ 2 f LM ρ N . (A31) where we have neglected the A 2 -term contribution, assumed a dense lattice, which for example results in Bragg scattering being neglected, and performed the rotating wave approximation by setting ω 2 − ω 2 nuc ≈ 2ω nuc (ω − ω nuc ) as well as ω 2 nuc ω 2 ≈ 1 [69]. We note that all of these approximations are already implicit in the nuclear refractive index formula Eq. (A30), with more general treatments including lattice effects and other scattering processes being available [60]. The Lamb-Mössbauer factor f LM encodes the fraction of nuclei that participate in the recoil free scattering process [11,60] and therefore modifies the nuclei number density.
A further approximation that is commonly performed in nuclear resonant scattering and is implicit in Eq. (A30) is the small response approximation n = √ ε = √ 1 + χ ≈ 1 + χ/2. Performing this approximation also in the linear dispersion theory case we obtain n − n electronic ≈ − |d| 2 ω − ω nuc + i γ 2 f LM ρ N 2 . (A32) In order to obtain the effective nuclear transition matrix element d, we can therefore straightforwardly compare Eq. (A32) to Eq. (A30), which gives |d57 Fe | 2 ≈ 2πγ k 3 0 1 2(1 + α) 2I e + 1 2I g + 1 . (A33) We recall that as our previous calculations this formula is given in natural units with = c = ε 0 = 1. We can confirm and understand the result by a comparison to the radiative spontaneous emission formula for a single dipole transition in free space [106], γ rad = ω 3 nuc |d57 Fe | 2 /3πε 0 c 3 , giving Noting that the radiative decay rate γ rad is given by the full decay rate γ and the internal conversion fraction, γ rad = γ/2(1 + α) ≈ 0.05γ, that the spin factor reduces to (2I e +1)/(2I g +1) = 3/2, and using consistent units, we see that Eq. (A34) essentially corresponds to a rearrangement of Eq. (A33). This second treatment also generalizes to more general transitions beyond dipoles, such that the transition moment vectors necessary for the Green's function approach in Sec. III can be extracted from the well-known properties of the nuclei [60]. This calculation also clarifies the connection to the semi-classical layer formalism [11], which is a standard tool for describing resonant x-ray scattering experiments [11,64] and can be used to include various experimental imperfections in the description [11,64]. In this context, we note that our quantum optical approach is based on nuclear transitions, whose properties are assumed to be known, as the starting point. In the established perturbative scattering theory [60], the transition structure of the nuclei is investigated in detail and effects such as interaction with the lattice are included. Related approaches such as Shvyd'ko's time and space picture [66] allow for the inclusion of additional effects such as inelastic scattering [67]. For further details on the relation between the formalisms refer to Sec. III B and Fig. 2. Appendix B: Detailed comparison of the ab initio and phenomenological few-mode approaches In Sec. II G 4, we presented a comparison of the phenomenological and ab initio few-mode approaches with regards to their capability to model nuclear spectra of the example cavity. However, this comparison is not unique, because the phenomenological models are based on fits of their parameters to predictions from semi-classical theories, and the fits can be obtained using different fit objectives, i.e., different approaches to quantify the differences between the model and the references. Note that these differences often are not of practical relevance close to the resonance at which the parameters are fitted, but deviations between different parameter sets are expected further away from the resonance. In contrast, the parameters of the ab initio model are unique, as they are calculated from the cavity structure without a fitting procedure.
Because of the ambiguity in the phenomenological parameters, in this appendix, we provide additional detail on this comparison to support our conclusions. In particular, we consider multiple ways of fitting the phenomenological model parameters. Fig. 18(a) shows the layer formalism spectrum used as a reference for this analysis. In order to obtain such a spectrum from the phenomenological model, in all cases, the empty cavity parameters are first fitted to the rocking curve as described in Sec. II G 4, yielding the best fit shown in Fig. 18(b). Next, we employ three fitting procedures for the mode-ensemble interaction parameters.
a. Fit method 1. In this method, we only fit one global interaction parameter. This method was already employed successfully in [45] for the EIT-scenario and uses an extraction of the individual mode parameters. The latter is achieved by decomposing the coupling constants as [45] N l g λl =Ẽ λl · ( N l g l ) , whereẼ λl is taken as the mode amplitude at the location of ensemble l when illuminated at the rocking minimum corresponding to mode λ [45]. For our single ensemble case, √ N l g l is then a single global scale that is fitted to the two-dimensional spectrum.
b. Fit method 2. In this approach we replace the mode parameter extraction prescription of fit method 1 and instead fit the √ N l g λl parameters for each mode. As a fit objective, we use the one-dimensional spectrum at resonance with the first cavity mode, that is at θ = θ 0 . c. Fit method 3. Here, we proceed as in fit method 2, only that agreement to the full two-dimensional spectrum (see Fig. 18) is used as a fit objective.
The results from the three methods are shown in Fig. 18. Panels (c,1)-(g,1) show the two-dimensional spectra calculated using the phenomenological model with the model parameters extracted from the respective fitting methods. The remaining panels show the corresponding residual deviation to the layer formalism result in panel (a), both on a linear scale in panels (c,2)-(g,2) and on a log scale in panels (c,3)-(g, 3).
We find that all models capture the behavior well, in particular at low incidence angles. At higher incidence angles, that is going towards the edge of the fitting range for the cavity parameters, the deviations become more significant. The agreement at higher incidence angles can be improved by including a larger set of modes, resulting in more fitting parameters. We further see that the various fitting procedures result in deviations due to Reference spectrum obtained using a semi-classical calculation. Panel (b) shows the best phenomenological model fit (five modes, yellow dashed) of the rocking curve (solid blue), which is used to constrain the cavity parameters (κ λ , κ R,λ , θ λ ) in the phenomenological models. Panels (c)-(g) show results for different few-mode approaches, with (c)-(e) comprising the phenomenological model with different parameter fitting procedures (fit method 1, 2 and 3, respectively, definition see text) and (f), (g) being the ab initio few-mode theory with 5 and 20 modes, respectively. Panels (c-g,1) show the model spectrum in each case, while panels (c-g,2) and (c-g, 3) show the residual deviation on a linear and logarithmic scale, respectively. the emphasis on different regions of the spectrum, but yield the same qualitative spectral features and the same overall quantitative level of agreement.
In comparison, Fig. 18(f) shows analogous results for the ab initio few-mode theory with 5 cavity modes. We see that the performance is similar to the phenomenological cases around the first resonance, but better at higher incidence angles. Panel (g) shows the corresponding 20 modes result, illustrating the systematic convergence when including more cavity modes. In the latter case, the residual deviations are mainly due to the thin layer approximation, with higher modes only having a tiny contribution.
In Fig. 19, we focus on the slice at the first mode reso- nance (θ = θ 0 ). To illustrate the excellent level of agreement of all approaches, the nuclear spectra are plotted on top of the semi-classical reference (panels a,b). The residual deviations are shown in panel (c). We find that the ab initio result at 20 modes is well converged. Indeed, the phenomenological result gives a comparable level of agreement with only 5 modes for the case of fit method 2, illustrating that the phenomenological model can yield very good agreement if one focuses on a particular spectral region. However, it is not known a priori which phenomenological fit works best for a given problem. The remaining fitting methods and the 5-mode ab initio theory show similar deviations, which are structurally comparable, with already a good quantitative agreement on the level of a few percent deviation. reflection coefficients, respectively, from layer i into layer j [121]. They can be calculated via a recursion formula [121] analogous to Parratt's formalism [65]. Specifically, the i/k coefficients can be expressed in terms of i/j, j/k coefficients for any in-between j by [121] r q i/k = r q i/j/k = r q i/j + (t q i/j t q j/i − r q i/j r q j/k )r q j/k e 2iβj dj D (ik) qj , (C3a) If the coefficient i/k is required, the conversion is conveniently started by choosing j = i + 1 for i < k or j = i − 1 for i > k, and terminated by the Fresnel coefficients for adjacent layers [121] r q ij = β i − γ q ij β j β j + γ q ij β j , (C4a) where the absence of a slash indicates that the layers are adjacent and γ p ij = ε i (ω)/ε j (ω) as well as γ s ij = 1. Together, the results from [121] summarized here provide a recursively analytic and hence numerically efficient way to calculate the Green's function at a given parallel wave vector for the layer cavities, which is the central quantity appearing in the effective nuclear level scheme described in the main text.