# Non-Hermitian quantum impurity systems in and out of equilibrium: Noninteracting case

Takato Yoshimura, Kemal Bidzhiev, Hubert Saleur

<figure><img src="image_1.png" /></figure>

► To cite this version:

Takato Yoshimura, Kemal Bidzhiev, Hubert Saleur. Non-Hermitian quantum impurity systems in and out of equilibrium: Noninteracting case. Physical Review B, 2020, 102 (12), pp.125124. 10.1103/Phys-RevB.102.125124. hal-03017010

HAL Id: hal-03017010

https://hal.science/hal-03017010

Submitted on 12 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Non-Hermitian quantum impurity systems in and out of equilibrium: Noninteracting case

Takato Yoshimura ◦,¹,² Kemal Bidzhiev,³ and Hubert Saleur⁴

¹Institut de Physique Théorique Philippe Meyer, École Normale Supérieure, PSL University, Sorbonne Universités, CNRS, 75005 Paris, France

²Department of Mathematics, King's College London, Strand, London WC2R 2LS, United Kingdom

³Université Paris-Saclay, CNRS, LPTMS, 91405 Orsay, France

⁴Institut de Physique Théorique, Université Paris Saclay, CEA, CNRS, F-91191 Gif-sur-Yvette, France

<figure><img src="image_3.png" /></figure>

(Received 12 March 2020; revised 27 August 2020; accepted 31 August 2020; published 15 September 2020)

We provide systematic analysis on a non-Hermitian *PT*-symmetric quantum impurity system both in and out of equilibrium, based on exact computations. In order to understand the interplay between non-Hermiticity and Kondo physics, we focus on a prototypical noninteracting impurity system, the resonant level model, with complex coupling constants. Explicitly constructing biorthogonal basis, we study its thermodynamic properties as well as the Loschmidt echo starting from the initially disconnected two free fermion chains. Remarkably, we observe the universal crossover physics in the Loschmidt echo, both in the *PT* broken and unbroken regimes. We also find that the ground state quantities we compute in the *PT* broken regime can be obtained by analytic continuation. It turns out that Kondo screening ceases to exist in the *PT* broken regime, which was also previously predicted in the non-Hermitian Kondo model. All the analytical results are corroborated against biorthogonal free fermion numerics.

DOI: 10.1103/PhysRevB.102.125124

I. INTRODUCTION

Quantum impurity systems, such as the Kondo model,
serve as representative examples where nonperturbative
quantum many-body effects give rise to unconventional low-
temperature behaviors that are in stark contrast with those of
Fermi liquids [1]. One of the cornerstones of such impurity
models is the screening of the impurity spin at sufficiently low
temperature (lower than the typical temperature scale called
Kondo temperature $T_K$), which is characterized by the spin
singlet made of the impurity spin and low-lying excitations
[2]. The impurity spin is therefore deactivated by healing with
the lead baths, thereby the typical low temperature behav-
ior of the spin susceptibility $\chi(T) \sim 1/T$ is replaced by the
approach towards a constant value $\chi(T) \to \chi_0$ as $T \to 0$.
A great deal of effort has been made to understand the na-
ture of Kondo screening in the past decades making use of
renormalization group [3], Bethe ansatz [4,5], and numerical
methods such as the numerical renormalization group (NRG)
[6] and density-matrix renormalization group (DMRG) [7].
Recently, experimental investigations of the Kondo physics
have been put forward using alkaline-earth atoms, opening up
an unprecedented avenue to detect the Kondo cloud, which
has remained experimentally elusive [8,9].

There recently has been an interest in exploring what
happens to Kondo physics in the non-Hermitian case [10,11].
This is partly motivated by setups such as the one in Ref. [10],
where a non-Hermitian Kondo model is obtained as an
effective model that characterizes ultracold mobile and
immobile atoms that undergo inelastic scatterings (resulting
in two-body losses). Another motivation comes from the old

observation [12–14] that making the couplings complex in
the Kosterlitz-Thouless (KT) renormalization group (RG)
flow can profoundly change the low-energy physics and lead
sometimes to “circular” behavior. This RG flow is relevant
for the sine-Gordon model with complex coupling constant, a
model that found early applications in 2d statistical mechanics
(e.g., polymers), and was also studied more recently in
Ref. [15] in the context of dissipative quantum mechanics.
The KT RG flow is of course also relevant to the (anisotropic)
Kondo model, and the circular behavior was explored in this
context in Ref. [11], where it was concluded that the Kondo
effect might disappear in the non-Hermitian case.

Non-Hermitian “Hamiltonians” of physical relevance are
usually *PT* symmetric and can exhibit both *PT*-broken and
unbroken regimes. The spectrum of a system with unbroken
*PT* symmetry is actually real, and such systems can essen-
tially be treated like Hermitian ones. The situation is different
in the *PT*-broken regime where generically eigenstates with
complex eigenvalues are present. The new low-energy behav-
ior in non-Hermitian versions of the sine-Gordon model or the
Kondo model is expected to occur in such broken symmetry
regimes.

Since the physics of interest in this context involves typ-
ically many-body effects, available methods to understand
what is going on are few, apart from perturbative RG cal-
culations. Attempts using the Bethe ansatz [12–14] have in
particular been much less successful than in the usual Her-
mitian case, due in part to the appearance of bound states
of complex or even purely imaginary energy (“monstrons”
[13]). Reference [10] does partly explore the issue in the (non-
Hermitian but not *PT*-symmetric) Kondo case, concluding
that the Kondo effect might disappear when the strong non-Hermiticity is imposed; the calculations in this reference are however only a starting point and leave many questions open.

In the Hermitian case, a great deal of the physics of quantum impurity problems both in and out of equilibrium can be learned without the Bethe ansatz, by studying the resonant level model (RLM) [16], which occurs as a simpler (Toulouse) limit of interacting impurity systems such as the anisotropic Kondo model and the interacting resonant level model. It is thus natural, in order to clarify what happens to quantum-impurity problems in the non-Hermitian case, to start by studying a PT-symmetric version of the RLM (PTRLM) model. The fact that the model is essentially free will allow us to carefully analyze the issues of biorthogonal bases, bound states, and PT-symmetry breaking, paving the way for further studies of the more delicate, interacting case. We will see that the conclusions of earlier studies must be taken with a grain of salt and that there is probably a lot of physics left to be understood in this problem. We also point out that a fermionic chain with localized losses, which is somewhat similar to the model we consider in this paper, have also been investigated with the aim to better understand the interplay between coherence and dissipation in quantum systems [17,18]. We finally mention that our model can be in principle realized using alkaline-earth atoms with a tunable PT symmetric optical lattice [8,19,20].

The paper is organized as follows. The spectrum of the model is obtained in Sec. II, where we identify the different regimes of the model, and the PT-symmetry breaking. The dot occupancy is studied in Sec. III, and the boundary free-energy in Sec. IV. These two equilibrium quantities are known to exhibit characteristic features of Kondo physics in the Hermitian case. While we observe similar features in the PT-symmetric regime, very different properties occur when that symmetry is broken. Section V deals with nonequilibrium physics, where we look for another angle to understand the physics of the broken PT-symmetric phase by studying the Loschmidt echo. Consideration of analytic continuation are discussed in Sec. VI, where we suggest another, more direct way to handle the problem that could be generalized to the interacting case. Further observations and conclusions are discussed in the conclusion.

## II. SPECTRUM OF THE PT-SYMMETRIC RESONANT LEVEL MODEL

### A. The model

The Hamiltonian of the PT-symmetric RLM model (PTRLM) we study takes the following form (see also Fig. 1):

$$ 
\begin{aligned}
H &= H_A + H_B + H_d \\
H_A &= -t \sum_{x=-N/2}^{-2} (c_x^\dagger c_{x+1} + \text{H.c.}) \\
H_B &= -t \sum_{x=1}^{N/2-1} (c_x^\dagger c_{x+1} + \text{H.c.}) \\
H_d &= -\gamma(c_{-1}^\dagger c_0 + \text{H.c.}) - \gamma^*(c_1^\dagger c_0 + \text{H.c.}) \\
    &= -Je^{i\varphi}(c_{-1}^\dagger c_0 + \text{H.c.}) - Je^{-i\varphi}(c_1^\dagger c_0 + \text{H.c.}),
\end{aligned} 
\quad (1) $$

<figure><img src="image_4.png" /><figcaption>FIG. 1. A schematics of the PT resonant level model (1). The impurity couples to the leads with hopping amplitudes $\gamma$ and $\gamma^*$ (or tunneling strength).</figcaption></figure>

where $H_{A/B}$ describes the kinetic energy of free spinless fermions in the left and right wires, $H_d$ defines tunneling from the wires to the dot impurity with a hopping amplitude $\gamma = \gamma_1 + i\gamma_2 = Je^{i\varphi}$. In what follows, we set tunneling strength between wires $t=1$ for brevity. That this model is PT symmetric follows from the following definitions [11,21] of how P and T act on operators

$$ \hat{P}c_x\hat{P}^{-1} = c_{-x}, \quad \hat{T}c_x\hat{T}^{-1} = c_x, \quad \hat{T}i\hat{T}^{-1} = -i. \quad (2) $$

$$ (\hat{P}\hat{T})H(\hat{P}\hat{T})^{-1} = H \quad (3) $$

### B. The eigenstates

In order to study the spectrum of PT-RLM (1), let us start with solving the eigenvalue problems of the model by constructing the biorthogonal basis [22]. We write a generic one-particle right eigenstate with the eigenvalue E as

$$ |R\rangle = \sum_x R_E(x)|x\rangle, \quad (4) $$

where $|x\rangle = c_x^\dagger|0\rangle$. First, for $1 < |x| < N/2$, $R_E(x)$ solves

$$ R_E(x-1) + R_E(x+1) = -E_R(x). \quad (5) $$

Let us take the following ansatz [23] for $x \neq 0$,

$$ R(k; x) = \begin{cases} \mathcal{N}_+(k)(e^{ikx} + \psi_+(k)e^{-ikx}) & x > 0 \\ \mathcal{N}_-(k)(e^{ikx} + \psi_-(k)e^{-ikx}) & x < 0 \end{cases}, \quad (6) $$

where $\psi_+ = e^{2ik\delta(k)}$, $\psi_- = e^{-2ik\delta(k)}$, and $\mathcal{N}_\pm(k)$ are normalizations. This solves (5) where the quasimomentum $k$ is related to $E_k$ by the dispersion relation $E_k = -2\cos k$. The phase shift $\delta(k)$ is determined by

$$ R(k; 2) + \gamma^* R(k; 0) = -E_k R(k; 1) \quad (7) $$

$$ \gamma R(k; -1) + \gamma^* R(k; 1) = -E_k R(k; 0) \quad (8) $$

$$ R(k; -2) + \gamma R(k; 0) = -E_k R(k; -1). \quad (9) $$

Solving these self-consistently, we find that, in fact, there are two types of solutions that are labeled by $\sigma = \pm 1$. The phase shifts corresponding to these solutions accordingly depend on $\sigma$ and satisfy

$$ e^{2ik\delta_\sigma(k)} = \frac{1 - \Gamma_\sigma e^{2ik}}{\Gamma_\sigma - e^{2ik}}, \quad (10) $$

where $\Gamma_\sigma = (1+\sigma)\tilde{\Gamma} - 1$ with $\tilde{\Gamma} = (\gamma^2 + (\gamma^*)^2)/2 = J^2 \cos 2\varphi$. Note that $\tilde{\Gamma} = J^2$ in the Hermitian case. Defining $\mathcal{N}_{R,+} = 2\mathcal{N}_+ e^{ik\delta_+(k)}$ and $\mathcal{N}_{R,-} = 2i\mathcal{N}_+ e^{ik\delta_-(k)}$, it is then a simple matter to observe that the wave function $R_\sigma(k; x)$ takes
<figure><img src="image_5.png" /><figcaption>FIG. 2. The phase diagram of the PTRLM for different coordinates $\gamma = \gamma_1 + i\gamma_2 = Je^{i\phi}$. Phases I and II are the PT-unbroken phases that are distinguished by the presence (phase II) and the absence (phase I) of bound states. Phase III corresponds to PT-broken phase with complex-valued single particle energies.</figcaption></figure>

the following form

$$R_{\sigma}(k; x) = \begin{cases} \mathcal{N}_{R,\sigma}(k) \cos[k(x - \delta_{\sigma}(k))] & x > 0 \\ \mathcal{M}_{\sigma}\mathcal{N}_{R,\sigma}(k) \cos[k(x + \delta_{\sigma}(k))] & x < 0 \end{cases}, \quad (11)$$

where $M_+ = \gamma^2/|\gamma|^2$ and $M_- = -(\gamma^*)^2/|\gamma|^2$. Note that $R_\sigma(k; 0)$ is given by (8). Finally, the quantization condition is obtained from the equation for $x = N/2$

$$\cos k(N/2 + 1 - \delta_{\sigma}(k)) = 0, \quad (12)$$

which, together with (10), gives

$$e^{ikN} = \frac{\Gamma_{\sigma} - e^{-2ik}}{\Gamma_{\sigma} - e^{2ik}}. \quad (13)$$

The construction of the left eigenstate $|L\rangle = \sum_x L_E(x)|x\rangle$ can also be carried out in the analogous way. For real $k$, it reads

$$L_{\sigma}(k; x) = \begin{cases} \mathcal{N}_{L,\sigma}(k) \cos[k(x - \delta_{\sigma}(k))] & x > 0 \\ \mathcal{M}_{\sigma}^{*}\mathcal{N}_{L,\sigma}(k) \cos[k(x + \delta_{\sigma}(k))] & x < 0 \end{cases} \quad (14)$$

(Fig. 1).

### C. Bound states

So far we have focused on the case where $k$ is real, but we can also find solutions of equation (13) with complex $k$, corresponding to bound (localized) states. In general it is hard to determine the precise form of their eigenfunctions except when thermodynamic limit ($N \to \infty$) is taken, which is the case we shall consider below. It turns out that there are three distinct regimes which we present in Fig. 2. In what follows, without loss of generality, we shall focus on $\gamma_1, \gamma_2 > 0$.

(i) Phase I ($\gamma_2 < \gamma_1 < \sqrt{1 + \gamma_2^2}$): PT-unbroken regime without bound states.

No complex $k$ is allowed in this phase, hence there is no bound state and the excitation is solely characterized by the real $k$. The right and left eigenfunctions are given by (11) and (14). The normalization condition for the left and right wave functions reads

$$\begin{aligned} \mathcal{N}_{L,+}(k)^* \mathcal{N}_{R,+}(k) &= \frac{(\gamma^*)^2}{2\tilde{\Gamma}} \left[ \frac{N}{4} + \frac{1}{2} + \frac{2\tilde{\Gamma}(1-\tilde{\Gamma})}{\Delta_+(k)} \right]^{-1} \\ \mathcal{N}_{L,-}(k)^* \mathcal{N}_{R,-}(k) &= \frac{\gamma^2}{2\tilde{\Gamma}} \left[ \frac{N}{4} + \frac{1}{2} \right]^{-1}, \end{aligned} \quad (15)$$

where $\Delta_+(k) = 1 - 2\tilde{\Gamma}_+ \cos 2k + \Gamma_+^2 = 4(\tilde{\Gamma}^2 + (1 - 2\tilde{\Gamma}) \cos^2 k)$. To determine $\mathcal{N}_{R,\pm}$ and $\mathcal{N}_{L,\pm}$, it is useful to remind ourselves of a simple fact $\hat{T}|R\rangle = |L\rangle$, which is a consequence of $\hat{T}H\hat{T}^{-1} = H^\dagger$. This relation then implies that $\mathcal{N}_{R,\pm} = \mathcal{N}_{L,\pm}^*$ with

$$\begin{aligned} \mathcal{N}_{R,+}(k) &= \frac{\gamma^*}{\sqrt{2\tilde{\Gamma}}} \left[ \frac{N}{4} + \frac{1}{2} + \frac{2\tilde{\Gamma}(1-\tilde{\Gamma})}{\Delta_+(k)} \right]^{-1/2} \\ \mathcal{N}_{R,-}(k) &= \frac{\gamma}{\sqrt{2\tilde{\Gamma}}} \left[ \frac{N}{4} + \frac{1}{2} \right]^{-1/2}. \end{aligned} \quad (16)$$

(ii) Phase II ($\gamma_1 > \sqrt{1 + \gamma_2^2}$): PT-unbroken regime with bound states.

In this phase, on top of the real eigenvalues as in phase I, bound states being localized around the dot can also be formed. The corresponding right and left eigenfunctions are

$$R_{b,\pm}^{\text{re}}(x) = \begin{cases} \mathcal{N}_{R,b}^{\text{re}}(\pm 1)^x e^{-x/\xi} & x > 0 \\ \mathcal{N}_{R,b}^{\text{re}} \frac{1}{\gamma^*} & x = 0, \\ \frac{\gamma^2}{|\gamma|^2} \mathcal{N}_{R,b}^{\text{re}}(\pm 1)^x e^{x/\xi} & x < 0 \end{cases} \quad (17)$$

$$L_{b,\pm}^{\text{re}}(x) = \begin{cases} \mathcal{N}_{L,b}^{\text{re}}(\mp 1)^x e^{-x/\xi} & x > 0 \\ \mathcal{N}_{L,b}^{\text{re}} \frac{1}{\gamma} & x = 0, \\ \frac{(\gamma^*)^2}{|\gamma|^2} \mathcal{N}_{L,b}^{\text{re}}(\mp 1)^x e^{x/\xi} & x < 0 \end{cases} \quad (18)$$

where $\xi = 2/\log|\Gamma_+| = 2/\log|\gamma^2 + (\gamma^*)^2 - 1|$ is the localization length of the bound states. Associated eigenvalues and the normalization condition for $\mathcal{N}_{R,b}^{\text{re}}$ and $\mathcal{N}_{L,b}^{\text{re}}$ are provided by

$$E_{b,\pm}^{\text{re}} = \pm 2 \frac{\gamma_1^2 - \gamma_2^2}{\sqrt{2\gamma_1^2 - 2\gamma_2^2 - 1}} = \pm \frac{\gamma^2 + (\gamma^*)^2}{\sqrt{\gamma^2 + (\gamma^*)^2 - 1}} \quad (19)$$

and

$$\mathcal{N}_{R,b}^{\text{re}} = \gamma^* \sqrt{\frac{\tilde{\Gamma} - 1}{2\tilde{\Gamma} - 1}}, \quad \mathcal{N}_{L,b}^{\text{re}} = (\mathcal{N}_{R,b}^{\text{re}})^*. \quad (20)$$

(iii) Phase III ($\gamma_2 > \gamma_1$): PT-broken regime with bound states.

In this regime, PT symmetry is spontaneously broken, and the bound states with imaginary eigenvalues emerge. Eigenmodes that give these states are $k_b = \frac{i}{2}\log(-1) + \frac{i}{2}\log|\Gamma| + n\pi$ for $n \in \mathbb{Z}$. Due to the two possible branches of $\log(-1)$, we have two right eigenfunctions and associated left eigenfunctions

$$R_{b,\pm}^{\text{im}}(x) = \begin{cases} \mathcal{N}_{R,b}^{\text{im}}(\pm i)^x e^{-x/\xi} & x > 0 \\ \mathcal{N}_{R,b}^{\text{im}} \frac{1}{\gamma^*} & x = 0, \\ \frac{\gamma^2}{|\gamma|^2} \mathcal{N}_{R,b}^{\text{im}}(\mp i)^x e^{x/\xi} & x < 0 \end{cases} \quad (21)$$

$$L_{b,\pm}^{\text{im}}(x) = \begin{cases} \mathcal{N}_{L,b}^{\text{im}}(\mp i)^x e^{-x/\xi} & x > 0 \\ \mathcal{N}_{L,b}^{\text{im}} \frac{1}{\gamma} & x = 0, \\ \frac{(\gamma^*)^2}{|\gamma|^2} \mathcal{N}_{L,b}^{\text{im}}(\pm i)^x e^{x/\xi} & x < 0 \end{cases} \quad (22)$$

Eigenvalues of these imaginary bound states are (see Fig. 3)

$$E_{b,\pm}^{\text{im}} = \mp i \frac{\gamma^2 + (\gamma^*)^2}{\sqrt{1 - \gamma^2 - (\gamma^*)^2}} = \mp 2i \frac{\tilde{\Gamma}}{\sqrt{1 - 2\tilde{\Gamma}}}. \quad (23)$$
<figure><img src="image_6.png" /><figcaption>FIG. 3. Imaginary bound state energy $E_{b, \pm}^{\text{im}}$ Eq. (23) [panel (a)] and a localization length $\xi$ [panel (b)] as a function of a coupling argument $\varphi(\gamma = Je^{i\varphi})$ for different values of a coupling magnitude $J$. Inset in panel (b): $|\psi|^2$ of an imaginary bound state for $\varphi = \pi/2$. Lines represent Eq. (23) (a) and Eqs. (22) and (21) (b), while points label numerical results. Localization length $\xi$ diverges and becomes larger than system size $N$ in the vicinity of $\varphi = \pi/4$ and, therefore, imaginary bound states with $\xi > N$ are absent (here $N = 1025$).</figcaption></figure>

Later, we shall also denote the absolute value of the bound state eigenenergies, whether real or imaginary, as $E_b = 2|\tilde{\Gamma}|/\sqrt{|1-2\tilde{\Gamma}|}$. It is readily seen that the right and left eigenfunctions are indeed orthogonal with the normalization constants

$$ \mathcal{N}_{R,b}^{\text{im}} = \gamma^* \sqrt{\frac{1-\tilde{\Gamma}}{1-2\tilde{\Gamma}}}, \quad \mathcal{N}_{L,b}^{\text{im}} = (\mathcal{N}_{R,b}^{\text{im}})^*, \quad (24) $$

which is consistent with the expected symmetry $(R_{b,+}^{\text{im}}(x))^* = R_{b,-}^{\text{im}}(-x)$.

The observations made above imply that, in our PT-symmetric RLM (1), the spontaneous breaking of PT-symmetry is solely caused by the presence of bound states that are localized across the impurity. Furthermore, as $\gamma_1$ and $\gamma_2$ approach, the two states characterized by $R_{b, \pm}^{\text{im}}(x)$ coalesce and form exponential points with diverging localization length $\xi$. In general eigenenergies with positive imaginary values inevitably give rise to exponential growth of physical quantities, such as correlation functions, in time evolution. This is obviously incompatible with unitarity and causality that are required for any physically sensible systems. The existence of such anomalous excitations in non-Hermitian systems was actually proposed decades ago in the imaginary Sine-Gordon model [12–14] and coined *monstrons*. Such

exponential divergency could be, however, canceled out by normalizing the time-evolved state $|\psi(t)\rangle$ by its biorthogonal norm $\langle\tilde{\psi}(t)|\psi(t)\rangle$, where $\langle\tilde{\psi}(t)|$ is the associated state of $|\psi(t)\rangle$ (see Appendix A for the discussion about the associated state and the biorthogonal norm).

We finally note that our approach does not work on the line of exceptional points $\gamma_1 = \gamma_2$ where the Hamiltonian (1) is not diagonalizable, and accordingly the wave function is not normalizable [24]. One therefore needs to devise an alternative method to study the spectrum at these particular points.

### III. DOT DENSITY

#### A. Diagonalized Hamiltonian and the ground state

Let us first define mode operators $c_{R,k}^{\dagger}$ and $anc_{L,k}^{\dagger}$ for real $k$ as

$$ c_{R,k\sigma}^{\dagger} = \sum_{x} R_{\sigma}(k; x) c_{x}^{\dagger}, \quad c_{L,k\sigma}^{\dagger} = \sum_{x} L_{\sigma}(k; x) c_{x}^{\dagger}, \quad (25) $$

and for complex $k$ as

$$ c_{R,\epsilon}^{\dagger} = \sum_{x} R_{b,\epsilon}(x) c_{x}^{\dagger}, \quad c_{L,\epsilon}^{\dagger} = \sum_{x} L_{b,\epsilon}(x) c_{x}^{\dagger}, \quad (26) $$

so that $|R_{\sigma}(k)\rangle = c_{R,k}^{\dagger}|0\rangle$ and $|L_{\sigma}(k)\rangle = c_{L,k}^{\dagger}|0\rangle$, and $|R_{b,\epsilon}\rangle = c_{R,\epsilon}^{\dagger}|0\rangle$ and $|L_{b,\epsilon}\rangle = c_{L,\epsilon}^{\dagger}|0\rangle$. Using the orthogonality relation $\int dx L_{\sigma}^{*}(k; x) R_{\sigma'}(k'; x) = \delta_{k,k'} \delta_{\sigma,\sigma'}$, these are equivalently expressed as

$$ c_x^\dagger = \sum_{k,\sigma} L_\sigma^*(k; x) c_{R,k\sigma}^\dagger, \quad c_x = \sum_{k,\sigma} L_\sigma(k; x) c_{R,k\sigma} \quad (27) $$

or

$$ c_x^\dagger = \sum_{k,\sigma} R_\sigma^*(k; x) c_{L,k\sigma}^\dagger, \quad c_x = \sum_{k,\sigma} R_\sigma(k; x) c_{L,k\sigma}, \quad (28) $$

where summation is taken over both real and complex $k$.<fnref n="1" /> Notice that $c_{R,k\sigma}^{\dagger}$ and $c_{L,k\sigma}^{\dagger}$ satisfy the anticommutation relation $\{c_{R,k\sigma}^{\dagger}, c_{L,k'\sigma'}^{\dagger}\} = \delta_{k,k'} \delta_{\sigma,\sigma'}$. In terms of these mode operators, (1) can be diagonalized as

$$ H = \sum_{k,\sigma} E_k c_{R,k\sigma}^{\dagger} c_{L,k\sigma} + \sum_{\epsilon=\pm} E_{b,\epsilon} c_{R,\epsilon}^{\dagger} c_{L,\epsilon}, \quad (29) $$

where the second term accounts for the possible bound states that appear in phases II and III (in phase I this term vanishes) with $E_{b,\epsilon}$ being either $E_{b,\epsilon}^{\text{re}}$ or $E_{b,\epsilon}^{\text{im}}$. We shall construct the ground state of PTRLM as in the Hermitian case, i.e., by filling modes that lower the real part of the energies. When PT symmetry is not broken, there is no ambiguity in defining the ground state as per the above principle, hence the left and the right ground state can be built as

$$ |\text{GS}\rangle_R = \prod_{\sigma} \prod_{k \in \text{Fermi sea}} c_{R,k\sigma}^{\dagger} |0\rangle, \quad (30) $$

$$ |\text{GS}\rangle_L = \prod_{\sigma} \prod_{k \in \text{Fermi sea}} c_{L,k\sigma}^{\dagger} |0\rangle. \quad (31) $$

<footnote n="1">¹It is important to note that actually $k$ depends on $\sigma$ through (13), hence the summation over $k$ and $\sigma$ are not independent. Here, for ease of notation we decided not to explicitly show this condition in writing the summations.</footnote>
Notice that the ground state can be uniquely constructed even if PT symmetry is spontaneously broken (i.e., phase III), except when the system is at half filling. At half filling ($\mu = 0$), the real part of the bound states lies on the boundary of the Fermi sea, hence it could be either filled or not filled. In this paper, we define the ground state in phase III at half filling by filling those two bound states so as to avoid complex ground state energies. We stress that there is no fundamental reason why we have to do so; we might as well fill just one of them at the price of having the ground state with complex energy. Our claim here is that our way of constructing the ground state is physically sensible (no complex ground state energy) and also gives rise to the behavior of the dot density that is consistent with previous literature [10]. Furthermore, the analytic continuation to the dot density, which we will elaborate on in Sec. VI, works only with this choice of the ground state.

## B. Dot density, susceptibility, and correlation function

One of the simplest yet nontrivial objects in quantum impurity systems is the dot density. Let us first compute it when PT symmetry is not spontaneously broken. It is defined by

$$d = \mathrm{Tr}(\rho c_0^\dagger c_0) = \sum_k \frac{L_+^*(k; 0)R_+(k; 0)}{1 + e^{\beta(E_k - \mu)}} + \sum_\epsilon \frac{L_{b, \epsilon}^*(0)R_{b, \epsilon}(0)}{1 + e^{\beta(E_{b, \epsilon} - \mu)}}. \quad (32)$$

The second term in (32) vanishes in phase I, and d reads simply, in the thermodynamic limit,

$$\int_{-\pi}^{\pi} \frac{dk}{2\pi} \frac{1}{1 + e^{\beta(E_k - \mu)}} \frac{\tilde{\Gamma}^2}{\tilde{\Gamma}^2 + (1 - 2\tilde{\Gamma})\cos^2 k} \quad (33)$$

If we take the real value limit $\gamma_2 = 0$, this reproduces the known result in the real RLM. In the phase II, we have contributions from the bound states, which can be calculated as

$$d_{\text{bound}}^{\text{re}} = \sum_{\epsilon} \frac{L_{b, \epsilon}^{*}(0) R_{b, \epsilon}(0)}{1 + e^{\beta(E_{b, \epsilon} - \mu)}} = \frac{2 - \gamma^2 - (\gamma^*)^2}{1 - \gamma^2 - (\gamma^*)^2} \frac{1 + e^{-\beta\mu} \cosh \beta E_b}{1 + 2e^{-\beta\mu} \cosh \beta E_b + e^{-2\beta\mu}}, \quad (34)$$

where we recall $E_b = |E_{b, \pm}^{\text{re}}|$. The dot density in phase II is then given by

$$d = \int_{-\pi}^{\pi} \frac{dk}{2\pi} \frac{1}{1 + e^{\beta(E_k - \mu)}} \frac{\tilde{\Gamma}^2 \sin^2 k}{\tilde{\Gamma}^2 + (1 - 2\tilde{\Gamma}) \cos^2 k} + d_{\text{bound}}^{\text{re}}. \quad (35)$$

Note that the dot densities in both phase I and II are identical to 1/2 when $\mu = 0$ due to the particle-hole symmetry (see Appendix B for the proof).

Next, we turn our attention to phase III where PT symmetry is broken by the presence of complex eigenstates. It turns out that thermodynamics is ill defined in this phase: Any thermal state characterized by a Gibbs ensemble is dynamically unstable. Concretely speaking, a Gibbs measure $\rho = \exp(-\beta(H - \mu Q))/Z$ is not invariant under time evolution

due to the breakdown of Heisenberg picture (see the discussion in Appendix A). Typically it is expected that stationary states in this phase are not universal and strongly depending on the initial state. Nonetheless, the average of an arbitrary observable in the ground state is still stationary and therefore makes physical sense. Then, we can calculate the dot density as

$$d = \int_{-k_F}^{k_F} \frac{dk}{2\pi} \frac{\tilde{\Gamma} \sin^2 k}{\tilde{\Gamma}^2 + (1 - 2\tilde{\Gamma}) \cos^2 k} + \frac{2 - \gamma^2 - (\gamma^*)^2}{1 - \gamma^2 - (\gamma^*)^2}. \quad (36)$$

Note that the appearance of the second term can amount to the dot density larger than 1. This is because we generically fill both of two bound states with pure imaginary eigenvalues when constructing the ground state. This is somewhat pathological in the sense that it contradicts with the Pauli's principle if one cares only about the real part of the energies. Nevertheless, because the imaginary parts of the two bound states differ, these modes are allowed to be occupied simultaneously, resulting in the anomalous behavior of the dot density.

In fact, in non-Hermitian systems, there is no restriction for the value of the dot density in the ground state. Indeed, in the usual Hermitian case, one can show $0 \le d = \langle \text{GS} | c_0^\dagger c_0 | \text{GS} \rangle \le 1$ as follows: First we note that d can be written as $d = ||c_0|\text{GS}\rangle||^2$, which is clearly non-negative. Likewise, $1 - d = \langle \text{GS} | c_0 c_0^\dagger | \text{GS} \rangle = ||c_0^\dagger | \text{GS} \rangle||^2$, which again is non-negative. Merging them together, it then immediately follows $0 \le d \le 1$. This line of argument, however, fails to be valid in the non-Hermitian case due to the defining feature of the biorthogonal basis: $(|\text{GS}\rangle_R)^\dagger \neq L \langle \text{GS}|$. To be more precise, $L \langle \text{GS} | c_0^\dagger$ is not the associated state of $c_0^\dagger | \text{GS} \rangle_R$, hence $L \langle \text{GS} | c_0^\dagger c_0 | \text{GS} \rangle_R$ cannot be cast into the form of the biorthogonal norm. On the physical ground, however, an interpretation of the behavior of the dot density in the PT-broken regime is still missing and calls for further investigations.

That said, we can alternatively look at the dot density that is made solely from the right ground state, i.e., $|\text{GS}\rangle_R$. Such a dot density $d_R = R \langle \text{GS} | c_0^\dagger c_0 | \text{GS} \rangle_R / R \langle \text{GS} | \text{GS} \rangle_R$ is real and smaller than 1 by construction. Although this is not a right quantity to look at from the biorthogonal quantum mechanical perspective, we provide comparison between two approaches in Fig. 10 in Appendix D.

Having pointed out the peculiar behavior of the dot density in phase III, let us address a natural question that arises from the above observation: what is actually happening when PT symmetry is broken, and in particular, does Kondo screening disappear in the phase III or not? To answer this question, let us recall the phenomenology of Kondo screening in phase I; Kondo screening in general refers to the fact that the impurity is healed with the leads in the IR, forming a spin singlet. In the scaling limit, this is typically signalled by the behavior of the dot density in the IR, i.e., if the dot density goes to 1/2 in the IR, we interpret it as the evidence that Kondo screening is taking place. This is in contrast to what happens in the UV, i.e., $d \to 1$, which indicates that the spin is effectively free. Now, when PT symmetry is spontaneously broken, the IR and UV behaviors of the dot density largely differ from those in the phase I. Namely, along the flow from the UV to the
IR, the value of the dot density changes from 1 to 3/2, which we shall see in the next subsection. We then interpret from this that the Kondo screening ceases to exist in phase III. We, however, emphasise that this is not exactly what we expect from the previous studies on ISG. For instance, in Ref. [12], it was discovered that the starting and end points of the RG flow in ISG are characterized by almost the same fixed points, the c = 1 CFT, with the only difference being the compactification radii of them. Translating this situation into PTRLM, we are led to expect that the boundary RG flow in phase III should start from the Neumann boundary condition (impurity being disconnected from the leads) and return to the same boundary condition again in the IR. This will be then reflected to the behavior of d, namely d should approach to 1 in the IR, which is not quite what we observe. For now we have no convincing explanation to reconcile this dichotomy; one possibility is that the IR behavior in phase III is a peculiarity in the free theory, and the expected IR physics (i.e., return to the Neumann boundary condition) might be recovered by the inclusion of interaction.

As a final remark for the dot density, it is also worth mentioning what would happen to the dot density in phase III if we were to fill just one of the bound states when constructing the ground state. It turns out that the value of the dot density is the same regardless of which bound state is filled, and, rather curiously, d = 1/2 at half filling. This implies that particle-hole symmetry is not broken in phase III with this choice of the ground state, or at least the breakdown of it, if any, cannot be deduced from the dot density at half filling.

Another quantity that characterizes the equilibrium property of the impurity system is the dot susceptibility defined by $\chi(T) = \partial d/\partial \mu|_{\mu=0}$. From (33) and (35), it can be easily seen that indeed the susceptibility follows Curie's law $\chi(T) \sim 1/T$, indicating that the impurity behaves like a free spin when temperature is high. Interesting physics emerges at low temperature; in any phase, we can readily see that the susceptibility reaches to the stationary value in the ground state $\chi(T=0) = 1/(2\pi\tilde{\Gamma})$, which is a hallmark of the Kondo screening. Therefore the existence of the PT phase transition cannot actually be inferred from the behavior of the susceptibility. Notice also that the susceptibility blows up as the system approaches the PT critical line defined by $\gamma_1 = \gamma_2$. This is a genuine non-Hermitian effect, caused by the competition between the gain and loss terms.

Finally, to further strengthen our conclusion that there is no Kondo screening in phase III, we numerically compute the correlation function $C(x) = \langle GS|c_0^\dagger c_x|GS\rangle_R$. In the scaling regime, it has been known that the presence of Kondo screening manifests itself as a crossover behavior of $C(x)$: across the characteristic scale $1/T_K$, where $T_K = 2\tilde{\Gamma}$ is the (generalized) Kondo temperature, the way $C(x)$ decays change [16]. Indeed, in Fig. 4(a), there is a crossover from the $-\log x$ decay to the power law decay $1/x$ across $x \sim 1/T_K$ when PT symmetry is unbroken. The situation significantly alters in the PT-broken case: There is no crossover nor a well-defined $1/x$ decay anymore, although the curve still has scaling properties, see Fig. 4(b). Therefore this observation is consistent with our conclusion that the Kondo screening is absent in the PT-broken case.

<figure><img src="image_7.png" /><figcaption>FIG. 4. Rescaled correlation function $\langle c_0^\dagger c_x \rangle$ vs rescaled distance $xT_K$ from the impurity. Panel (a): PT-unbroken regime with the coupling argument $\varphi = \pi/8$ and different coupling constants $J = 0.02, \dots, 0.5$. For short distances the correlation function behaves as $\sim \log(x)$, whereas far away it algebraically decays $\sim 1/x$. Panel (b): PT-broken regime with $\varphi = 3\pi/8$ and $J = 0.05, \dots, 0.5$. The curve preserves scaling properties, however there is no clear algebraic or log-like decay. System size in both panels $N = 2007$.</figcaption></figure>

### C. UV and IR expansions of the dot density at $T=0$

At $T=0$, it is in fact possible to carry out the integrations explicitly (see Appendix B for the details). For instance, in phase I, we have the following UV expansion (small coupling)

$$d = 1 + \frac{\tilde{\Gamma}}{\pi(1 - 2\tilde{\Gamma})} \arctan \frac{\sqrt{4 - \mu^2}}{\mu} + \frac{1 - \tilde{\Gamma}}{\pi u(\tilde{\Gamma})(1 - 2\tilde{\Gamma})} \sum_{n=1}^{\infty} \frac{(-1)^n}{2n-1} u(\tilde{\Gamma})^{2n}, \quad (37)$$

where $u(\tilde{\Gamma}) = -\tilde{\Gamma} \tan k_F/(1 - \tilde{\Gamma}) = \tilde{\Gamma} \sqrt{4 - \mu^2}/(\mu(1 - \tilde{\Gamma}))$ which is a monotonically increasing function of $\tilde{\Gamma}$ in phase I (i.e., $\tilde{\Gamma} < 1$). Note that the convergence of radius of the series is $u(\tilde{\Gamma}) < 1$. The IR expansion reads, on the other hand,

$$d = \frac{1}{2} - \frac{\tilde{\Gamma}}{\pi(1 - 2\tilde{\Gamma})} \arctan \frac{\mu}{\sqrt{4 - \mu^2}} - \frac{u(\tilde{\Gamma})(1 - \tilde{\Gamma})}{\pi(1 - 2\tilde{\Gamma})} \sum_{n=1}^{\infty} \frac{(-1)^n}{2n-1} u(\tilde{\Gamma})^{-2n}, \quad (38)$$
which converges only when $u(\tilde{\Gamma}) < 1$. Similarly, in phase III, the UV and IR expansions are

$$d = 1 + \frac{\tilde{\Gamma}}{\pi(1 - 2\tilde{\Gamma})} \arctan \frac{\sqrt{4 - \mu^2}}{\mu} - \frac{1 - \tilde{\Gamma}}{\pi v(\tilde{\Gamma})(1 - 2\tilde{\Gamma})} \sum_{n=1}^{\infty} \frac{(-1)^n}{2n - 1} v(\tilde{\Gamma})^{2n} \quad (39)$$

and

$$d = \frac{1}{2} + \frac{1 - \tilde{\Gamma}}{1 - 2\tilde{\Gamma}} - \frac{\tilde{\Gamma}}{\pi(1 - 2\tilde{\Gamma})} \arctan \frac{\mu}{\sqrt{4 - \mu^2}} + \frac{v(\tilde{\Gamma})(1 - \tilde{\Gamma})}{\pi(1 - 2\tilde{\Gamma})} \sum_{n=1}^{\infty} \frac{(-1)^n}{2n - 1} v(\tilde{\Gamma})^{-2n}, \quad (40)$$

where $v(\tilde{\Gamma}) = \tilde{\Gamma} \tan k_F/(1 - \tilde{\Gamma}) = -\tilde{\Gamma}\sqrt{4 - \mu^2}/(\mu(1 - \tilde{\Gamma}))$, and the UV and IR expansions are valid when $v(\tilde{\Gamma}) < 1$ and $v(\tilde{\Gamma}) > 1$, respectively. We note that, in the scaling limit $\gamma_1, \gamma_2, \mu \ll 1$, (37) and (38) reproduce the known expression [26,27] in the Toulouse limit of the anisotropic Kondo model. Similar formulas can also be obtained by the same manipulations in phase II and presented in Appendix B. Profiles of the dot density at $T = 0$ as a function of $\mu$ in different phases are depicted in panel (a) of Fig. 5.

These exact UV and IR expansions of the dot density allow us to observe the aforementioned unusual phenomenon: The dot density no longer equals 1/2 at half filling, indicating that particle-hole symmetry is spontaneously broken. This observation further implies a peculiar behavior of the magnetization $m = 1/2 - d$. For the sake of simplicity, let us focus on the scaling limit $\gamma_1, \gamma_2, \mu \ll 1$. When $\gamma_2$ crosses the value of $\gamma_1$ from below, the dot density jumps from 1/2 to 3/2, see Fig. 5. This implies the PT phase transition with the concomitant breakdown of particle-hole symmetry triggers the jump of magnetization from 0 to -1. A similar behavior of the magnetization was also observed in Ref. [10] where a jump of the magnetization occurred as one increases the non-Hermiticity of the non-Hermitian Kondo model. Furthermore, the emergence of edge modes with pure imaginary eigenvalues, which is associated to spontaneous broken particle-hole symmetry, was also previously observed in photonic graphene with gain and loss [25].

# IV. g FUNCTION AND BOUNDARY FREE ENERGY

## A. g function

It is always instructive to study the g function to see the UV and IR behavior of a given impurity system. It is known that, in quantum impurity systems in the scaling regime, the g function undergoes a monotonic change under the boundary renormalization group flow, satisfying $g_{UV} < g_{IR}$ [28]. In PTRLM, we can readily observe the crossover physics of the g function as T increases. The g function is defined in terms of the impurity entropy

$$\log g = S_{\text{imp}} = S - S_0, \quad (41)$$

<figure><img src="image_8.png" /><figcaption>FIG. 5. Panel (a): dot density as a function of the chemical potential for different values of the coupling magnitude J and argument $\varphi$ ($\gamma = Je^{i\varphi}$). Panel (b): dot density as a function of the coupling argument $\varphi$ for different values of the coupling magnitude J and the fixed chemical potential $\mu = 0$. A nonanalytic behavior appears as a jump between the PT-unbroken and the PT-broken phases in the vicinity of the line of exceptional points $\varphi = \pi/4$. Numerical predictions depart from analytic curves near $\varphi = \pi/4$ but slowly converges to it with increasing a system size N, here N = 2019. Panel (c): the dot density as a function of temperature T for $\mu = 0.2$ and $J = 0.3$, the color encodes fractions of $\varphi$ between real line $\varphi = 0$ and the critical line $\varphi = \pi/4$.</figcaption></figure>

where S is the thermodynamic entropy of the RLM, and $S_0$ is the same quantity without the impurity (i.e., two decoupled leads). Recall that the partition function is given by $Z = \prod_{\sigma, k_n} (1 + e^{-\beta(E_k - \mu)})$, where $k_n$ satisfies

$$k_n = \frac{(2n + 1)\pi}{N + 2} + \frac{k\delta_{\sigma}(k)}{N + 2}, \quad n = 1, \dots, \frac{N + 2}{2}. \quad (42)$$
<figure><img src="image_9.png" /><figcaption>FIG. 6. Impurity entropy $log g$ as a function of temperature $T$ for different values of chemical potential $\mu$. The color encodes $\mu$, while the line types encodes value of a coupling constant argument $\varphi$ ($\gamma = Je^{i\varphi}$). Points represent numerics, while lines represent analytical result Eq. (45).</figcaption></figure>

Once again, we note that $k$ implicitly depends on $\sigma$. Then, the root density $\rho(k)$ is

$$ \rho_{\sigma}(k_n) = \lim_{N \to \infty} \frac{1}{(N+2)(k_{n+1} - k_n)} \\ = \frac{1}{2\pi} \left( 1 + \frac{2}{N+2} \frac{1 - \Gamma_{\sigma}^2}{\Delta_{\sigma}(k)} + \mathcal{O}(N^{-2}) \right). \quad (43) $$

Therefore, under the thermodynamic limit, the partition function becomes

$$ \log Z = (N + 2) \sum_{\sigma, k_n} \frac{1}{N + 2} \rho_{\sigma}(k_n) \log(1 + e^{-\beta(E_{k_n} - \mu)}) \\ = 2(N + 2) \int_{0}^{\pi} \frac{dk}{2\pi} \log(1 + e^{-\beta(E_k - \mu)}) \\ + \sum_{\sigma} \int_{0}^{\pi} \frac{dk}{2\pi} \frac{1 - \Gamma_{\sigma}^2}{\Delta_{\sigma}(k)} \log(1 + e^{-\beta(E_k - \mu)}), \quad (44) $$

where the nonextensive part represents the impurity contribution. Since, without the impurity, $\log Z$ is simply $\log Z = (N + 2) \int_{0}^{\pi} \frac{dk}{\pi} \log(1 + e^{-\beta(E_k - \mu)})$, using the formula $S = \log Z - \beta \partial_{\beta} \log Z$, the $g$ function can be obtained as

$$ \log g = -2 \int_{0}^{\pi} \frac{dk}{\pi} \left[ \frac{(2\tilde{\Gamma} - 1) \cos 2k - 1}{4(\tilde{\Gamma}^2 + (1 - 2\tilde{\Gamma}) \cos^2 k)} + \frac{1}{2} \right] \\ \times \left[ \log(1 + e^{-\beta(E_k - \mu)}) + \frac{\beta(E_k - \mu)}{1 + e^{\beta(E_k - \mu)}} \right]. \quad (45) $$

One can readily confirm that when $\beta \gg 1$, $\log g \to 0$, indicating that at low temperature, the impurity is effectively hybridized with the wires. On the other hand, for $\beta \ll 1$, $\log g \to \log 2$, which suggests that at high temperature, the impurity is decoupled from the wires, attaining two degrees of freedom (occupied or unoccupied dot state), see Fig. 6.

Next we treat phase II. The additional term $S_{\text{imp,b}}$ to (41) comes from two bound states, reading

$$ S_{\text{imp,b}} = \log[1 + e^{2\beta\mu} + 2e^{\beta\mu} \cosh(\beta E_b)] \\ - 2\beta e^{\beta\mu} \frac{\mu e^{\beta\mu} + \mu \cosh \beta E_b + E_b \sinh \beta E_b}{1 + e^{2\beta\mu} + 2e^{\beta\mu} \cosh(\beta E_b)}. \quad (46) $$

Therefore the $g$ function in phase II is the sum of (45) and (46)

$$ \log g = -2 \int_{0}^{\pi} \frac{dk}{\pi} \left[ \frac{(2\tilde{\Gamma} - 1) \cos 2k - 1}{4(\tilde{\Gamma}^2 + (1 - 2\tilde{\Gamma}) \cos^2 k)} + \frac{1}{2} \right] \\ \times \left[ \log(1 + e^{-\beta(E_k - \mu)}) + \frac{\beta(E_k - \mu)}{1 + e^{\beta(E_k - \mu)}} \right] + S_{\text{imp,b}}. \quad (47) $$

The asymptotic behavior in phase II turns out to be the same as in phase I. This is because the bound state energy satisfies $E_b > 2$ in phase II ($\gamma_1^2 - \gamma_2^2 > 1$), and as such $E_b > \mu$ always holds when the ground state is nontrivial (the ground state becomes nontrivial only when $\mu < E_k$ for some $k \in [-\pi, \pi]$, i.e., $0 < \mu < 2$).

## B. Boundary free energy

Now, we turn our attention to phase III. As argued in the previous section, we do not have any invariant statistical measure when PT symmetry is spontaneously broken, and hence, $g$ function is not suitable for analyzing the UV and IR behaviors of the system. That being said, there exists an alternative quantity, the boundary free energy, from which we can infer the two asymptotic behaviors of the system. The boundary free energy $f$ can be obtained from the definition $m = \frac{\partial f}{\partial \mu}$ [note that $\mu$ in the RLM language is the opposite of the natural magnetic field in the corresponding (anisotropic) Kondo problem]: We will compare the behaviors of the boundary free energy in phase I and III in the rest of this section. To avoid unnecessary complication, let us focus on the real coupling case $\gamma = \gamma_1$ (phase I) and the pure imaginary coupling case $\gamma = i\gamma_2$ (phase III) in the scaling limit. Then, according to (37) and (39), the boundary free energy in phase I and III is given by, respectively,

$$ f_1 = -\frac{\mu}{2} + \frac{2\gamma_1^2}{\pi} \left( \frac{\mu}{2\gamma_1^2} \arctan \frac{2\gamma_1^2}{\mu} + \frac{1}{2} \log \left( 1 + \frac{\mu^2}{4\gamma_1^4} \right) \right) $$

$$ f_{\text{III}} = -\frac{\mu}{2} - \frac{2\gamma_2^2}{\pi} \left( \frac{\mu}{2\gamma_2^2} \arctan \frac{2\gamma_2^2}{\mu} + \frac{1}{2} \log \left( 1 + \frac{\mu^2}{4\gamma_2^4} \right) \right). \quad (48) $$

Accordingly their asymptotic behaviors are

$$ f_1 \to \begin{cases} -\frac{\mu}{2} & \gamma_1 \to 0 \\ 0 & \gamma_1 \to \infty \end{cases} \quad (49) $$

$$ f_{\text{III}} \to \begin{cases} -\frac{\mu}{2} & \gamma_2 \to 0 \\ -\mu & \gamma_2 \to \infty \end{cases}. \quad (50) $$

Recall that the boundary free energy $f$ encodes information on the impurity degrees of freedom at the UV (resp. IR) fixed point when $\gamma \to 0$ (resp. $\gamma \to \infty$). In the spin $j$ anisotropic Kondo model, it is known that $f$ behaves as $f \to -j\mu$ and $f \to -(j - \frac{1}{2})\frac{\mu}{\tilde{\epsilon}}$ near the UV and IR fixed
points, respectively. Note that near the IR fixed point, the spin is renormalized by the Luttinger coupling constant $\ell$, which comes from the fact that the underlying model is anisotropic, i.e., SU(2) symmetry is spoiled.

The choice $j = \frac{1}{2}$, $\ell = \frac{1}{2}$ corresponds to the scaling limit of the resonant level model for which we recover (49). This is essentially what we saw in (45): Near the UV fixed point, the impurity spin is effectively decoupled while near the IR fixed point, the impurity is screened and accordingly $f \to 0$. The UV behavior of $f$ in phase III is the same as that in phase I, as the impurity is decoupled anyway even if coupling constants are complex. The IR behavior, however, comes as a surprise; the impurity remains active even when the couplings are strong enough. This is in accordance with the anomalous behavior of the dot density, and again we understand this as the signal of the lack of Kondo screening. Remarkably, we also observe [29] that this peculiar IR behavior of $f$ persists even away from the free limit $\ell = \frac{1}{2}$ and for arbitrary spin $j$ in the anisotropic Kondo model, with $f \to -j \frac{\mu}{\ell}$, suggesting that this is a universal phenomenon in quantum impurity systems when PT symmetry is broken.

## V. LOSCHMIDT ECHO IN THE PT-SYMMETRIC RLM

One of the simplest objects that encodes nontrivial information on dynamics is the *Loschmidt echo*, which is defined by

$$ \bar{\mathcal{L}}(t) = {}_0\langle \text{GS} | e^{itH_0} | \Psi(t) \rangle, \quad (51) $$

where $|GS\rangle_0$ denotes the ground state of two disconnected chains $H_0 = H_L + H_R$. Here $|\Psi(t)\rangle$ denotes the normalized time-evolved state $|\Psi(t)\rangle = |\psi(t)\rangle/|\langle\bar{\psi}(t)|\psi(t)\rangle|$, where $|\psi(t)\rangle = e^{-iHt}|\text{GS}\rangle_0$, and $\langle\bar{\psi}(t)|$ is its associated state (A8). Note that the norm is 1 as long as PT symmetry is not broken, in which case dynamics is unitary. Practically the norm plays a role of taming the exponential growth of the Loschmidt echo in phase III. We note that a similar normalization is necessary whenever Loschmidt echoes for Hermitian Hamiltonians are calculated in imaginary time using CFT [30] or integrability methods [31]: In that case, phase terms $e^{iEt}$ indeed cannot be neglected as they acquire exponential growth. Loschmidt echo with this protocol was also studied previously in noninteracting impurity systems using the form factor approach [32] and the finite volume free fermion technique which is akin to ours [33].

The Loschmidt echo simply characterizes the return probability of the initial state. Each chain can be diagonalized easily as

$$ H_L = \sum_{k>0} E_k d_{1k}^\dagger d_{1k}, \quad H_R = \sum_{k>0} E_k d_{2k}^\dagger d_{2k}, \quad (52) $$

where canonical fermionic operators are given by

$$ d_{1k} = \sum_{x>0} \phi(k; x) c_x, \quad d_{2k} = \sum_{x<0} \phi(k; x) c_x \quad (53) $$

with $\phi(k; x) = 2 \sin kx / \sqrt{N+2}$ and $k$ is quantized as $k_n = 2\pi n/(N+2)$, $n = 1, \dots, (N+2)/2$. Therefore $H_0$ can be written as $H_0 = \sum_{k>0} E_k (d_{1k}^\dagger d_{1k} + d_{2k}^\dagger d_{2k})$ with the ground

state

$$ |\text{GS}\rangle_0 = \prod_{n \in \text{FS}} c_{k_n}^\dagger, \quad c_{k_n} = \begin{cases} d_{1k_n} & n > 0 \\ d_{2k_n} & n < 0 \end{cases}, \quad (54) $$

where the integer set FS = $\{-M, \dots, -1, 1, \dots, M\}$ corresponds to the modes comprising Fermi seas of each lead with $M = (N+2)/2$. Invoking that the Loschmidt echo (51) has a structure of the Slater determinant, we immediately realize that $\mathcal{L}(t)$ can be written as a determinant of a matrix involving overlap matrices

$$ \mathcal{L}(t) = e^{itE_{\text{GS}}} \det_{i,j \in \text{FS}} \mathbf{M}(k_i, k_j; t), \quad (55) $$

where $\mathbf{M}(k_i, k_j; t) = \sum_m e^{-iE_{\text{mt}}t} \mathbf{L}^*(k_i; m)\mathbf{R}(k_j; m)$. Here the overlap matrices are defined as $\mathbf{L}^*(k; m) = \sum_x \vartheta(kx)\phi(|k|; x)\mathbf{L}^*(m; x)$ and $\mathbf{R}(k; m) = \sum_x \vartheta(kx)\phi(|k|; x)\mathbf{R}(m; x)$, where $\vartheta(x)$ is the Heaviside step function. After carefully evaluating the matrix $\mathbf{M}(k_i, k_j; t)$ (see Appendix C for the details), we end up with the following Fredholm determinant

$$ \bar{\mathcal{L}}(t) = \begin{cases} \det_{i,j \in \text{FS}} \left(1 + \frac{2\pi}{N+2} \mathbf{B}(t)\right)_{ij} & \text{Phase I} \\ \det_{i,j \in \text{FS}} \left(1 + \frac{2\pi}{N+2} \mathbf{B}(t) + \mathbf{B}_b^{\text{re}}(t)\right)_{ij} & \text{Phase II,} \\ \frac{\det_{i,j \in \text{FS}} \left(1 + \frac{2\pi}{N+2} \mathbf{B}(t) + \mathbf{B}_b^{\text{im}}(t)\right)_{ij}}{\sqrt{|\det_{i,j \in \text{FS}} \left(1 + \frac{2\pi}{N+2} \mathbf{N}(t)\right)|}} & \text{Phase III} \end{cases} \quad (56) $$

where the matrix $\mathbf{B}(t)$, $\mathbf{B}_b^{\text{re}}(t)$, $\mathbf{B}_b^{\text{im}}(t)$, and $\mathbf{N}(t)$ is defined in (C30), (C43), (C52), and (C54), respectively. We are primarily interested in the behavior of $\sqrt{t}|\bar{\mathcal{L}}(t)|^2$ as a function of $t$. The reason for having the additional factor $\sqrt{t}$ is the following: The large $t$ behavior of the echo is dominated by the low-energy process, under which case the dot is fully healed with the leads. As such, for large $t$, $\bar{\mathcal{L}}(t)$ is expected to match with that describing the quench from the two disconnected chains to the fully homogeneous chain, which, for large enough $t$, goes as $t^{-1/4}$ [32], where the exponent is related to the scaling dimension of the boundary condition changing operator from free to fixed boundary condition in the critical boundary Ising model [32]. We therefore multiply $\sqrt{t}$ to $|\bar{\mathcal{L}}(t)|^2$ in order to cancel out such decay.

Let us first look at the behavior in phase I. We find that the echo exhibits a crossover behavior whose time scale is controlled by the generalized Kondo temperature $T_K$. The temperature is in fact nothing but the absolute value of the eigenenergies of the bound states: $T_K = |E_b|$, appearing in phase II and III. This phenomenon, which can be regarded as a dynamical manifestation of the Kondo screening, was already observed in Refs. [32,34], but a new observation here is that such dynamical Kondo effect persists even away from the Hermitian case, as long as the PT symmetry of the system is unbroken. Furthermore, we also observe that the echo possesses the universal scaling form $\sqrt{t}T_K f(tT_K)$, which was also shown in Ref. [32] in the Hermitian case. We depict the behavior of the echo in phase I in Fig. 7 where we compare the analytic results (56) with independent free fermion numerics. We not only observe a beautiful agreement but also clearly see that the echo follows the universal curve, which has a crossover around $t \sim 1/T_K$. Such a crossover is a simple consequence of interpolating the UV and IR behavior of the echo: For a short time (UV), the dot is effectively screened from the
<figure><img src="image_10.png" /><figcaption>FIG. 7. Panel (a): rescaled Loschmidt echo $\sqrt{x}|\tilde{L}(x)|^2$ as a function of rescaled time $x = tT_K$ in phases I and III for different values of $\varphi$ (coupling constant $\gamma = 0.2e^{i\pi\varphi}$). Panel (b): rescaled Loschmidt echo $|\tilde{L}(x)|^2/T_K^2$ (scaling is $1/T_K^2$) as a function of rescaled time in phase II ($\gamma = 2e^{i\varphi}$). Scaling factor is chosen to be $T_K = |E_b|$. The colors label the values of $\varphi$, lines represent Fredholm determinant Eq. (56), while dots are obtained by numerics.</figcaption></figure>

chains, hence the system behaves as two disconnected chains, giving rise to a sharp buildup of the echo. For a large time (IR), however, the echo decays as $t^{-1/2}$ for the reason above, thus presents a bump when undergoing the transition from the UV to IR time scale. It is also noteworthy that, since $E_b \sim \tilde{\Gamma}$, the time needed for the crossover can be arbitrarily long near the PT critical line $\gamma_1 = \gamma_2$.

We next turn to phase II. In this phase the Loschmidt echo displays persistent oscillations whose frequency is controlled by the Kondo temperature $T_K$. These oscillations are caused by the bound states, thereby spoiling the dynamical Kondo screening which was observed in phase I. With the appropriate rescaling of the Loschmidt echo, which can be found numerically, we observe that the curves for different $\varphi$ can collapse onto a single universal curve for large $t$ [see panel (b) in Fig. 7].

Having analyzed the Loschmidt echo in phase I and II, let us finally deal with the most controversial phase: phase III. In this regime the echo is dominated by the contribution from the eigenmode with positive imaginary energy, hence $\mathcal{L}(t)$ exponentially grows. The norm then serves to cancel out such

explosion, yielding a finite result. Curiously, the normalized echo seems to show a similar crossover physics as time increases, and rather surprisingly, it turns out to decay in the same fashion as in phase I, i.e., decays as $\mathcal{L}(t) \sim t^{-1/4}$. Such crossover behavior in the PT-broken regime can be expected by invoking the similar argument as in the PT-unbroken regime. What is surprising here however is that, with the division by the nontrivial norm, the asymptotic behavior of the normalized echo follows the same pattern as in phase I, and the position of the bump is still controlled by the absolute value of the Kondo temperature. We emphasise that whether the mechanism of this decay is the same as in phase I is not at all obvious and necessitates a further understanding of phase III. The collapse of curves onto the single universal curve also occurs in this phase in the scaling regime, see Fig. 7.

## VI. APPROACHING THE DIFFERENT PHASES VIA ANALYTIC CONTINUATION

First, we note that in all regimes, the dot density is an analytical function of $\tilde{\Gamma}$ at fixed $\mu$. $d$ admits two different expansions depending on whether the coupling $\tilde{\Gamma}$ is small (UV) or large (IR), but these are two expansions of the same function, whose singularities lie in general elsewhere (more on this below). For instance, we have in regime I

$$d = 1 + \frac{\tilde{\Gamma}}{\pi(1 - 2\tilde{\Gamma})} \arctan \frac{\sqrt{4 - \mu^2}}{\mu} - \frac{1 - \tilde{\Gamma}}{\pi(1 - 2\tilde{\Gamma})} \arctan u(\tilde{\Gamma}) \quad \text{UV} \quad (57)$$

$$d = 1 + \frac{\tilde{\Gamma}}{\pi(1 - 2\tilde{\Gamma})} \arctan \frac{\sqrt{4 - \mu^2}}{\mu} - \frac{1 - \tilde{\Gamma}}{\pi(1 - 2\tilde{\Gamma})} \left( \frac{\pi}{2} - \arctan \frac{1}{u(\tilde{\Gamma})} \right) \quad \text{IR,} \quad (58)$$

where recall $u(\tilde{\Gamma}) = \frac{\tilde{\Gamma}\sqrt{4-\mu^2}}{\mu(1-\tilde{\Gamma})}$. We see that the IR expansion is obtained from the UV one by using the continuation of the function arctan on the positive real axis when its argument gets larger than one.

The same feature holds in regimes II and III. For instance in regime III we now have

$$d = 1 + \frac{\tilde{\Gamma}}{\pi(1 - 2\tilde{\Gamma})} \arctan \frac{\sqrt{4 - \mu^2}}{\mu} + \frac{1 - \tilde{\Gamma}}{\pi(1 - 2\tilde{\Gamma})} \arctan v(\tilde{\Gamma}) \quad \text{UV} \quad (59)$$

$$d = 1 + \frac{\tilde{\Gamma}}{\pi(1 - 2\tilde{\Gamma})} \arctan \frac{\sqrt{4 - \mu^2}}{\mu} - \frac{1 - \tilde{\Gamma}}{\pi(1 - 2\tilde{\Gamma})} \left( -\frac{\pi}{2} + \arctan \frac{1}{v(\tilde{\Gamma})} \right) \quad \text{IR,} \quad (60)$$

where $v(\tilde{\Gamma}) = -\tilde{\Gamma}\sqrt{4-\mu^2}/(\mu(1-\tilde{\Gamma})) = -u(\tilde{\Gamma})$. Here again we see that the IR expansion is obtained by “straightforward” continuation of the UV one.

This feature is common to other problems of this type in the Hermitian case: It has been observed in the Kondo and anisotropic Kondo models, and the problem of tunneling
between edges in the fractional quantum Hall effect, at least in the scaling limit (where, in particular, the FQHE tunneling problem is described by the boundary sine-Gordon model [35]). It is expected whenever the short coupling expansion of a physical quantity admits a finite radius of convergence where there is however no physical singularity.

More surprising maybe is the fact that the UV expansion of $d$ in regimes II and III can be obtained as well by analytical continuation of the UV expansion in regime I. This is clear if one considers for instance the UV expansions in regime I and III: The result in equation (59) is in fact the same as the one in (57) when expressed in terms of $\tilde{\Gamma}$. The result is in fact the same in regime II as well, and completely independent of the phase of $\Gamma$. This result is in fact quite reasonable if one imagines calculating a quantity such as $d$ perturbatively in $\tilde{\Gamma}$: This will produce a series integer powers of $\tilde{\Gamma}$, with amplitudes given by integrals of various correlations functions over the vacuum at $\tilde{\Gamma} = 0$. Since this vacuum does not depend on $\tilde{\Gamma}$, there should be in fact a unique UV expansion for $d$ (and other such perturbative physical quantities).

On the other hand, one has to be very careful with what happens in the IR regime. To obtain the results at large coupling in regime III from those in regime I, one needs to pay attention to the determination of the arctan involved in the problem. Specifically, starting with equation (58), the corresponding result in regime III is obtained not only by observing that $u(\tilde{\Gamma}) = -v(\tilde{\Gamma})$: Instead one must set

$$ \arctan \frac{1}{u(\tilde{\Gamma})} \rightarrow \pi - \arctan \frac{1}{v(\tilde{\Gamma})} \quad (61) $$

even though we recall that $v(\tilde{\Gamma}) = -u(\tilde{\Gamma})$. What happens of course is that the arctan function admits multiple determinations, and to see what happens to analytical continuations requires following trajectories on the corresponding Riemann surface. In other words, sending $z$ into $e^{i\pi z}$ and $z$ into $1/z$ are two operations that do not commute for the functions arctan $z = \frac{1}{2i} \ln \frac{1+i z}{1-i z}$, whose Riemann surface is branched along the imaginary axis for $|z| \ge 1$.

Extra factors of $\pi$ such as in (61) will be needed only after $\arg u(\tilde{\Gamma}) > \pi/2$. In our model where $u$ is real,<fnref n="2" /> this corresponds to $\varphi > \pi/2$, i.e., regime III indeed. We can predict from this that the results in regime II at large coupling are obtained by the same formulas as those in regime I not only in the UV (as argued at the beginning) but also in the IR—this is detailed in the Appendix.

Finally, contrary to what we expect, we find that analytic continuation does not seem to work for out-of-equilibrium quantities. To test the idea, we numerically expand two Loschmidt echoes and compare coefficients of Taylor series, one with $\gamma \in \mathbb{R}$ (Hermitian, phase I) and another one with $\gamma \in i\mathbb{R}$ (pure imaginary coupling, phase III) in terms of $\gamma$. We expect that, for analytic continuation to work, the expansion

<figure><img src="image_11.png" /><figcaption>FIG. 8. Taylor coefficient $a_2$, $b_2$ and $a_4$, $b_4$ as a function of time for phase I (real coupling $\gamma = J$) and phase III (imaginary coupling $\gamma = iJ$). Lower order coefficients are in agreement with analytical continuation picture [panel (a)], whereas the fourth order coefficient $a_4$ behaves differently for real and imaginary couplings.</figcaption></figure>

of the latter case be obtained by performing a transformation $\gamma = J \mapsto iJ$ to the former one. We study the real and imaginary parts of the Loschmidt echo $\tilde{\mathcal{L}}$ separately, each of which is given by

$$ \begin{aligned} \mathrm{Re}\tilde{\mathcal{L}}(t) &= a_0 + a_1\gamma + a_2\gamma^2 + a_3\gamma^3 + a_4\gamma^4 + \mathcal{O}(\gamma^5) \\ \mathrm{Im}\tilde{\mathcal{L}}(t) &= b_0 + b_1\gamma + b_2\gamma^2 + b_3\gamma^3 + b_4\gamma^4 + \mathcal{O}(\gamma^5). \end{aligned} \quad (62) $$

The coefficient $a_0 = 1$ corresponds to the case $\gamma = 0$, i.e., the initial and the time evolved states are equal, so $\langle\psi_0|\psi(t)\rangle = 1$. Note also that coefficients corresponding to the odd powers of $\gamma$ are zeros $a_1 = a_3 = b_0 = b_1 = b_3 = 0$, indicating that the Loschmidt echo is a function of $\tilde{\Gamma}$ in both phase I and III at least in the scaling regime. Nontrivial coefficients appear with the even powers of $\gamma$, and are presented in Fig. 8. We first notice that $a_2(\gamma = J) = -a_2(\gamma = iJ)$ and $b_2(\gamma = J) = -b_2(\gamma = iJ)$ always hold at any time, which indicates that the analytic continuation works at this order. The next order, the coefficients of $\gamma^4$, however, clearly invalidate the continuation. Namely, in Fig. 8(b), we observe a firm agreement for imaginary parts of the Loschmidt echo in phases I and III, while the curves representing the real parts start departing as soon as $t > 0$. For now we have no explanation as to why the continuation works for the imaginary part only and not for the real part. We also confirm that the discrepancy cannot be removed by increasing the system size.

<footnote n="2">²There are variants of the model which are not PT symmetric, such as the RLM with a complex coupling $\gamma$. In this case, $d$ can be expanded in powers of a variable similar to $u$ and proportional to $\gamma^2$. For $\gamma = Je^{i\psi}$, the transition between the equivalent of regime I and regime III occurs when $\arg\gamma^2 = \pi/2$, i.e., $\varphi = \pi/4$. This will be discussed elsewhere.</footnote>
## VII. CONCLUSIONS

We find that the behaviors of the magnetization and free energy in the PT-broken regime are markedly different from those in the case where PT symmetry is unbroken. Our results suggest that, when PT symmetry is spontaneously broken, the impurity is not screened in the infrared. This lack of Kondo screening is further suggested by a jump of the magnetization associated with the PT phase transition at half filling, as in the non-Hermitian Kondo model studied in Ref. [10].

Meanwhile, our study of the Loschmidt echo shows that the overlap at large times decays as $t^{-1/2}$ whatever the regime. This decay is usually associated with healing, i.e., a situation where the two wires are fully connected [30]. In other words, we seem to find that, at low energy, the impurity is not screened and the wires are healed. We do not really understand how this can be made compatible with the expected circular behavior of the RG explored in particular in Refs. [12–14]. Of course, this circular behavior is under control only in the perturbative regime (close to the isotropic Kondo problem), and it may be that it disappears in the RLM limit: We hope to get back to this question soon. In particular, we find the fact that equilibrium quantities in the PT-broken regime can be calculated by analytical continuation from those in the unbroken regime very encouraging, since a lot is known about (UV) perturbative expansions of physical quantities in the (anisotropic) Kondo problem and boundary sine-Gordon problem.

We also studied the Loschmidt echo starting from the initially disconnected chains. This requires further study: To our knowledge, it was not even known that this echo would obey universal behavior in the different regimes, with different universal curves.

Finally, we note that the quench protocol we examined in this paper was in fact already studied in the context of quantum impurity systems [36,37]. It therefore might also be possible to study a local quench in PTRLM experimentally by performing a quench to a PT-symmetric optical trap.

## ACKNOWLEDGMENTS

T.Y. thanks B. Doyon and R. Pereira, and O. Shpielberg for useful discussions. T.Y. acknowledges the financial support from Takenaka Scholarship Foundation, and hospitality at the Institut de Physique Théorique, CEA, where the present work was initiated. K.B. is grateful to G. Misguich, A. Michailidis, and O. Gamayun for valuable suggestions. K.B. acknowledges the support from the European Research Council under the Starting Grant No. 805252 LoCoMacro. H.S. wishes to thank N. Andrei for discussions. This research was supported in part by the ERC Advanced Grant NuQFT.

## APPENDIX A: NON-HERMITIAN QUANTUM MECHANICS AND BIORTHOGONAL BASIS

In general, when a given Hamiltonian $H$ is not Hermitian, its eigenstates are not orthogonal. This in turn results in the loss of ordinary properties that Hermitian quantum systems have, such as the orthogonality of wave functions and positivity of inner products. This problem can be in fact

circumvented by properly extending the notion of basis to the biorthogonal basis $\{|R_n\rangle, |L_n\rangle\} \in \mathcal{H} \otimes \mathcal{H}^*$, where $\mathcal{H}$ and $\mathcal{H}^*$ are the Hilbert space and its dual space, respectively, satisfying

$$H|R_n\rangle = E_n|R_n\rangle, \quad H^\dagger|L_n\rangle = E_n^*|L_n\rangle. \quad (A1)$$

Biorthogonality further entails

$$\langle R_n|L_m\rangle = \delta_{n,m}, \quad \sum_n |R_n\rangle\langle L_n| = \mathbb{I}. \quad (A2)$$

The trace of any observable over $\mathcal{H}$ can be expanded with respect to the biorthogonal basis. This can be best seen by first expanding $\text{Tr} \mathcal{O}$ with respect to a generic orthonormal basis $|\xi\rangle$ (which is always guaranteed to exist in a Hilbert space) as $\text{Tr} \mathcal{O} = \sum_{\xi} \langle\xi|\mathcal{O}|\xi\rangle$. Inserting the spectral decomposition $\sum_n |R_n\rangle\langle L_n| = \mathbb{I}$, we then have

$$\text{Tr} \mathcal{O} = \sum_{\xi} \sum_{n} \langle L_n | \mathcal{O} | \xi \rangle \langle \xi | R_n \rangle = \sum_{n} \langle L_n | \mathcal{O} | R_n \rangle. \quad (A3)$$

Now suppose that there is a parity operator $\mathcal{P}$ such that

$$H^\dagger = \mathcal{P}\mathcal{H}\mathcal{P}, \quad (A4)$$

and

$$\mathcal{P}|R_n\rangle = p_n|L_n\rangle \quad (A5)$$

with $\mathcal{P}^2 = \mathbb{I}$ and $p_n = \pm 1$. It is then customary to define an operator called $\mathcal{C}$ operator

$$C = \sum_n p_n |R_n\rangle\langle L_n| \quad (A6)$$

so that $\rho = \mathcal{P}\mathcal{C} = \sum_n |L_n\rangle\langle L_n|$ (not to be confused with density operator) which relates $H$ and $H^\dagger$ as

$$H^\dagger \rho = \rho H. \quad (A7)$$

With this operator, we can provide an alternative way of interpreting the biorthogonal basis. Namely, the use of the biorthogonal basis is equivalent to working with the Hilbert space endowed with the metric $\rho$. For a given state $|\alpha\rangle = \sum_n \alpha_n |R_n\rangle \in \mathcal{H}$, let us define the associate state as

$$\langle \tilde{\alpha} | = \sum_n \alpha_n^* \langle L_n | \in \mathcal{H}^*. \quad (A8)$$

The inner product of two states $|\alpha\rangle$ and $|\beta\rangle$ for a biorthogonal system is then defined by

$$\langle \beta | \alpha \rangle_{\text{PT}} := \langle \tilde{\beta} | \alpha \rangle = \sum_n \beta_n^* \alpha_n. \quad (A9)$$

This inner product can be viewed as a conventional quantum mechanical inner product with the metric $\rho$ defined as follows

$$\langle \beta | \alpha \rangle_{\rho} := \langle \beta | \rho | \alpha \rangle = \sum_n \beta_n^* \alpha_n, \quad (A10)$$

hence the equivalence of two point of views. It is readily seen that this inner product is positive definite having discussed the states in a biorthogonal basis, now let us turn to observables. In terms of a given biorthogonal basis $\{|R_n\rangle, |L_n\rangle\}$, we can expand a generic operator $A$ as

$$A = \sum_{n,m} a_{nm} |R_n\rangle\langle L_m|. \quad (A11)$$
The expectation value of A in a pure state $|\psi\rangle = \sum_n c_n|R_n\rangle$, which we denote $\langle A\rangle_\psi$, is then defined by

$$ \langle A\rangle_\psi = \frac{\langle\tilde{\psi}|A|\psi\rangle}{\langle\tilde{\psi}|\psi\rangle} = \sum_{n,m} \frac{c_n^* c_m a_{nm}}{\sum_n c_n^* c_n}, \quad (A12) $$

where we recall that $\langle\tilde{\psi}| = \sum_n c_n^* \langle L_n|$. Likewise, an arbitrary density operator $\rho$ can be decomposed as

$$ \rho = \sum_{n,m} \rho_{nm} |R_n\rangle\langle L_m|. \quad (A13) $$

The statistical average of an observable A with respect to $\rho$ then reads

$$ \langle A\rangle = \mathrm{Tr}(\rho A) = \sum_n \langle L_n|\rho A|R_n\rangle = \sum_{n,m} a_{nm} \rho_{mn}. \quad (A14) $$

Next, we turn to time evolution. A pure state $|\psi\rangle$ is expected to time evolve according to the following Schrodinger equation as in Hermitian systems:

$$ i \frac{\partial}{\partial t} |\psi\rangle = H |\psi\rangle. \quad (A15) $$

Suppose that we start with an initial pure state $|\psi\rangle_0 = \sum_n c_n|R_n\rangle$. The time-evolved state $|\psi(t)\rangle$ is given by

$$ |\psi(t)\rangle = \sum_n c_n e^{-iE_n t} |R_n\rangle \quad (A16) $$

with the corresponding bra state

$$ \langle\tilde{\psi}(t)| = \sum_n c_n^* e^{iE_n^* t} \langle L_n|. \quad (A17) $$

The biorthogonal norm of the state $|\psi(t)\rangle$ then reads

$$ \langle\tilde{\psi}(t)|\psi(t)\rangle = \sum_n c_n^* c_n e^{-i(E_n - E_n^*)t}. \quad (A18) $$

We observe that the norm is time independent when the spectrum is entirely real, i.e., $E_n \in \mathbb{R}$ for all $n$. Furthermore, if the spectrum is real, then the associated state $|\tilde{\psi}(t)\rangle = \sum_n c_n e^{-iE_n t} |L_n\rangle$, which is the conjugate of (A17), satisfies another Schrodinger equation $i \frac{\partial}{\partial t} |\tilde{\psi}\rangle = H^\dagger |\tilde{\psi}\rangle$, implying that the Heisenberg picture is still valid when PT symmetry is unbroken:

$$ \langle\tilde{\psi}(t)|A|\psi(t)\rangle = \langle\tilde{\psi}|e^{iHt}Ae^{-iHt}|\psi\rangle = \langle\tilde{\psi}|A(t)|\psi\rangle \quad (A19) $$

with $A(t) = e^{iHt}Ae^{-iHt}$. We however note that the validity of the Heisenberg picture breaks down when at least one eigenvalue is complex, as the associated state $|\tilde{\psi}(t)\rangle$ no longer follows the Schrodinger equation, i.e., $i \frac{\partial}{\partial t} |\tilde{\psi}\rangle \neq H^\dagger |\tilde{\psi}\rangle$ [22].

# APPENDIX B: DETAILS ON THE DOT DENSITY

## 1. Dot density at $T = 0$

### a. Phase I

We first note that $k_F = \arccos(-\frac{\mu}{2}) \ge \frac{\pi}{2}$. Then we can safely decompose the integral as

$$ \begin{aligned} d &= \int_0^{k_F} \frac{dk}{\pi} \frac{\tilde{\Gamma} \sin^2 k}{\tilde{\Gamma}^2 + (1-2\tilde{\Gamma})\cos^2 k} = \int_0^{\pi/2} \frac{dk}{\pi} \frac{\tilde{\Gamma} \sin^2 k}{\tilde{\Gamma}^2 + (1-2\tilde{\Gamma})\cos^2 k} + \lim_{\epsilon \to 0^+} \int_{\pi/2+\epsilon}^{k_F} \frac{dk}{\pi} \frac{\tilde{\Gamma} \sin^2 k}{\tilde{\Gamma}^2 + (1-2\tilde{\Gamma})\cos^2 k} \\ &= \frac{1}{2} + \frac{\tilde{\Gamma}}{1-2\tilde{\Gamma}} \int_{-\infty}^{\tan k_F} \frac{dx}{\pi} \left( \frac{1}{1+\frac{\tilde{\Gamma}^2}{(1-\tilde{\Gamma})^2}x^2} - \frac{1}{1+x^2} \right) \\ &= \frac{1}{2} + \frac{\tilde{\Gamma}}{1-2\tilde{\Gamma}} \left( \int_0^\infty \frac{dx}{\pi} - \int_0^{-\tan k_F} \frac{dx}{\pi} \right) \left( \frac{1}{1+\frac{\tilde{\Gamma}^2}{(1-\tilde{\Gamma})^2}x^2} - \frac{1}{1+x^2} \right) \\ &= \frac{1}{2} + \frac{1-\tilde{\Gamma}}{1-2\tilde{\Gamma}} \int_{u(\tilde{\Gamma})}^\infty \frac{dx}{\pi} \frac{1}{1+x^2} - \frac{\tilde{\Gamma}}{1-2\tilde{\Gamma}} \int_{-\tan k_F}^\infty \frac{dx}{\pi} \frac{1}{1+x^2} \\ &= 1 + \frac{1-\tilde{\Gamma}}{\pi u(\tilde{\Gamma})(1-2\tilde{\Gamma})} \sum_{n=1}^\infty \frac{(-1)^n}{2n-1} u(\tilde{\Gamma})^{2n} + \frac{\tilde{\Gamma}}{\pi \tan k_F (1-2\tilde{\Gamma})} \sum_{n=1}^\infty \frac{(-1)^n}{2n-1} \tan^{2n} k_F \quad (B1) \end{aligned} $$

where we remarked $\tan k_F < 0$ and defined $u(\tilde{\Gamma}) = -\tilde{\Gamma} \tan k_F/(1-\tilde{\Gamma}) = \tilde{\Gamma} \sqrt{4-\mu^2}/(\mu(1-\tilde{\Gamma}))$. We also assumed that $u(\tilde{\Gamma}) < 1$ for the time being. The point here is that we cannot naively change the integration variable $x = \tan k$ in the second equality as the integration region $[0, k_F]$ does not one-to-one correspond to $[0, \tan k_F]$. We then have the UV expansion

$$ d = 1 - \frac{1}{\pi(1-2\tilde{\Gamma})} \left[ (1-\tilde{\Gamma}) \arctan u(\tilde{\Gamma}) - \tilde{\Gamma} \arctan \frac{\sqrt{4-\mu^2}}{\mu} \right], \quad (B2) $$

where the function arctan is defined by its series expansion (this point has to be made clear since multiple definitions of the function arctan $z = \frac{1}{2i} \ln \frac{1+i\epsilon}{1-i\epsilon}$ play a role in the discussion of analytical continuations in the body of this paper):

$$ \arctan x = \sum_{n=0}^{\infty} (-1)^{2k+1} \frac{x^{2k+1}}{2k+1}, \quad |x| < 1. \quad (B3) $$
For the opposite case $u(\tilde{\Gamma}) > 1$, one can perform the dual transformation $x \mapsto 1/x$ in (B1) and get the IR expansion

$$d = \frac{1}{2} + \frac{1 - \tilde{\Gamma}}{1 - 2\tilde{\Gamma}} \int_0^{1/u(\tilde{\Gamma})} \frac{dx}{\pi} \frac{1}{1+x^2} - \frac{\tilde{\Gamma}}{1 - 2\tilde{\Gamma}} \int_0^{-1/\tan k_F} \frac{dx}{\pi} \frac{1}{1+x^2} \\ = \frac{1}{2} - \frac{u(\tilde{\Gamma})(1 - \tilde{\Gamma})}{\pi(1 - 2\tilde{\Gamma})} \sum_{n=1}^{\infty} \frac{(-1)^n}{2n-1} u(\tilde{\Gamma})^{-2n} - \frac{\tan k_F \tilde{\Gamma}}{\pi(1 - 2\tilde{\Gamma})} \sum_{n=1}^{\infty} \frac{(-1)^n}{2n-1} \left( \frac{1}{\tan k_F} \right)^{2n} \\ = \frac{1}{2} + \frac{1}{\pi(1 - 2\tilde{\Gamma})} \left[ (1 - \tilde{\Gamma}) \arctan \frac{1}{u(\tilde{\Gamma})} - \tilde{\Gamma} \arctan \frac{\mu}{\sqrt{4 - \mu^2}} \right]. \qquad (B4)$$

Therefore, in summary,

$$d = \begin{cases} 1 - \frac{1}{\pi(1-2\tilde{\Gamma})} \left[ (1-\tilde{\Gamma}) \arctan u(\tilde{\Gamma}) - \tilde{\Gamma} \arctan \frac{\sqrt{4-\mu^2}}{\mu} \right] & u(\tilde{\Gamma}) < 1 \\ \frac{1}{2} + \frac{1}{\pi(1-2\tilde{\Gamma})} \left[ (1-\tilde{\Gamma}) \arctan \frac{1}{u(\tilde{\Gamma})} - \gamma_1^2 \arctan \frac{\mu}{\sqrt{4-\mu^2}} \right] & u(\tilde{\Gamma}) > 1 \end{cases}, \qquad (B5)$$

and in particular, under the scaling limit $\tilde{\Gamma}, \mu \ll 1$, we have $u(\tilde{\Gamma}) \simeq \tilde{\Gamma}/\mu$, and

$$d = \begin{cases} 1 - \frac{1}{\pi} \arctan u(\tilde{\Gamma}) & u(\tilde{\Gamma}) < 1 \\ \frac{1}{2} + \frac{1}{\pi} \arctan \left( \frac{1}{u(\tilde{\Gamma})} \right) & u(\tilde{\Gamma}) > 1 \end{cases}. \qquad (B6)$$

b. Phase II

We can calculate the dot density in phase II in the same way for $\mu > 0$.

$$\begin{aligned} d &= \int_0^{k_F} \frac{dk}{\pi} \frac{\tilde{\Gamma} \sin^2 k}{\tilde{\Gamma}^2 + (1-2\tilde{\Gamma})\cos^2 k} + \frac{1-\tilde{\Gamma}}{1-2\tilde{\Gamma}} \\ &= \int_0^{\pi/2} \frac{dk}{\pi} \frac{\tilde{\Gamma} \sin^2 k}{\tilde{\Gamma}^2 + (1-2\tilde{\Gamma})\cos^2 k} + \frac{1-\tilde{\Gamma}}{1-2\tilde{\Gamma}} + \lim_{\epsilon \to 0^+} \int_{\pi/2+\epsilon}^{k_F} \frac{dk}{\pi} \frac{\tilde{\Gamma} \sin^2 k}{\tilde{\Gamma}^2 + (1-2\tilde{\Gamma})\cos^2 k} \\ &= \frac{1}{2} + \lim_{\epsilon \to 0^+} \int_{\pi/2+\epsilon}^{k_F} \frac{dk}{\pi} \frac{\tilde{\Gamma} \sin^2 k}{\tilde{\Gamma}^2 + (1-2\tilde{\Gamma})\cos^2 k} \\ &= \frac{1}{2} - \frac{1-\tilde{\Gamma}}{1-2\tilde{\Gamma}} \int_{v(\tilde{\Gamma})}^{\infty} \frac{dx}{\pi} \frac{1}{1+x^2} - \frac{\tilde{\Gamma}}{1-2\tilde{\Gamma}} \int_{-\tan k_F}^{\infty} \frac{dx}{\pi} \frac{1}{1+x^2} \\ &= \frac{\tilde{\Gamma}}{2\tilde{\Gamma}-1} - \frac{1-\tilde{\Gamma}}{\pi v(\tilde{\Gamma})(1-2\tilde{\Gamma})} \sum_{n=1}^{\infty} \frac{(-1)^n}{2n-1} v(\tilde{\Gamma})^{2n} - \frac{\tilde{\Gamma}}{\pi \tan k_F (1-2\tilde{\Gamma})} \sum_{n=1}^{\infty} \frac{(-1)^n}{2n-1} \tan^{2n} k_F, \end{aligned} \qquad (B7)$$

where we assumed $v(\tilde{\Gamma}) := -u(\tilde{\Gamma}) < 1$. Note that when passing from the second line to the third line, we used

$$\int_0^{\pi/2} \frac{dk}{\pi} \frac{\tilde{\Gamma} \sin^2 k}{\tilde{\Gamma} + (1 - 2\tilde{\Gamma}) \cos^2 k} + \frac{1 - \tilde{\Gamma}}{1 - 2\tilde{\Gamma}} = \frac{1}{2}. \qquad (B8)$$

The UV expansion in phase II then reads

$$d = \frac{\tilde{\Gamma}}{2\tilde{\Gamma} - 1} + \frac{1}{\pi(2\tilde{\Gamma} - 1)} \left[ (\tilde{\Gamma} - 1) \arctan v(\tilde{\Gamma}) - \tilde{\Gamma} \arctan \frac{\sqrt{4 - \mu^2}}{\mu} \right]. \qquad (B9)$$

As in phase I, the dual transformation in the integrations instead gives the IR expansion

$$\begin{aligned} d &= \frac{1}{2} + \frac{v(\tilde{\Gamma})(\tilde{\Gamma} - 1)}{\pi(2\tilde{\Gamma} - 1)} \sum_{n=1}^{\infty} \frac{(-1)^n}{2n-1} v(\tilde{\Gamma})^{-2n} + \frac{\tilde{\Gamma}}{\pi(2\tilde{\Gamma} - 1)} \arctan \frac{\mu}{\sqrt{4 - \mu^2}} \\ &= \frac{1}{2} - \frac{1}{\pi(2\tilde{\Gamma} - 1)} \left[ (\tilde{\Gamma} - 1) \arctan \frac{1}{v(\tilde{\Gamma})} - \tilde{\Gamma} \arctan \frac{\mu}{\sqrt{4 - \mu^2}} \right]. \end{aligned} \qquad (B10)$$
### c. Phase III

A similar analysis can be carried out in phase III and for $\mu > 0$

$$ 
\begin{aligned}
d &= \int_0^{k_F} \frac{dk}{\pi} \frac{\tilde{\Gamma} \sin^2 k}{\tilde{\Gamma}^2 + (1-2\tilde{\Gamma})\cos^2 k} + \frac{2-2\tilde{\Gamma}}{1-2\tilde{\Gamma}} \\
&= \int_0^{\pi/2} \frac{dk}{\pi} \frac{\tilde{\Gamma} \sin^2 k}{\tilde{\Gamma} + (1-2\tilde{\Gamma})\cos^2 k} + \frac{2-2\tilde{\Gamma}}{1-2\tilde{\Gamma}} + \lim_{\epsilon \to 0^+} \int_{\pi/2+\epsilon}^{k_F} \frac{dk}{\pi} \frac{\tilde{\Gamma} \sin^2 k}{\tilde{\Gamma}^2 + (1-2\tilde{\Gamma})\cos^2 k} \\
&= \frac{1}{2} + \frac{1-\tilde{\Gamma}}{1-2\tilde{\Gamma}} + \lim_{\epsilon \to 0^+} \int_{\pi/2+\epsilon}^{k_F} \frac{dk}{\pi} \frac{\tilde{\Gamma} \sin^2 k}{\tilde{\Gamma}^2 + (1-2\tilde{\Gamma})\cos^2 k} \\
&= \frac{1}{2} + \frac{1-\tilde{\Gamma}}{1-2\tilde{\Gamma}} - \frac{1-\tilde{\Gamma}}{1-2\tilde{\Gamma}} \int_{-u(\tilde{\Gamma})}^{\infty} \frac{dx}{\pi} \frac{1}{1+x^2} - \frac{\tilde{\Gamma}}{1-2\tilde{\Gamma}} \int_{-\tan k_F}^{\infty} \frac{dx}{\pi} \frac{1}{1+x^2} \\
&= 1 - \frac{1-\tilde{\Gamma}}{\pi v(\tilde{\Gamma})(1-2\tilde{\Gamma})} \sum_{n=1}^{\infty} \frac{(-1)^n}{2n-1} v(\tilde{\Gamma})^{2n} - \frac{\tilde{\Gamma}}{\pi \tan k_F(1-2\tilde{\Gamma})} \sum_{n=1}^{\infty} \frac{(-1)^n}{2n-1} \tan^{2n} k_F,
\end{aligned} 
\quad (B11) $$

where we again used (B8) in the third line, which is still valid in phase III. We therefore have the UV expansion in phase III

$$ d = 1 + \frac{1}{\pi(1 - 2\tilde{\Gamma})} \left[ (1 - \tilde{\Gamma}) \arctan v(\tilde{\Gamma}) + \tilde{\Gamma} \arctan \frac{\sqrt{4 - \mu^2}}{\mu} \right]. \quad (B12) $$

Likewise, the IR expansion ($v(\tilde{\Gamma}) > 1$) reads

$$ 
\begin{aligned}
d &= \frac{1}{2} + \frac{1 - \tilde{\Gamma}}{1 - 2\tilde{\Gamma}} - \frac{1 - \tilde{\Gamma}}{1 - 2\tilde{\Gamma}} \int_0^{1/v(\tilde{\Gamma})} \frac{dx}{\pi} \frac{1}{1 + x^2} - \frac{\tilde{\Gamma}}{1 - 2\tilde{\Gamma}} \int_0^{-1/\tan k_F} \frac{dx}{\pi} \frac{1}{1 + x^2} \\
&= \frac{1}{2} + \frac{1 - \tilde{\Gamma}}{1 - 2\tilde{\Gamma}} - \frac{1}{\pi(1 - 2\tilde{\Gamma})} \left[ (1 - \tilde{\Gamma}) \arctan \frac{1}{v(\tilde{\Gamma})} + \tilde{\Gamma} \arctan \frac{\mu}{\sqrt{4 - \mu^2}} \right].
\end{aligned} 
\quad (B13) $$

### 2. Proof of $d = \frac{1}{2}$ at half filling

Defining $D(k) = \tilde{\Gamma} \sin^2 k/(\tilde{\Gamma}^2 + (1-2\tilde{\Gamma})\cos^2 k)$, we rewrite (33) as

$$ d = \int_{0}^{\pi/2} \frac{dk}{\pi} \frac{1}{1 + e^{-2\beta \cos k}} D(k) + \int_{\pi/2}^{\pi} \frac{dk}{\pi} \frac{e^{2\beta \cos k}}{1 + e^{2\beta \cos k}} D(k). \quad (B14) $$

In each integration region, the Fermi distribution can be expanded as

$$ 
\begin{aligned}
\frac{1}{1 + e^{-2\beta \cos k}} &= \sum_{n=0}^{\infty} (-e^{-2\beta \cos k})^n = \sum_{n=0}^{\infty} (-1)^n \sum_{m=0}^{\infty} \frac{(-2n\beta)^m}{m!} \cos^m k \\
\frac{e^{2\beta \cos k}}{1 + e^{2\beta \cos k}} &= \sum_{n=0}^{\infty} (-1)^n e^{2\beta(n+1)\cos k} = \sum_{n=0}^{\infty} (-1)^n \sum_{m=0}^{\infty} \frac{[2(n+1)\beta]^m}{m!} \cos^m k.
\end{aligned} 
\quad (B15) $$

As we will see later, it is useful to divide the case $0 < \tilde{\Gamma} \le 1/2$ and $1/2 \le \tilde{\Gamma} < 1$ and analyze separately. Let us deal with the former case first in which case we can expand $D(k)$ as

$$ D(k) = \frac{1}{(\tilde{\Gamma} - 1)^2 + (2\tilde{\Gamma} - 1)\sin^2 k} = \frac{1}{(\tilde{\Gamma} - 1)^2} \sum_{l=0}^{\infty} \left( -\frac{2\tilde{\Gamma} - 1}{(\tilde{\Gamma} - 1)^2} \right)^l \sin^{2l} k, \quad (B16) $$

where we noted that $-1 < -(2\tilde{\Gamma} - 1)/(\tilde{\Gamma} - 1)^2 \le 0$. We then perform the integrations each of which reads

$$ 
\begin{aligned}
\int_0^{\pi/2} dk \cos^m k \sin^{2(l+1)} k &= \frac{1}{2} B\left(\frac{m+1}{2}, \frac{2l+3}{2}\right) = \frac{1}{2} \frac{\Gamma\left(\frac{m+1}{2}\right) \Gamma\left(l + \frac{3}{2}\right)}{\Gamma\left(l + \frac{m}{2} + 2\right)} \\
\int_{\pi/2}^{\pi} dk \cos^m k \sin^{2(l+1)} k &= \frac{(-1)^m}{2} \frac{\Gamma\left(\frac{m+1}{2}\right) \Gamma\left(l + \frac{3}{2}\right)}{\Gamma\left(l + \frac{m}{2} + 2\right)}.
\end{aligned} 
\quad (B17) $$
Using them, we can express the dot density d as

$$d = \frac{\tilde{\Gamma}}{2\pi(\tilde{\Gamma}-1)^2} \sum_{n=0}^{\infty} (-1)^n \sum_{m=0}^{\infty} \frac{(-2n\beta)^m + (-2(n+1)\beta)^m}{m!} \sum_{l=0}^{\infty} \left(-\frac{2\tilde{\Gamma}-1}{(\tilde{\Gamma}-1)^2}\right)^l \frac{\Gamma\left(\frac{m+1}{2}\right)\Gamma\left(l+\frac{3}{2}\right)}{\Gamma\left(l+\frac{m}{2}+2\right)} \\ = \frac{\tilde{\Gamma}}{2\pi(\tilde{\Gamma}-1)^2} \sum_{l=0}^{\infty} \left(-\frac{2\tilde{\Gamma}-1}{(\tilde{\Gamma}-1)^2}\right)^l \frac{\Gamma\left(\frac{1}{2}\right)\Gamma\left(l+\frac{1}{2}\right)}{\Gamma(l+2)} = \frac{\tilde{\Gamma}}{4(\tilde{\Gamma}-1)^2} {}_2F_1\left(1, \frac{3}{2}; 2; -\frac{2\tilde{\Gamma}-1}{(\tilde{\Gamma}-1)^2}\right). \quad (B18)$$

Since ${}_2F_1(1, \frac{3}{2}; 2; z) = (\frac{1}{2} + \frac{1}{2}\sqrt{1-z})^{-1}/\sqrt{1-z}$ (note the positive branch of square root), it follows that $d = \frac{1}{2}$. For the other case $1/2 \le \tilde{\Gamma} < 1$, noticing that $-1 < (1-2\tilde{\Gamma})/\tilde{\Gamma}^2 \le 0$, we can use another expansion

$$D(k) = \frac{1}{\tilde{\Gamma}^2} \sum_{l=0}^{\infty} \left( \frac{1 - 2\tilde{\Gamma}}{\tilde{\Gamma}^2} \right)^l \cos^{2l} k, \quad (B19)$$

yielding the same result.

## APPENDIX C: LOSCHMIDT ECHO

Here we elucidate the computations of the Loschmidt echo, which is defined by

$$\mathcal{L}(t) = {}_0\langle \text{GS} | e^{itH_0} e^{-itH} | \text{GS} \rangle_0 = e^{itE_{\text{GS}}} {}_0\langle \text{GS} | e^{-itH} | \text{GS} \rangle_0, \quad (C1)$$

where $|\text{GS}\rangle_0$ denotes the ground state of two disconnected chains. We consider the case where an arbitrary chemical potential is imposed. Let us first expand $|\text{GS}\rangle_0$ in terms of the mode operator in the PTRLM. Defining the cutoff mode $M = \lfloor \frac{N+2}{2\pi} k_F \rfloor$ with the Fermi momentum $k_F = \arccos(-\frac{\mu}{2})$,

$$|\text{GS}\rangle_0 = \prod_{i \in \text{FS}} c_{k_i}^{\dagger} |0\rangle = \prod_{i \in \text{FS}} \sum_{x_i} \phi(k_i; x_i) c_{x_i}^{\dagger} |0\rangle = \prod_{i \in \text{FS}} \sum_{x_i} \sum_{l_i} \phi(k_i; x_i) L^*(l_i; x_i) c_{R, l_i}^{\dagger} |0\rangle = \prod_{i \in \text{FS}} \sum_{l_i} L^*(k_i, l_i) c_{R, l_i}^{\dagger} |0\rangle, \quad (C2)$$

where we defined $L^*(k_i, l_i) := \sum_x \phi(k_i; x)L^*(l_i; x)$. In the same way, we also expand ${}_0\langle \text{GS} |$ as

$${}_0\langle \text{GS} | = \langle 0 | \prod_{i \in \text{FS}} \sum_{l_i} \mathbb{R}(k_i, l_i) c_{L, l_i}. \quad (C3)$$

Now, inserting the biorthogonal resolution of identity $1 = \sum_K \frac{1}{K!} \prod_{j=1}^K \sum_{m_j} c_{R,m_1}^\dagger \cdots c_{R,m_K}^\dagger |0\rangle \langle 0| c_{L,m_K} \cdots c_{L,m_1}$ into (C1), we obtain

$$\mathcal{L}(t) = \frac{e^{itE_{\text{GS}}}}{(2M+1)!} \prod_{i \in \text{FS}} \sum_{m_i} e^{-iE_{m_i}t} \det_{i,j \in \text{FS}} L^*(k_i, m_j) \det_{i',j' \in \text{FS}} \mathbb{R}(k_{i'}, m_{j'}) = e^{itE_{\text{GS}}} \det_{i,j \in \text{FS}} \mathbb{M}(k_i; k_j), \quad (C4)$$

where $\mathbb{M}(k_i, k_j; t) = \sum_m e^{-iE_{m}t} L^*(k_i, m) \mathbb{R}(k_j, m)$. To proceed, we divide $\mathbb{M}(k_i, k_j; t)$ into $\mathbb{M}(k_i, k_j) = \mathbb{M}_+(k_i, k_j) + \mathbb{M}_-(k_i, k_j)$ each of which is defined as $\mathbb{M}_\pm(k_i, k_j; t) = \sum_{m_\pm} e^{-iE_{m_\pm}t} L^*(k_i, m_\pm) \mathbb{R}(k_j, m_\pm)$.

### 1. Phase I

Let us start with $\mathbb{M}_+(k_i, k_j; t)$. Note that $L^*(k, m_+) $ and $\mathbb{R}(k, m_+) $ can be explicitly evaluated as

$$\begin{aligned} L^*(k; m_+) &= \sum_x \phi(k; x) L_+^*(m_+; x) = 2\mathcal{N}_{L,+}^*(m_+) \sqrt{\frac{1}{N+2}} \sum_{x>0} (\vartheta(k) + \mathcal{M}_+\vartheta(-k)) \sin kx \cos[m_+(x - \delta_+(m_+))] \\ &= \frac{\mathcal{N}_{L,+}^*(m_+)}{\sqrt{N+2}} (\vartheta(k) + \mathcal{M}_+\vartheta(-k)) \frac{\sin k \cos[m_+\delta_+(m_+)]}{\cos m_+ - \cos k}, \end{aligned} \quad (C5)$$

and

$$\mathbb{R}(k; m_+) = \frac{\mathcal{N}_{R,+}(m_+)}{\sqrt{N+2}} (\vartheta(k) + \mathcal{M}_+\vartheta(-k)) \frac{\sin k \cos[m_+\delta_+(m_+)]}{\cos m_+ - \cos k}. \quad (C6)$$

Recalling that $\cos^2[m_+\delta_+(m_+)] = 4\tilde{\Gamma}^2 \sin^2 m/\Delta_+(m)$, $\mathbb{M}_+(k_i, k_j; t)$ can therefore be expressed as

$$\mathbb{M}_+(k_i, k_j; t) = \frac{4\tilde{\Gamma}^2 \sin k_i \sin k_j}{N+2} (\vartheta(k_i) + \mathcal{M}_+\vartheta(-k_i))(\vartheta(k_j) + \mathcal{M}_+\vartheta(-k_j)) \sum_{m_+} \frac{\mathcal{N}_+(m_+) \sin^2 m_+ e^{-itE_{m_+}}}{\Delta_+(m_+)(\cos m_+ - \cos k_i)(\cos m_+ - \cos k_j)}. \quad (C7)$$
We focus on computing the building block

$$Q_+(k_i, k_j; t) := \sum_{m_+ > 0} w_1(m_+; k_i, k_j; t), \quad w_1(m_+; k_i, k_j; t) = \frac{\mathcal{N}_+(m_+) \sin^2 m_+ e^{-itE_{m_+}}}{\Delta_+(m_+)(\cos m_+ - \cos k_i)(\cos m_+ - \cos k_j)}. \quad (C8)$$

In the thermodynamic limit, the calculation of (C8) is apparently plagued by poles, thereby necessitating a systematic way to deal with them. To this end, we introduce a counting function $Q_σ(k)$ that satisfies

$$Q_{\sigma}(k) = kN + i \log S_{\sigma}(k), \quad (C9)$$

where the $S$ matrix is $S_σ(k) = (\Gamma_σ - e^{-2ik})/(\Gamma_σ - e^{2ik})$. Only $k$ such that $Q_σ(k) = 2\pi n$ for some quantum number $n \in \mathbb{Z}$ is allowed by the quantization condition (12). Using $Q(k)$, we can rewrite the sum over $k$ (suppose that $k$ satisfies the quantization condition with respect to the sign $\sigma$) as

$$\sum_k \rightarrow \sum_n \oint_{C_n} \frac{dk}{2\pi} \frac{-Q'_\sigma(k)}{e^{-iQ_\sigma(k)} - 1}, \quad (C10)$$

where $C_n$ is a closed contour encircling each pole. We then rewrite $Q_+(k_i, k_j; t)$ as

$$Q_+(k_i, k_j; t) = \left( \int_{0-i\epsilon}^{\pi-i\epsilon} - \int_{0+i\epsilon}^{\pi+i\epsilon} \right) \frac{dm}{2\pi} \frac{-Q'_+(m)}{e^{-iQ_+(m)} - 1} w_1(m; k_i, k_j; t) - i \sum_l \mathrm{Res}|_{m_+=\bar{k}_i} \frac{-Q'_+(m_+)}{e^{-iQ_+(m_+)} - 1} w_1(m_+; k_i, k_j; t), \quad (C11)$$

where $\bar{k}_j > 0$ are the possible poles in $w_1(m_+; k_i, k_j; t)$. Noticing that $\mathcal{N}_+(m)Q'_+(m) = 2(\gamma^*)^2/\tilde{\Gamma}$, we can further recast it into the following form

$$Q_+(k_i, k_j; t) = \frac{2i(\gamma^*)^2}{\tilde{\Gamma}} (\mathcal{I}_{ij}(t) + \mathcal{W}(k_i, k_j; t)), \quad (C12)$$

where

$$\mathcal{I}_{ij}(t) = i \left( \int_{0-i\epsilon}^{\pi-i\epsilon} - \int_{0+i\epsilon}^{\pi+i\epsilon} \right) \frac{dm}{2\pi} w_2(m; k_i, k_j; t), \quad \mathcal{W}(k_i, k_j; t) = \sum_l \mathrm{Res}|_{m_+=\bar{k}_l} w_2(m_+; k_i, k_j; t) \quad (C13)$$

with

$$w_2(m_+; k_i, k_j; t) = \frac{\sin^2 m_+ e^{-itE_{m_+}}}{\Delta_+(m_+)(e^{-iQ_+(m_+)} - 1)(\cos m_+ - \cos k_i)(\cos m_+ - \cos k_j)}. \quad (C14)$$

Let us first analyze how $\mathcal{I}_{ij}(t)$ behaves for large $N$. Notice that the second term in $\mathcal{I}_{ij}(t)$ goes to zero when $N \to \infty$. Therefore only the first term remains finite in the thermodynamic limit. Its behavior distinctively differs depending on the value of $\tilde{\Gamma}$ and goes as

$$\lim_{\epsilon \to 0} \lim_{N \to \infty} \mathcal{I}_{ij}(t) = i \left( \mathcal{J}(k_i, k_j; t) + \mathcal{J}\left(\frac{\pi}{2} - k_i, \frac{\pi}{2} - k_j; t\right)^* \right) \\ + \vartheta\left(\frac{1}{2} - \tilde{\Gamma}\right) \frac{i(1 - \tilde{\Gamma})}{8|\tilde{\Gamma}|} \frac{e^{-tE_b}}{(|\tilde{\Gamma}| + i\sqrt{1 - 2\tilde{\Gamma}} \cos k_i)(|\tilde{\Gamma}| + i\sqrt{1 - 2\tilde{\Gamma}} \cos k_j)}, \quad (C15)$$

where

$$\mathcal{J}(k_i, k_j; t) = ie^{2it} \int_0^\infty \frac{dx}{2\pi} \frac{e^{-tx}\sqrt{-x^2 + 4ix}}{(4\tilde{\Gamma}^2 + (1 - 2\tilde{\Gamma})(2 + ix)^2)(2 + ix - 2\cos k_i)(2 + ix - 2\cos k_j)}. \quad (C16)$$

Numerically it can be easily seen that the second term is dominant when $\tilde{\Gamma} < 1/2$, and the first term, which is decaying algebraically, characterizes small ripples around the exponential decay. However the second term vanishes when $\tilde{\Gamma} < 1/2$, and the dominant contribution is taken over by the first term, which for large $t$ behaves as

$$i \left( \mathcal{J}(k_i, k_j; t) + \mathcal{J}\left(\frac{\pi}{2} - k_i, \frac{\pi}{2} - k_j; t\right)^* \right) = \frac{i}{128\sqrt{\pi}(1 - \tilde{\Gamma})^2} \left( \frac{e^{i(2t+\pi/4)}}{\sin^2 k_i \sin^2 k_j} + \frac{e^{-i(2t+\pi/4)}}{\cos^2 k_i \cos^2 k_j} \right) t^{-3/2} + \mathcal{O}(t^{-2}). \quad (C17)$$

(i) $k_i \neq k_j > 0$: (the usual trick) $\mathcal{W}(k_i, k_j; t)$ becomes

$$\mathcal{W}(k_i, k_j; t) = (\mathrm{Res}_{m_+=k_i} + (\mathrm{Res}_{m_+=k_j}) w_2(m_+; k_i, k_j; t) = \mathbf{X}_{ij}(t) + \mathbf{X}_{ji}(t), \quad (C18)$$

where

$$\mathbf{X}_{ij}(t) = \frac{1}{2} \left( \sin k_i - i \frac{(1 - \tilde{\Gamma}) \cos k_i}{\tilde{\Gamma}} \right) \frac{e^{-iE_{k_i}t}}{\Delta_+(k_i)(\cos k_i - \cos k_j)}. \quad (C19)$$
(ii) $k_j < 0 < k_i, k_i \neq -k_j$: In this case $W(k_i, k_j; t)$ can be calculated as

$$W(k_i, k_j; t) = (\mathrm{Res}_{|m_+|=k_i} + (\mathrm{Res}_{|m_+|=-k_j}) w_2(m_+; k_i, k_j; t)) = \mathbf{X}_{ij}(t) + \mathbf{Y}_{ji}(t), \qquad (\mathrm{C}20)$$

where

$$\mathbf{Y}_{ji}(t) = -\frac{1}{2} \left( \sin k_j + i \frac{(1 - \tilde{\Gamma}) \cos k_j}{\tilde{\Gamma}} \right) \frac{e^{-iE_{kj}t}}{\Delta_+(k_j)(\cos k_j - \cos k_i)}. \qquad (\mathrm{C}21)$$

(iii) $k_i < 0 < k_j, k_i \neq -k_j$: Following the same manipulation as above, $W(k_i, k_j; t)$

$$W(k_i, k_j; t) = (\mathrm{Res}_{|m_+|=-k_i} + (\mathrm{Res}_{|m_+|=k_j}) w_2(m_+; k_i, k_j; t)) = \mathbf{Y}_{ij}(t) + \mathbf{X}_{ji}(t). \qquad (\mathrm{C}22)$$

(iv) $k_i \neq k_j < 0$: Likewise,

$$W(k_i, k_j; t) = \mathbf{Y}_{ij}(t) + \mathbf{Y}_{ji}(t). \qquad (\mathrm{C}23)$$

The remaining case is $|k_i| = |k_j|$, for which we have a double pole that can be treated by the same trick above.

(v) $k_i = k_j > 0$: We can evaluate $W(k_i, k_i; t)$ as

$$W(k_i, k_i; t) = \mathrm{Res}_{|m_+|=k_i} w_2(m_+; k_i, k_i; t) = -\frac{ie^{-iE_{kj}t}}{16\tilde{\Gamma}^2 \sin^2 k_i} (N+2) + \mathbf{D}_i(t), \qquad (\mathrm{C}24)$$

where

$$\mathbf{D}_i(t) = \frac{e^{-iE_{kj}t}}{\Delta_+(k_i)} \left( \sin k_i - i \frac{(1 - \tilde{\Gamma}) \cos k_i}{\tilde{\Gamma}} \right) \left( it - \frac{4(1 - 2\tilde{\Gamma}) \cos k_i}{\Delta_+(k_i)} \right) - \frac{e^{-iE_{kj}t}}{2\Delta_+(k_i)} \left( \frac{\cos k_i}{\sin k_i} + i \frac{1 - \tilde{\Gamma}}{\tilde{\Gamma}} \right). \qquad (\mathrm{C}25)$$

(vi) $0 < k_i = -k_j$: In this case $W(k_i, -k_i; t)$ is the same as in the previous case, so $W(k_i, -k_i; t) = \mathbf{D}_i(t)$.

(vii) $0 < k_j = -k_i$: Likewise

$$W(k_i, -k_i; t) = \mathrm{Res}_{|m_+|=-k_i} w_2(m_+; k_i, -k_i; t) = -\frac{ie^{-iE_{kj}t}}{16\tilde{\Gamma}^2 \sin^2 k_i} (N+2) + \mathbf{E}_i(t), \qquad (\mathrm{C}26)$$

where

$$\mathbf{E}_i(t) = \frac{e^{-iE_{kj}t}}{\Delta_+(k_i)} \left( \sin k_i + i \frac{(1 - \tilde{\Gamma}) \cos k_i}{\tilde{\Gamma}} \right) \left( -it + \frac{4(1 - 2\tilde{\Gamma}) \cos k_i}{\Delta_+(k_i)} \right) + \frac{e^{-iE_{kj}t}}{2\Delta_+(k_i)} \left( \frac{\cos k_i}{\sin k_i} - i \frac{1 - \tilde{\Gamma}}{\tilde{\Gamma}} \right). \qquad (\mathrm{C}27)$$

(viii) $k_i = k_j < 0$: In this case $W(k_i, k_i; t)$ is the same as in the previous case, therefore $W(k_i, k_i; t) = \mathbf{E}_i(t)$.

To summarize, we have

$$W(k_i, k_j; t) = \begin{cases} \mathbf{X}_{ij}(t) + \mathbf{X}_{ji}(t) & k_i \neq k_j > 0 \\ \mathbf{X}_{ij}(t) + \mathbf{Y}_{ji}(t) & k_i > 0 > k_j, k_i \neq -k_j \\ \mathbf{Y}_{ij}(t) + \mathbf{X}_{ji}(t) & k_j > 0 > k_i, k_i \neq -k_j \\ \mathbf{Y}_{ij}(t) + \mathbf{Y}_{ji}(t) & k_i \neq k_j < 0 \\ -\frac{ie^{-iE_{kj}t}}{16\tilde{\Gamma}^2 \sin^2 k_i} (N+2) + \mathbf{D}_i(t) & k_i = \pm k_j, k_i > 0 \\ -\frac{ie^{-iE_{kj}t}}{16\tilde{\Gamma}^2 \sin^2 k_i} (N+2) + \mathbf{E}_i(t) & k_i = \pm k_j, k_i < 0 \end{cases}. \qquad (\mathrm{C}28)$$

Therefore

$$\mathbf{M}_+(k_i, k_j; t) = \frac{(\gamma^*)^2}{2\tilde{\Gamma}} \begin{cases} 1 & k_i = k_j > 0 \\ -\mathcal{M}_+ & k_i = -k_j \\ \mathcal{M}_+^2 & k_i = k_j < 0 + \frac{2\pi}{N+2} \mathbf{A}_{ij}(t), \\ 0 & \text{otherwise} \end{cases} \qquad (\mathrm{C}29)$$

where

$$\mathbf{A}_{ij}(t) = \frac{8i(\gamma^*)^2 \tilde{\Gamma} \sin k_i \sin k_j}{2\pi} \begin{cases} \mathbf{I}_{ij}(t) + \mathbf{X}_{ij}(t) + \mathbf{X}_{ji}(t) & k_i \neq k_j > 0 \\ \mathcal{M}_+(\mathbf{I}_{ij}(t) + \mathbf{X}_{ij}(t) + \mathbf{Y}_{ji}(t)) & k_i > 0 > k_j, k_i \neq -k_j \\ \mathcal{M}_+(\mathbf{I}_{ij}(t) + \mathbf{Y}_{ij}(t) + \mathbf{X}_{ji}(t)) & k_j > 0 > k_i, k_i \neq -k_j \\ \mathcal{M}_+^2(\mathbf{I}_{ij}(t) + \mathbf{Y}_{ij}(t) + \mathbf{Y}_{ji}(t)) & k_i \neq k_j < 0 \\ \mathbf{I}_{ii}(t) + \mathbf{D}_i(t) & k_i = k_j > 0 \\ \mathcal{M}_+(\mathbf{I}_{ii}(t) + \mathbf{D}_i(t)) & k_i = -k_j > 0 \\ \mathcal{M}_+(\mathbf{I}_{ii}(t) + \mathbf{E}_i(t)) & k_i = -k_j < 0 \\ \mathcal{M}_+^2(\mathbf{I}_{ii}(t) + \mathbf{E}_i(t)) & k_i = k_j < 0 \end{cases}. \qquad (\mathrm{C}30)$$
Next we turn to $M_-(k_i, k_j; t)$, which is much simpler. Observe

$$L_-(k; m_-) = i \mathcal{N}_{L,-}^* \frac{\sqrt{N+2}}{2} (\vartheta(k) + \mathcal{M}_- \vartheta(-k)) e^{-im\delta_{-(m_-)}} (\delta_{k,m} - \delta_{k,-m}) \\ R_-(k; m_-) = i \mathcal{N}_{R,-}^* \frac{\sqrt{N+2}}{2} (\vartheta(k) + \mathcal{M}_- \vartheta(-k)) e^{-im\delta_{-(m_-)}} (\delta_{k,m} - \delta_{k,-m}). \qquad (C31)$$

Hence

$$M_-(k_i, k_j; t) = \frac{\gamma^2}{2\tilde{\Gamma}} \begin{cases} 1 & k_i = k_j > 0 \\ -\mathcal{M}_- & k_i = -k_j \\ \mathcal{M}_-^2 & k_i = k_j < 0 \\ 0 & \text{otherwise} \end{cases} \qquad (C32)$$

Combining everything, we finally obtain the full $M(k_i, k_j; t)$, which turns out to be of the following form

$$M(k_i, k_j; t) = \delta_{k_i, k_j} e^{-iE_{k_i}t} + \frac{2\pi}{N+2} A_{ij}(t). \qquad (C33)$$

This implies that the Loschmidt echo has the form of the Fredholm determinant as expected

$$\mathcal{L}(t) = \det_{i,j \in \mathrm{FS}} \left( 1 + \frac{2\pi}{N+2} \mathbf{B}(t) \right)_{ij} \underset{N \to \infty}{\longrightarrow} \det(1 + \mathbf{B}(t)), \qquad (C34)$$

where $\mathbf{B}_{ij}(t) = (\mathrm{diag}(e^{iE_{k_i}t})\mathbf{A}(t))_{ij}$, and the determinant in the last line is performed over the $L^2$ space.

## 2. Phase II

When $\gamma_1^2 > \sqrt{1 + \gamma_2^2}$, the system acquires a new type of excitations, which are bound states that are localized across the dot. Labelling them with $m = \pm$, the corresponding overlap matrices read

$$\begin{align*}
L^*(k, +) &= \frac{2(\mathcal{N}_{L,b}^{\mathrm{re}})^*}{\sqrt{N+2}} \sum_{x>0} \left( \vartheta(k) + \frac{\gamma^2}{|\gamma|^2} \vartheta(-k) \right) \sin kx e^{-x/\xi} \\
L^*(k, -) &= \frac{2(\mathcal{N}_{L,b}^{\mathrm{re}})^*}{\sqrt{N+2}} \sum_{x>0} \left( \vartheta(k) + \frac{\gamma^2}{|\gamma|^2} \vartheta(-k) \right) (-1)^x \sin kx e^{-x/\xi} \\
R(k, +) &= \frac{2\mathcal{N}_{R,b}^{\mathrm{re}}}{\sqrt{N+2}} \sum_{x>0} \left( \vartheta(k) + \frac{\gamma^2}{|\gamma|^2} \vartheta(-k) \right) \sin kx e^{-x/\xi} \\
R(k, -) &= \frac{2\mathcal{N}_{R,b}^{\mathrm{re}}}{\sqrt{N+2}} \sum_{x>0} \left( \vartheta(k) + \frac{\gamma^2}{|\gamma|^2} \vartheta(-k) \right) (-1)^x \sin kx e^{-x/\xi}. \tag{C35}
\end{align*}$$

Since

$$\sum_{x>0} \sin kx e^{-x/\xi} = \frac{\sin k + (e^{-1/\xi})^{N/2+1} \sin \frac{kN}{2}}{2(\cos 1/\xi - \sin k)} = \frac{\sqrt{2\tilde{\Gamma} - 1} \sin k}{2(\tilde{\Gamma} - \sqrt{2\tilde{\Gamma} - 1} \cos k)} + \mathcal{O}(e^{-N/(2\xi)}), \qquad (C36)$$

and

$$\sum_{x>0} (-1)^x \sin kx e^{-x/\xi} = - \frac{\sin k + (-e^{-1/\xi})^{N/2+1} \sin \frac{kN}{2}}{2(\cos 1/\xi + \sin k)} = - \frac{\sqrt{2\tilde{\Gamma} - 1} \sin k}{2(\tilde{\Gamma} + \sqrt{2\tilde{\Gamma} - 1} \cos k)} + \mathcal{O}(e^{-N/(2\xi)}), \qquad (C37)$$

we obtain

$$L^*(k, \pm) = \frac{(\mathcal{N}_{L,b}^{\mathrm{im}})^* \sqrt{2\tilde{\Gamma} - 1}}{\sqrt{N+2}} \mathbf{A}_{\pm}^{\mathrm{re}}(k), \quad R(k, \pm) = \frac{\mathcal{N}_{R,b}^{\mathrm{im}} \sqrt{2\tilde{\Gamma} - 1}}{\sqrt{N+2}} \mathbf{A}_{\pm}^{\mathrm{re}}(k) \qquad (C38)$$

with

$$\mathbf{A}_{\pm}^{\mathrm{re}}(k) = \pm \left( \vartheta(k) + \frac{\gamma^2}{|\gamma|^2} \vartheta(-k) \right) \frac{\sin k}{\tilde{\Gamma} \mp \sqrt{2\tilde{\Gamma} - 1} \cos k}. \qquad (C39)$$

The additional contribution to the matrix $M(k_i, k_j; t)$, which we denote $M_b^{\mathrm{re}}(k_i, k_j; t)$, reads

$$M_b^{\mathrm{re}}(k_i, k_j; t) = \frac{(\gamma^*)^2(\tilde{\Gamma} - 1)}{N+2} \left( e^{-iE_{b,t}} \mathbf{A}_+^{\mathrm{re}}(k_i) \mathbf{A}_+^{\mathrm{re}}(k_j) + e^{iE_{b,t}} \mathbf{A}_-^{\mathrm{re}}(k_i) \mathbf{A}_-^{\mathrm{re}}(k_j) \right) + \mathcal{O}(e^{-N/(2\xi)}), = \frac{2\pi}{N+2} (\mathbf{A}_b^{\mathrm{im}})_{ij}(t) + \mathcal{O}(e^{-N/(2\xi)}), \qquad (C40)$$
where

$$ (A_b^{\text{re}})_{ij}(t) = \frac{(\gamma^*)^2(\tilde{\Gamma} - 1)}{2\pi} (e^{-iE_b t} A_+^{\text{im}}(k_i) A_+^{\text{im}}(k_j) + e^{iE_b t} A_-^{\text{im}}(k_i) A_-^{\text{im}}(k_j)). \quad (C41) $$

The Loschmidt echo in phase II can be then obtained by

$$ \mathcal{L}(t) = \det_{i,j \in \text{FS}} \left( 1 + \frac{2\pi}{N+2} (\mathbf{B}(t) + \mathbf{B}_b^{\text{re}}(t)) \right)_{ij} \underset{N \to \infty}{\longrightarrow} \det(1 + \mathbf{B}(t) + \mathbf{B}_b^{\text{re}}(t)), \quad (C42) $$

where

$$ (\mathbf{B}_b^{\text{re}})_{ij}(t) = (\text{diag}(e^{iE_{b,t}}) \mathbf{A}_b^{\text{re}}(t))_{ij}. \quad (C43) $$

### 3. Phase III

In this regime, we have additional bound states with pure imaginary energy on top of the other particle-hole type modes. We can take them into account in exactly the same way as in phase II. What is different from phase II is that the biorthogonal norm, which we also compute in this section, becomes time dependent in this phase. Again labeling these two modes with $m = \pm$, the overlap matrices are given by

$$ \begin{aligned} \mathcal{L}^*(k, \pm) &= \frac{2(\mathcal{N}_{\text{L,b}}^{\text{im}})^*}{\sqrt{N+2}} \left( \vartheta(k) + \frac{\gamma^2}{|\gamma|^2} \vartheta(-k) \right) \sum_{x>0} (\pm i)^x \sin kx e^{-x/\xi} \\ \mathcal{R}(k, \pm) &= \frac{2\mathcal{N}_{\text{R,b}}^{\text{im}}}{\sqrt{N+2}} \left( \vartheta(k) + \frac{\gamma^2}{|\gamma|^2} \vartheta(-k) \right) \sum_{x>0} (\pm i)^x \sin kx e^{-x/\xi}. \end{aligned} \quad (C44) $$

Since

$$ \sum_{x>0} i^x \sin kx e^{-x/\xi} = \frac{\sin k + i(e^{-1/\xi})^{N/2+1} \sin \frac{kN}{2}}{2(\sin 1/\xi - \sin k)} = \frac{\sqrt{1-2\tilde{\Gamma}} \sin k}{2(i\tilde{\Gamma} - \sqrt{1-2\tilde{\Gamma}} \cos k)} + \mathcal{O}(e^{-N/(2\xi)}) \quad (C45) $$

and

$$ \sum_{x>0} (-i)^x \sin kx e^{-x/\xi} = - \frac{\sqrt{1-2\tilde{\Gamma}} \sin k}{2(i\tilde{\Gamma} + \sqrt{1-2\tilde{\Gamma}} \cos k)} + \mathcal{O}(e^{-N/(2\xi)}), \quad (C46) $$

we obtain

$$ \mathcal{L}^*(k, \pm) = \frac{(\mathcal{N}_{\text{L,b}}^{\text{im}})^* \sqrt{1-2\tilde{\Gamma}}}{\sqrt{N+2}} \mathbf{A}_{\pm}(k), \quad \mathcal{R}(k, \pm) = \frac{\mathcal{N}_{\text{R,b}}^{\text{im}} \sqrt{1-2\tilde{\Gamma}}}{\sqrt{N+2}} \mathbf{A}_{\pm}(k) \quad (C47) $$

with

$$ \mathbf{A}_{\pm}^{\text{im}}(k) = \pm \left( \vartheta(k) + \frac{\gamma^2}{|\gamma|^2} \vartheta(-k) \right) \frac{\sin k}{2(i\tilde{\Gamma} \mp \sqrt{1-2\tilde{\Gamma}} \cos k)}. \quad (C48) $$

Notice that, in terms of $\tilde{\Gamma}$, $\sin 1/\xi$ can also be expressed as $\sin 1/\xi = i(1 - \tilde{\Gamma})/\sqrt{1 - 2\tilde{\Gamma}}$. The additional contribution to the matrix $M(k_i, k_j; t)$, which we denote $M_b^{\text{im}}(k_i, k_j; t)$ then reads

$$ M_b^{\text{im}}(k_i, k_j; t) = \frac{(\gamma^*)^2(1 - \tilde{\Gamma})}{N+2} \left( e^{-E_b t} A_+^{\text{im}}(k_i) A_+^{\text{im}}(k_j) + e^{E_b t} A_-^{\text{im}}(k_i) A_-^{\text{im}}(k_j) \right) + \mathcal{O}(e^{-N/(2\xi)}) = \frac{2\pi}{N+2} (A_b^{\text{im}})_{ij}(t) + \mathcal{O}(e^{-N/(2\xi)}), \quad (C49) $$

where

$$ (A_b^{\text{im}})_{ij}(t) = \frac{(\gamma^*)^2(1 - \tilde{\Gamma})}{2\pi} \left( e^{-E_b t} A_+^{\text{im}}(k_i) A_+^{\text{im}}(k_j) + e^{E_b t} A_-^{\text{im}}(k_i) A_-^{\text{im}}(k_j) \right). \quad (C50) $$

The Loschmidt echo in phase III can be computed as

$$ \mathcal{L}(t) = \det_{i,j \in \text{FS}} \left( 1 + \frac{2\pi}{N+2} (\mathbf{B}(t) + \mathbf{B}_b^{\text{im}}(t)) \right)_{ij} \underset{N \to \infty}{\longrightarrow} \det(1 + \mathbf{B}(t) + \mathbf{B}_b^{\text{im}}(t)), \quad (C51) $$

where

$$ (\mathbf{B}_b^{\text{im}})_{ij}(t) = (\text{diag}(e^{iE_{b,t}}) \mathbf{A}_b^{\text{im}}(t))_{ij}. \quad (C52) $$

Next we compute the biorthogonal norm of $|\psi(t)\rangle = e^{-iHt}|\text{GS}\rangle$. Thanks to the fact that $|\mathcal{N}_L|^2 = |\mathcal{N}_R|^2$ [22], we can express the norm $\langle\tilde{\psi}(t)|\psi(t)\rangle$ as

$$ \langle\tilde{\psi}(t)|\psi(t)\rangle = \det_{i,j \in \text{FS}} N(k_i; k_j), \quad (C53) $$
where

$$
\begin{align}
N(k_i; k_j) &= \sum_m e^{-i(E_m - E_m^*)t} L^*(k_i, m) R(k_j, m) = \delta_{k_i, k_j} + \sum_{\epsilon=\pm} (e^{-2iE_{b,\epsilon}^*t} - 1) L^*(k_i, \epsilon) R(k_j, \epsilon) \nonumber \\
&= \delta_{k_i, k_j} + \frac{(\gamma^*)^2(1-\tilde{\Gamma})}{N+2} \sum_{\epsilon=\pm} (e^{-2\epsilon E_{b,\epsilon}} - 1) A_{\epsilon}^{\text{im}}(k_i) A_{\epsilon}^{\text{im}}(k_j) =: \delta_{k_i, k_j} + \frac{2\pi}{N+2} \tilde{N}(t) \tag{C54}
\end{align}
$$

which is again of the form of the Fredholm determinant. Therefore the norm is given by

$$
\langle \tilde{\psi}(t) | \psi(t) \rangle = \det_{i,j \in \text{FS}} \left( 1 + \frac{2\pi}{N+2} \tilde{N}(t) \right). \qquad (\text{C55})
$$

APPENDIX D: DETAILS ABOUT THE NUMERICAL
SIMULATIONS

1. Single particle energies and states

A Hamiltonian Eq. (1) of two leads coupled to the impurity
can be diagonalized using a Bogoliubov transformation

$$
H = \sum_{i,j} c_i^\dagger J_{ij} c_j = \sum_k \mathcal{E}_k d_{R,k}^\dagger d_{L,k}, \quad (\text{D1})
$$

$$
d_{R,k}^{\dagger} = \sum_{j} c_{j}^{\dagger} R_{jk}, \quad d_{L,k} = \sum_{j} L_{kj} c_{j} \quad (\text{D2})
$$

$$
\sum_{ij} L_{ki} J_{ij} R_{jp} = \delta_{kp} \mathcal{E}_p, \quad (\text{D3})
$$

where $\mathcal{E}_p$ are single excitation energies, and $L$ and $R = L^{-1}$
are corresponding left and right eigenvectors of $J$. In the
standard Hermitian case $J = J^\dagger$ a Bogoliubov transformation
simplifies to a unitary one

$$
H_0 = \sum_{i,j} c_i^\dagger J_{ij} c_j = \sum_k \mathcal{E}_k^0 f_k^\dagger f_k, \quad (\text{D4})
$$

$$
f_k^\dagger = \sum_j c_j^\dagger U_{jk}^\dagger, \quad f_k = \sum_j U_{kj} c_j \qquad (\text{D5})
$$

$$
\sum_{ij} U_{ki} J_{ij} U_{jp}^{\dagger} = \delta_{kp} \mathcal{E}_k^0. \qquad (\text{D6})
$$

Single excitation energies $\mathcal{E}_k^0$ are real, whereas two of $\mathcal{E}_k$ in Eq. (D1) are purely imaginary ($E_{b,\pm}^{\text{im}} = \mathcal{E}_{\pm}$) and correspond to the bound states in phase III, see Fig. 3. It turns out that for finite systems sizes $N = 5$ (mod 4) the imaginary energies vanish $E_{b,\pm}^{\text{im}}$ (num) = 0, as in Fig. 9. They arise for a system size larger than the localization length $N > 2\xi$ and $E_{b,\pm}^{\text{im}}$ (num) converge to the value given by Eq. (23) in the thermodynamic limit $N \to \infty$.

2. Two-point correlation functions

A correlation function, e.g., the dot density Eq. (32), can
be defined in two inequivalent ways. The first way is to define
it in a biorthogonal basis through an associated state $\langle \tilde{\psi} |$, as
$\langle \tilde{\psi} | c_0^\dagger c_0 | \psi \rangle$. In the second way, the observable can be calcu-
lated in the standard Hermitian approach, i.e., $\langle \psi | c_0^\dagger c_0 | \psi \rangle$. We
discuss both approaches below.

a. *Correlations in the biorthogonal basis*

A correlation matrix C in the ground state has the following
form

$$
C_{ij} = \langle \tilde{\psi} | c_i^\dagger c_j | \psi \rangle, \qquad (\text{D7})
$$

where a state $|\psi\rangle$ is constructed as

$$
|\psi\rangle = d_{R,k_1}^{\dagger} \dots d_{R,k_{N/2}}^{\dagger} |0\rangle, \quad \langle \tilde{\psi} | = \langle 0 | d_{L,k_{N/2}} \dots d_{L,k_1} |, \quad (\text{D8})
$$

and therefore $\langle \tilde{\psi} | \psi \rangle = 1$. Momentum indices $k_i$ belong to a Fermi sea $\mathcal{E}_{k_i} < 0$ or two imaginary bound states $\text{Im}(\mathcal{E}_{k_i}) \neq 0$.

Creation or annihilation operators can be expressed as a
linear combination of $d_R^\dagger$ or $d_L$ operators

$$
C_{ij} = \sum_{pq} L_{ip} R_{jq} \langle \tilde{\psi} | d_{R,p}^{\dagger} d_{L,q} | \psi \rangle = \sum_{pq} L_{ip} R_{jq} P_{pq} = (RPL)_{ji}, \qquad (\text{D9})
$$

where *P* is a projection operator

$$
P_{pq} = \langle d_{R,p}^{\dagger} d_{L,q} \rangle = \begin{cases} 1, & \text{if } p = q, \mathcal{E}_p < 0 \\ \text{Im}\mathcal{E}_p \neq 0 & \text{or} \\ 0, & \text{otherwise,} \end{cases} \quad (\text{D10})
$$

where we suppressed the parity index $\sigma$.

<figure><img src="image_12.png" /><figcaption>FIG. 9. Bound state energy $E_{b,\pm}^{\text{im}}$ in phase III as a function of system size $L$ for different values of a coupling magnitude $J$ [panel (a)] and of a coupling argument $\varphi$ [panel (b)]. For particular system sizes $L = 5 \, \text{mod} \, 4$ and $L < 2\xi$ the single particle energies are real, so $E_{b,\pm}^{\text{im}} = 0$.</figcaption></figure>
We use Mathematica for numerical calculations. To construct the projection matrix $P_{pq}$, one needs to diagonalize the $J_{ij}$ matrix by using the eigensystem function. Eigenvalues of $J_{ij}$ correspond to single particle excitations $\mathcal{E}_p$ and eigenvectors correspond to matrix $L$. So, $P_{pq}$ can be defined in a for-loop with the condition Eq. (D10).

### b. Correlations in the Hermitian basis

The dot density $d(J, \varphi)$, or any other correlation functions, can be obtained by using Wick theorem, which leads to a determinant form. To demonstrate that, we consider a case with two excitations:

$$ 
\begin{aligned}
d(J, \varphi) &= \frac{\langle \psi | c_0^\dagger c_0 | \psi \rangle}{\langle \psi | \psi \rangle} \\
&= \frac{\langle d_{R,k_1} d_{R,k_2} | c_0^\dagger c_0 | d_{R,k_2}^\dagger d_{R,k_1}^\dagger \rangle}{\langle d_{R,k_1} d_{R,k_2} | d_{R,k_2}^\dagger d_{R,k_1}^\dagger \rangle}.
\end{aligned}
\quad (D11) $$

The denominator can be expanded in $c_p^\dagger, c_q$ fermions:

$$ 
\begin{aligned}
|d_{R,k_2}^\dagger d_{R,k_1}^\dagger\rangle &= \sum_{p,q} R_{pk_2} R_{qk_1} c_p^\dagger c_q^\dagger |0\rangle, \quad (D12) \\
\langle d_{R,k_1} d_{R,k_2} | d_{R,k_2}^\dagger d_{R,k_1}^\dagger \rangle &= \sum_{p,p',q,q'} R_{k_2 p'}^\dagger R_{k_1 q'}^\dagger R_{q k_1} R_{p k_2} R_{q k_1} \langle c_{q'} c_{p'} | c_p^\dagger c_q^\dagger \rangle \\
&= (R^\dagger R)_{k_1 k_1} (R^\dagger R)_{k_2 k_2} - (R^\dagger R)_{k_1 k_2} (R^\dagger R)_{k_2 k_1} \\
&= \det (R^\dagger R)_{i,j}.
\end{aligned}
\quad (D13) $$

We do similar manipulations with the numerator:

$$ 
\begin{aligned}
\langle d_{R,k_1} d_{R,k_2} | c_0^\dagger c_0 | d_{R,k_2}^\dagger d_{R,k_1}^\dagger \rangle \\
&= R_{k_1 0}^\dagger R_{0,k_1} (R^\dagger R)_{k_2 k_2} - R_{k_1 0}^\dagger R_{0,k_2} (R^\dagger R)_{k_2 k_1} \\
&\quad - R_{k_2 0}^\dagger R_{0,k_1} (R^\dagger R)_{k_1 k_2} + R_{k_2 0}^\dagger R_{0,k_2} (R^\dagger R)_{k_1 k_1},
\end{aligned}
\quad (D14) $$

where $(R^\dagger R)_{ij}$ are minors of the matrix $R^\dagger R$. The two-point correlation function for the system of generic size $N$ has the following expression

$$ \frac{\langle \psi | c_i^\dagger c_j | \psi \rangle}{\langle \psi | \psi \rangle} = \frac{1}{\det (R^\dagger R)} \sum_{pq} (-1)^{p+q} R_{pi}^\dagger R_{jq} \tilde{M}_{pq}, \quad (D15) $$

where $p, q, p', q'$ belongs to a Fermi sea or two imaginary bound states. $\tilde{M}_{pq}$ are minors of the matrix $R^\dagger R$, i.e., determinants of the matrix obtained by deleting the $p$th row and the $q$th column of $R^\dagger R$. To obtain the results as for the biorthogonal basis, the matrix $R^\dagger$ should be replaced by $L$ in Eq. (D15).

We compare these two approaches of Eqs. (D9) and (D15) for the dot density on Fig. 10. Although a discontinuity appears in a vicinity of the phase transition line $\varphi = \pi/4$ in both cases, we expect that the jump in the Hermitian basis disappears in the thermodynamic limit.

<figure><img src="image_13.png" /><figcaption>FIG. 10. Dot density as a function of a coupling argument $\varphi(\gamma = 0.2e^{i\varphi}, N = 519)$ for different approaches Eqs. (D9) and (D15). The red curve with $\langle R | c_0^\dagger c_0 | R \rangle$ is not exceeding 1, however it is not a right quantity.</figcaption></figure>

### 3. Loschmidt echo

#### a. Unnormalized echo $\mathcal{L}(t)$

Recall that we define a normalized returning amplitude, or normalized Loschmidt echo, as

$$ \tilde{\mathcal{L}}(t) = \frac{\mathcal{L}(t)}{\sqrt{|\langle \tilde{\psi}(t) | \psi(t) \rangle|}}, \quad (D16) $$

where $\mathcal{L} = \langle \psi_0 | e^{-iHt} | \psi_0 \rangle$. The initial state $|\psi_0\rangle$ is a ground state of a system $H_0(\gamma = 0)$ of the impurity and two disjoint wires, a Hamiltonian $H$ describes wires coupled ($\gamma \neq 0$) to the dot. To calculate the time evolution, we notice that

$$ 
\begin{aligned}
e^{-iHt} |\psi_0\rangle &= e^{-iHt} f_{k_1}^\dagger \dots f_{k_{N/2}}^\dagger |0\rangle \\
&= e^{-iHt} f_{k_1}^\dagger e^{iHt} \dots e^{-iHt} f_{k_{N/2}}^\dagger e^{iHt} |0\rangle \\
&= f_{k_1}^\dagger(t) \dots f_{k_{N/2}}^\dagger(t) |0\rangle = |\psi(t)\rangle,
\end{aligned}
\quad (D17) $$

where each creation operator evolves according to

$$ 
\begin{aligned}
f_k^\dagger(t) &= e^{-iHt} f_k^\dagger e^{iHt} = \sum_p e^{-iHt} d_p^\dagger (LU^\dagger)_{pk} e^{iHt} \\
&= \sum_p d_p^\dagger e^{-it\mathcal{E}_p} (LU^\dagger)_{pk} = \sum_q f_q^\dagger (URe^{-it\mathcal{E}}LU^\dagger)_{qk} \\
&= \sum_q f_q^\dagger A(t)_{qk}.
\end{aligned}
\quad (D18) $$

Then the numerator of the Loschmidt echo is given by

$$ 
\begin{aligned}
\mathcal{L}(t) &= \langle \psi_0 | \psi(t) \rangle \\
&= \langle 0 | f_{k_{N/2}} \dots f_{k_1} f_{k_1}^\dagger(t) \dots f_{k_{N/2}}^\dagger(t) | 0 \rangle \\
&= \det \langle 0 | f_k f_p^\dagger(t) | 0 \rangle = \det A(t)_{kp},
\end{aligned}
\quad (D19) $$

indices run a Fermi sea $\mathcal{E}_{k,p}^0 < 0$, or, simply, we just eliminate rows in $U$, which correspond to positive singe excitation energies.

#### b. Norm of a time evolving state

According to the recipe of constructing a biorthogonal basis, in the conjugated state $\langle \tilde{\psi}(t) |$ we replace $R^\dagger$ by $L = R^{-1}$
<figure><img src="image_14.png" /><figcaption>FIG. 11. Rescaled Loschmidt echo as a function of rescaled time for imaginary coupling $\gamma = 0.2e^{i\pi/2}$. Scaling factor $T_K = |E_b^{im}|$ Eq. (23). Color labels different values of a system size.</figcaption></figure>

and correspondingly $L^\dagger$ by $R = L^{-1}$. Then the PT norm of the time evolving state $|\psi(t)\rangle$ is

$$ \langle \tilde{\psi}(t) | \psi(t) \rangle = \det(URe^{-it(\mathcal{E}-\mathcal{E}^*)}LU^\dagger)_{kp}, \quad (\text{D20}) $$

where indices run a Fermi sea $\mathcal{E}_{k,p}^0 < 0$, as in Eq. (D19). For numerical values for normalized Loschmidt echo converges for the large system sizes, see Fig. 11.

[1] A. C. Hewson, *The Kondo Problem to Heavy Fermions*, Cambridge Studies in Magnetism (Cambridge University Press, Cambridge, 1993).

[2] E. S. Sørensen and I. Affleck, Phys. Rev. B **53**, 9153 (1996).

[3] P. W. Anderson, J. Phys. C: Solid State Phys. **3**, 2436 (1970).

[4] N. Andrei, Phys. Rev. Lett. **45**, 379 (1980).

[5] A. Tsvelick and P. Wiegmann, Adv. Phys. **32**, 453 (1983).

[6] R. Bulla, T. A. Costi, and T. Pruschke, Rev. Mod. Phys. **80**, 395 (2008).

[7] E. Boulat, H. Saleur, and P. Schmitteckert, Phys. Rev. Lett. **101**, 140601 (2008).

[8] M. Kanász-Nagy, Y. Ashida, T. Shi, C. P. Moca, T. N. Ikeda, S. Fölling, J. I. Cirac, G. Zaránd, and E. A. Demler, Phys. Rev. B **97**, 155156 (2018).

[9] M. Ferrier, T. Arakawa, T. Hata, R. Fujiwara, R. Delagrange, R. Deblock, Y. Teratani, R. Sakano, A. Oguri, and K. Kobayashi, Phys. Rev. Lett. **118**, 196803 (2017).

[10] M. Nakagawa, N. Kawakami, and M. Ueda, Phys. Rev. Lett. **121**, 203001 (2018).

[11] J. A. S. Lourenço, R. L. Eneias, and R. G. Pereira, Phys. Rev. B **98**, 085126 (2018).

[12] P. Fendley, H. Saleur, and A. B. Zamolodchikov, Int. J. Mod. Phys. A **08**, 5717 (1993).

[13] A. Zamolodchikov, Phys. Lett. B **335**, 436 (1994).

[14] H. Saleur and S. Skorik, Phys. Lett. B **336**, 205 (1994).

[15] Y. Ashida, S. Furukawa, and M. Ueda, Nat. Commun. **8**, 15791 (2017).

[16] S. Ghosh, P. Ribeiro, and M. Haque, J. Stat. Mech. (2014) P04011.

[17] H. Fröml, A. Chiocchetta, C. Kollath, and S. Diehl, Phys. Rev. Lett. **122**, 040402 (2019).

[18] H. Fröml, C. Muckel, C. Kollath, A. Chiocchetta, and S. Diehl, Phys. Rev. B **101**, 144301 (2020).

### 4. Analytical continuation of Loschmidt echo

We expand Eq. (D16) in Taylor series of the following form:

$$ \tilde{\mathcal{L}}(t) = \sum_{n=0} a_n \gamma^n. \quad (\text{D21}) $$

In order to obtain coefficients $a_n$ from numerical data, we implement central finite difference operators (discrete derivative) of higher orders [38] in the vicinity of $\gamma = 0$, which help to increase precision of calculations. For instance

$$ \begin{aligned} h f'(x)_{x=0} = & \frac{1}{280} f(x - 4h) + \frac{-4}{105} f(x - 3h) + \frac{1}{5} f(x - 2h) \\ & + \frac{-4}{5} f(x - h) + \frac{4}{5} f(x + h) + \frac{-1}{5} f(x + 2h) \\ & + \frac{4}{105} f(x + 3h) + \frac{-1}{280} f(x + 4h) + \mathcal{O}(h^8), \end{aligned} \quad (\text{D22}) $$

where $h$ is a distance between neighboring points, which we chose $h = \Delta\gamma = 0.01$ and $\gamma \in (-0.2, 0.2)$. We want to pay attention that decreasing $h$ may lead to increasing errors, so one should choose an optimal $h$.

[19] J. Bauer, C. Salomon, and E. Demler, Phys. Rev. Lett. **111**, 215304 (2013).

[20] L. Riegger, N. Darkwah Oppong, M. Höfer, D. R. Fernandes, I. Bloch, and S. Fölling, Phys. Rev. Lett. **120**, 143601 (2018).

[21] K. Kawabata, Y. Ashida, H. Katsura, and M. Ueda, Phys. Rev. B **98**, 085116 (2018).

[22] D. C. Brody, J. Phys. A: Math. Theor. **47**, 035305 (2013).

[23] L. D. Landau and L. M. Lifshitz, *Quantum Mechanics Non-Relativistic Theory* (Pergamon Press, London, 1981).

[24] M.-A. Miri and A. Alù, Science **363**, 67709 (2019).

[25] Z. Oztas and C. Yuce, Phys. Rev. A **98**, 042104 (2018).

[26] E. Boulat and H. Saleur, Phys. Rev. B **77**, 033409 (2008).

[27] G. Camacho, P. Schmitteckert, and S. T. Carr, Phys. Rev. B **99**, 085122 (2019).

[28] I. Affleck and A. W. W. Ludwig, Phys. Rev. Lett. **67**, 161 (1991).

[29] T. Yoshimura, K. Bidzhiev, and H. Saleur (unpublished).

[30] J.-M. Stéphan and J. Dubail, J. Stat. Mech. (2011) P08019.

[31] L. Piroli, B. Pozsgay, and E. Vernier, Nucl. Phys. B **933**, 454 (2018).

[32] R. Vasseur, K. Trinh, S. Haas, and H. Saleur, Phys. Rev. Lett. **110**, 240601 (2013).

[33] O. Gamayun, O. Lychkovskiy, and J.-S. Caux, SciPost Phys. **8**, 036 (2020).

[34] D. M. Kennes, V. Meden, and R. Vasseur, Phys. Rev. B **90**, 115101 (2014).

[35] P. Fendley, A. W. W. Ludwig, and H. Saleur, Phys. Rev. Lett. **74**, 3005 (1995).

[36] C. Latta, F. Haupt, M. Hanl, A. Weichselbaum, M. Claassen, W. Wuester, P. Fallahi, S. Faelt, L. Glazman, J. von Delft, H. E. Türeci, and A. Imamoglu, Nature (London) **474**, 627 (2011).

[37] B. Urbaszek, X. Marie, T. Amand, O. Krebs, P. Voisin, P. Malentinsky, A. Högele, and A. Imamoglu, Rev. Mod. Phys. **85**, 79 (2013).

[38] B. Fornberg, Mathematics of Computation **51**, 699 (1988).