Non-equilibrium entanglement asymmetry for discrete groups: the example of the XY spin chain

In the last few years, the interplay between entanglement and symmetries has become the centre of intense research activity, yielding up a more refined vision of the behaviour of many-body quantum systems [1–6] and creating a new framework to study entanglement [7–21]. For example, novel quantities, such as the symmetry-resolved entanglement entropy [7–9], have been conceived to analyse how entanglement is distributed in the symmetry sectors of a theory. On the other hand, the idea of using entanglement tools to study symmetry breaking has recently arisen. To this end, entanglement asymmetry was proposed in [22] as a measure of how much a symmetry is broken in a part of an extended quantum system. So far, the entanglement asymmetry has been applied to examine the time evolution of an initially broken U(1) symmetry after a quench with a Hamiltonian that preserves it, both in free [22, 23] and interacting integrable systems [24]. In general, it is expected that the symmetry is dynamically restored in a subsystem after the quench. However, unexpectedly, the entanglement asymmetry reveals that the symmetry can be restored earlier when it is initially more broken [22]. This surprising phenomenon can be seen as a quantum version of the very counter-intuitive Mpemba effect: the more a system is initially out of equilibrium, the faster it relaxes. Entanglement asymmetry was also employed in [23] to analyse the special case in which symmetry is not restored and the subsystem relaxes to a non-Abelian generalised Gibbs ensemble (GGE).

The previous works are restricted to the continuous U(1) group. The goal of the present manuscript is to extend the study of the entanglement asymmetry to finite, discrete cyclic $\mathbb{Z}_N$ groups. In particular, we consider the $\mathbb{Z}_2$ spin flip parity symmetry of the XY spin chain which, in the thermodynamic limit, is spontaneously broken by the ground state in the ferromagnetic phase [25]. Usually, this symmetry breaking is detected using the one-site longitudinal magnetisation as an order parameter. This order parameter presents the inconvenience that, despite a non-zero value indicating that the symmetry is broken, the converse is not always true [26]. On the contrary, the entanglement asymmetry unambiguously determines whether the symmetry is respected or not. Moreover, although for a single spin the order parameter can be used to measure symmetry breaking, for larger subsystems there are longer range correlations that can break the symmetry and are only taken into account by the entanglement asymmetry. Using this, we study the non-equilibrium time evolution of the $\mathbb{Z}_2$ spin parity symmetry after a global quantum quench in the XY spin chain Hamiltonian, extending the study performed in [27–29] employing the order parameter. Global quenches in this system have been studied from many perspectives; see e.g. [30–44]. We distinguish between chains with open and periodic boundary conditions (OBC and PBC), since the approach to calculate the entanglement asymmetry and the results obtained differ: while in the periodic case this $\mathbb{Z}_2$ symmetry is always restored, the same does not always occur in the open chain due to the presence of boundary modes. Moreover, in the periodic case, by comparing with the analytic results obtained in [27, 28] for the order parameter in the thermodynamic limit, we conjecture an effective expression that describes the evolution of the entanglement asymmetry after the quench.

The paper is organised as follows: in section 2, we define the entanglement asymmetry, we review its main known features for the U(1) group and we extend it to cyclic $\mathbb{Z}_N$ groups, describing some of its principal properties in this case. In section 3, we introduce the XY spin chain and we present the behaviour of the entanglement asymmetry of the spin parity in the ground state. Section 4 is devoted to calculating and discussing the entanglement asymmetry after a quantum quench of the XY spin chain with OBC, while in section 5 we carry out a similar analysis but using PBC. In section 6, we draw our conclusions. The manuscript also includes several appendices that describe in detail the more technical points of the main text.

Let us consider an extended quantum system that can be divided into two spatial parts A and B. If we assume that the total system is in a pure state $\ket{\Psi}$, then the state of one of the regions, for example A, is described by the reduced density matrix $\rho_A = {\mathrm{Tr}}_B (\ket{\Psi} \bra{\Psi})$, which is generally a mixed state. We further consider a local charge operator Q given by the sum of the charge in A and B, $Q = Q_A+Q_B$. We assume that Q generates an Abelian group with elements of the form $\mathrm{e}^{\mathrm{i}\alpha Q}$, where α is a real number. This includes the unitary group U(1) or the cyclic groups $\mathbb{Z}_N$.

It is well known that, when $\ket{\Psi}$ respects the symmetry associated with Q, i.e. it is an eigenstate of Q, the reduced density matrix of any subsystem A displays a block-diagonal structure in the eigenbasis of Q_A , and the blocks correspond to the sectors of a given charge. In this case, the entanglement entropy of subsystem A can be decomposed into the contributions of each charge sector, known as symmetry-resolved entanglement entropies [7–9].

In the case of symmetry-breaking states, entanglement entropy can be used to quantify how much the symmetry generated by Q is broken in a subsystem. In fact, if the state $\ket{\Psi}$ locally breaks the symmetry in the region A, then ρ_A does not commute with Q_A . This means that, in the eigenbasis of Q_A , ρ_A presents non-zero entries outside the sectors of fixed charged $\rho_{q_j}$,

Zoom In Zoom Out Reset image size

here represented by $*$. From this state, it is always possible to define another density matrix $\rho_{A, Q}$ that respects the symmetry, i.e. $[\rho_{A, Q}, Q_A] = 0$. If we take the projection of ρ_A on the different charge sectors of Q_A and we sum them all, we obtain a block diagonal matrix

Zoom In Zoom Out Reset image size

where $\Pi_{q_j}$ is the projector on the subspace of eigenvalue q_j . The density matrix $\rho_{A, Q}$ can be physically interpreted as the result of performing a non-selective measurement of Q_A in ρ_A . From ρ_A and its projection $\rho_{A, Q}$, we can define the entanglement asymmetry [22]

$$\begin{align} \Delta S_A = S\left(\rho_{A, Q}\right)-S\left(\rho_A\right), \end{align} \tag{ 1 }$$

in terms of their von Neumann entropies

$$\begin{align} S\left(\rho\right) = -{\mathrm{Tr}}\left(\rho\log\rho\right). \end{align} \tag{ 2 }$$

The entanglement asymmetry measures how much the symmetry generated by Q is broken in the subsystem A. In fact, $\Delta S_A$ satisfies two crucial properties: it is always non-negative, $\Delta S_A\unicode{x2A7E} 0$, and it vanishes if and only if $[\rho_A, Q_A] = 0$ [45]. These properties are satisfied whether the total system is in a pure or mixed state. Therefore, the entanglement asymmetry can also be applied to the study of symmetry breaking in mixed states and, in particular, to a system at a certain temperature.

Since the direct calculation of the von Neumann entropy is normally a difficult task, the usual strategy [46, 47] is to replace it in equation (1) by the Rényi entropies

$$\begin{align} S^{\left(n\right)}\left(\rho\right) = \frac{1}{1-n}\log {\mathrm{Tr}}\left(\rho^n\right). \end{align} \tag{ 3 }$$

Then we define the Rényi entanglement asymmetry as

$$\begin{align} \Delta S_A^{\left(n\right)} = S^{\left(n\right)}\left(\rho_{A, Q}\right)-S^{\left(n\right)}\left(\rho_A\right). \end{align} \tag{ 4 }$$

Note that, if we take in this expression the limit n → 1, we recover equation (1). Moreover, Rényi entanglement asymmetries are also non-negative quantities, $\Delta S_A^{(n)}\unicode{x2A7E} 0$, and they cancel if and only if the symmetry is respected in subsystem A, i.e. $[\rho_A, Q_A] = 0$ [48]. It is also important to remark that the Rényi entanglement asymmetries for an integer $n\unicode{x2A7E} 2$ are experimentally measurable in ion traps through randomised measurements [49–53].

2.1. Review of entanglement asymmetry for the U(1) group

In the literature, the entanglement asymmetry has only been analysed for the Abelian unitary group U(1). In this case, the paradigmatic example to understand how it works is the tilted ferromagnetic state,

$$\begin{align} \ket{\theta} = \mathrm{e}^{\mathrm{i}\frac \theta 2 \sum_{j} \sigma^y_j} \ket{\uparrow \uparrow \uparrow \cdots}. \end{align} \tag{ 5 }$$

For θ ≠ 0, this state breaks the U(1) symmetry corresponding to rotations around the z axis, which are generated by the transverse magnetisation $Q = \frac12 \sum_j \sigma^z_j$. It has been shown that, if the subsystem A is an interval of contiguous spins of length $\ell$, the entanglement asymmetry behaves as [22]

$$\begin{align}\Delta S_A^{(n)}=\frac 1 2 \log \ell +\frac 1 2\log \frac{\pi n^{\frac{1}{n-1}}\sin^2\theta }{2}+\mathcal{O}(\ell^{-1})\end{align} \tag{ 6 }$$

for large $\ell$. To better understand this result, let us focus on the n = 2 Rényi asymmetry $\Delta S_A^{(2)}$. Using the Fourier representation of the projectors $\Pi_{q_j}$, the projected matrix $\rho_{A, Q}$ can be written as

$$\begin{align}\rho_{A, Q} = \int_{0}^{2\pi} \frac{d\alpha}{2\pi} e^{i\alpha Q_A}\rho_A e^{-i\alpha Q_A}.\end{align} \tag{ 7 }$$

Since $\ket{\theta}$ is a separable state, $\rho_A = \ket{\theta, \ell}\bra{\theta,\ell}$, and $\ket{\theta,\ell}$ is the restriction of $\ket{\theta}$ to $\ell$ contiguous spins. We can see equation (7) as the state of A after applying a random element of the U(1) group to ρ_A with a uniform probability density. Consequently, $\Delta S_A$ can be interpreted as the quantum information loss of the subsystem A after such an operation. Despite the infinite number of elements in the U(1) group, this information loss is finite because the two states $\ket{\theta, \ell}$ and $\mathrm{e}^{\mathrm{i}\alpha Q_A}\ket{\theta, \ell}$ are not orthogonal for $\alpha\neq 0$ and finite $\ell$ but have a non-zero overlap. In fact, the n = 2 Rényi entanglement asymmetry can be written as

$$\begin{align}\Delta S_A^{(2)} = -\log \left(\int_{0}^{2\pi}\frac{d\alpha}{2\pi}|\bra{\theta, \ell}e^{i\alpha Q_A}\ket{\theta,\ell}|^2 \right).\end{align} \tag{ 8 }$$

This integral is the average overlap between the states $\ket{\theta, \ell}$ and $\mathrm{e}^{\mathrm{i}\alpha Q_A}\ket{\theta, \ell}$ with respect to α. Hence, according to equation (8), $\Delta S_A^{(2)}$ grows as the overlap between $\ket{\theta,\ell}$ and $\mathrm{e}^{\mathrm{i}\alpha Q_A}\ket{\theta,\ell}$ decreases and, consequently, the maximal entanglement asymmetry is obtained when $\ket{\theta, \ell}$ and $\mathrm{e}^{\mathrm{i}\alpha Q_A}\ket{\theta,\ell}$ are orthogonal for all $\alpha\in(0,2\pi)$. A simple calculation gives that

$$\begin{align}|\bra{\theta,\ell}e^{i\alpha Q_A}\ket{\theta,\ell}|^2=\left(\cos^2(\alpha/2)+\cos^2\theta\sin^2(\alpha/2)\right)^\ell,\end{align} \tag{ 9 }$$

which, for $\alpha,\theta\in (0,2\pi)$, only vanishes in the limit $\ell\to \infty$. Therefore, for a given value of the tilting angle θ, the entanglement asymmetry of the U(1) group generated by the transverse magnetisation is unbounded. Moreover, using equation (9), we obtain that the average overlap between $\ket{\theta,\ell}$ and $\mathrm{e}^{\mathrm{i}\alpha Q_A}\ket{\theta,\ell}$ decreases as $1/\sqrt{\ell}$ for large $\ell$, giving rise to the logarithmic divergence of $\Delta S_A^{(2)}$ in equation (6).

2.2. Entanglement asymmetry for cyclic groups $\mathbb{Z}_N$

Let us now consider the finite cyclic group $\mathbb{Z}_N$ of rotations around the z axis by multiple angles of $2\pi/N$, which is a subgroup of the U(1) symmetry generated by the transverse magnetisation Q. The density matrix $\rho_{A, Q_N}$ of this group is obtained by projecting on the different charge sectors of the transverse magnetisation Q modulo N. Using the Fourier representation of the projectors, it can be written as

$$\begin{align}\rho_{A, Q_N} =\frac{1}{N}\sum_{j=1}^{N} e^{i\alpha_j Q_A}\rho_A e^{-i\alpha_j Q_A},\quad \alpha_j =\frac{2\pi j}{N}.\end{align} \tag{ 10 }$$

This matrix can be interpreted either as the state of subsystem A after a non-selective measurement of the transverse magnetisation Q modulo N or after randomly applying to ρ_A one of the elements of the $\mathbb{Z}_N$ group with equal probabilities. The second interpretation gives the intuition that the maximal information loss after such an operation is $\log N$, since N is the dimension of the group, and one can always recover all the information on the state of A by learning which element of the group has been applied. In the left panel of figure 1, we study the entanglement asymmetry of the tilted ferromagnetic state $\ket{\theta}$ as a function of the subsystem length $\ell$ for different subgroups $\mathbb{Z}_N$ of the U(1) group of rotations around the z axis seen before. For a small subsystem size, the entanglement asymmetry presents the same behaviour as in the U(1) case. However, for large $\ell$, it precisely saturates to $\log N$. This saturation occurs at larger values of $\ell$ as the dimension of the subgroup increases, recovering the result of equation (6) for the U(1) group in the limit $N\to \infty$. In particular, for N = 2, we have

$$\begin{align} \Delta S_A^{\left(n\right)} = H^{\left(n\right)}\left(\cos^\ell \theta\right), \end{align} \tag{ 11 }$$

where

$$\begin{align} H^{\left(n\right)}\left(x\right) = \frac{1}{1-n}\log \left[\left(\frac{1+x}{2}\right)^n + \left(\frac{1-x}{2}\right)^n\right] \end{align} \tag{ 12 }$$

are the functions for the Rényi-n entropies of one bit of information. Their limit n → 1 reads

$$\begin{align} H^{\left(1\right)} \left(x\right) : = \lim_{n\to 1} H^{\left(n\right)}\left(x\right) = -\frac{1+x}{2}\log\left(\frac{1+x}{2}\right) - \frac{1-x}{2}\log\left(\frac{1-x}{2}\right). \end{align} \tag{ 13 }$$

In the right panel of figure 1, we plot equation (11) (i.e. for N = 2) as a function of the tilting angle θ for n = 1 and different subsystem sizes. We can see that $\Delta S_A$ grows monotonically with θ until $\theta = \pi/2$, at which value the U(1) symmetry is maximally broken since all the spins point in the x direction. We observe that, as the subsystem size increases, the asymmetry saturates faster to its maximum value, $\log 2$, when the tilted angle θ turns on. In figure 2, we study the Rényi entanglement asymmetry for the subgroups $\mathbb{Z}_2$ (left panel) and $\mathbb{Z}_3$ (right panel) as a function of the tilting angle θ for different Rényi indices n and a fixed subsystem length. As we can see, for the subgroup $\mathbb{Z}_2$, $\Delta S_A^{(n)}$ is monotonic in θ for $n\to\infty$. On the other hand, for $\mathbb{Z}_3$, it is non-monotonic in θ when $n\to\infty$, presenting the same peaked behaviour as in the U(1) case [22]. This is expected since, as we show in the left panel of figure 1, $\Delta S_A^{(n)}$ for the subgroups $\mathbb{Z}_N$ tends to the U(1) asymmetry in the limit $N\to \infty$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We can generally prove that, for $\mathbb{Z}_N$ groups, the entanglement asymmetry $\Delta S_A$ is bounded by the dimension of the group. This is actually a consequence of the following inequality for the von Neumann entropy. If we consider a set of density matrices $\{\rho_j\}$ and of positive coefficients $\{\lambda_j\}$ with $\sum_j \lambda_j = 1$, then [54]

$$\begin{align} S\left(\sum_j \lambda_j \rho_j\right)\unicode{x2A7D} \sum_j \lambda_j S\left(\rho_j\right) -\sum_j \lambda_j \log\lambda_j. \end{align} \tag{ 14 }$$

Therefore, if in the definition (1) of the entanglement asymmetry we take into account that for a $\mathbb{Z}_N$ group the projected density matrix $\rho_{A, Q_N}$ can be written in the form of equation (10), and we apply the previous inequality identifying $\lambda_j = 1/N$ and $\rho_j = \mathrm{e}^{\mathrm{i}\alpha_j Q_A} \rho_A \mathrm{e}^{-\mathrm{i}\alpha_j Q_A}$, then we can easily conclude that

$$\begin{align} \Delta S_A\unicode{x2A7D} \log N. \end{align} \tag{ 15 }$$

The equality in equation (14) is attained if and only if the states ρ_j have support in orthogonal subspaces. This implies that the entanglement asymmetry for a $\mathbb{Z}_N$ group saturates to its maximal value, $\Delta S_A = \log N$, if and only if all the states $\mathrm{e}^{\mathrm{i}\alpha_j Q_A}\rho_A \mathrm{e}^{-\mathrm{i}\alpha_j Q_A}$ have orthogonal support. Therefore, the $\mathbb{Z}_N$ symmetry is maximally broken if the action of any element of the group, except the identity, maps ρ_A to different orthogonal subspaces. In the tilted ferromagnetic state, this is equivalent to requiring that $\ket{\theta,\ell}$ and $\mathrm{e}^{\mathrm{i}\alpha_j Q_A}\ket{\theta,\ell}$ are orthogonal, similarly to the discussion done in section 2.1 for the U(1) case, but now the average (8) is over a finite number of states and, therefore, is bounded. In appendix A, we prove that these properties are also satisfied by the Rényi entanglement asymmetry (4) for any index n. Interestingly, the equality in equation (14) corresponds to a particular instance of Holevo theorem [54], a result of quantum information theory that states that the maximum information that one can extract from a mixed state such as $\rho_{A, Q_N}$ is given by the difference between $S(\sum_j\lambda_j\rho_j)$ and $\sum_j\lambda_j S(\rho_j)$, which in our case is precisely $\log N$, as we already stated.

In the rest of the paper, we study the non-equilibrium dynamics of a $\mathbb{Z}_2$ broken symmetry in the XY spin chain. This section is devoted to introducing this model and discussing some of its general features that will be useful afterwards. The Hamiltonian of the XY spin chain reads

$$\begin{align} H = -\sum_{j = 1}^{L-1}\left[ \frac{1+\gamma}{4}\sigma^x_j \sigma^x_{j+1} + \frac{1-\gamma}{4}\sigma^y_j \sigma^y_{j+1}\right] - \frac{1}{2}\sum_{j = 1}^{L} h\sigma^z_j\,\; \left(+ H_\mathrm{b}\right), \end{align} \tag{ 16 }$$

where

$$\begin{align} H_\mathrm{b} = - \left[ \frac{\left(1+\gamma\right)}{4}\sigma^x_L \sigma^x_{1} + \frac{\left(1-\gamma\right)}{4}\sigma^y_L \sigma^y_{1} \right] \end{align} \tag{ 17 }$$

is a boundary term closing the chain with PBCs. If this term is not present, then the chain has OBCs. We observe that the Hamiltonian contains nearest-neighbour interactions between the x and y components of the spins. The parameter γ > 0 tunes the anisotropy between the couplings in these two directions. There is also an external magnetic field of intensity h > 0 in the z direction.

The reason why we are interested in this system is that it exhibits a $\mathbb{Z}_2$ symmetry corresponding to rotations of π rad around the z axis, a subgroup of the U(1) symmetry discussed in section 2.1 and in [22, 23]. It is associated with the parity operator

$$\begin{align} P = \prod_{j = 1}^L \sigma^z_j = \left(-i\right)^L \mathrm{e}^{\mathrm{i}\pi Q},\quad \text{where}\,\, Q = \frac{1}{2}\sum_{j = 1}^L\sigma_j^z, \end{align} \tag{ 18 }$$

which maps the spins $\sigma_j^{x, y}$ to $-\sigma_j^{x,y}$ leaving $\sigma_j^z$ invariant. This symmetry is spontaneously broken by the ground state of equation (16) in the thermodynamic limit $L\to\infty$ when h < 1 (ferromagnetic phase) while, for h > 1 (paramagnetic phase), the symmetry is respected. Both phases are separated by the critical line h = 1.

These two phases can be distinguished using the one-site longitudinal magnetisation $\langle \sigma_j^x\rangle$ as the order parameter. In fact, for h < 1, we have that $\langle \sigma_j^x \rangle \neq 0$ while, for h > 1, $\langle \sigma_j^x\rangle = 0$. Since a non-zero magnetisation in the x direction obviously breaks the $\mathbb{Z}_2$ symmetry of π rad rotations about the z axis, then the condition $\langle \sigma_j^x\rangle\neq 0$ implies that such symmetry is broken in the ground state.

We can also analyse this spontaneous breaking of the $\mathbb{Z}_2$ parity symmetry with the Rényi entanglement asymmetry (4). In figure 3, we represent the $(h, \gamma)$-parameter plane of the XY spin chain, indicating the large subsystem value of the asymmetry in each phase. We find that, in the ferromagnetic phase h < 1, $\Delta S_A^{(n)}$ quickly saturates for any Rényi index n to its maximal value, $\log 2$, as the length of the subsystem $\ell$ increases. On the other hand, in the paramagnetic region h > 1, the asymmetry vanishes, signalling that the ground state respects the symmetry. Therefore, for infinitely large subsystems $\Delta S_A^{(n)}$ acts somehow as a topological indicator of symmetry breaking. It is also worth mentioning that the tilted ferromagnetic state $\ket{\theta}$, for which we obtained its exact $\mathbb{Z}_2$ entanglement asymmetry in section 2.2, is the ground state of the XY spin chain along the circle $\gamma^2+h^2 = 1$ [55, 56].

Zoom In Zoom Out Reset image size

In this work, we are interested in studying the time evolution of this $\mathbb{Z}_2$ symmetry in a global quantum quench from the ground state of the Hamiltonian with couplings $(h_0, \gamma_0)$ in the ferromagnetic phase, i.e. $h_0\lt1$, to another XY Hamiltonian with different couplings $(h_\mathrm{f}, \gamma_\mathrm{f})$. One has $[H, P] = 0$ for any post-quench XY Hamiltonian H. Hence, after this quench, the state that describes a subsystem A is expected to relax to a GGE that respects the $\mathbb{Z}_2$ symmetry [28, 29, 38]. This symmetry restoration has actually been studied using the order parameter $\langle \sigma_j^x\rangle$ in [27, 28]. It was found to exhibit an exponential decay in time, both for quenches in the ferromagnetic and paramagnetic phases. In our analysis, we distinguish two cases. In section 4, we consider the entanglement asymmetry of an interval attached to one of the extremes of an open chain, a situation that in the thermodynamic limit $L\to \infty$ corresponds to a semi-infinite chain. In section 5, we examine the entanglement asymmetry in the thermodynamic limit of a chain with PBCs. In each case, we apply a different approach and we find that, in some situations, the evolution of the asymmetry after the quench differs.

3.1. Jordan–Wigner transformation

One of the most important features of the XY spin chain (16) is that it is free, i.e. it can be mapped to a Hamiltonian describing non-interacting spinless fermions [57]. This can be done by means of the Jordan–Wigner transformation

$$\begin{align} c_j = \left[\prod_{l = 1}^{j-1}\sigma^z_l\right] \frac{\sigma^x_j + i\sigma^y_j}{2}, \quad c_j^\dagger = \left[\prod_{l = 1}^{j-1}\sigma^z_l\right] \frac{\sigma^x_j - i\sigma^y_j}{2}, \quad \sigma_j^z = 1-2c_j^\dagger c_j, \end{align} \tag{ 19 }$$

where c_j and $c_j^\dagger$ are creation and annihilation fermionic operators, satisfying the anticommutation relations

$$\begin{align} \left\{c_j, c_{j^{\prime}}\right\} = \left\{c_j^\dagger, c_{j^{\prime}}^\dagger\right\} = 0\text{ and } \left\{c_j^\dagger, c_{j^{\prime}}\right\} = \delta_{jj^{\prime}}. \end{align} \tag{ 20 }$$

Under this mapping, the parity operator (18) becomes

$$\begin{align} P = \prod_{j = 1}^L \left(1-2c_j^\dagger c_j\right) = \mathrm{e}^{\mathrm{i}\pi \mathcal N},\quad \text{where } \mathcal N = \sum_{j = 1}^L c_j^\dagger c_j, \end{align} \tag{ 21 }$$

which corresponds to the parity of the number of fermions $\mathcal N$. The Hamiltonian (16) is mapped to the Kitaev fermionic chain

$$\begin{align} H = -\frac{1}{2}\left[\sum_{j = 1}^{L-1}\left(c_j^\dagger c_{j+1}+c_{j+1}^\dagger c_j\right)+\gamma \left(c_j^\dagger c_{j+1}^\dagger+c_{j+1}c_j\right)+ \sum_{j = 1}^{L} h\left(1-2 c_j^\dagger c_j\right)\right] \left(+ H_\mathrm{b}\right), \end{align} \tag{ 22 }$$

which is local in terms of the fermionic operators in spite of the non-locality of the Jordan–Wigner transformation. Observe that the transverse magnetic field h in equation (16) translates into a chemical potential and the anisotropy parameter γ controls the superconducting pairing terms $c_j^\dagger c_{j+1}^\dagger + \textrm{h.c.}$ Importantly, these pairing terms break the fermion number symmetry, but not their parity since they create/destroy particles in pairs.

After the Jordan–Wigner transformation, the boundary term $H_\mathrm{b}$ reads as

$$\begin{align} H_\mathrm{b} = \frac{1}{2}P\left[\left(c_L^\dagger c_{1}+c_{1}^\dagger c_L\right)+\gamma \left(c_L^\dagger c_{1}^\dagger+c_{1}c_L\right)\right]. \end{align} \tag{ 23 }$$

We observe that it includes the parity operator P and, therefore, depending on the parity sector, yields periodic or antiperiodic boundary conditions in the Kitaev chain.

For OBC, there is no boundary term $H_\mathrm{b}$ in the XY spin chain Hamiltonian (16), which is a quadratic form in terms of the fermionic operators c_j and $c_j^\dagger$, as we have seen in equation (22). In order to simplify the calculations and to better understand some physical results later, it is very useful to introduce the Majorana fermionic operators [58]

$$\begin{align} \check a_{2j-1} = c_j^\dagger + c_j \text{ and } \check a_{2j} = i\left(c_j^\dagger - c_j\right), \end{align} \tag{ 24 }$$

by doubling the sites of the lattice. These operators are real-valued, $\check a_m^\dagger = \check a_m$, and satisfy the anticommutation relation

$$\begin{align} \left\{\check a_m, \check a_{m^{\prime}}\right\} = 2\delta_{mm^{\prime}}.\end{align} \tag{ 25 }$$

Since the Hamiltonian (22) is quadratic, following [57], it can be diagonalised by performing a Bogoliubov rotation to a new set of fermionic operators η_k and $\eta_k^\dagger$. If we organise the Majorana $\check a_m$ and Bogoliubov operators η_k , $\eta_k^\dagger$ in single vectors,

$$\begin{align} \boldsymbol{\eta} = \begin{pmatrix} \eta_{k_0}^\dagger \\ \eta_{k_0} \\ \substack{\displaystyle\cdot\\ \displaystyle\cdot\\ \displaystyle\cdot\\ \displaystyle\cdot\\ \displaystyle\cdot\\ } \\ \eta_{k_{L-1}}^\dagger \\ \eta_{k_{L-1}} \end{pmatrix},\quad \boldsymbol{\check a} = \begin{pmatrix} \check a_{1}\\ \check a_{2} \\ \substack{\displaystyle\cdot\\ \displaystyle\cdot\\ \displaystyle\cdot\\ \displaystyle\cdot\\ \displaystyle\cdot\\ } \\ \check a_{2L-1} \\ \check a_{2L} \end{pmatrix}, \end{align} \tag{ 26 }$$

then there exists a matrix $\mathcal{R}_{h,\gamma}$ of dimension $2L\times 2L$ that relates η and $\boldsymbol{\check{a}}$

$$\begin{align} \boldsymbol{\eta} = \mathcal{R}_{h,\gamma} \boldsymbol{\check a} \end{align} \tag{ 27 }$$

such that the Hamiltonian (22) is diagonal in terms of the Bogoliubov modes

$$\begin{align} H = \sum_{j = 0}^{L-1} \epsilon_{k_j} \left(\eta_{k_j}^\dagger \eta_{k_j} -\frac{1}{2}\right), \end{align} \tag{ 28 }$$

where

$$\begin{align} \epsilon_k = \sqrt{\left(h-\cos k\right)^2+\gamma^2\sin^2k} \end{align} \tag{ 29 }$$

is the single-particle dispersion relation and the momenta k_j satisfy a particular quantisation condition. In actual applications, we will restrict for simplicity to the line γ = 1 (quantum Ising chain). In appendix B, we report in detail how we obtain the matrix $\mathcal{R}_{h, \gamma}$ in that case.

Therefore, the eigenstates of H are the different configurations $\mathbb{K}$ of occupied Bogoliubov modes

$$\begin{align} \ket{\mathbb{K}} = \prod_{k_j\in\mathbb{K}} \eta_{k_j}^\dagger \ket{\emptyset}, \end{align} \tag{ 30 }$$

where $\ket{\emptyset}$ is the Bogoliubov vacuum defined as $\eta_k\ket{\emptyset} = 0$ for all k satisfying the quantisation condition. In particular, for h < 1, $\eta_{k_0}^\dagger\ket{\emptyset}$ is a boundary bound state; see also appendix B. It can be re-expressed as a Majorana fermionic excitation $\check b \ket{\emptyset}$ localised at one of the edges of the chain, as we illustrate in figure 4, with

$$\begin{align} \check b = \eta_{k_0}+\eta_{k_0}^\dagger = \sum_{m = 1}^{2L} \beta_m\check a_m \end{align} \tag{ 31 }$$

and the coefficients β_m are determined by the Bogoliubov rotation $\mathcal{R}_{h,\gamma}$ of equation (27) that relates η and $\boldsymbol{\check a}$. One can check that the energy of this boundary mode $\check b$ tends to zero as $L\to\infty$. Therefore, in the thermodynamic limit, the states $\ket{\emptyset}$ and $\check b\ket{\emptyset}$ are degenerate and the ground state

$$\begin{align} \ket{\mathrm{GS}} = \frac{\ket{\emptyset}\pm \check b\ket{\emptyset}}{\sqrt{2}} \end{align} \tag{ 32 }$$

spontaneously breaks the $\mathbb{Z}_2$ parity symmetry since $\langle \sigma_j^x\rangle\neq 0$.

Zoom In Zoom Out Reset image size

The state (32) is a linear combination of Slater determinants that in general does not satisfy Wick theorem. This means that its reduced density matrix ρ_A is not Gaussian and the usual methods to univocally determine it in terms of the two-point fermionic correlations [59], see also appendix C, cannot be applied in this case. For example, let us take h = 0 and γ = 1. One of the two degenerate symmetry-broken ground states (32) is

$$\begin{align} \ket{\mathrm{GS}; h = 0, \gamma = 1} = \frac{1 + \check a_1}{\sqrt 2}\ket{\emptyset} = \ket{\rightarrow \rightarrow \rightarrow \cdots}. \end{align} \tag{ 33 }$$

In this configuration, it is easy to see that a non-zero order parameter implies a non-zero odd-point fermionic correlation function

$$\begin{align} \langle \sigma_j^x \rangle_{h = 0, \gamma = 1} = \langle \check a_1 \rangle_{h = 0, \gamma = 1} = 1, \end{align} \tag{ 34 }$$

and, therefore, the Wick theorem is not satisfied.

4.1. Generalised Wick theorem

For symmetry-breaking states of the form (32), it is actually possible to obtain a generalised version of the Wick theorem that allows any correlator to be expressed in terms of the one-point and two-point fermionic correlation functions. The rules that we must apply depend on the behaviour of the operators under the parity transformation (18).

For even operators $[\mathcal{O}_e, P] = 0$, which correspond to even-point fermionic correlators, we have

$$\begin{align} \langle \mathcal{O}_e \rangle = \frac{1}{2}\bra{\emptyset}\mathcal{O}_e\ket{\emptyset}+\frac{1}{2}\bra{\emptyset}\check b \mathcal{O}_e \check b \ket{\emptyset}, \end{align} \tag{ 35 }$$

where we take into account that the odd-point functions $\bra{\emptyset} \check b \mathcal{O}_e\ket{\emptyset} = \bra{\emptyset} \mathcal{O}_e \check b \ket{\emptyset} = 0$ since $\ket{\emptyset}$ is a Slater determinant and we can apply the Wick theorem. Using the definition (31) of the boundary mode $\check b$, we find

$$\begin{align} \langle \mathcal{O}_e\rangle = \frac{1}{2}\bra{\emptyset}\mathcal{O}_e\ket{\emptyset}+\frac{1}{2}\sum_{m, m^{\prime}}\beta_m\beta_{m^{\prime}}\bra{\emptyset} \check a_m \mathcal{O}_e \check a_{m^{\prime}}\ket{\emptyset}. \end{align} \tag{ 36 }$$

We observe that, in the correlators that appear in the right-hand side of this equation, we can apply the usual Wick theorem and they can be expressed in terms of the two-point fermionic correlation matrix

$$\begin{align} \left(I+i\Gamma_0\right)_{mm^{\prime}} = \bra{\emptyset}\check a_m \check a_{m^{\prime}}\ket{\emptyset}. \end{align} \tag{ 37 }$$

The only extra ingredients are the coefficients β_m , which determine the boundary mode $\check b$ of equation (31) and actually correspond to the one-point correlation functions

$$\begin{align} \langle \check a_m\rangle = \frac{1}{2}\bra{\emptyset}\left\{\check a_m, \check b\right\}\ket{\emptyset} = \beta_m. \end{align} \tag{ 38 }$$

For odd operators, $\{\mathcal{O}_o, P\} = 0$, which correspond to odd-point fermionic correlations, one can check that the presence of the boundary mode makes them non-zero,

$$\begin{align} \langle \mathcal{O}_o \rangle = \frac12 \bra{\emptyset}\left\{\mathcal{O}_o, \check b\right\}\ket{\emptyset} = \frac12 \sum_m \beta_m \bra{\emptyset} \left\{\mathcal{O}_o, \check a_m\right\}\ket{\emptyset}, \end{align} \tag{ 39 }$$

where $\{\mathcal{O}_o, \check a_m\}$ is an even operator and its expectation value in the vacuum $\ket{\emptyset}$ can be computed from the usual Wick theorem.

4.2. Reduced density matrix

In the light of the previous discussion, for states of the form (32), it is in principle possible to use the generalised Wick theorem of equations (36) and (39) for even and odd correlators to reconstruct the reduced density matrix of any subsystem A from the observables with support on it through the expression

$$\begin{align} \rho_A = \sum_{\mathcal{O}\in A}\frac{1}{2^\ell}\langle \mathcal{O}\rangle \mathcal{O}. \end{align} \tag{ 40 }$$

The problem with this formula is that, since the number of operators in a subsystem grows exponentially with its size, one cannot use it to study large subsystems. However, in [60], an effective expression of ρ_A is proposed for certain states that break parity symmetry. This ansatz can be generalised to the ground state of equation (32), for which we find

$$\begin{align} \rho_A = \rho_{A, e}+\frac{1}{2}\langle \check b_A\rangle \left(\check b_A \rho_{A,e} +\rho_{A,e} \check b_A\right), \end{align} \tag{ 41 }$$

provided A is an interval of contiguous spins starting from the edge as in figure 5. In this ansatz, $\rho_{A,e}$ is the even part of ρ_A under parity P, which is given by the sum of two Gaussian density matrices

$$\begin{align} \rho_{A, e} = \frac{1}{2}\left(\rho_0+\tilde \rho_{0}\right), \end{align} \tag{ 42 }$$

where ρ₀ and $\tilde \rho_{0}$ are the reduced density matrices to A of the states $\ket{\emptyset}$ and $\check b\ket{\emptyset}$, respectively. We have numerically checked that, when we take the thermodynamic limit $L\to\infty$ (semi-infinite chain), $\rho_0\approx \tilde{\rho}_0$; therefore, in that case, we assume that both ρ₀ and $\tilde{\rho}_0$ are given by the restriction to A, $\Gamma_0^A$, of the correlation matrix Γ₀ in equation (37). The operator $\check b_A$ in equation (41) is the (renormalised) restriction of the boundary mode $\check b$ of equation (31) to the subsystem A,

$$\begin{align} \check b = \langle \check b_A\rangle \check b_A+ \langle \check b_B\rangle \check b_B. \end{align} \tag{ 43 }$$

One can check that equation (41) for ρ_A correctly reproduces all the expectation values of even and odd operators predicted by the generalised Wick theorem (36) and (39) and, importantly, it allows us to efficiently compute the entanglement asymmetry of the parity, as we will see in section 4.4.

Zoom In Zoom Out Reset image size

4.3. One- and two-point correlation functions

In the ground state: In order to compute the one- and two-point correlation functions of equations (37) and (38) in the ground state (32), the knowledge of the matrix $\mathcal{R}_{h,\gamma}$ of equation (27) that diagonalises the fermionic Hamiltonian (22) is sufficient. In fact, the vacuum $\ket{\emptyset}$ is characterised by the absence of excitations, i.e. $\langle \eta_k^\dagger \eta_k\rangle = 0$, and therefore

$$\begin{align} \langle \boldsymbol{\eta}^\dagger \boldsymbol{\eta} \rangle = \mathcal{C}_0 \text{ with } \mathcal{C}_0 = \mathrm{diag}\left(1, 0, \dots, 1, 0\right). \end{align} \tag{ 44 }$$

Consequently, by expressing the Majorana fermions $\check a_m$ in terms of the Bogoliubov operators η_k , $\eta_k^\dagger$ using the matrix $\mathcal{R}_{h,\gamma}$, one can compute the total system two-point correlation matrix (37) as

$$\begin{align} I + i\Gamma_0 = \langle \boldsymbol{\check a}^\dagger \boldsymbol{\check a} \rangle = \mathcal{R}_{h,\gamma}^{-1} \mathcal{C}_0 \mathcal{R}_{h,\gamma} . \end{align} \tag{ 45 }$$

On the other hand, since $\check b = \eta_{k_0}^\dagger + \eta_{k_0}$, the one-point fermionic functions β_m of equation (38) are given by the sum of the two first lines of the matrix $\mathcal{R}_{h,\gamma}$ which correspond to the boundary mode according to equation (27) or, alternatively,

$$\begin{align} \beta_m = 2\textrm{Re}\left[\left(\mathcal{R}_{h,\gamma}\right)_{1,m}\right]. \end{align} \tag{ 46 }$$

After a quench: Once we let the ground state (32) for $(h_0, \gamma_0)$, with $h_0\lt1$, unitarily evolve with the XY Hamiltonian (16) of another set of couplings $(h_\mathrm{f}, \gamma_\mathrm{f})$, both the one- and two-point correlation functions $\langle \check a_m(t)\rangle$ and $\langle \check a_m(t) \check a_{m^{\prime}}(t)\rangle$ can be computed numerically using the Heisenberg picture. The Majorana operators $\check a_m$ can be evolved in time using the change of basis that relates the Bogoliubov modes that diagonalise the pre- and post-quench Hamiltonians. Since the time evolution in the Heisenberg picture of the post-quench Bogoliubov operators $\tilde{\eta}_k$ and $\tilde{\eta}_k^\dagger$ is determined by

$$\begin{align} \boldsymbol{\tilde\eta}\left(t\right) = \mathcal{U}\left(t\right) \boldsymbol{\tilde\eta}\left(0\right) \text{with }\mathcal{U}\left(t\right) = \mathrm{diag}\left(\mathrm{e}^{-\mathrm{i}\tilde{\epsilon}_{k_0}t}, \mathrm{e}^{\mathrm{i}\tilde{\epsilon}_{k_0}t}, \dots, \mathrm{e}^{-\mathrm{i}\tilde{\epsilon}_{k_{L-1}}t}, \mathrm{e}^{\mathrm{i}\tilde{\epsilon}_{k_{L-1}}t}\right), \end{align} \tag{ 47 }$$

we can then apply successive changes of basis to compute the dynamics of the correlation matrix (37),

$$\begin{align} I +i\Gamma_0\left(t\right) = \mathcal{R}_{h_\mathrm{f},\gamma_\mathrm{f}}^{-1} \mathcal{U}^{-1} \left(t\right)\mathcal{R}_{h_0, \gamma_0; h_\mathrm{f},\gamma_\mathrm{f}}^{-1} \mathcal{C}_0 \mathcal{R}_{h_0\gamma_0; h_\mathrm{f},\gamma_\mathrm{f}}\mathcal{U}\left(t\right)\mathcal{R}_{h_\mathrm{f}, \gamma_\mathrm{f}} , \end{align} \tag{ 48 }$$

where $\mathcal{R}_{h_0, \gamma0; h_\mathrm{f}, \gamma_\mathrm{f}}$ is the matrix that expresses the modes of the pre-quench Hamiltonian in terms of the new ones $\mathcal{R}_{h_0, \gamma_0; h_\mathrm{f}, \gamma_\mathrm{f}} = \mathcal{R}_{h_0, \gamma_0}\mathcal{R}_{h_\mathrm{f}, \gamma_\mathrm{f}}^{-1}$. Similarly, we can get the one-point functions $\beta_m(t)$ after the quench with the equation

$$\begin{align} \beta_m\left(t\right) = \sum_{m^{\prime}}2\textrm{Re}\left[\left(\mathcal{R}_{h_0,\gamma_0}\right)_{1,m^{\prime}} \left(\mathcal{R}_{h_\mathrm{f}, \gamma_\mathrm{f}}^{-1} \mathcal{U}\left(t\right)\mathcal{R}_{h_\mathrm{f}, \gamma_\mathrm{f}}\right)_{m^{\prime},m}\right]. \end{align} \tag{ 49 }$$

In this way, we have access to all the information on the dynamics of the symmetry-broken state (32) after the quench.

4.4. Entanglement asymmetry

Using the results of the previous subsections, we can compute the entanglement asymmetry of the spin parity after a quench from the state (32) for a subsystem attached to the boundary of the chain, as in figure 5.

In general, according to its definition (4), the calculation of the Rényi entanglement asymmetry requires us to determine the moments of ρ_A and its projection $\rho_{A, P}$ over the charged sectors of the parity operator P in equation (18). In the eigenbasis of the parity in the subsystem A, P_A , the reduced density matrix ρ_A takes the form

Zoom In Zoom Out Reset image size

where the charge sectors $\rho_{++}$ and $\rho_{-}$ are even, i.e. they commute with the parity operator P_A , while the off-diagonal blocks $\rho_{+-}$ and $\rho_{-+}$ are odd, i.e. they anticommute with it. Thus, we can always write

$$\begin{align} \rho_A = \rho_{A, e} + \rho_{A, o} \text{ and } \rho_{A, P} = \rho_{A, e}. \end{align} \tag{ 51 }$$

Therefore, for the parity operator, the moments of $\rho_{A, P}$ are equal to the ones of the even part of ρ_A , ${\mathrm{Tr}}(\rho_{A, P}^n) = {\mathrm{Tr}}(\rho_{A, e}^n)$, which is a Gaussian operator determined by the restriction to A of the correlation matrix Γ₀. On the other hand, since the trace of odd operators is zero, the moments of ρ_A can be obtained from those of the even part of $\rho_A^n$, which we denote as $\rho_A^n|_e$, such as ${\mathrm{Tr}}(\rho_A^n) = {\mathrm{Tr}}(\rho_A^n|_e)$. Therefore, the definition of the Rényi entanglement asymmetry in equation (4) can be re-expressed in the form

$$\begin{align} \Delta S_A^{\left(n\right)} = \frac{1}{1-n}\left(\log {\mathrm{Tr}}\left(\rho_{A, e}^n\right)-\log {\mathrm{Tr}}\left(\rho_{A}^n|_e\right)\right). \end{align} \tag{ 52 }$$

In appendix D.1, we employ this result together with the ansatz of equation (41) for ρ_A to obtain an expression for the entanglement asymmetry of an interval attached to one of the boundaries in terms of the one- and two-point correlation functions β_m and Γ₀ seen before. The final formula reads

$$\begin{align} \Delta S_A^{\left(n\right)} = \frac{1}{n-1}\log\left[\frac{1}{2^n}\sum_{\left\{\tau_i = 0, 1\right\}}\zeta\left[\left\{\tau_i\right\}\right]\frac{\left\{\Gamma_{\tau_1}^A, \dots, \Gamma_{\tau_n}^A\right\}}{\left\{\Gamma_0^A, \overset{n}{\dots}, \Gamma_0^A\right\}}\right], \end{align} \tag{ 53 }$$

where $\Gamma_0^A$ denotes the restriction to the subsystem A of the matrix Γ₀, and $\Gamma_1^A$ is a $2\ell\times 2\ell$ matrix with entries

$$\begin{align} \left(\Gamma_1^A\right)_{mm^{\prime}} = \left(\Gamma_0^A\right)_{mm^{\prime}}+\frac{2}{\langle \check b_A\rangle^2}\sum_{l = 1}^{2\ell}\beta_l\left[\left(\Gamma_0^A\right)_{lm}\beta_{m^{\prime}}-\left(\Gamma_0^A\right)_{lm^{\prime}}\beta_m\right], \end{align} \tag{ 54 }$$

and

$$\begin{align} \zeta \left[\left\{\tau_i\right\}\right] = \sum_{\substack{\left\{w_i = 0, 1\right\} \\ w_0 = w_n = 0}}\prod_{i = 1}^n \left(\langle\check b_A\rangle^{|w_i -w_{i-1}|} \left(1 + \left(-1\right)^{\tau_i}\left(1-w_i-w_{i-1}\right)\right) \right). \end{align} \tag{ 55 }$$

The brackets $\{\Gamma_{\tau_1}, \dots, \Gamma_{\tau_n}\}$ stand for the composition of the correlation matrices of Gaussian operators; see appendix C for its precise definition. While equation (53) becomes very cumbersome as the Rényi index n increases, it simplifies to a rather simple expression for the small ones; for example, for n = 2, we find

$$\begin{align} \Delta S_A^{\left(2\right)} = \log\left[1+\frac{\langle \check b_A\rangle^2}{2}+\frac{\langle \check b_A\rangle^2}{2}\frac{\left\{\Gamma_0^A, \Gamma_1^A\right\}}{\left\{\Gamma_0^A, \Gamma_0^A\right\}}\right].\end{align} \tag{ 56 }$$

We can check that equations (53) and (56) are correct by comparing them with the results obtained using exact diagonalisation, as we do in figure 6 for Rényi indexes n = 2 and n = 5 in two different quenches. The points are the exact numerical values of $\Delta S_A^{(n)}$ calculated using the reduced density matrix ρ_A determined from the exact diagonalisation of the XY spin chain Hamiltonian (16) with OBC. The continuous lines correspond to equation (53) for n = 5 and equation (56) in the case n = 2.

Zoom In Zoom Out Reset image size

Of course, the exact diagonalisation approach is only efficient for small systems, such as the ones considered in figure 6. For larger sizes, one must resort to the expression obtained in equation (53) in terms of the one- and two-point correlations. In what follows, we discuss the results that we find for the entanglement asymmetry using it and taking $L\to\infty$ (semi-infinite line) in quenches from the symmetry-breaking ground state (32) in the ferromagnetic region to different XY Hamiltonians both in the paramagnetic and in the ferromagnetic phases.

Before proceeding, it is interesting to analyse with equation (53) the behaviour of the entanglement asymmetry in the initial ground state (32) as a function of the subsystem size. As figure 5 illustrates, A is an interval of contiguous spins of length $\ell$ attached to the edge of the chain where the boundary mode $\check b$ of equation (31) is localised. This mode has characteristic length $\ell_{\check b} = -1/(\log h)$. Therefore, we expect that the asymmetry of the interval depends on which fraction of the boundary mode is included in the interval: if $\ell\lt\ell_{\check b}$, then $0\lt \Delta S_A^{(n)}\lt\log 2$ while $\Delta S_A^{(n)}\to \log 2$ when $\ell\gg \ell_{\check b}$. This is confirmed by figure 7, where we plot the entanglement asymmetry (53) as a function of $\ell/\ell_{\check b}$ for the ground state (32) for different values of h, obtaining an excellent overlap between them, signalling that the entanglement asymmetry is a scaling function of $\ell/\ell_{\check b}$.

Zoom In Zoom Out Reset image size

Quenches to the paramagnetic phase: When we quench to the critical point h = 1 (left panel of figure 8) or to the paramagnetic phase h > 1 (right panel of figure 8), we observe that the entanglement asymmetry $\Delta S_A^{(n)}$ tends to zero at large times. This implies that the $\mathbb{Z}_2$ spin parity symmetry, initially broken, is dynamically restored in subsystem A due to the melting of the boundary mode. As evident in figure 8, the main feature is that the entanglement asymmetry does not instantaneously decrease after the quench, but it describes a plateau until times proportional to the subsystem length $\ell$.

Zoom In Zoom Out Reset image size

This result has a very simple interpretation: when quenching to the paramagnetic phase, the boundary mode of the pre-quench Hamiltonian can be expressed as a sum of the quasi-particle excitations of the post-quench Hamiltonian, which all have a non-zero speed of propagation. Consequently, all these modes end up escaping the subsystem. Therefore, the time at which $\Delta S_A^{(n)}$ starts to drop is roughly the time necessary for the fastest excitations to reach the end of the subsystem; that is, $t\sim \ell/v_\textrm{max}$, where $v_\textrm{max} = \max_k|\tilde{\epsilon}_k^{^{\prime}}|$ is the maximal velocity of propagation of the post-quench Bogoliubov modes, which for γ = 1 and $h\unicode{x2A7E} 1$ is one.

In figure 9, we study the dependence of the results on the Rényi index n. From that plot, it is clear that it does not have great relevance as the entanglement asymmetry essentially shows the same behaviour when varying its value.

Zoom In Zoom Out Reset image size

Quenches to the ferromagnetic phase: If we quench to the ferromagnetic phase, h < 1, the spin parity symmetry is only partially restored, as shown in figure 10. In the left panel, we observe that the farther the initial magnetic field h₀ is from the final one $h_\mathrm{f}$, the smaller the stationary value of $\Delta S_A^{(n)}$ is, i.e. the more the symmetry is restored. This phenomenon is due to the existence of a boundary mode in the post-quench Hamiltonian. In fact, the pre-quench boundary mode $\check b$, defined in equation (31), can be expressed in terms of the post-quench Bogoliubov excitations as

$$\begin{align} \check b = \alpha_0 \check b_\mathrm{f}+\sum_{j\neq 0}\left(\alpha_j\eta_{k_j}+\alpha_j^*\eta_{k_j}^\dagger\right). \end{align} \tag{ 57 }$$

Therefore, it has a non-zero overlap α₀ with the boundary mode $\check b_\mathrm{f}$ of the post-quench Hamiltonian, which has zero velocity. Consequently, part of the initial boundary mode, which is responsible for the parity symmetry breaking at t = 0, stays in the subsystem and the symmetry is only partially restored. In the right panel of figure 10, we consider different quenches from the same initial state. We observe that the symmetry is more restored in A as we quench closer to the critical point, where it is fully restored. This can be explained by the fact that the length of the post-quench boundary mode diverges as we quench closer to $h_\mathrm{f} = 1$ and, moreover, the overlap α₀ diminishes too.

Zoom In Zoom Out Reset image size

In this section, we study the time evolution after a quench of the entanglement asymmetry of spin parity in a XY spin chain (16) with PBCs. In this case, the presence of the boundary term H_b , given by equation (23), in the fermionic Hamiltonian (22) makes it non-quadratic. However, using the fact that the Hamiltonian (16) is invariant under parity P transformations (18), we can solve this issue by separating the Hilbert space into two sectors, called Neveu–Schwarz (NS) for positive parity and Ramond (R) for the negative one. Then, the XY Hamiltonian (16) can be written as

$$\begin{align} H = \frac{1+P}{2}H_\mathrm{NS}\frac{1+P}{2} + \frac{1-P}{2}H_\mathrm{R}\frac{1-P}{2}, \end{align} \tag{ 58 }$$

where $H_\mathrm{NS}$ and $H_\mathrm{R}$ are the fermionic Hamiltonian (22) with antiperiodic and periodic boundary conditions, respectively, according to equation (23). Note that both commute with the parity P. Therefore, if we are restricted to one of these sectors, the Hamiltonian is quadratic and, by performing a Bogoliubov rotation (27) as in the case of OBCs, we can diagonalise them. For the NS sector, we find

$$\begin{align} H_\mathrm{NS} = \sum_{j}\epsilon_{k_j}\left(\eta_{k_j}^\dagger \eta_{k_j}-\frac{1}{2}\right),\end{align} \tag{ 59 }$$

where ε_k is the one-particle dispersion relation introduced in equation (29). In this case, the momenta k_j satisfy the quantisation condition

$$\begin{align} k_j = \frac{2\pi\left(j-\frac{1}{2}\right)}{L},\quad j = 0, \dots, L-1. \end{align} \tag{ 60 }$$

For the R sector, we have

$$\begin{align} H_\mathrm{R} = \sum_{j\neq 0}\epsilon_{p_j}\left(\eta_{p_j}^\dagger\eta_{p_j}-\frac{1}{2}\right)-2\left(1-h\right)\left(\eta_{0}^\dagger\eta_0-\frac{1}{2}\right) \end{align} \tag{ 61 }$$

with the momenta p_j satisfying

$$\begin{align} p_j = \frac{2\pi j}{L},\quad j = 0, \dots, L-1.\end{align} \tag{ 62 }$$

These two sectors define two different vacuums, $\ket{\emptyset}_\mathrm{NS}$ and $\ket{\emptyset}_\mathrm{R}$, such that

$$\begin{align} \eta_{k_j} \ket{\emptyset}_\mathrm{NS} = 0 \text{ and } \eta_{p_j} \ket{\emptyset}_\mathrm{R} = 0, \forall j. \end{align} \tag{ 63 }$$

If the chain has a finite size, then its ground state is the vacuum $\ket{\emptyset}_\mathrm{NS}$. In the thermodynamic limit $L\to\infty$, as we already mentioned, the states $\ket{\emptyset}_\mathrm{NS}$ and $\eta_0^\dagger \ket{\emptyset}_\mathrm{R}$ become degenerate in the ferromagnetic phase h < 1 and the ground state is the linear superposition

$$\begin{align} \ket{\Psi} = \frac{\ket{\emptyset}_\mathrm{NS} + \eta_0^\dagger \ket{\emptyset}_\mathrm{R}}{\sqrt 2}, \end{align} \tag{ 64 }$$

which spontaneously breaks the $\mathbb{Z}_2$ parity symmetry (21). As happens for the ground state (32) with OBC, the state (64) does not satisfy the Wick theorem. However, the main difference is that the odd part in equation (64) consists of an excitation over a vacuum (R vacuum) different from the vacuum of the even part (NS vacuum). Thus, one cannot use the ansatz (41) for ρ_A in OBC. We have to take an alternative route based on cluster decomposition to calculate the entanglement asymmetry in a quench from the state (64).

In the thermodynamic limit, the correlation functions of even operators $\mathcal{O}_e$ are equal for $\ket{\emptyset}_\mathrm{NS}$ and $\eta_0^\dagger\ket{\emptyset}_\mathrm{R}$ [28]. Therefore, the even-point fermionic correlation functions of the state (64) can be calculated using the Wick theorem and, consequently, are fully characterised by the two-point correlation matrix Γ₀ defined in equation (37), which can be determined using the same procedure as in OBC. In the thermodynamic limit, it takes the following form after the quench [28]:

$$\begin{align} \left(\Gamma_0\left(t\right)\right)_{jj^{\prime}} = \int_{-\pi}^{\pi}\frac{\mathrm{d}k}{2\pi} \mathcal{G}\left(k, t\right) \mathrm{e}^{\mathrm{i}k\left(j-j^{\prime}\right)}, \quad j, j^{\prime} = 1, 2,\dots, \end{align} \tag{ 65 }$$

where $\mathcal{G}(k, t)$ is the $2\times 2$ matrix

$$\begin{align} \mathcal{G}\left(k, t\right) = \left(\begin{array}{cc} f\left(k\right) & -g\left(k\right)\\ g\left(-k\right) & -f\left(k\right) \end{array}\right) \end{align} \tag{ 66 }$$

with

$$\begin{align} f\left(k\right) = i\sin\Delta_k\sin\left(2\tilde{\epsilon}_kt\right),\quad g\left(k\right) = \mathrm{e}^{-\mathrm{i}\theta_k^\mathrm{f}}\left(\cos\Delta_k+i\sin\Delta_k\cos\left(2\tilde{\epsilon}_kt\right)\right). \end{align} \tag{ 67 }$$

Here, $\Delta_k$ is the difference $\Delta_k = \theta_k^\mathrm{f}- \theta_k^0$ between the Bogoliubov angles of the post- and pre-quench Hamiltonians with couplings $(h_\mathrm{f}, \gamma_\mathrm{f})$ and $(h_0, \gamma_0)$, respectively, and $\tilde{\epsilon}_k$ is the dispersion relation of the post-quench Hamiltonian. The Bogoliubov angle for arbitrary couplings h, γ is given by

$$\begin{align} \mathrm{e}^{\mathrm{i}\theta_k} = \frac{\cos k - h +i\gamma\sin k}{\sqrt{\left(h-\cos k\right)^2+\gamma^2\sin^2 k}}.\end{align} \tag{ 68 }$$

On the other hand, the calculation of odd-point fermionic correlations in the state (64) is a rather difficult task. However, following the ideas in [61], we can use the cluster decomposition principle to determine them from the even-point ones, exploiting the fact that the latter are equal in both parity sectors in the thermodynamic limit and are given by the matrix Γ₀. In fact, let us take an arbitrary odd operator $\mathcal{O}_o$ with support in the subsystem A. We further consider the auxiliary spin operator $\sigma_r^x$ on the site r that does not belong to the interval A. If we take the even correlation function $\langle \mathcal{O}_o\sigma_r^x\rangle$, then we expect that in the infinite distance separation, $r\to\infty$, the spin at the site r is uncorrelated from the spins in the subsystem A and

$$\begin{align} \langle \mathcal{O}_o \sigma^x_r \rangle \underset{r\to \infty}{ = }\langle \mathcal{O}_o \rangle \langle \sigma_r^x \rangle. \end{align} \tag{ 69 }$$

Therefore,

$$\begin{align} \langle \mathcal{O}_o \rangle = \pm \lim_{r\to \infty} \frac{\langle \mathcal{O}_o \sigma^x_r \rangle}{\sqrt{\langle \sigma^x_1 \sigma^x_r \rangle}}. \end{align} \tag{ 70 }$$

Using this idea, in appendix E, we show that it is possible to obtain the odd part of ρ_A for the state (64). Applying some of the methods to compute the entanglement entropy of disjoint intervals in the XY spin chain [61], we find that ρ_A can be expressed in terms of two Gaussian operators ρ₀ and ρ₁

$$\begin{align} \rho_{A} = \rho_0 + \langle \sigma_1^x \rangle \sigma_1^x \rho_1, \end{align} \tag{ 71 }$$

which are respectively determined by the restriction to A, $\Gamma_0^A$, of the correlation matrix (65) and by $\Gamma_1^A$, defined as

$$\begin{align} \left(I + i\Gamma_1^A\right)_{mm^{\prime}} = \frac{\langle \check a_m \check a_{m^{\prime}} \sigma^x_1\rangle}{\langle \sigma^x_1 \rangle}, \quad m,m^{\prime}\in A. \end{align} \tag{ 72 }$$

As we show in appendix E, the latter can be obtained from Γ₀ using cluster decomposition as a Schur complement.

5.1. Entanglement asymmetry

From the expression (71) of the reduced density matrix ρ_A for the ground state (64), we can derive using equation (52) an efficient formula to calculate the entanglement asymmetry of the spin parity symmetry for a large subsystem. Its obtention involves some tedious algebraic manipulations, which are described in detail in section D.2. The final result reads

$$\begin{align} \Delta S_A^{\left(n\right)} = \frac{1}{n-1}\log \left[\sum_{\left\{\chi_i = 0,1\right\}}\frac{1+\left(-1\right)^{s_n}}{2}\langle \sigma^x \rangle^{s_n} \frac{\left\{P_1^{s_1}\Gamma_{\chi_1}^AP_1^{s_1}, \ldots, P_{1}^{s_n}\Gamma_{\chi_i}^AP_{1}^{s_n}\right\}}{\left\{\Gamma_0^A, \overset{n}{\dots}, \Gamma_0^A\right\}}\right], \end{align} \tag{ 73 }$$

where $\Gamma_0^A$ and $\Gamma_1^A$ are the restrictions to A of the correlation matrices (65) and (72), respectively, $s_n = \sum_{j = 1}^n\chi_j$ and P₁ is the diagonal matrix $P_1 = \textrm{diag}(-1, 1, \ldots, 1)$. The Rényi-2 asymmetry has a simpler expression

$$\begin{align} \Delta S_A^{\left(2\right)} = \log\left[1+\frac{\left\{\Gamma_1^A, P_1\Gamma_1^AP_1\right\}}{\left\{\Gamma_0^A, \Gamma_0^A\right\}}\right].\end{align} \tag{ 74 }$$

As we did for OBC, we check that the formula of equation (73) produces the correct entanglement asymmetry by comparing with calculations done using exact diagonalisation for small systems. In figure 11, we analyse the time evolution of $\Delta S_A^{(n)}$ for n = 7 in a particular quench from the state (64). The continuous line corresponds to equation (73), while the points have been obtained with exact diagonalisation for two different total sizes L of the chain. We observe that, while the curve and the points match just after the quench, they deviate at large times. The reason is that equation (73) gives $\Delta S_A^{(n)}$ in the thermodynamic limit $L\to\infty$. In fact, recall that this expression relies on the cluster decomposition between the spins of subsystem A and an auxiliary spin. Therefore, it is expected that our method only agrees with the exact diagonalisation results at early times and small correlation lengths, before finite-size effects, as the revivals of figure 11, appear.

Zoom In Zoom Out Reset image size

The rest of this section is devoted to discussing the results that we obtain for the entanglement asymmetry after quenches from the symmetry-breaking ground state (64) in the ferromagnetic phase to different points of the $(h,\gamma)$-plane of the XY spin chain. We first consider the case of a subsystem made of a single site in order to compare the entanglement asymmetry with the usual order parameter $\langle \sigma_1^x\rangle$. We will then take larger subsystem sizes and obtain an analytic expression that describes the time evolution of the asymmetry at leading order.

5.1.1. One-site subsystems—comparison with the order parameter.

As we already pointed out in section 3, in the thermodynamic limit, the breaking of the $\mathbb{Z}_2$ symmetry associated with spin parity can be detected with the standard order parameter $\langle \sigma_1^x\rangle$. It is then interesting to compare the results for $\langle \sigma_1^x\rangle$ with those of the asymmetry for a subsystem of one site. By cluster decomposition, we can indeed get $\langle \sigma_1^x \rangle^2$ by computing the expectation value of a string of Majorana fermions as described in [28].

According to equation (40), the reduced density matrix of a one-site subsystem is simply

$$\begin{align} \rho_A = \frac12 \left[ I + \langle \sigma_1^z \rangle \sigma_1^z + \langle \sigma_1^x \rangle \sigma_1^x + \langle \sigma_1^y \rangle \sigma_1^y\right] \end{align} \tag{ 75 }$$

and its projection on the $\mathbb{Z}_2$ parity sectors is

$$\begin{align} \rho_{A, P} = \frac12 \left[ I + \langle \sigma_1^z \rangle \sigma_1^z\right]. \end{align} \tag{ 76 }$$

Using now the anticommutation relations of Pauli matrices, one can derive

$$\begin{align} \Delta S^{\left(n\right)}_A = H^{\left(n\right)}\left(\langle \sigma_1^z \rangle\right) - H^{\left(n\right)}\left(\sqrt{\langle \sigma_1^x \rangle^2+ \langle \sigma_1^y \rangle^2+ \langle \sigma_1^z \rangle^2}\right), \end{align} \tag{ 77 }$$

where the function $H^{(n)}(x)$ is defined in equation (12). In this expression, we can take the limit n → 1 according to equation (13). We observe that the asymmetry of one site is not exactly equivalent to the order parameter $\langle \sigma_1^x\rangle$ since it takes into account the symmetry breaking in both directions x and y. In fact, $\Delta S_A^{(n)} = 0$ if and only if both $\langle \sigma_1^x \rangle$ and $\langle \sigma_1^y \rangle$ are zero. In the left panel of figure 12, we compare both quantities for a particular quench. If we take logarithmic scale in the y axis, as we do in the right panel of figure 12, their differences are more obvious. We can see that, while the order parameter $\langle \sigma_1^x\rangle^2$ periodically vanishes at certain times due to the oscillations between the x and y components of the spin, the asymmetry does not, since it takes both directions of symmetry breaking into account. This is a physically important result that shows the power of the asymmetry compared to traditional observables: the symmetry is never restored at intermediate times, although the order parameter vanishes.

Zoom In Zoom Out Reset image size

Despite the differences between the order parameter $\langle \sigma_1^x\rangle^2$ and the asymmetry, a remarkable result is that the leading term at late times is given by the same exponential decay for both quantities. In [27, 28], it was found that

$$\begin{align} \langle \sigma_1^x \left(t\right) \rangle^2 \underset{t\to \infty} {\sim}\exp\left(2t\int_{0}^\pi \frac{\mathrm{d}k}{\pi} |\tilde{\epsilon}_k^{\prime}| \log |\cos \Delta_k |\right), \end{align} \tag{ 78 }$$

which corresponds to the solid line in the right panel of figure 12.

5.1.2. Larger subsystems.

In figure 13, we study the time evolution of the entanglement asymmetry for larger subsystems in different quenches starting from an initial state of the form (64) that spontaneously breaks the $\mathbb{Z}_2$ spin parity symmetry and letting it evolve with an XY spin chain Hamiltonian with different couplings h and γ both in the paramagnetic and in the ferromagnetic phases. In all cases, we observe a plateau in the asymmetry at early times after the quench, followed by an exponential decay. In fact, in the upper left panel of figure 13, we study a particular quench to the paramagnetic phase for two different subsystem lengths; in the upper right panel, we consider the same quench, but taking logarithmic scale in the y axis to confirm the exponential relaxation to zero of the entanglement asymmetry. In the case of OBC, we explain the existence of the plateau in terms of the boundary mode present in the subsystem A. In PBC, even if we do not have the boundary mode picture, by analogy, we can consider that the quasi-particle excitations generated after the quench emerge from the two end-points of the subsystem. This would explain why the length of the plateau is twice shorter, $\sim \ell/(2v_\textrm{max})$, than in OBC (recall that, in that case, the interval is attached to the boundary and the quasi-particles can only emerge from one point).

Zoom In Zoom Out Reset image size

We find that the exponential decay is rather universal for all the quenches. The change of the initial and final couplings $(h_0, \gamma_0)$ and $(h_\mathrm{f},\gamma_\mathrm{f})$ only affects the decay rate, which is independent of the subsystem size $\ell$ and the Rényi index n. This observation is crucial since, as we have seen above, both the asymmetry of a single site and the order parameter $\langle \sigma_1^x\rangle^2$ show the same exponential decay at late times, which is described by equation (78). We can use this fact to conjecture an effective analytic expression for the asymmetry of a subsystem of length $\ell$ after the quench. We find that it is actually well reproduced by the formula

$$\begin{align} \Delta S_A^{\left(n\right)}\left(t\right)\approx \frac{\log 2}{\log C_n}\log\left[1\!+\!\left(C_n-1\right)\mathrm{exp}\left(\int_{0}^\pi \!\!\frac{\mathrm{d}k}{\pi}\max\left(2|\tilde{\epsilon}_k^{\prime}|t\!-\!\ell \!+\!\mathcal{O}\left(1\right), 0\right)\log|\cos\Delta_k|\right)\right], \end{align} \tag{ 79 }$$

where $C_n = 2+(n-2)C_{h_0, \gamma_0; h_\mathrm{f}, \gamma_\mathrm{f}}$ and $C_{h_0, \gamma_0; h_\mathrm{f}, \gamma_\mathrm{f}}$ is an unknown coefficient that depends on the couplings of the pre- and post-quench Hamiltonian. Observe that equation (79) satisfies the required properties, for $t\lt l/2v_{\max}$ it simplifies to $\Delta S_A^{(n)}(0) = \log 2$, while at large times it exponentially decays at leading order as the order parameter (78),

$$\begin{align} \Delta S_A^{\left(n\right)}\left(t\right)\underset{t\to \infty} {\sim}\mathrm{exp}\left(2t\int_{0}^\pi \frac{\mathrm{d}k}{\pi}|\epsilon_k^{\prime}|\log|\cos\Delta_k|\right), \end{align} \tag{ 80 }$$

which is independent of the subsystem length $\ell$ and Rényi index n, and in particular is valid in the limit n → 1. The solid curves in all the plots of figure 13 correspond to the analytic expression (79). For each quench, the value of the unknown coefficient $C_{h_0, \gamma_0;h_\mathrm{f},\gamma_\mathrm{f}}$ has been adjusted by hand to obtain the best agreement with the numerical points calculated using equation (73).

Quantum Mpemba effect? In the case of the U(1) symmetry associated with the transverse magnetisation, a quantum version of the Mpemba effect has been observed, in which the symmetry is restored faster after a quench for the initial states that break it more [22]. For the $\mathbb{Z}_2$ spin parity symmetry, we cannot in principle observe the same effect since the asymmetry at time t = 0 converges to $\log 2$ very fast as the length $\ell$ of the subsystem increases. Therefore, for large enough intervals, the spin parity symmetry is maximally broken for any initial state of the form (64) that we can consider. However, when studying a series of quenches to the same point in the paramagnetic region from different initial states such as the ones plotted in the left panel of figure 14, it is very tempting to conclude that the symmetry is restored faster for the ground state corresponding to $h_0 = 0$ than for $h_0 = 0.2$ or $h_0 = 0.4$, which are parameters closer to the paramagnetic region.

Zoom In Zoom Out Reset image size

On the other hand, it is important to remember that a value $\log 2$ for the asymmetry of a $\ell$-site subsystem means that the symmetry is completely broken at that length scale. If we analyse the entanglement asymmetry of the initial states (64) as a function of the subsystem size, as we do in the middle panel of figure 14, we can see that, for any subsystem length, the entanglement asymmetry for $h_0 = 0.8$ is smaller than for $h_0 = 0.5$ and $h_0 = 0.2$. In the right panel of figure 14, we plot the value of entanglement asymmetry at time t = 10 after performing a quench to the same point $(h_\mathrm{f}, \gamma_\mathrm{f}) = (1.2, 1)$ from the initial states considered in the middle panel. As we see, at t = 10, the points are inverted with respect to the ones at the initial time, and the symmetry is now more broken for the state corresponding to $h_0 = 0.8$ than for $h_0 = 0.5$ and $h_0 = 0.2$. This is a phenomenon somewhat reminiscent of the quantum Mpemba effect found in the U(1) case in the sense that the more the spin parity symmetry is broken at any length scale, the faster it is restored at any length scale.

In this manuscript, we extend to finite discrete cyclic groups $\mathbb{Z}_N$ the study of the recently introduced entanglement asymmetry. Contrary to the U(1) symmetry considered in [22], the entanglement asymmetry is bound by the dimension of the group. As an application, we have examined the entanglement asymmetry of the $\mathbb{Z}_2$ group generated by the spin parity in the XY spin chain. In particular, we have analysed its time evolution after a quench from the ferromagnetic phase, where this symmetry is spontaneously broken by the ground state in the thermodynamic limit. This extends the analysis carried out in [27, 28] using the standard order parameter. While the order parameter can be employed as a measure of symmetry breaking for one-site subsystems, the longer range correlations that break the symmetry in extended subsystems are only taken into account in the asymmetry, which is then a proper quantifier of symmetry breaking for larger subsystems. Furthermore, we show that the asymmetry is more suited to detecting symmetry breaking, since we found examples in which the order parameter vanishes when the symmetry is still broken.

We were faced with the technical issue that the states that break this $\mathbb{Z}_2$ symmetry do not satisfy the Wick theorem, and that the powerful techniques for Gaussian reduced density matrices employed in the case of the U(1) symmetry in [22, 23] cannot be directly applied here. We solved the problem of obtaining the reduced density matrix by using a generalised version of the Wick theorem that takes into account the one-point functions in OBCs and by applying the cluster decomposition principle when the chain has PBCs.

According to our results, the entanglement asymmetry rapidly saturates for large subsystems to its maximal value, $\log 2$, in all of the ferromagnetic phase, being zero in the paramagnetic region; therefore, it acts as a topological indicator of symmetry breaking for large subsystems. The entanglement asymmetry also allows us to determine when and how the $\mathbb{Z}_2$ parity symmetry is restored in a subsystem after a quench to another set of couplings of the XY spin chain. In particular, we have found that, just after the quench, the entanglement asymmetry describes a plateau of length proportional to the subsystem size and then it suddenly decreases its value. For OBCs, if the quench is to the paramagnetic region, the asymmetry of an interval starting from the boundary drops to zero at large times and the symmetry is fully restored, while in a quench to the ferromagnetic phase, it tends to a non-zero value and the symmetry is only partially restored due to the boundary mode of the post-quench Hamiltonian. For PBCs, the symmetry is fully restored in all cases. Remarkably, we have found that, at large times and any Rényi index, the entanglement asymmetry shows the same exponential decay to zero as the order parameter, for which an analytic expression was derived in [27, 28]. From this, we have conjectured an analytical expression for the entanglement asymmetry that effectively reproduces its time evolution after the quench.

We thank Sara Murciano for useful discussions. F F acknowledges the financial support from the Erasmus+ program and the excellence IDEX grant of Université Paris-Saclay. P C and F A acknowledge support from ERC under Consolidator Grant Number 771536 (NEMO).

In this appendix, we prove that, for a finite discrete group G, the Rényi entanglement asymmetry $\Delta S_A^{(n)}$ is bounded by the dimension N of the group. Let us write the projected reduced density matrix $\rho_{A, Q_G}$ of a finite group G of elements $g_1, \dots, g_N$ as

$$\begin{align} \rho_{A, Q_G} = \frac{1}{N}\sum_{j = 1}^N \rho_{j}\ \text{where } \rho_{j} = g_j \rho_A g_j^{-1}. \end{align} \tag{ 81 }$$

The trace of $\rho_{A, Q_G}^n$ is given by

$$\begin{align} {\mathrm{Tr}}\left(\rho_{A, Q_G}^n\right) = \frac{1}{N^n}\sum_{j_1, \ldots, j_n = 1}^N {\mathrm{Tr}}\left(\rho_{j_1}\cdots \rho_{j_n}\right). \end{align} \tag{ 82 }$$

Since the ρ_j are all positive semi-definite,

$$\begin{align} {\textrm{Re}}\left[{\mathrm{Tr}}\left(\rho_{j_1}\cdots \rho_{j_n}\right)\right]\unicode{x2A7E} 0, \quad \sum_{j_1, \ldots, j_n = 1}^N \mathrm{Im}\left[{\mathrm{Tr}}\left(\rho_{j_1}\cdots \rho_{j_n}\right)\right] = 0, \end{align} \tag{ 83 }$$

and they have the same spectrum ${\mathrm{Tr}}(\rho_{j}^n) = {\mathrm{Tr}}(\rho_A^n)$, one consequently has

$$\begin{align} \sum_{j_1, \ldots, j_n = 1}^N {\mathrm{Tr}}\left(\rho_{j_1}\cdots \rho_{j_n}\right)\unicode{x2A7E} N {\mathrm{Tr}}\left(\rho_A^n\right). \end{align} \tag{ 84 }$$

Applying this result in the definition (4) of the Rényi entanglement asymmetry, one finds that

$$\begin{align} \Delta S_A^{\left(n\right)}\unicode{x2A7D} \log N. \end{align} \tag{ 85 }$$

If the state ρ_A maximally breaks the symmetry, i.e. the ρ_j have orthogonal support,

$$\begin{align} {\mathrm{Tr}}\left(\rho_{j} \rho_{j^{\prime}}\right) = \delta_{jj^{\prime}} {\mathrm{Tr}}\left(\rho_A^2\right), \end{align} \tag{ 86 }$$

then if two indices differ in the product (84) of ρ_j , by the cyclic invariance of the trace we can always write using the Cauchy–Schwarz inequality,

$$\begin{align} 0\unicode{x2A7D} {\mathrm{Tr}}\left(\rho_{j_1}\cdots \rho_{j_n}\right) \unicode{x2A7D} {\mathrm{Tr}}\left(\rho_{j_1} \cdots \rho_{j_{n-2}}\right){\mathrm{Tr}}\left(\rho_{j_{n-1}}\rho_{j_n}\right) = 0, \quad j_{n-1}\neq j_n, \end{align} \tag{ 87 }$$

and the bound of equation (85) is then saturated.

As we discuss in section 4.4, in order to diagonalise the Hamiltonian (16) for OBCs, it is sufficient to know the rotation matrix $\mathcal{R}_{h,\gamma}$ between the Majorana and Bogoliubov fermions introduced in equation (27). In the quantum Ising line, γ = 1, we obtain it following the detailed discussion presented in [62]. First, it is necessary to solve the quantisation condition

$$\begin{align} k = \frac{\pi}{L+1}\left(j+\frac{1}{\pi} \arctan\left(\frac{\sin k}{\cos k-h}\right)\right)\text{with }k \in \left(0, \pi\right) \text{ and } j\in \left\{1, 2, \dots, L \right\}, \end{align} \tag{ 88 }$$

which has L − 1 real solutions for h < 1 and L solutions for $h\unicode{x2A7E} 1$. For h < 1, the last mode is the boundary mode, given by the complex momentum $k_0 = \pi + i\nu$, where ν is obtained by solving

$$\begin{align} \tanh \left(\nu \left(L+1\right)\right) = \frac{\sinh \nu}{h + \cosh \nu}\quad \text{with }\nu>0. \end{align} \tag{ 89 }$$

Once we have numerically obtained the momenta with equation (88), we can construct the modes η_k that diagonalise the Hamiltonian (22) as

$$\begin{align} \eta_{k}^\dagger & = \frac{1}{2}\sum_{j = 1}^{2L} A_k \left(\sin\left(kj-\theta_k\right) \check a_{2j-1} +i\sin\left(kj\right) \check a_{2j} \right), \end{align} \tag{ 90 }$$

$$\begin{align} \eta_{k} & = \frac{1}{2}\sum_{j = 1}^{2L} A_k \left(\sin\left(kj-\theta_k\right) \check a_{2j-1} -i\sin\left(kj\right) \check a_{2j} \right), \end{align} \tag{ 91 }$$

where A_k is a normalisation constant,

$$\begin{align} A_k = \left(\frac{L}{2} - \frac{h\left(\cos k-h\right)}{2\epsilon_k^2}\right)^{-\frac 1 2}. \end{align} \tag{ 92 }$$

Similarly, the boundary modes are given by

$$\begin{align} \eta_{k_0}^\dagger = \frac{1}{2}\sum_{j = 1}^{2L} A_{k_0} \left(-1\right)^j \left(\sinh\left(\nu \left(L+1-j\right)\right) \check a_{2j-1} +i\sinh\left(\nu j\right) \check a_{2j} \right), \end{align} \tag{ 93 }$$

$$\begin{align} \eta_{k_0} = \frac{1}{2}\sum_{j = 1}^{2L} A_{k_0} \left(-1\right)^j \left(\sinh\left(\nu \left(L+1-j\right)\right) \check a_{2j-1} -i\sinh\left(\nu j\right) \check a_{2j} \right), \end{align} \tag{ 94 }$$

where

$$\begin{align} A_{k_0} = \left(\sum_{j = 1}^{L}\sinh^2\left(\nu \left(L+1-j\right)\right)+\sinh^2\left(\nu j\right) \right)^{-\frac 1 2}. \end{align} \tag{ 95 }$$

In general, a Gaussian density matrix is generated from a quadratic form, such that

$$\begin{align} \rho_W = \frac{1}{\mathcal Z\left(W\right)} \exp\left(\frac{1}{4}\sum_{m,m^{\prime}}W_{mm^{\prime}} \check a_m \check a_{m^{\prime}}\right), \end{align} \tag{ 96 }$$

where the factor $\mathcal{Z}(W)$ ensures that ${\mathrm{Tr}}(\rho_W) = 1$. A Gaussian density matrix ρ_W satisfies the Wick theorem and, therefore, it is univocally determined by the two-point correlation matrix [59]

$$\begin{align} \delta_{mm^{\prime}} + i\Gamma_{mm^{\prime}} = {\mathrm{Tr}}\left(\rho_W \check a_m \check a_{m^{\prime}}\right), \end{align} \tag{ 97 }$$

which can be written as

$$\begin{align} \Gamma = i\tanh\left(\frac{W}{2}\right). \end{align} \tag{ 98 }$$

We can thus use the notation $\rho_W\equiv \rho[\Gamma]$. For large subsystems, it is in practice impossible to directly use full reduced density matrices because their dimension grows exponentially as $2^\ell$ with the subsystem size $\ell$, whereas correlation matrices have a size $2\ell\times 2\ell$.

The computation of the Rényi entanglement asymmetry (4) involves traces of products of different Gaussian matrices ${\mathrm{Tr}}\left(\rho[\Gamma_1] \cdots \rho[\Gamma_n]\right)$. To that end, let us first consider two Gaussian density matrices $\rho[\Gamma]$ and $\rho[\Gamma^{\prime}]$ associated with the correlation matrices Γ and $\Gamma^{\prime}$. The special properties of Gaussian operators allow us to write their product in the form [61, 63]

$$\begin{align} \rho\left[\Gamma\right] \rho\left[\Gamma^{\prime}\right] = {\mathrm{Tr}}\left[\rho\left[\Gamma\right] \rho\left[\Gamma^{\prime}\right]\right] \rho\left[\Gamma \times \Gamma^{\prime}\right], \end{align} \tag{ 99 }$$

where the composition rule × for correlation matrices is defined as

$$\begin{align} \Gamma \times \Gamma^{\prime} : = i\left[I - \left(I+ i \Gamma^{\prime}\right)\left(I - \Gamma \Gamma^{\prime}\right)^{-1} \left(I + i \Gamma\right) \right] \end{align} \tag{ 100 }$$

and the trace of two Gaussian density matrices can be obtained from their respective correlation matrices as

$$\begin{align} \left\{\Gamma, \Gamma^{\prime}\right\}: = {\mathrm{Tr}}\left(\rho\left[\Gamma\right]\rho\left[\Gamma^{\prime}\right]\right) = \prod_{j} \frac{1-\mu_j}{2} = \pm \sqrt{\det \left(\frac{I - \Gamma \Gamma^{\prime}}{2}\right)} \end{align} \tag{ 101 }$$

where µ_j are the eigenvalues of $\Gamma \Gamma^{\prime}$ with halved degeneracy. Observe that, in principle, there is an ambiguity in the global sign in the last equality of the expression above, which is fixed in our case since the traces we compute always have a positive sign.

Using these rules, any trace of a product of Gaussian operators can be expressed as

$$\begin{align} {\mathrm{Tr}}\left[\rho_1 \cdots \rho_n\right] = \left\{\Gamma_1, \ldots, \Gamma_n\right\}, \end{align} \tag{ 102 }$$

where

$$\begin{align} \left\{\Gamma_1, \ldots, \Gamma_n\right\} = \left\{\Gamma_1, \Gamma_2\right\} \left\{\Gamma_1 \times \Gamma_2, \Gamma_3, \ldots, \Gamma_n\right\} \end{align} \tag{ 103 }$$

can be computed by recurrence.

In this appendix, we describe in detail how we obtain the formulas of equations (53) and (73) that we use in the main text to efficiently compute the time evolution of entanglement asymmetry after a quench in OBC and PBC, respectively.

D.1. In OBC

According to equation (52), in order to calculate the Rényi entanglement asymmetry we must determine the traces ${\mathrm{Tr}}(\rho_{A, e}^n)$ and ${\mathrm{Tr}}(\rho_{A}^n|_{e})$. In section 4, we have seen that in OBC the reduced density matrix ρ_A is given by the ansatz of equation (41). Therefore, in this case, the even and odd parts of ρ_A are

$$\begin{align} \rho_{A, e} = \rho_0, \quad \rho_{A, o} = \frac12 \langle \check b_A \rangle \left(\check b_A \rho_0 + \rho_0 \check b_A \right). \end{align} \tag{ 104 }$$

The nth power of ρ_A can be then obtained as

$$\begin{align} \rho_A^n = \sum_{\left\{\pi_j = e, o\right\}}\prod_{j = 1}^n \rho_{A, \pi_j}. \end{align} \tag{ 105 }$$

Therefore, its projection on the even parity sector, $\rho_A^n|_e$, only includes the terms in the sum that have an even number of times the factor $\rho_{A, o}$.

Let us now introduce the variable χ_j corresponding to the parity of each factor, i.e. $\chi_j = 0$ for $\pi_j = e$ and $\chi_j = 1$ for $\pi_j = o$, such that the even and odd parts (104) of ρ_A can be expressed in the same foot,

$$\begin{align} \rho_{A, \pi_j} = \frac{1}{2} \langle \check b_A \rangle^{\chi_j} \check b_A^{\chi_j} \left(\rho_0 + \check b_A^{\chi_j} \rho_0 \check b_A^{\chi_j}\right), \end{align} \tag{ 106 }$$

where

$$\begin{align} \check b_A^0 \equiv I, \quad \check b_A^1 \equiv \check b_A. \end{align} \tag{ 107 }$$

Observe that, in a product of several matrices $\rho_{A, \pi_j}$, each factor can be written as $(\rho_0 + \check b_A^{\chi_j} \rho_0 \check b_A^{\chi_j})$ with an extra Majorana fermion $\langle \check b_A\rangle \check b_A$ inserted on the left if it corresponds to $\rho_{A, o}$. Using now the fact that $\check b_A^2 \equiv I$, we can rewrite the product, e.g.

$$\begin{align} \rho_{A, o} \rho_{A, e} \rho_{A, o} = \langle \check b_A\rangle^2 \check b_A\left( \frac12 \rho_0 +\frac12 \check b_A \rho_0 \check b_A\right)\rho_0 \check b_A \left(\frac12 \rho_0 + \frac12 \check b_A \rho_0 \check b_A\right) \end{align} \tag{ 108 }$$

in the form

$$\begin{align} \rho_{A, o} \rho_{A, e} \rho_{A, o} = \langle \check b_A\rangle^2 \left[\check b_A\left( \frac12 \rho_0 +\frac12 \check b_A \rho_0 \check b_A\right) \check b_A\right] \left[\check b_A\rho_0 \check b_A\right] \left(\frac12 \rho_0 + \frac12 \check b_A \rho_0 \check b_A\right). \end{align} \tag{ 109 }$$

That is, whenever we have a product of the matrices $\rho_{A, e}$ and $\rho_{A, o}$ that contains an even number of terms $\rho_{A, o}$, we can rewrite it as

$$\begin{align} \prod_{j = 1}^n \rho_{A, \pi_j} = \frac{1}{2^n}\prod_{j = 1}^n \check b_A^{s_j} \left(\rho_0 + \check b_A^{\chi_j} \rho_0 \check b_A^{\chi_j}\right) b_A^{s_j}, \end{align} \tag{ 110 }$$

where $s_j = \sum_{j^{\prime} = 1}^j \chi_j$. Observe that we insert the operator $\check b_A$ on the left and the right of the j factor $(\rho_0 + \check b_A^{\chi_j} \rho_0 \check b_A^{\chi_j})$ if until that position there is an odd number of terms $\rho_{A,o}$. Consequently, it is possible to write the even part of $\rho_A^n$ as

$$\begin{align} \rho_A^n|_e = \frac{1}{2^n} \sum_{\left\{\chi_i = 0, 1\right\}} \frac{1+\left(-1\right)^{s_n}}{2} \langle \check b_A \rangle^{s_n} \prod_{j = 1}^n \check b^{s_j}\left(\rho_0 + \check b^{\chi_j} \rho_{0} \check b^{\chi_j}\right) \check b^{s_j}, \end{align} \tag{ 111 }$$

where the coefficient $(1+(-1)^{s_n})/2$ ensures that we only take products containing an even number of $\rho_{A, o}$. Since $\check b_A^2\equiv I$, it is now useful to perform in equation (122) the change of variables $w_i = s_i(\!\!\!\!\!\!\mod 2)$,

$$\begin{align} \rho_A^n|_e = \frac{1}{2^n} \sum_{\substack{\left\{w_i = 0, 1\right\} \\ w_0 = w_n = 0}} \prod_{j = 1}^n \langle \check b_A \rangle^{|w_j-w_{j-1}|}\left( \check b_A^{w_j}\rho_0 \check b_A^{w_j} + \check b_A^{w_{j-1}} \rho_{0} \check b_A^{w_{j-1}}\right), \end{align} \tag{ 112 }$$

where the condition $w_0 = 0$ is due to the fact that the empty sum is 0 and $w_n = 0$ ensures that we take an even number of factors $\rho_{A, o}$. Equation (112) can be rewritten as

$$\begin{align} \rho_A^n|_e = \frac{1}{2^n} \sum_{\substack{\left\{w_i = 0, 1\right\} \\ w_0 = w_n = 0}} \prod_{j = 1}^n \langle \check b_A \rangle^{|w_j-w_{j-1}|}\left(\left(w_j+w_{j-1}\right)\rho_0 + \left(2-w_j-w_{j-1}\right) \check b_A\rho_0 \check b_A\right). \end{align} \tag{ 113 }$$

Once expanded, it finally takes the form

$$\begin{align} \rho_A^n|_e = \frac{1}{2^n}\sum_{\left\{\tau_i = 0, 1\right\}}\zeta\left[\left\{\tau_i\right\}\right]\prod_{j = 1}^n \check b_A^{\tau_j}\rho_0 \check b_A^{\tau_j}, \end{align} \tag{ 114 }$$

where

$$\begin{align} \zeta \left[\left\{\tau_i\right\}\right] = \sum_{\substack{\left\{w_i = 0, 1\right\} \\ w_0 = w_n = 0}}\prod_{j = 1}^n \left(\langle\check b_A\rangle^{|w_j -w_{j-1}|} \left(1 + \left(-1\right)^{\tau_j}\left(1-w_j-w_{j-1}\right)\right) \right). \end{align} \tag{ 115 }$$

According to equation (113), the calculation of the trace ${\mathrm{Tr}}(\rho_A^n|_e)$ boils down to computing the traces of products ${\mathrm{Tr}}(\prod_{j = 1} \check b_A^{\chi_j} \rho_0 b_a^{\chi_j})$. One can show that the matrix $\check b_A \rho_0 \check b_A$ is Gaussian,

$$\begin{align} \check b_A \rho\left[\Gamma_0^A\right] \check b_A = \rho\left[\Gamma_1^A\right], \end{align} \tag{ 116 }$$

and is univocally determined by the correlation matrix $(\Gamma_1^A)_{mm^{\prime}} = {\mathrm{Tr}}(\check b_A \rho_0 \check b_A a_m a_{m^{\prime}})$, $m, m^{\prime} \in A$, whose entries are given by the expression reported in equation (54). Therefore, we can apply the formula of equation (102) to compute the trace of a product of Gaussian operators in terms of their correlation matrices. We then find

$$\begin{align} {\mathrm{Tr}}\left(\rho_A^n|_e\right) = \frac{1}{2^n}\sum_{\left\{\tau_i = 0, 1\right\}}\zeta\left[\left\{\tau_i\right\}\right]\left\{\Gamma_{\tau_1}^A, \ldots, \Gamma_{\tau_n}^A\right\}.\end{align} \tag{ 117 }$$

The other ingredient left is the trace ${\mathrm{Tr}}(\rho_{A, e}^n)$, since $\rho_{A, e} = \rho_0$ is a Gaussian operator with correlation matrix $\Gamma_0^A$, then the straightforward application of equation (102) gives

$$\begin{align} {\mathrm{Tr}}\left(\rho_{A, e}^n\right) = \left\{\Gamma_0^A, \dots, \Gamma_0^A\right\}. \end{align} \tag{ 118 }$$

Plugging the results of equations (117) and (118) into equation (52), we obtain expression (53).

D.2. In PBC

In PBC, the reduced density matrix ρ_A takes the form of equation (71), and its even and odd parts correspond to

$$\begin{align} \rho_{A, e} = \rho_0, \quad \rho_{A, o} = \langle \sigma_1^x\rangle \sigma_1^x \rho_1. \end{align} \tag{ 119 }$$

As in OBC, the nth power of ρ_A takes the form of equation (105). To compute the entanglement asymmetry (52), we need to obtain the trace of the even part of $\rho_A^n$, which consists of the terms of the sum (105) that contain an even number of factors $\rho_{A, o}$.

Let us consider a particular product of operators $\rho_{A, e}$ and $\rho_{A, o}$, similar to the ones that appear in equation (105); for example,

$$\begin{align} \rho_{A, e} \rho_{A, o} \rho_{A, e} \rho_{A, e} \rho_{A, o} = \langle \sigma^x_1 \rangle^2 \rho_{0} \sigma^x_1 \rho_1 \rho_0 \rho_0 \sigma^x_1 \rho_1. \end{align} \tag{ 120 }$$

We can now use the fact that $\sigma^x_1 = \check a_1$, and that the product on both sides of a Gaussian operator by a Majorana fermionic excitation, e.g. $\check a_1 \rho_0 \check a_1$ and $\check a_1 \rho_1 \check a_1$, is also Gaussian. Since $\check a_1^2 = 1$, equation (120) can then be re-expressed as the product of Gaussian operators

$$\begin{align} \rho_{A, e} \rho_{A, o} \rho_{A, e} \rho_{A, e} \rho_{A, o} = \langle \sigma^x_1 \rangle^2 \rho_0 \left(\check a_1 \rho_1 \check a_1\right)\left( \check a_1 \rho_0 \check a_1\right)\left( \check a_1 \rho_0 \check a_1\right) \rho_1. \end{align} \tag{ 121 }$$

Observe that, whether or not the Majorana fermion $\check a_1$ multiplies on both sides, either ρ₀ or ρ₁ at position j depends on the number of $\rho_{A, o}$ factors in the product up to that position. If the number of $\rho_{A, o}$ is odd in the total string, there is an additional $\check a_1$ term at the end and the trace vanishes.

Using the previous observation, we can write the even part of $\rho_A^n$ as

$$\begin{align} \rho_A^n|_e = \sum_{\left\{\chi_i = 0, 1\right\}} \frac{1+\left(-1\right)^{s_n}}{2} \langle \sigma^x_1 \rangle^{s_n} \prod_{j = 1}^n \left(\check a_1^{s_j} \rho_{\chi_j} \check a_1^{s_j}\right), \end{align} \tag{ 122 }$$

where we defined $s_j = \sum_{i = 1}^j \chi_i$ and the term $(1+(-1)^{s_n})/2$ ensures that we only take products containing an even number of $\rho_{A, o}$. Thanks to equation (122), the trace of $\rho_A^n|_e$ can be obtained as a sum of traces of products of Gaussian operators ${\mathrm{Tr}}(\prod_{j = 1}^n\check a_1^{s_j}\rho_{\chi_j}\check a_1^{s_j})$, which can be computed in terms of their correlation matrices using equation (102). The only missing step is thus to determine the correlation matrices of $\check a_1 \rho_0 \check a_1$ and $\check a_1 \rho_1 \check a_1$. Employing the anticommutation relations of Majorana fermions, we can obtain the identity

$$\begin{align} \check a_1 \rho\left[\Gamma\right]\check a_1 = \rho\left[P_1 \Gamma P_1\right] \end{align} \tag{ 123 }$$

for any Gaussian matrix $\rho[\Gamma]$, where $P_1 = \textrm{diag}(-1, 1, \ldots, 1)$. Applying this equality, we conclude that the Gaussian operator $\check a_1 \rho_{\chi_j} \check a_1$ is determined by the correlation matrix $P_1 \Gamma_{\chi_j}^A P_1$. Therefore, using the identity (102) when we take the trace in equation (122), we have

$$\begin{align} {\mathrm{Tr}}\left(\rho_{A}^n|_e\right) = \sum_{\left\{\chi_i = 0, 1\right\}}\frac{1+\left(-1\right)^{s_n}}{2}\langle \sigma_1^x\rangle^{s_n}\left\{P_1^{s_1}\Gamma_{\chi_1}^AP_1^{s_1}, \ldots, P_{1}^{s_n}\Gamma_{\chi_i}^AP_{1}^{s_n}\right\}. \end{align} \tag{ 124 }$$

On the other hand, since $\rho_{A, e} = \rho_0$, we can directly apply equation (102) to calculate ${\mathrm{Tr}}(\rho_{A, e}^n)$ in terms of the correlation matrix $\Gamma_0^A$,

$$\begin{align} {\mathrm{Tr}}\left(\rho_{A, e}^n\right) = \left\{\Gamma_0^A, \overset{n}{\dots}, \Gamma_0^A\right\}. \end{align} \tag{ 125 }$$

Plugging equations (124) and (125) into equation (52), we finally find the formula (73) to study the entanglement asymmetry of spin parity in PBC.

In this appendix, we show how to write the reduced density matrix ρ_A for the ground state (64) in the form of equation (71). Our goal is to reconstruct the odd part of it, $\rho_{A, {o}}$, in terms of a Gaussian operator. Since we are in PBC, without loss of generality and to simplify the notation, we will take as subsystem A an interval from the sites j = 1 to $j = \ell$; that is, attached to the beginning of the Jordan–Wigner string (19). First, let us take the density matrix of the total system in the ground state (64), $\rho = \ket{\Psi}\bra{\Psi}$. It can be decomposed in even and odd parts

$$\begin{align} \rho = \rho_{e}+\rho_{o}, \end{align} \tag{ 126 }$$

where

$$\begin{align} \rho_{e} = \frac{1}{2}\ket{\emptyset}_\mathrm{NS~NS}{\bra{\emptyset}} + \frac{1}{2}\ket{\emptyset}_\mathrm{R~R}{\bra{\emptyset}} \end{align} \tag{ 127 }$$

and

$$\begin{align} \rho_{o} = \frac{1}{2}\ket{\emptyset}_\mathrm{R~NS}{\bra{\emptyset}}+ \frac{1}{2}\ket{\emptyset}_\mathrm{NS~R}{\bra{\emptyset}}. \end{align} \tag{ 128 }$$

In principle, $\rho_{A, o} = {\mathrm{Tr}}_B(\rho_o)$, but this is an odd operator, while Gaussian ones are even. Inspired by [61], the strategy is to construct a fake density matrix ρ₁ multiplying $\rho_{A, o}$ by an odd operator, such as $\sigma_1^x$,

$$\begin{align} \rho_1 = \frac{\sigma_1^x}{\langle \sigma_1^x\rangle } \rho_{A,o} ,\end{align} \tag{ 129 }$$

where the factor $1/\langle \sigma_1^x\rangle$ ensures that ${\mathrm{Tr}} \rho_1 = 1$. Observe that $\rho_1 = {\mathrm{Tr}}_B(\sigma_1^x \rho_{o})/\langle \sigma_1^x\rangle$ since for any operator $\mathcal{O}_A$ with support on A,

$$\begin{align} {\mathrm{Tr}}_A\left({\mathrm{Tr}}_B\left(\sigma_1^x \rho_{o}\right)\mathcal{O}_A\right) = {\mathrm{Tr}}\left(\sigma_1^x \rho_{o} \mathcal{O_A}\right) = {\mathrm{Tr}}_A\left(\sigma_1^x \rho_{A, o} \mathcal{O}_A\right) \end{align} \tag{ 130 }$$

and $\sigma_1^x$ acts on $\mathcal{H}_A$.

Although ρ₁ is an even operator, as such it cannot be computed in any practical way. This is where cluster decomposition enters into the game. We consider an auxiliary spin $\sigma_r^x$ at the site r and let us introduce

$$\begin{align} \rho_{1, r} = \frac{{\mathrm{Tr}}_B\left(\sigma_1^x \sigma_r^x \rho_{e}\right)}{\langle \sigma_1^x \sigma_r^x \rangle}. \end{align} \tag{ 131 }$$

If we consider an arbitrary operator $\mathcal{O}_A$ with support in A, then

$$\begin{align} {\mathrm{Tr}}_{A}\left(\rho_{1, r}\mathcal{O}_A\right) = \frac{{}_\mathrm{NS}{\bra{\emptyset}} \mathcal{O}_A \sigma_1^x \sigma_r^x \ket{\emptyset}_\mathrm{NS}}{\langle \sigma_1^x \sigma_r^x \rangle} = \begin{cases} 0~\text{if}~\mathcal{O}_A~\text{is odd,}\\[10pt] \dfrac{\langle \mathcal{O}_A \sigma_1^x \sigma_r^x \rangle}{\langle \sigma_1^x \sigma_r^x \rangle}~\text{if}~\mathcal{O}_A~\text{is even.} \end{cases} \end{align} \tag{ 132 }$$

If we now take the limit $r\to \pm \infty$ in the expression above, we can use the cluster decomposition principle so that

$$\begin{align} \lim_{r\to \pm \infty} {\mathrm{Tr}}_A\left(\rho_{1, r} \mathcal{O}_A\right) = \left.\begin{cases} 0~\text{if}~\mathcal{O}_A~\text{is odd}\\[10pt] \dfrac{\langle \mathcal{O}_A \sigma_1^x \rangle}{\langle \sigma_1^x \rangle}~\text{if}~\mathcal{O}_A~\text{is even} \end{cases}\right\} = {\mathrm{Tr}}\left(\rho_1 \mathcal{O}_A\right). \end{align} \tag{ 133 }$$

Therefore, $\rho_{1, r}$ is Gaussian since it is built from the even part of ρ, which is Gaussian itself and, consequently, it is univocally determined by a two-point correlation matrix. Moreover, it correctly reproduces the same correlations as ρ₁ when $r\to \pm \infty$. For practical reasons, we will now take $r\to -\infty$ so that the Jordan–Wigner string between the sites 1 and r does not belong to the subsystem A. We can now write

$$\begin{align} \rho_1 = \lim_{r\to - \infty} \rho_{1, r}. \end{align} \tag{ 134 }$$

The correlation matrix $\Gamma_{1, r}^A$ that determines $\rho_{1, r}$ has entries $\delta_{mm^{\prime}}+i(\Gamma_{1, r}^A)_{mm^{\prime}} = \langle \check a_m \check a_{m^{\prime}}\sigma_1^x \sigma_r^x \rangle$, with $m, m^{\prime}\in A$, which are expectation values of strings of Majorana fermions since

$$\begin{align} P_\mathrm{R}: = \sigma^x_1 \sigma^x_r = \sigma^x_r \sigma^x_1 = \prod_{j = r}^0 \left(-i\check a_{2j} \check a_{2j+1}\right). \end{align} \tag{ 135 }$$

Consequently, the correlation matrix of $\rho_{1,r}$ has the form

$$\begin{align} \left(\delta_{mm^{\prime}}+i\Gamma_{1, r}^A\right)_{mm^{\prime}} = \frac{\langle \check a_m \check a_{m^{\prime}} P_\mathrm{R} \rangle}{\langle P_\mathrm{R} \rangle}, \text{with }m, m^{\prime}\in A. \end{align} \tag{ 136 }$$

These correlations can be computed as a Pfaffian of a submatrix of Γ₀, defined in equation (65), as shown in [28]. However, a more computationally efficient way to obtain them is to use the results on reduced density matrices of disjoint intervals found in [61]. Indeed, correlation matrices of the form (136) can be computed as Schur complements. Since the string of equation (135) overlaps with subsystem A at m = 1, we need to consider the subsystem Aʹ as defined in figure 15. Then we introduce the correlation matrix $\Gamma_{1, r}^{A^{\prime}}$ with entries $\delta_{mm^{\prime}}+i(\Gamma_{1, r}^{A^{\prime}})_{mm^{\prime}} = \langle \check a_m \check a_{m^{\prime}} P_\mathrm{R}\rangle$, with $m, m^{\prime}\in A^{\prime}$.

Zoom In Zoom Out Reset image size

Let us now consider the restriction of the two-point correlation matrix Γ₀ to the subsystem $R\cup A^{\prime}$,

Zoom In Zoom Out Reset image size

The correlation matrix $\Gamma_{1, r}^{A^{\prime}}$ can then be obtained as a Schur complement [61]

$$\begin{align} \Gamma_{1, r}^{A^{\prime}} = \Gamma_0^{A^{\prime}} - \Gamma_0^{A^{\prime}, R}\left(\Gamma_0^{R}\right)^{-1}\Gamma_0^{R, A^{\prime}} \end{align} \tag{ 137 }$$

and it has the following structure:

Zoom In Zoom Out Reset image size

By definition, $(\Gamma_{1, r}^{A^{\prime}})_{mn} = (\Gamma_1^A)_{mn}$ for the sites $m, n\in A\cap A^{\prime}$. Using now the translational invariance of the chain and the fact that when $r\to -\infty$ we can take a separation of one less site with no consequence, we have

$$\begin{align} \langle \check a_1 \check a_j P_\mathrm{R} \rangle = \langle \check a_1 \check a_j \left(-i\check a_0 \check a_1\right) \prod_{j = r}^{-1} \left(-i \check a_{2j} \check a_{2j+1}\right) \rangle = i\langle \check a_0 \check a_j \prod_{j = r}^{-1} \left(-i \check a_{2j} \check a_{2j+1}\right) \rangle \underset{r\to \infty}{ = } i\langle \check a_2 \check a_{j+2} P_\mathrm{R} \rangle. \end{align} \tag{ 138 }$$

Therefore, we can finally reconstruct $\Gamma_1^A$ as

Zoom In Zoom Out Reset image size

Note that all this construction is based on the assumption that $\langle \sigma_1^x \rangle \neq 0$. It is a priori possible to derive $\rho_{A, o}$ from other odd correlators, but we found that it is always sufficient numerically to use this expression. When $\langle \sigma^x_1\rangle$ vanishes because of oscillations after the quench, we simply used the same reasoning applied to $\sigma^y_1$ and $\sigma^y_r$ to cure the divergence.