Kubo-Martin-Schwinger relation for an interacting mobile impurity

In this work we study the Kubo-Martin-Schwinger (KMS) relation in the Yang-Gaudin model of an interacting mobile impurity. We use the integrability of the model to compute the dynamic injection and ejection Green’s functions at finite temperatures. We show that due to separability of the Hilbert space with an impurity, the ejection Green’s in a canonical ensemble cannot be reduced to a single expectation value as per microcanonical picture. Instead, it involves a thermal average over contributions from different subspaces of the Hilbert space which, due to the integrability, are resolved using the so-called spin rapidity. It is then natural to consider the injection and ejection Green’s functions within each subspace. We rigorously prove by reformulating the refined KMS condition as a Riemann-Hilbert problem, and then we verify numerically, that such Green’s functions obey a refined KMS relation from which the original one naturally follows.

The notion of thermal equilibrium is one of the fundamental concepts of physics.In quantum many-body systems, thermal states are described by the density matrix ρ = e −β Ĥ /Z with trρ = 1, where Ĥ is the Hamiltonian of the system.The same operator Ĥ is responsible for the time evolution of the system.This double role of the Hamiltonian is crucially important in establishing the Kubo-Martin-Schwinger [1][2][3] (KMS) condition between the Green's functions tr (ρA(t − iβ)B(0)) = tr (ρB(0)A(t)) , where A, B are operators and A(t) follows the Heisenberg evolution A(t) = e i Ĥt Ae −i Ĥt .The KMS condition involves analytic continuation of the Green's function tr(ρA(t)B(0)) to imaginary times in a strip 0 < Im t < β.
Such a relation is central to the concept of thermal equilibrium.For example, the thermal density matrix ρ can be actually defined as the one for which the relation holds [3][4][5].Moreover, the KMS condition is an example of a detailed balance relation which guarantees stability of thermal equilibrium under fluctuations [6][7][8] and is also a foundation for fluctuation-dissipation relations as established in the original works [1,2].Finally, the relation can also be promoted to higher point functions [9].Whereas the KMS condition simply follows from cyclicity of the trace, it is generally very difficult to show this relation explicitly by evaluating thermal two-point functions in interacting many-body systems and comparing both sides of the equality.The aim of our work is to demonstrate such computations for an impurity Green's functions.
The thermal expectation values appearing in (1) involve averaging over different eigenstates.In the thermodynamically large system and under the equivalence between the grand canonical (GCE) and microcanonical (MCE) ensembles, they can be computed over a single eigenstate, tr (ρA(t)B(0)) = ⟨ρ|A(t)B(0)|ρ⟩, where the eigenstate |ρ⟩ is chosen such that its energy equals the average energy in the canonical ensemble.Our results show that unexpectedly the thermal average reduces to a single expectation value on one side of the KMS relation but not on the other.This makes the relation itself highly non-trivial.Without the trace, there is no simple derivation of it and it must rely on a more subtle structure of the thermal expectation values which we unravel.
The two important correlation functions for the impurity problems are injection and ejection Green's function (we define them precisely below).The two Green's functions are related by the KMS condition, see Fig. 1.However, as we show, the equivalence between GCE and MCE holds only for the injection Green's function and it is broken for the ejection.We attribute its breaking to the Hilbert space separability -the Hilbert space of states with a single impurity separates into different sectors -with each sector contributing a different value to the thermal correlator.Then the latter cannot be reconstructed by a single expectation value.Interestingly, we can define the injection and ejection Green's function between individual sectors of the Hilbert space.As we will show later, this refined Green's function obey a refined KMS relation.
The separability of the Hilbert space in the Yang-Gaudin model is ultimately related to its integrability.However, a similar mechanism called Hilbert space fragmentation [31,32] appears more generically and leads to breaking of the Eigenstate Thermalization Hypothesis describing the thermalization of closed quantum manybody systems [33][34][35][36][37][38][39].
The model and its correlation functions: The Yang-Gaudin model [40] describes a 1d gas of spinful fermions with spin polarized interactions, The Hamiltonian density is (3) and we consider repulsive interactions, g > 0. The total number of particles of each spin, N σ = ´dxψ † σ (x)ψ σ (x) is a constant of motion and the Hamiltonian can be diagonalized in a subspace with fixed numbers of particles of the two kinds.The relevant for the impurity problem are then subspaces with N spin up particles and 0 or 1 spin down particle.We denote them by H 0 and H 1 respectively.The dynamics in H 0 is then this of a free spinless Fermi gas, in H 1 it is of the free spinless Fermi gas with a single impurity.
We define two equilibrium dynamic Green's functions of the impurity where the trace is over either H 0 or H 1 .We note that these are not normalized Green's functions [41].To re-store proper normalization one needs to divide them by the partition function Z i = tr i e −β(H−µN ↑ ) .The two (normalized) Green's functions describe response of the system to injection (in) and ejection (ej) of the impurity respectively and can be measured in spectroscopy experiments [42][43][44][45][46].The KMS relation between the two Green's function was derived in [47,48] and in our notation takes the following form The extra factor can be seen as a ratio of Z 0 /Z 1 defined for the Green's functions (4), (5).
After this introduction we can now state our results.These confirm the microcanonical picture for the injection function and disprove it for the ejection function.In the former case we find with |ρ⟩ a representative state of the thermal equilibrium for free fermionic gas.Instead, for the ejection function there is a single remaining degree of freedom, the rapidity Λ, which is related to the way the impurity and the gas share the total momentum in the system, with the free energy F(Λ) associated with varying Λ, the corresponding partition function Z ej , and the Λ-resolved ejection function G ej (x, t, Λ).Whereas the microcanonical picture does not work for the whole correlator, once the value of Λ is fixed, the thermal average reduces to a single expectation value, with |ρ Λ ⟩ denoting a representative thermal eigenstate with fixed value of Λ.This is similar to the generalized ETH [49,50] appearing in the equilibration processes of integrable models.There, the expectation value can be represented by a single eigenstate once all the conserved charges are fixed.Here, it is sufficient to fix a single additional degree of freedom Λ.This observation is formalized by realizing the Hilbert space with a single impurity through a direct sum of subspaces.In a finite system, values of Λ are quantized and infinite.Denoting the possible values by Λ m with m ∈ Z we can formally write (10) This structure allows us to formally define the projection operator in different subspaces denoted H 1 (Λ), The projection operator can be now used to define a Λresolved injection function such that The two Λ-resolved correlation functions obey a refined KMS condition Prove of this relation is the main result of work.This new KMS relation implies that Λ acquires a thermodynamic meaning, its refines the concept of thermal equilibrium to be specified not only by the temperature and the chemical potential µ but also by the rapidity Λ.The Λ-resolved KMS relation implies then the stability of this generalized equilibrium state.Integrating ( 14) over Λ provides an alternative to [47,48] derivation of the KMS relation (6).

Bethe ansatz solution to the McGuire model:
We present now the relevant ingredients of the Bethe ansatz solution to the Yang-Gaudin model with a single spin down and refer to [51] for details.This special case is known as McGuire model [52,53].The system's eigenstates |{k j }, Λ⟩ are specified by a set of (N+1) rapidities {k j } together with an extra rapidity Λ.For a system of length L with periodic boundary conditions the rapidities k j are solutions to the Bethe equations where quantum numbers n j are integers and obey the Pauli principle and the phase shift is The allowed values for the rapidity Λ are obtained by requiring that the total momentum is quantized Therefore, the state of the system is characterized by a choice of quantum numbers {n j } and I and the rapidities Λ and {k j } follow from the Bethe equations.Because the choice of the quantum number I is independent of the other quantum numbers the whole Hilbert space separates into subspaces with fixed value of I, or equivalently, with fixed value of Λ as anticipated in (10).In a large system there is then a density of states a(Λ) associated to a given Λ, where σ(k) = 1/(1 + e β(k 2 /2−µ) ) is the Fermi-Dirac distribution.The density of states which is expected to appear in (8) and ( 13) is absorbed in the definition of the Λ-resolved Green's functions.Finally, the impurity free energy is [54,55] Impurity Green's functions: The program of computations of impurity Green's functions was initiated in [56] where the zero temperature static injection Green's function was computed.This was subsequently generalized to dynamic correlators [57] and to finite temperatures [58].Similar techniques were later applied to determine the momentum distribution function of the impurity in a zero-temperature polaron state [59].On the other hand, the static finite temperature ejection Green's function was computed in [55].This can be then generalized to time-dependent Green's function as we show in the Supplementary Materials [60].This approach culminates in the following expressions for the two Green's functions: Here det σ (1 + K) denotes a Fredholm determinant of the kernel K(q, q ′ ) with a measure given by the Fermi-Dirac distribution σ(k) [61].The kernels appearing above take the following universal form V (q, q ′ ) = 1 π e + (q)e − (q ′ ) − e + (q ′ )e − (q) q − q ′ , ( W ± (q, q ′ ) = ± 1 π e ± (q)e ± (q ′ ), (23) where e + (q) = e(q)e − (q).Functions e − (q) and e(q) have different expressions for the injection and ejection.Namely, for the injection we have e in − (q) = e itq 2 /4−ixq/2 (24) and additionally For a real argument we specify how we deform a contour to pass over or under a real line by ±i0 shifts.
For the ejection the formulas read The expressions for the two Green's functions, however sharing similar structures, in the end involve different functions and a priori any simple relation between G in (x, t, Λ) and G ej (x, t, Λ) is unexpected.It is clear from the formula for G in (x, t, Λ) that it depends on the spin rapidity Λ in a non-trivial way which is a sufficient condition for breaking of the equivalence between GCE and MCE.
Refined KMS relation and the Riemann-Hilbert problem: We now sketch a derivation of the refined detailed balance relation [62].The idea of the proof relies on interpreting the two Green's function as solutions to the Riemann-Hilbert problem (RHP).The RHP in its simplest formulation is a problem of determining a function which is analytic everywhere but along the real line and which asymptotically approaches 1.The non-analytic behaviour along the real line is characterized by the jump condition: χ + (x) = χ − (x)G(x), where χ ± (x) = lim ϵ→0 χ(x ± iϵ).For our application we need a correspondent generalization to matrix-valued functions [63,64].
Representing the Green's functions through a solution to the RHP problem corresponds to identifying the jump matrices G in,ej .We then perform a series of manipulations of jump matrices which demonstrates that they can be made equal after an appropriate change of the coordinate and time.This allows us to conclude that with a space-time-independent function Therefore W (β, Λ) can be found by considering t = 0 and taking x → ∞ limit, where the x-asymptote of the Fredholm determinants can be evaluated with the effective form-factors approach [65][66][67].The result is W (β, Λ) = exp(βF(Λ)), thus finishing the proof of the refined KMS relation.[68] The Green's function can be evaluated numerically [69].This requires a numerical evaluation of the Fredholm determinant and can be done by the quadrature method [70].In Fig. 2 we show the KMS relation between the injection and ejection Green's function.This constitutes a numerical proof of this relation.
In the limit of vanishing impurity-gas interactions, the Green's functions become those of a free particle with its dispersion controlled by Λ [71] and additionally for the ejection weighted by the thermal gas distribution In the momentum space the correlators are then proportional to δ(k ± Λ) -in the non-interacting limit the spin rapidity Λ can be identified with momentum k.This shows that the Λ degree of freedom emerges from the interacting nature of the system.It is an open question how to observe this effect for non-integrable models, for example using the variational method [47,48].

Conclusions:
In this work we studied exact thermal Green's functions of an interacting impurity model.We have shown that due to the interactions the theory acquires a new thermodynamic parameter, the spin rapidity Λ of the nested Bethe Ansatz.As a result, the ejection Green's function involves a sum over an ensemble of systems weighted with Λ-dependent free energies.The injection Green's function can also be resolved into contributions from different Λ's.This appears as a result of quantum-mechanical averaging rather then thermal av-eraging.Despite that, the two Λ-resolved Green's functions obey the KMS relation.The existence of such a relation effectively promotes Λ to a thermodynamic parameter describing the equilibrium state of the system.On the other hand the KMS relation is a statement about analytic continuation of correlation functions.In order to prove it, it requires non-perturbative control over the correlators, which is provided within our mathematical framework.We then employed the Riemann-Hilbert problem to prove the KMS relation and we further verified it numerically.
The methods developed here apply beyond the thermal equilibrium.For example, the expressions for the Λ-resolved correlation functions are valid for any distribution σ(k), not necessarily thermal.An interesting question in this direction is about the existence of KMS relation beyond the thermal equilibrium.Closely related statements of detailed balance were found in similar circumstances for the 1d interacting Bose gas [72,73].

S1 Yang-Gaudin model and the impurity problem
In this section we review the Bethe ansatz solution of the Yang-Gaudin model.We focus on the sectors of the theory with zero or one spin up particle in a thermodynamically large sea of spin down particles.These sectors are relevant for the impurity Green's function.In our presentation we follow [1,2].We also discuss the excitation spectrum of the theory, specifically we show a degeneracy excitations over finite temperature state.This degeneracy can be traced back as a microscopic foundation for the refined detailed balance.
The Hamiltonian of the Yang-Gaudin model is The Hamiltonian commutes with the number operators of spin up and down particles and with the total momentum operator, and particle number of each type is a conserved quantity.The eigenstates of the system can be then characterized by total number of particles N = N ↑ + N ↓ and M = N ↓ number of spin down particles.We call tuple (N, M ) a sector.
In each sector (N, M ) the Hamiltonian can be written in the first quantized form In the following we set ℏ = 1 and m = 1.

The Bethe ansatz solution
The eigenstates in the sector (N, M ) are characterized by two sets of rapidities (quasi-momenta), {k} = k 1 , . . ., k N and {λ} = λ 1 , . . ., λ M which are real numbers and obey Pauli principle in each set separately.The basis of the Fock space of the model is constructed in the standard way by defining a vacuum state |0⟩ such that The eigenstates take then the following form with wave function Ψ {s} N,M ({k}, {λ}|{x}) and set of positions {x} = x 1 , . . ., x N .For the wave function in sector (N, M ) to be non-zero the set of spins {σ} = σ 1 , . . ., σ N must be such that exactly M elements among them are spin downs.Finally, the rapidities are solutions to the nested Bethe equations The momentum and energy of an eigenstate are (S1.9) The wave function Ψ {s} N,M ({k}, {λ}|{x}) takes the typical structure of Bethe ansatz solvable models, it consists of a superposition of plane waves with rapidity-dependent amplitudes.In the following we focus on sectors (N, 0) and (N + 1, 1) which are relevant for the impurity problem.
In sector (N, 0) the system is a free Fermi gas of N spin up particles.The second set of rapidities is empty, {λ} = 0 and the Bethe equations reduce to the standard quantization conditions with the wave function given by a Slater determinant of an N × N matrix, (S1.11) In sector (N + 1, 1) the wave-function can also be represented by a determinant.To achieve this one changes the coordinate system to the impurity's (the spin down particle) rest frame, a transformation known as the Lee-Low-Pines transformation [3], applied to the McGuire model in [4,5].Here we follow [2].The resulting wave-function is . . . . . . . . .
with the normalization factor and with and is related to the original wave function through two relations where by ↓ j we abbreviate a set of spins {σ} with a single spin down at position j in the sea of spin up particles.The rapidities {k j } and the single spin rapidity Λ obey The second equation can be understood as the quantization condition of the total momentum.Indeed taking the product over all l of the first equation and then using the second equation we find

Bethe equations and excited states
The Bethe equations in the logarithmic form are where quantum numbers n j are integers and obey the Pauli principle.The phase shift is The spin rapidity Λ can be fixed by specifying values of other integrals of motions in this model.Traditionally, we require that the total momentum given by is fixed, i.e.P ({k j }, Λ) = Q.The Λ dependence in (S1.23) is implicit through k j as solutions to the Bethe equations.Therefore, for given Q and the set of integers one can resolve condition (S1.23) and thus solve the Bethe equations.The Bethe equations for rapidities {k j } can be rewritten in terms of a single function obeying a transcendental equation, where Function f depends only on the parameters of the system α and L but not on the rapidities.The expression for k j gives a bound on the total momentum of the state.Denoting and using that arctan(x) is a bounded function we find The maximal and minimal values of rapidity correspond to Λ = ±∞.In that case the rapidities can be computed analytically.The solution is For the other cases Λ is finite and equations have to be solved numerically.Degeneration of states: Consider a set of quantum numbers {n j } and corresponding Λ such that the Bethe equations are fulfilled and denote rapidities by k j .We can now build another state related to it by a parity operation.Choose nj = −n j + 1 and Λ = −Λ.The corresponding rapidities are kj = −k j and obey the Bethe equations.The two states have the same energy and opposite momenta.Derivation where in the last step we used the symmetry property of f (x) and that arctan(x) is an odd function.This can be also derived directly from the original formulation.
Excitations far from the ground state: In the thermodynamically large system, a state can be described by the distribution of rapidities σ(k) and value of the rapidity Λ.With respect to such state we consider excitations.This take form of modifying some subset of quantum numbers n j and/or adding/removing the impurity.
Consider a state without the impurity described by quantum numbers {n j } N j=1 and a state with the impurity described by quantum numbers {n ′ j } N +1 j=1 and Λ.We assume that the two sets share most of the quantum numbers with m quantum numbers different such that m/L → 0 in the thermodynamic limit.Then the two states are described by the same distribution n(k).Consider subleading in the system size corrections to the values of the rapidities.For j = 1, . . ., N , Unless n ′ j − n j ∼ L, the difference between the rapidities is of the order 1/L.For the quantum numbers that were not modified we then have where we replaced the quantum number 2πn ′ j /L by the rapidity k j as the error is of the order 1/L and therefore has subleading contribution to k j − k ′ j .Consider now difference in momenta between the two states Similar computations for the energy difference give These excitation can be summarized by We can also consider an opposite type of excitation which involves annihilating the impurity, The momentum and energy carried by this excitation is We now show that there exist a symmetry between the two types of excitations.This symmetry underlies the refined detailed balance.Choosing Λ ′ = −Λ, p ′ j = h j and h ′ j = p j the two excited states carry exactly opposite momenta and energies with respect to their background states if n(k) is a symmetric function of k.Reformulating this, for every excited state in which we create an impurity there exist an excited state in which we annihilate the impurity such that the two excitations have exactly opposite energies, momenta and Λ's.

S2 Impurity Green's functions at finite temperatures
In this Section we derive the expression for the ejection Green's function at finite temperature and show that it breaks the Eigenstate Thermalization Hypothesis.We also derive explicit expression for the finite temperature ejection and injection Green's function.For the ejection Green's function we generalize the static finite temperature result of [6].For the injection the finite tempearture dynamic Green's function was computed in [2].
We consider now evaluation of the ejection Green's function in a finite system.In thermal equilibrium it is given by with thermal density matrix Denoting eigenstates in finite system by |{k}, Λ⟩ we write with the partition function Z having analogous representation The building block of the finite temperature correlation function is an expectation value in a single eigenstate, which we denote The normalized form-factors ⟨{q}|Ψ(0)|{k}, Λ⟩ were computed in [1] and take the following form Here the derivative of k j over Λ is formal and explicitly equals to Notice that in the sum in the denominator, one can ignore 1/L corrections and obtain The result of [6], here generalized to dynamic correlation function, is that in the thermodynamic limit the expectation value ⟨{k}, Λ|Ψ † (x, t)Ψ(0, 0)|{k}, Λ⟩ depends on the underlying distribution of {k j } and on the spin rapidity Λ ⟨σ(k), Λ|Ψ † (x, t)Ψ(0, 0)|σ(k), Λ⟩ = lim td ⟨{k}, Λ|Ψ † (x, t)Ψ(0, 0)|{k}, Λ⟩. (S2.9) This implies that when performing the thermal averaging one can use the saddle-point argument to localize the sum over {k j } at configuration that maximizes the free energy.The free energy of the system takes the following form [6] where F th (σ(k)) corresponds to the free fermions energy, thence not depending on the rapidity Λ.However, the subleading contribution F(σ(k), Λ) does depend on it.Therefore, the thermal sum in the numerator of the ejection Green's function evaluates in the thermodynamic limit to with constant factors coming from the saddle-point evaluation.The same factors will appear in the evaluation of the partition function and therefore they cancel in the final expression for the Green's function.This takes then the following form , (S2.12) Rewriting the sum over Λ through integrals we obtain where we defined with a(Λ) defined in (S2.8).Finally, because the saddle point configuration comes from the extremum of the free energy of the non-interacting Fermi gas, the distribution σ(k) takes the Fermi-Dirac form The final formula for G ej (x, t, Λ) can be obtained from generalization of the static Green's function.In [6] it was derived The kernels are V (q, q ′ ) = e + (q)e − (q ′ ) − e − (q)e + (q ′ ) q − q ′ , Ŵ− (q, q ′ ) = − 1 π e − (q)e − (q ′ ), (S2.17 with functions e ± (q) given by e + (q) = 1 π e iqx/2+iδ(q) , e − (q) = e −iqx/2 sin δ(q).(S2.18) We also note that G ej (x, −Λ) is a complex function such that G ej (x, −Λ) = G * ej (x, Λ).Therefore the resulting one-body function G ej (x) is a real function.
To include the time dependence we include the energy contribution to functions e ± (q), Notice that for x > 0 and t = 0 the integral in e + (q) vanishes.In this way we obtain the formula for the ejection Green's function G ej (x, t; σ, Λ) reported in the main text.Injection Green's function: Similar analysis can be performed for the injection Green's function.We find there the ETH works directly at the level of the full correlator reducing the thermal averaging to a single expectation value.However, as first investigated in [2], it is possible to truncate the internal sum to states with fixed value of Λ thus in practice realizing the projection operator P Λ .The formula for the Λ-resolved correlation function was derived in [2].The expressions are Here det [−Q,Q] (1 + K) denotes a Fredholm determinant of the kernel K(q, q ′ ) with a measure given by the Fermi-Dirac distribution σ(k) with Q the Fermi rapidity.The kernels appearing above are where e + (q) = e(q)e − (q) and functions ϵ − (q) and ϵ(q) given by e − (q) = e itq 2 /4−ixq/2 , (S2.22) and additionally H = 1 − F (κ + )/α + F (κ − )/α.Here κ ± = (Λ ± i)/α and This function can be expressed in terms of error function.For a real argument we specify how we deform a contour to pass over or under a real line by ±i0 shifts.

S3 Riemann-Hilbert problem and the KMS relation
In this section, we reformulate the computation of Fredholm determinants in terms of a Riemann-Hilbert problem (RHP).We perform a transformation on the RHP for the injection and ejection cases to derive the refined detailed balance and mostly follow Ref. [7].Let us start by recalling the statement of the matrix Riemann-Hilbert problem.
Let Σ be an oriented contour in a complex plane and let G be a matrix valued function defined on Σ.The task is to find a matrix valued function χ(z) which is holomorphic in the complement of Σ.Additionally, let us denote by χ ± (z) the value of χ(q) when it approaches Σ from one of its two sides.We require χ(z) to fulfill two conditions • the limiting values are related through the jump matrix G, i.e. χ + = χ − (1 + 2πiG), • the function χ(z) approaches the identity matrix when |z| → ∞ in the complement of the contour Σ.
In our case, the Σ contour will be simply the real line.
Our strategy to prove the refined the detailed balance is the following.We will reformulate the computation of the two Green's as Riemann-Hilbert problems.This amounts to specifying the contour Σ (which in both cases is real line) and the jump matrix.We will then perform some transformations of the RHP to arrive at the same jump matrix.From this, assuming the uniqueness of the solution to the RHP, we derive the following relation between the two Green's functions Importantly, from the RHP we can infer that function the ratio of the determinants is independent of x and t.Therefore it can be evaluated by taking a suitable limit of t = 0 and x → ∞.Computing the ratio amounts then to extracting the asymptotic behavior of the two determinants.This can be achieved with the effective form-factors approach.The result is thus proving the refined detailed balance relation.
We start with the injection case and rewrite the kernel V in in the vector notations where we have introduced "bra" and "ket" vectors | = (e + (q), e − (q)).(S3.4) In these notations, we can immediately recognize the kernel as a generalized sine-kernel [8].Notice also that in the original kernels for both injection and ejection Eqs.(S2.17),(S2.21)we can replace factors responsible for time dynamics e itq 2 /2 → e itε(q) , with the shifted energy ε(q) = q 2 /2 − µ, which results in the rescalings of the components e ± (q) → e ± (q)e ∓iµt/2 and the total correlation functions Because of the special (integrable) form of the operator V its resolvent R defined via also have an integrable form Here, similar to (S3.4), we have put The components satisfy linear integral equations where we remind that the action of the operator is given by the convolution with its kernel.This system of linear integral equations can be reformulated as a Riemann-Hilbert problem.To do so we formally introduce 2 × 2 matrices and verify that |F in (q)⟩ = χ in (q)|E in (q)⟩, ⟨F in (q)| = ⟨E in (q)| χin (q), χ in (q) χin (q) = 1. (S3.11) The matrix function χ is analytic everywhere except the real line, where it experiences a jump with G in (q) = |E in (q)⟩⟨E in (q)| = e − (q)e + (q) −(e + (q)) 2 (e − (q)) 2 −e − (q)e + (q), (S3.12) Therefore χ in (q) is a solution to the Riemann-Hilbert problem.Formulation of χ(q) as a solution to the RHP allows us to infer about its asymptotics.Asymptotically, the solution to the RHP takes the following form where using notations (S3.4) and (S3.8) we explicitly write These expressions are also usually referred to as potentials.The symmetry of the kernel V (q, q ′ ) = V (q ′ , q) is reflected in the relation B +− = B −+ .One can now easily express Green's function via potentials [7].In particular, using relation (S3.9), for the injection Green's function we obtain where we used that W (q, q ′ ) = e + (q)e + (q ′ ) is rank one and introduced the τ -function as the determinant of the operator Vin , Now let us consider different RH problem for the matrix ϕ in The corresponding jump Ḡin in ϕ in (q) across the real line is Ḡin = 1 π σ(q)s − σ(q)e −i(tε(q)t−qx) s − s + (1 − σ(q))e i(tε(q)−qx) −σ(q)s + .(S3.18) To derive this expression we have extensively used that s − − s + = 2is − s + and φ + (q) − φ − (q) 2πi = e −itε(q)+iqx π s − s + , (e + − φ ± e − )(q) = σ(q)s ± (q) e itε(q)/2−iqx/2 π .
(S3. 19) The jump matrix is valid for any distribution σ(q).For the thermal distribution it simplifies to Ḡin = 1 π(e βε(q) + 1) For q → ∞, ϕ(q) has the same expansion as χ(q) but with modified potentials b ij and c ij .Moreover, from Eq. (S3.17) we conclude the relation between the old and new potentials For the Green's function we then have This expresses the injection Green's function through the potential b ++ and the τ -function.The potential b ++ appears from a solution to the RHP with specific jump matrix Ḡin .We will find now a similar representation for the ejection Green's function.Now let us consider the ejection.We denote the corresponding RHP matrix as χ ej and the potentials with the subscripts ej.Following the same procedure as for the injection we obtain We perform again the conjugation The jump matrix ϕ ej for reads Gej = 1 π −s + σ(q) s − s + σ(q)e itε(q)−iqx (1 − σ(q))e itε(q)−iqx σ(q)s − .(S3.26) This jump matrix is structurally similar to Ḡin from eq. (S3.18) but with s + and s − replaced.This can be fixed by performing a further conjugation, ϕ ej = σ 1 ψ ej σ 1 . 1 The resulting jump matrix is Finally specifying σ(q) to be a thermal distribution we conclude the following jump matrix for Ψ Comparing this expression with the jump matrix (S3.20) we observe that which leads to the following identity By comparing the asymptotic expansion of both sides of this equality we conclude that For the ratio of the Green's function we find We will now prove that the ratio of the determinants is independent of x and t.To this end consider derivatives of the kernels V in and V ej .For the injection, taking into account specific dependence of the kernel, we find 1 Matrix σ 1 appears here on both sides to ensure the same asymptotic behavior of φ ej and ψ ej .

S4 Fredholm determinant representation of the Green's functions
In this Section we recall the definition of the Fredholm determinant and provide details on the numerical evaluation of it that we use in our work.
The kernels are of a Sine-type and are of the form (f (x) − f (x ′ ))/(x − x ′ ).To avoid the problematic x = x ′ point we discretize x and x ′ on grids shifted with respect to each other such that x j ̸ = x ′ k for any j and k.To establish the convergence of the results we evaluate the determinants on grids with different number of points.To quantify the convergence we define conv = G (N2) (x, t, Λ) − G (N1) (x, t, Λ) G (N1) (x, t, Λ) , (S4.5) where G (N ) (x, t, Λ) is the Green's function (either injection or ejection) computed on a grid consisting of N points.
The results for the grids of N = 100, 200, 500, 1000 points and for both correlators are shown in fig.S1.We have verified that the fast convergence rate holds also for other values of the parameters of the system.
FIG. 1. a) The standard KMS relation infers the detailed balance condition Gin(k, ω) = e −βω Gej(k, ω), between the Fourier transforms of injection and ejection Green's function, which implies that probabilities for fluctuations creating or destroying impurity (denoted by a red circle) are the same in the thermal equilibrium.b) The thermal state with an impurity is an ensemble of different states enumerated by Λn. c) The refined KMS relation guarantees the stability of each sub-ensemble under the fluctuations creating or destroying the impurity described by Gin,ej(k, ω, Λ).
Kubo-Martin-Schwinger relation for an interacting mobile impurity Contents S1.Yang-Gaudin model and the impurity problem 1 S2.Impurity Green's functions at finite temperatures 5 S3.Riemann-Hilbert problem and the KMS relation 7 S4.Fredholm determinant representation of the Green's functions 12 References 13 FIG.S1: Plots of the convergence measure (S4.5).We plot real and imaginary parts of the two Green's function for the parameters of the system shown above the plots.The results show that increasing the size of the grid from N = 500 to N = 1000 points changes the values of the two functions by at most 0.1%.