Driven Spin-Boson Luttinger Liquids

We introduce a lattice model of interacting spins and bosons that leads to Luttinger-liquid physics, and allows for quantitative tests of the theory of bosonization by means of trapped-ion or superconducting-circuit experiments. By using a variational bosonization ansatz, we calculate the power-law decay of spin and boson correlation functions, and study their dependence on a single tunable parameter, namely a bosonic driving. For small drivings, Matrix-Product-States (MPS) numerical methods are shown to be efficient and validate our ansatz. Conversely, even static MPS become inefficient for large-driving regimes, such that the experiment can potentially outperform classical numerics, achieving one of the goals of quantum simulations.

This system not only leads to various equilibrium phases, such as spin-boson Mott [18] and Anderson [19] insulators, or spinboson Ising magnets [20], but also to non-equilibrium phenomena, such as emerging causality [21]. So far, LLs have not been explored in this context.
To fill in this gap, we study a driven spin-boson Hamiltonian H = H 0 +V , with H 0 describing the uncoupled dynamics where J s (J b ) represents the spin (bosonic) tunnelling strengths, ω b is the bosonic on-site energy (h = 1), and Ω d , q d are the driving amplitude and momentum, respectively. We also introduce a spin-boson interaction of strength U, namely Instances of this spin-boson model can be implemented with trapped ions or superconducting circuits, but we shall postpone the experimental details, focusing first on its physical content. In Eq. (1), the spin part describes the well-known isotropic XY model, whose groundstate corresponds to a halffilled Fermi sea of spinless fermions [22]. In the absence of driving, the bosons minimize their energy in the global vacuum ω b > 2J b , such that the spin-boson interaction (2) does not induce any spin-boson LL behaviour. As shown below, as the driving is switched on, the bosonic mode of momentum q d gets macroscopically populated. As a consequence, the spinboson interaction induces non-trivial correlations, which will be shown to correspond to a new phase: a spin-boson LL. Variational ansatz.-According to the above discussion, we introduce spinless fermions through a Jordan-Wigner map- [23], and build our ansatz gradually. First, we note that the macroscopic population of the driven mode resembles the Bogoliubov theory of superfluids, where a q = ∑ i e −iqdi a i / √ N → √ N d δ q,q d + δ a q , such that δ a q are the Gaussian fluctuations for quasi-momentum q ∈ [−π/d, π/d), and δ q,q d is the Kronecker delta fixing the mode with macroscopic occupation N d [24]. This can be accounted variationally by the ansatz , where {α q , ψ f } are the variational parameters, and |Ψ G is a fiducial state. Second, in order to account for the bosonfermion correlations caused by the interaction (2), we use a polaron transformation [25] for the fiducial state where { f iq } are such that the boson-fermion correlations are minimized in the unitarily-transformed picture, |0 a is the boson vacuum, and |ψ f is a general fermionic variational state. The variational minimization yields α q = √ N(Ω d /ω q )δ q,q d to lowest order in the interaction strength [26], where we have introduced the bosonic dispersion ω q = ω b −(U +2J b cos qd). This is the analogue of the macroscopic population in the Bogoliubov theory [24], which motivates neglecting contributions of quartic terms. Within this approximation, one can see that the interaction (2) tries to create/annihilate Gaussian excitations conditioned on the fermion number n i = c † i c i through a Holstein-type coupling [27]. This will lead to second-order processes where distant fermions exchange Gaussian bosons virtually, resulting in a fermion-fermion long-range interaction [28]. Additionally, as expected from the polaron transformation, a bosonic cloud will dress the fermions modifying their tunnelling (i.e. band narrowing).
To make this description more quantitative, we specify the polaron parameters Then, the variational minimization sets |ψ f as the groundstate of the fermionic Hamiltonian As announced above, we obtain boson-mediated interactions To guarantee the absence of negative frequencies and thus the stability of bosons, we focus on λ = 2J b /(ω b − U) < 1. By using the displacement-operator algebra, the binomial theorem, and some Taylor series [26], we find where the explicit dependences of η 1 , and V 0 , ξ 0 on the experimentally tunable parameters are listed in [30,31]. As announced, the dressed tunnelling gets exponentially suppressed as the bosonic driving increases, and the interactions are longranged. For instance, letting ω b ≈ 2J b , the length scale ξ 0 diverges, such that the interaction does not decay with distance. Conversely, for ω b 2J b , interactions decay very rapidly. In contrast to the dressed tunnelling, the interaction strength increases with the driving. Additionally, its attractive/repulsive character oscillates along the chain, which corresponds to frustration effects in the original spin model.    To test the correctness of Eq. (4), we compare it to the corresponding exact expressions evaluated numerically. In Figs. 1(a)-(b), we see that (i) the renormalization of the tunnelling and the exponential decay of the interactions (4) are very accurate. Therefore, the interaction range can be tuned over ξ 0 ∈ (0, ∞) by controlling λ ∈ (0, 1). (ii) The dependence of the degree of frustration on the driving momentum (4) is also very accurate: while Fig. 1(a) corresponds to unfrustrated attractive interactions, Fig. 1(d) shows that alternating attractive/repulsive interactions occur for q d = π/d. Finally, (iii) the ratio of the interactions to the dressed tunnelling can be tuned across |V 0 |/J s = 1 by controlling a single parameter, the driving strength Ω d ( Fig. 1(c)). This is quite remarkable as we started from the constraint U J s , J b imposed by the Bogoliubov theory of the bosons. Nonetheless, the role of fermion interactions is enhanced by increasing the driving strength.
Bosonization and MPS.-The nearest-neighbor limit of the model (3) is a paradigm of LLs [32]. Thermodynamic quantities given by the Bethe ansatz can be combined with LLtheory to obtain various correlation functions [33]. Additionally, the numerical Density-Matrix-Renormalization-Group (DMRG) gives accurate predictions in this limit [34,35] used to benchmark the LL [36]. The situation gets more involved for the full model (3), since Bethe-ansatz integrability is lost, and DMRG with long-range interactions is more intricate (i.e. typically, models with only a few neighbors are studied [36]).
An analytical LL-theory of the long-range model (3) can be obtained by phenomenological bosonization [37]. We use instead the constructive approach [38], which allows us to find a constraint on the interaction range [26]. Provided that (i) V i− j ≤ const/|i − j| at large distances, which is fulfilled by (4) except for λ → 1, and (ii) V 0 J s , the low-energy excitations are described by two bosonic branches characterized by the sound velocity The new operators d q± are related to particle-hole excitations of the original fermionic system by a squeezing transformation [26] that depends on the Luttinger parameter In the absence of driving, we get V 0 = 0 and recover the noninteracting value K = 1. When we switch it on, it is possible to tune K ≶ 1 over a wide range of values as displayed in Fig. 2. Following this discussion, the fermionic part of our variational groundstate |ψ f is the vacuum of the new bosonized modes |0 d , and we thus obtain a spin-boson LL groundstate We can now calculate any connected two-point correla- which coincide exactly with those of the Heisenberg-Ising or XXZ model [32], albeit with a different Luttinger parameter (7) due to the long range and possible frustration in (3). The off-diagonal spin correlators are The symbols are numerical data, and the red solid lines represent the power-law decay obtained from fitting the LL-predictions (9)-(11) at distances up to |i − j| = 100.
which display the power law of the bare XXZ model [32] with an additional distance-dependent renormalization e −η due to the bosonic cloud dressing the spins [30]. At long-distances, this polaron effect is constant C σ + σ − ( ) ∼ e −η 0 /2 −1/2K , and does not modify the power-law exponent. The off-diagonal bosonic correlators can also be understood from a polaron perspective which is a sum of diagonal spin correlations (9) with a polaron weight C ,˜ exp(−| −˜ |d/ξ 0 ) exponentially suppressed at large distances, where C ,˜ is listed in [39]. For ξ 0 d, all terms except˜ = are negligible, and we thus find a powerlaw decay only determined by the Luttinger parameter.
Finally, we have concrete power-law predictions (9)-(11) that can be numerically benchmarked. We use DMRG-type methods based on matrix product states (MPS) for the thermodynamic limit [40]. By implementing the original short-range spin-boson model (1)-(2), we avoid the intricacies associated to the long-range model (3) mentioned above. The results displayed in Fig. 3 agree with the aforementioned power-law decay at intermediate distances, and depart at longer distances due to technical limitations in the MPS dimension. These results confirm the validity of our ansatz.
Experimental LLs.-Experiments on this versatile LL can access different regimes by controlling a single parameter, the driving (Fig. 2). As a bonus, regimes that differ markedly from K = 1 require very large drivings, which increases both the bosonic population |α q d | 2 ∝ Ω 2 d and the dimension of the MPS. Therefore, MPS simulations, limited by available computing resources, will eventually cease to be trustworthy. Instead, the experiment would act as a reliable quantum simulator [5] capable of beating its classical MPS counterpart. Trapped-ion (TI) and superconducting-circuit (SC) architectures meet the requirements for this LL quantum simulator. We first focus on the bosonic degrees of freedom. For laser-cooled linear chains of TIs [41], the bosons are the transverse vibrational excitations of each ion, and display The driving is due to an oscillating potential at one of the electrodes, which has a frequency ω d = ω t + ∆ with detuning ∆ from the transverse trap frequency ω t /2π ∼1MHz. This leads to Ω d /2π ∼1-100kHz, and q d = 0. To obtain q d = 0, one should use instead the ac-Stark shift of a pair of lasers with beatnote ω 1 − ω 2 = ω t + ∆, such that q d = (k 1 − k 2 ) · e z depends on the laser wavevectors projected along the chain, and Ω d /2π ∼1-10kHz. Note that the crossed-beam ac-Stark shifts must coincide for each atomic level forming the spin. For cryogenically-cooled SCs [3], the bosons are the photonic excitations of superconducting resonators of frequency ω r /2π ∼1-10GHz, which are capacitively-coupled yielding nearest-neighbor tunnelings A microwave drive, detuned from the resonator frequency ω d = ω r + ∆, is injected in each resonator, and its amplitude/phase is individually controlled by quadrature mixers providing Ω d /2π ∼1-100 GHz and a site-dependent phase ϕ i = q d di.
We now introduce the spin-1/2 degrees of freedom. For TIs, two levels are selected from either the hyperfine groundstate or a dipole-forbidden optical transition [2]. Using lasers, a Jaynes-Cummings coupling gσ + i a i e −iδt can be introduced, where g/2π ∼1-50kHz, and δ /2π ∼0-0.1MHz is the redsideband detuning. For SCs, two states with different values of charge/flux variables are separated from the rest by exploiting the Josephson effect [42]. These spins can be coupled to the resonator photons via the same Jaynes-Cummings terms, albeit reaching g/2π ∼1-100MHz and δ /2π ∼0.1-1GHz. In the dispersive regime g δ , and setting J b ω t , ω r , these Jaynes-Cummings couplings are highly off-resonant, and lead to spin-spin J s = g 2 J b /δ 2 and spin-boson U = −g 2 /δ interactions to second order, with the peculiarity that TIs yield a dipolar decay J s /|i − j| 3 [19], while nearest-neighbor couplings are the leading contribution in SCs.
Let us note that the equivalence of the TI or SC Hamiltonian to (1)-(2) occurs in a rotating frame, where ω b = ∆ is the detuning of the driving, and not the trap or resonator frequencies ω t , ω r . In order to avoid spurious cross terms of the driving with the Jaynes-Cummings terms, we should also impose g |δ − ∆|. We remark that these realistic experimental parameters allow for a wide tunability of the LL parameters.
Once the driven spin-boson Hamiltonian is implemented, and the groundstate adiabatically prepared for a certain K of Fig. 2, one must probe the two-point correlators in Fig. 3. TIs excel at measuring any spin correlator through siteresolved spin-dependent fluorescence [2], whereas SCs are better suited to measure the photonic correlations by collecting the output from the cables used to drive the resonators. We thus conclude that either TIs or SCs are promising candidates to realize this spin-boson LL-liquid.
Conclusions.-We have presented a theoretical study for a new class of spin-boson LLs based on a variational bosonization approach benchmarked by MPS numerics, and proposed its implementation with TI and SC technologies. This model offers a flexible platform to test qualitatively LL-theory, and displays certain regimes where the LL quantum simulator can beat MPS numerics on any classical computer.
The authors acknowledge support from EU Project PROMISCE, Spanish MINECO Project FIS2012-33022, and CAM regional research consortium QUITEMAD S2009-ESP-1594. coupling expansion that treats the electron hopping as a perturbation J s ω b ,U, leads to fermion-fermion interactions of strength J 2 s /ω b U by a sort of super-exchange [29]. This limit is not consistent with our treatment of the fermion-boson coupling U J s , J b . Yet, the dispersive nature of our bosons allows for long-range interactions mediated by virtual boson exchange rather than by super-exchange.
[30] η = [32] A. Luther and I. Peschel, Phys. Rev. B 12, 3908 (1975 [39] C ,˜ = − In order to minimize over the variational family of groundstates introduced in the main text, and given by we first apply the displacement operator U B a q U † B = a q + α q , and proceed in the spirit of the Bogoliubov theory of superfluids [24]. Accordingly, we linearise the bosonic density in the Fermi-Bose interaction assuming sufficiently weak interactions, and obtain directly the variational parameter α q = √ N(Ω d /ω q )δ q,q d . For this particular choice, the linear terms in the Fermi-Bose Hamiltonian vanish, and we obtain the following variational problem As noted in the main text, the spin-boson interaction reduces to a Holstein-type coupling [27] within this variational formalism, which corresponds to the last term in the above equation. We note that, contrary to the Holstein model, the Fermi-Bose coupling strength here depends on momentum and the bosons are dispersive. We can now identify the parameters of the polaron unitary by following the premise that the Fermi-Bose coupling must vanish in the transformed picture. This leads to the polaron parameters Using the polaron transformation rules where D a (α) = exp αa † −α * a is the bosonic displacement operator, and taking the expectation value over the bosonic vacuum, one reduces the variational problem to a purely fermionic one Here, H f is an effective long-range XXZ Hamiltonian described in Eq. (3) of the main text, rewritten here for convenience where n i = c † i c i is the fermion number operator, and we have introduced the following microscopic parameters In order to obtain closed expressions for the parameters of this model, let us introduce λ = 2J b /(ω b − U) and assume that λ < 1 to guarantee the absence of negative boson frequencies. We can then use the geometric Taylor series, such that together with the binomial theorem after introducing the binomial coefficients and apply the identity ∑ q e iqx = Nδ x,0 . This allows us to express the interaction strengths as a single series, which can be exactly summed by considering a number of combinatorial identities. In this way, we find We now use the identity for the overlap of coherent states |α = D a (α) |0 a , namely α|β = e − 1 2 (α * β − β * α) e − 1 2 |α − β | 2 . This allows us to express the dressed tunnelling as an exponential renormalization of the bare tunnellingJ s = J s e −χ with This sum can be evaluated following a similar procedure as above. We first use the binomial Taylor series, such that This expression can be analytically summed by using again Eq. (19) , together with some combinatorial identities, such that which coincides with Eq. (4) of the main text with the parameters listed in [30]. The validity of this derivation has been confronted to the exact numerics of Fig. 1(b).

Constructive bosonization of the spin-boson lattice model
In this part of the Supplemental Material, we present the details for the derivation of variational bosonization ansatz of Eqs. (5)- (8) in the main text, and the corresponding two-point correlators in Eqs. (9)- (11).
We start by setting ω s = 0 in Eq. (3), which is reasonable for sufficiently weak interactions such thatJ s ω s . Assuming periodic boundary conditions, the kinetic part of the fermionic Hamiltonian H f = K f +V f can be written as where the phase of the dressed tunnelling e iarg(J s ) has been absorbed in the fermionic operators via a U(1) gauge transformation. This yields a band structure ε s (q) = 2|J s | cos(qd), where groundstate is obtained by filling all negative-energy levels −k F ≤ q ≤ k F , with k F = π/2d, and the low-energy excitations correspond to right-and left-moving fermions η ∈ {R, L} with momentum close to ±k F , respectively. In the continuum limit, we let d → 0 and N → ∞ such that the length L = Nd remains constant. The low-energy properties are described by a continuum field theory of slowly-varying fields for the left/right-moving fermions After Taylor expanding the slowly-varying fieldsΨ η (x ± d) ≈Ψ η (x) ± d∂ xΨη (x) + 1 2 d 2 ∂ 2 xΨη (x) ± . . . , the kinetic energy corresponds to a 1 + 1 Dirac quantum field theory for massless fermions where rapidly-oscillating terms e ±i2k F x = (−1) x/d average out under the integral, higher-order derivatives can be neglected for d → 0, and we have introduced the Fermi velocity that plays the role of the speed of light in the continuum limit 2|J s |d → v F . The first step of the constructive bosonization [38] is to extend this description, which is expected to be valid around the Fermi points, to all possible momenta k = 2π L n k where n k ∈ Z (i.e. the fermionic spectrum is linearized for all momenta), such that where we have introduced the normal ordering : ( ) : with respect to the groundstate of the free theory. In this way, we can focus on certain low-energy excitations above the Dirac sea of filled negative-energy states by defining operators associated to the particle-hole excitations of momentum q = 2π L n q > 0 for the left/right fermionic branches, namely which turn out to be bosonic. From the commutation rules of these operators with the Dirac Hamiltonian (26), it follows that in the thermodynamic limit L → ∞. Therefore, the low-energy excitations of the free fermionic model are described by two bosonic branches for the particle-hole excitations around each Fermi point (i.e. bosonization).
It is customary to express the Hamiltonian in terms of bosonic fields where a > 0 is a regularization constant that cuts off large momenta to ensure convergence. Since we are interested in low-energy properties, this cutoff does not change the physics, and we can set it to zero at the end of the calculations. The kinetic energy, and thus the 1+1 Dirac Hamiltonian, can then be expressed as a free bosonic field theory where we have also introduced the so-called dual fields Θ(x) = (φ R (x) + φ L (x))/ √ 2, and Φ(x) = (φ R (x) − φ L (x))/ √ 2. The equivalence of the fermionic (26) and bosonic (30) Hamiltonians suggests the existence of a direct mapping between fermionic and bosonic fields. To obtain the correct fermionic anticommutation relations, one has to introduce the so-called Klein The bosonization identity [38] then relates the fermionic and bosonic fields in the thermodynamic limit as follows Equipped with this operator identity, we can bosonize the interaction (3) which, in terms of the fermion fields (25), reads where we must define V d →Ṽ in the continuum limit in analogy to the Fermi velocity below Eq. (26). The bosonization of these interactions is more intricate, as one must avoid possible divergences by normal ordering. Besides, the long-range tail of the interaction allows for → ∞ in the continuum limit, such that special care must be taken in the truncation of the Taylor series of the fermionic fields for d → 0. This will impose restrictions on the interaction range tractable by bosonization.
Assuming that the interactions are small enoughṼ v F , such that the slowly-varying fields (25) are only slightly perturbed, we identify the interaction terms V f = V 0 f +V 2k F f +V 4k F f , where we have again neglected rapidly-oscillating contributions , , These terms are expressed in terms of the bosonic fields by means of the identity (31). To eliminate possible divergences in the first term, we use the so-called point-splitting regularisation : , where |0 is a reference state without particle-hole excitations, and ε → 0 such that we can Taylor expand the field operators. Accordingly, we find : We can now Taylor expand the bosonic fields φ η (x + d) ≈ φ η (x) + d∂ x φ η (x) + . . . and, provided that the interactions decay fast enough, namely the corresponding fermionic quartic term can be expressed as a quadratic bosonic one in the d → 0 limit ∑ η,η This means that only interactions that decay faster than, or equal to, the Coulomb interaction (35) can be treated via the constructive bosonization. This is crucial to neglect higher-order derivatives, which scale with C n−1 d n → 0 in the continuum limit provided that Eq. (35) is fulfilled, even for → ∞ in the thermodynamic limit. It is also clear that the bosonization predictions will be more accurate the faster the decay is, since the contributions of the higher-order derivatives will be less and less important. We can proceed similarly with the second term V 2k F f . Using the Baker-Campbell-Hausdorff formula with [ϕ † η (x), ϕ η (x )] = −δ η,η (δ η,R log 2π L (a + i(x − x ) + δ η,L log 2π L (a − i(x − x ) ), which can be checked from Eq. (29), we can write the second interaction term in the a → 0, and d → 0, limit as follows After Taylor expansion, and making use of the constraint on the interaction decay (35), we find The last term V 4k F f , which corresponds to the so-called Umklapp scattering, is also bosonized by using the Baker-Campbell-Hausdorff formula and a Taylor expansion, such that we find in the limit d → 0, we obtain ∑ η =η where the definition of the cosine operator in the last term must incorporate the combination of Klein factors. Using all these expressions, it is possible to rewrite the fermionic Hamiltonian as a sine-Gordon quantum field theory composed of the Hamiltonian of a Luttinger liquid H LL , and a non-linearity due to the Umklapp scattering H U , where u is a renormalized Fermi velocity, K is the so-called Luttinger parameter, and g U the Umklapp interacting strength. For the small interactionsṼ v F considered here, these parameters are obtained through the above expressions which coincide with a phenomenological bosonization [37] up to a different choice of fermionic algebra and dual fields. Using the exact expression of the mediated interactions in Eq. (4), we can find the analytical expressions for the Luttinger parameters by using the geometric Taylor series and some trigonometric identities, such that . (42) These expressions yield the Luttinger parameters in Eqs. (6) and (7) of the main text. For small-enough interactions, the sine-Gordon non-linearity is irrelevant [7], and the bosonized groundstate is solely determined by the Luttinger-liquid Hamiltonian H f ≈ H LL in Eq. (40). To obtain the corresponding groundstate, let us express the dual fields in terms of new bosonic fields Θ( are defined in terms of two species of bosonic creation-annihilation operators d † q± , d q± . By direct substitution, one finds that in analogy to the free theory, the spectrum of the full interacting theory can again be described by two linear bosonic branches with a renormalized Fermi velocity v F → u (41), which proves Eq. (5) in the main text. Therefore, the fermionic part of the variational groundstate is the vacuum of the new bosonic modes |ψ f = |0 d , and directly leads to the complete variational bosonization ansatz of Eq. (8), namely |Ψ G ({α q , ψ f }) = U † B U † P |0 a ⊗ |0 d . Once the variational bosonization ansatz is known, we can embark upon the calculation of the two-point correlators. We start from the simplest one, the diagonal spin-spin correlators C σ z σ z ( ) = ∑ i ( σ z +i σ z i − σ z +i σ z i )/N. Making use of the Jordan-Wigner mapping, one finds and since both unitaries commute with the fermion operators, it follows that C σ z σ z (x) = d 2 π 2 0 d | : Ψ † (x)Ψ (x) :: Ψ † (0)Ψ (0) : |0 d , where we have already taken the continuum limit in Eq. (25). This expectation value coincides exactly with that of the bosonized XXZ model [32], which has become a textbook example [7]. Using the above bosonization relations, and the Klein-factor algebra, one finds C σ z σ z (x) = 2d 2 ). The first term can be evaluated easily after expressing the dual field in terms of the new bosonic creation-annihilation operators (43), while the second term requires the additional identity for Gaussian states 0 Altogether, one obtains in the thermodynamic limit L → ∞ which corresponds to Eq. (9) after setting x = d, and 2k F d = π.
For the off-diagonal spin correlators C σ + σ − ( ) = ∑ i ( σ + +i σ − i − σ + +i σ − i )/N, we proceed analogously to find C σ + σ − ( ) = Π q,q 0 a | D a q ( f 1+ ,q )D a q (− f 1,q ) |0 a 0 d | σ + The bosonic contribution is due to the polaron cloud dressing the spins, and can be calculated in analogy to the renormalization of the tunnelling in Eq. (21) by using the coherent-state algebra. The spin contribution reduces once more to that of the bare XXZ model, which requires additional care due to the Jordan-Wigner string when expressed in terms of fermions [7]. Anyhow, it can be expressed again as a collection of expectation values of exponentials of the dual fields in the continuum limit, which can be calculated exactly using the same procedure as above. Altogether, this leads to whereχ(x) = ∑ q u 2 0 Nω 2 q 1 − e i(q−q d )x . This function can be evaluated by using once more the binomial series together with combinatorial relations, and after setting x = d, and 2k F d = π, we find precisely the expression in Eq. (10) of the main text.
Moving onto the bosons, let us address the off-diagonal correlators C a † a ( ) = ∑ i ( a † +i a i − a † +i a i )/N. Using the variational bosonization ansatz, we find that such correlators depend on the spin structure factors S σ z σ z (q) as follows To evaluate this expression, we define the inverse Fourier series 1/Ω 2 j = ∑ q e iqd j /ω 2 q √ N, which can be evaluated again making use of the binomial series and binomial theorem 1 Using this expression and the identity ∑ q e iqx = Nδ x,0 , one can express the bosonic correlators as which leads directly to Eq. (11) upon substitution of the diagonal spin-spin correlators in Eq. (9) already derived in this Supplemental Material.