Coupling function from bath density of states

Modelling of an open quantum system requires knowledge of parameters that specify how it couples to its environment. However, beyond relaxation rates, realistic parameters for specific environments and materials are rarely known. Here we present a method of inferring the coupling between a generic system and its bosonic (e.g., phononic) environment from the experimentally measurable density of states (DOS). With it we confirm that the DOS of the well-known Debye model for three-dimensional solids is physically equivalent to choosing an Ohmic bath. We further match a real phonon DOS to a series of Lorentzian coupling functions, allowing us to determine coupling parameters for gold, yttrium iron garnet (YIG) and iron as examples. The results illustrate how to obtain material-specific dynamical properties, such as memory kernels. The proposed method opens the door to more accurate modelling of relaxation dynamics, for example for phonon-dominated spin damping in magnetic materials.

Introduction. -Quantum technologies face many challenges, often arising due to the unavoidable coupling of any system to its environment. The prediction of their dynamics requires open quantum system methods that include such coupling effects, for example the Caldeira-Leggett model [1] and the spin-boson model [2]. These methods are successfully employed in many physical contexts, e.g., quantum optics [3][4][5], condensed matter [6][7][8][9][10][11], quantum computation [12][13][14], nuclear physics [15] and quantum chemistry [16]. For instance, modelling circuit quantum electrodynamics with the spin-boson model shows that the heat transport of a superconducting qubit within a hybrid environment changes significantly, depending on the qubit-resonator and resonator-reservoir couplings [6].
In the mathematical treatment of an open quantum system, a coupling function C ω is typically introduced that describes how strongly the system interacts with bath degrees of freedom (DoF). Its functional form determines the temporal memory of the bath and whether the noise is coloured or not [1,2,17], critically affecting the system dynamics [8,18,19]. A large body of theoretical results exist for various toy models that make specific assumptions on the coupling function C ω [1,2,20]. However, a major drawback is a somewhat lacking connection to system-or material-specific characteristics to which these methods could be applied: for a given DoF, in a given material, which coupling function C ω should one choose to model its dynamics?
An alternative approach is taken in the condensed matter literature, where open quantum systems are usually characterized by the density of states (DOS) of their environ-ment [21]. Measurement of, for example, the phonon DOS is well-established using different inelastic scattering techniques [22,23]. Modes in the environment typically couple to the system with a wave vector-dependent strength g k [2,24,25], which in many cases can be captured by a frequency-dependent g ω .
In this paper, we present a useful relation that translates the coupling function C ω of an open quantum system into an experimentally measurable DOS D ω , and vice versa. While a similar relation has previously been reported for onedimensional quantum spin impurities [26,27], the relation obtained here is valid for a generic system coupled to a bosonic bath, capturing dimensionality and anisotropy. It paves the way to parametrizing realistic coupling functions for a range of applications, for example, for spins in a magnetic material that experience damping through the coupling to the crystal lattice [17,28] or for nitrogen vacancy centers, a solid-state analogue of trapped atoms, whose coherence lifetime in optical transitions is also limited by interaction with phonons [29,30]. The link is explicitly established for a generic quantum system that couples locally to a bosonic environment. Extensions to other environments, such as fermionic environments, will be possible using similar arguments.
The paper is organised as follows: we first introduce the two approaches involving D ω and C ω , respectively. Setting up the dynamics of the environment, we evaluate its memory kernel and establish the link between D ω and C ω , allowing for general g ω . In the second part of the paper, we choose a flat g for simplicity, and illustrate the application of the relation with a few examples. We demonstrate that the widely p-1 arXiv:2112.04001v2 [quant-ph] 15 Dec 2022 used Debye approximation is equivalent to the well-known Ohmic coupling function. While this approximation suffices at low frequencies, experimental DOS show peaks at higher frequencies, leading to non-trivial dissipation regimes. We parametrize two measured phonon DOS, those of gold and iron (see Supplementary Material (SM)), and one theoretically computed phonon DOS of yttrium iron garnet (YIG) and extract key parameters for the corresponding coupling functions C ω . These give direct insight into the impact of memory for any phonon-damped dynamics in these materials.
Two approaches. -The Hamiltonian of a quantum system in contact with a bath iŝ where the bath HamiltonianĤ B and the system Hamiltonian H S may contain the internal interactions among their own components. The system-bath interaction is assumed to be of product form,Ĥ whereŜ is a (Hermitian) system operator andB is a bath operator, each with d s components. The form of the bath HamiltonianĤ B and of the bath operatorB depends on the context. We consider here a bosonic bath, i.e. an infinite set of harmonic oscillators. In the literature, one can broadly distinguish two representations of the bath, working either in wave vector (WV) or frequency (F) space, as illustrated in Fig. 1.
The wave vector approach is common in condensed matter physics [2,21] where the bath Hamiltonian is expressed as a sum over all possible modes k Here ω = ω k gives the dispersion relation of a normal mode with wave vector k andb k (b † k ) are bosonic annihilation (creation) operators of a mode excitation with commutation relations [b k ,b † k ] = δ kk . Usually one considers a three-dimensional (3D) structure with wave vectors k = (k x , k y , k z ). For example, in a cubic 3D lattice with number of lattice sites N , lattice constant a and volume V = N a 3 , each component of k runs through the range − , . . . , 0, . . . , For large N and V , and for any function f (ω k ) that only depends on the frequency ω k , one can approximate sums over the wave vectors as This equation defines D ω as the DOS per unit volume of bath modes at frequency ω [21]. For bosonic baths, we choose the standard interaction [2] where the bath operatorB is linear in the bosonic mode operators (single phonon processes), where ξ k = k hg 2 k /(2ω k ) 1/2 with k a d s -dimensional unit polarisation vector [1] and g k the wave vector-dependent coupling, see Fig. 1. Eq. (2) may be generalized to the situation that several system componentsŜ m are located at different positions R m , and sum over interaction terms, i.e. H SB = − mŜ m ·B(R m ). The field operators would then be R-dependent, i.e.B W V (R) = 1 √ V k ξ kbk e ik·R +h.c.. For simplicity, we will concentrate in the following on just one system site and drop summation over m again.
Another approach to setting up the bath HamiltonianĤ B and the interactionĤ SB is based on a frequency expansion often employed in the open quantum systems literature [1,2]. In contrast to Eq. (3), hereĤ B is written directly as a sum or integral over frequencies, whereP ω andX ω are 3D [in general, d-dimensional (dD)] momentum and position operators, respectively, for the bath oscillator with frequency ω. Their components obey [X ω,j ,P ω ,l ] = ih δ jl δ(ω − ω ). In this approach, the bath operator in Eq. (2) is often chosen as [17] where the coupling function C ω (in general a d s × d tensor) is weighting the system-bath coupling at frequency ω. The p-2 Coupling Function From Bath Density Of States system operators couple isotropically to the bath if C ω C T ω = 1 ds C 2 ω . The scalar coupling function C ω is related to the bath spectral density J ω , which alternatively quantifies the effect of the environment on the system as J ω ∝ C 2 ω /ω [1,2]. The bath dynamics can be categorised [2] based on the low-ω exponent of the spectral density, J ω ∝ ω s , into three different classes, called Ohmic (s = 1), sub-Ohmic (s < 1), and super-Ohmic (s > 1).
The difference between wave vector approach and frequency approach is that at a fixed frequency ω, there is in Eq. (7) just one bath operatorX ω that couples to the system, while according to Eq. (5), the interaction is distributed over several wave vector modes k with weighting factors ξ k , their number being set by the DOS D ω (see Fig. 1).
We now want to address the question of the connection between the DOS D ω and the coupling function C ω . To achieve this, we consider one relevant quantity in both approaches and equate the corresponding formulas. In the following, we choose the memory kernel K which encodes the response of the bath to the system operatorŜ. Note that the choice ofB in Eq. (5) restricts the discussion to the linear response of the bath, as is reasonable for a bath that is thermodynamically large [1,2].
Memory kernel in both approaches. -To find an explicit relation in the wave vector approach for the dynamics of the bath operatorB W V in Eq. (5), the starting point is the equation of motion forb k , whose retarded solution contains two termŝ Therefore, the time evolution of the bath operator can be written asB W V (t) =B W V induced (t) +B W V response (t). The first term represents the internally evolving bath which is contains information about the system's past trajectory, where K W V (t − t ) is the memory kernel (a tensor), Here, the ξ k have been expressed by the unit polarisation vectors k [see after Eq. (5)] and Θ(t − t ) is the Heaviside function, which ensures that the bath responds only to the past state of the system, i.e. t < t. For large volume V , the summation over k in Eq. (11) can be transformed into a frequency integration as in Eq. (4). The projection on polarization vectors, averaged over an isofrequency surface Ω, is taken into account by a (d s × d s ) positive Hermitian matrix g 2 where the matrix M ω is normalized to unit trace and g 2 ω is a scalar. With these notations, the memory tensor in the wave vector approach is Turning now to the frequency approach, the dynamics of the bath operatorX ω in Eq. (7) follows a driven oscillator equation Its exact solution iŝ where Green's function. Inserting this solution in Eq. (7) leads again to induced and response evolution parts given, respectively, Comparing with Eq. (10) one can identify the memory kernel tensor in the frequency approach as Coupling function C ω versus DOS D ω . -Since Eqs. (12) and (16) describe the same memory effects, we may set them equal, leading to This relation links the system-bath couplings in the two approaches, i.e. the DOS D ω is proportional to the Hermitian "square" of the coupling function C ω . This is the first result of the paper. The result (17) may be applied to any quantum system that interacts linearly with a bosonic bath. For instance, magnetic materials in which spinsŜ relax in contact with a phonon reservoir have been studied extensively [13,17,31]. The noise-affected occupation of fermionic modes in a double quantum dot [32], and the behaviour of an impurity in a Bose-Einstein condensate environment [33] are other examples.
Note that in Eq. (17) the dimension of the system is either smaller or equal to the dimension of the bath, i.e. d s ≤ d.
A rectangular (d s × d) coupling matrix C ω may model a graphene-on-substrate structure, where the electronic system (d s = 2) is in contact with a 3D phononic bath [34]. An example for equal dimensions is a 3D spin vector that couples to a 3D phononic environment [17].

p-3
Specific examples. -In this second part of the paper, we wish to use Eq. (17) to obtain coupling function estimates from experimentally measurable quantities. To do so we have to drop generality and make a number of simplifying assumptions.
First, we assume isotropic coupling to an isotropic bath, and set C ω = 1 ds C ω with scalar C ω , and M ω = 1 ds /d s .
Second, for simplicity, we assume a frequency-independent g so that the frequency-dependent impact of the coupling is captured by D ω alone. This is a common approximation for quantum optics systems [35,36], while examples of condensed matter systems where this approximation holds over a range of frequencies are limited [2]. Whenever a non-trivial g ω is known for a specific context, such as a power law behaviour ∝ ω p , this can be included in Eq. (17) separately from the DOS' ω-dependence. To establish such g ω requires microscopic models for specific physical situations, which make several approximation steps. For example, a derivation of the electron-phonon interactions in quantum dots [25] was given in [24]. These assumptions reduce Eq. (17) to the scalar equation with system dimension d s . We will base the following discussion of examples for the coupling functions C ω on this simpler scalar form.
Debye approximation. -In condensed matter physics, the Debye model is used to describe the phonon contribution to a crystal's thermodynamic properties. It assumes an acoustic dispersion, i.e. ω = c|k| with an averaged sound speed c, resulting in 3D in [21] D Deb Here ω D is the Debye frequency, i.e. the maximum bath frequency, which in practice is taken to be near the edge of the Brillouin zone. For example, for gold, see Fig. 2 The scaling of C Deb ω implies that the spectral density J(ω) ∝ C 2 ω /ω is Ohmic, i.e. J(ω) ∝ ω. Hence, the 3D Debye model with constant coupling g in the wave vector approach captures the same relaxation dynamics as an Ohmic bath in the frequency approach.
Beyond 3D cubic lattices, D ω will depend on the dimensionality and lattice symmetry. What happens if the lattice is effectively two-or one-dimensional? To answer this, let us imagine a dD isotropic lattice with volume V = N a d . The volume element of such a lattice in k-space corresponds to   )) and two-peak Lorentzian DOS (blue solid line, Eq. (22)) fitted to a measured phonon DOS for gold (red dots) reported as in Ref. [37]. The Debye frequency for gold is ωD/2π = 3.54 THz given in Ref. [21]. Fit specified peak frequencies ω0,j, widths Γj and peak ratios Aj/A1 are given in Table 1. The grey dashed lines separate three frequency regimes discussed in the main text. (b) Memory kernels K(t − t ) corresponding to Debye DOS and two-peak Lorentzian DOS.
Analogously to the 3D lattice, using the acoustic dispersion with an averaged sound speed c, one finds the dD Debye DOS Via Eq. (18) we obtain the power-law C ω ∝ ω (d−1)/2 for the corresponding coupling functions which implies spectral densities J(ω) ∝ ω d−2 . Thus, isotropic baths in 2D or 1D behave in a distinctly sub-Ohmic way.
Inferring coupling functions from DOS data. -Here we wish to go beyond the conceptually useful Debye model, and fully specify the functional form of C ω , given experimentally accessible DOS data that characterise the phononic environment.
A generic feature of real materials is a structured DOS, which shows several peaks [37,38]. Sums of Lorentzian or Gaussian functions are two convenient candidates to approximate such peaky shaped densities [39]. Here, we fit experimentally measured DOS for gold [37] (and iron [38] in SM) and theoretically computed DOS for YIG [40] to a function consisting of multiple Lorentzians, The fits, see Figs. 2 (a), 3 and figure in SM, reveal the material specific peak frequencies ω 0,j , peak widths Γ j and peak ratios A j /A 1 , see Table 1 and tables in SM, while the first peak amplitude A 1 remains undetermined. Fixing A 1 would require information additional to the DOS, such as the system's relaxation rate due to interaction with the phonon bath. Note that phonon DOS are generally slightly temperature dependent [38]. Hence the fit parameters in Eq. (22) will be (usually weak) functions of temperature, a dependence that only matters when a large range of temperatures is considered. The peak widths in Eq. (22) determine a characteristic memory time 1/Γ j . However, beyond this single timescale number, the functional dependence of the memory is fully determined by the kernel Eq. (12), which for multi-peak Lorentzians is proportional to with ω 1,j = ω 2 0,j − Γ 2 j /4. The degree of memory introduced by this kernel into a system's dynamics could be quantified in terms of several non-Markovianity measures, see e.g. [41][42][43][44].
For gold, Fig. 2 (a) shows the phonon DOS measured by Muñoz et al. [37], together with our two-peak Lorentzian fit. The fit gives good agreement in all frequency regimes, with a slightly slower decay in region III than the measured DOS. For a system coupled to phonons in gold, the peak widths (see Table 1) immediately imply a characteristic memory time in the picosecond range. The relevant kernel is shown (blue) in Fig. 2 (b) for the two-peak fitted DOS of gold shown in (a). Using the Debye model instead would give a qualitatively different behaviour: the pink curve shows a distinctly slower long-time tail, due to the sharp cutoff at the Debye frequency. Note also that without any cutoff, the kernel would be K(t − t ) ∝ ∂ t δ(t − t ), implying no memory [17]. In contrast, the Lorentzian fit (blue) provides a quantitatively accurate memory kernel.
Our approach may provide a more realistic picture of the magnetization dynamics based on actual material data. YIG Theo.  [40]. The grey dashed line shows a single-peak Lorentzian fit. The fitted peak frequencies ω0,j, widths Γj and amplitude ratios Aj/A1 can be found in Table 2 in the SM. [45,46] is a typical magnetic insulator in which the relaxation of a spin DoFŜ is dominated by the coupling to phonons [47], similar to magnetic alloys like Co-Fe [48], while in metallic materials, the coupling to electrons is more relevant [49]. Fig. 3 illustrates a theoretically computed DOS for YIG [40] with a fit that contains eighteen Lorentzians. (Parameters are displayed in Table 2 in the SM.) In this fit, a few negative amplitudes A j in Eq. (22) are needed to reproduce the gap near 16 THz, however, the total D ω remains positive. Using additional information of the typical Gilbert damping parameter for this material [50], also the peak amplitude A 1 can be determined (see the SM).

18-Lor Fit
More generally, via Eq. (18) the parameters of the multipeak DOS (22) immediately specify the functional form of the coupling C ω of a system to a phononic bath in real materials. This second result of the paper will be useful for modelling the Brownian motion of spins [17,51] and in applications such as quantum information processing with solid-state spin systems [52].
Conclusion. -We have derived the general relation (17) that translates the function C ω , determining the coupling of a generic system to a bosonic bath at various frequencies, into the density of states D ω of the latter. This was achieved by evaluating the memory kernel of dynamical bath variables in two equivalent approaches. Several applications of the relation were then discussed. We demonstrated how for systems damped by phonons in 3D with a frequency-independent g, Debye's quadratic DOS captures the same physics as a linear coupling function C ω which corresponds to an Ohmic spectral density. Secondly, we have established how to infer C ω from the measured DOS of a material, such that it reflects the specific properties of the material. Given that real materials have densities of states with multiple peaks, the typical picture which emerges from our general relation (17) is that the coupling function is non-Ohmic and memory effects in the system dynamics become important. The corresponding time scales (in the ps range, e.g., for gold in Fig. 2 (b)) can be conveniently determined by fitting multiple Lorentzians to the bath DOS.
Future work could address how to extend relation (17) to systems interacting with multiple independent baths. This should be suitable for non-equilibrium settings involving different temperatures [53], as used in heat transport [54]. The impact of memory may also change the behaviour of systems like superconducting qubits or two-level systems that are in contact with two baths [55,56].   [57]. Experimental data are shifted vertically for clarity. The fit parameters for the Lorentzian peaks are given in Table 2.
Supplementary Material. -Fig. 4 shows the phonon DOS for iron, a well-known magnetic conductor. In this figure, we use experimental data reported in Ref. [38] and microscopic theoretical results calculated in Ref. [57]. Fig. 4 (22) with a single peak]. The Debye cut-off frequency is taken from the Debye temperature (420 K) according tohω D = k B T D . Both models reproduce the quadratic scaling at low frequencies, but the Lorentzian also captures the first two peaks. It does not properly capture the data for large frequencies, however.
In Fig. 4 (b), a three-peak Lorentzian [Eq. (22)] fits the measured DOS spectrum remarkably well. It slightly deviates at higher frequencies because the Lorentzians decay less rapidly. Therefore, one can conclude that all three models are reliable fits to the measured data at low frequencies. However, if details are required in a broader frequency range (i.e. at shorter time scales), the three-peak Lorentzian DOS will be the better choice. Fig. 4(c) illustrates the additional structure in the DOS that becomes visible when fitting to the theoretical DOS of Ref. [57]. Here, we fitted a five-peak Lorentzian [see Table 2] to reproduce two additional shoulders. Compared to the experimental data, these are probably hidden by the finite resolution of the measurement. Taking experimental data with such a broadening at face value, one would therefore produce a fitted coupling function with a somewhat shorter memory time compared to the real material. Fit parameters are given in Table 2.
In Table 3, we give the fitting parameters for the DOS of the material Yttrium Iron garnet (YIG). This is a magnetic material which has phonon-damped spin dynamics. Here, an absolute magnitude for the peak amplitudes can be extracted using the known Gilbert damping parameter η 5 × 10 −4 [50] and the electron gyromagnetic ratio γ e 28×10 9 Hz/T [58]. We consider for the system operatorŜ in Eq. (2), as in Ref. [17], a spin vector s multiplied by the gyromagnetic ratio γ e (with |s| =h/2). For the eighteen-peak fit, in Fig. 3 of the main paper, we determine the absolute peak height to be A 1 71.14 (rad THz T) 2 /meV. The other fit parameters are given in Table 3. For the single-peak fit, which is easier to use in simulations, the absolute peak height is A 1 1194.19 (rad THz T) 2 /meV and the other fit parameters are given in Table 3. Table 3: Fit parameters of single-peak and eighteen-peak Lorentzians -shown in Fig. 3 of the paper -matched to the theoretical DOS for YIG reported in Ref. [40].