Loschmidt echo in XXZ Heisenberg spin chains

From the Quantum Transfer Matrix to the Quench Action: The Loschmidt echo in Heisenberg spin chains

Lorenzo Piroli, Balázs Pozsgay, Eric Vernier SISSA and INFN, via Bonomea 265, 34136 Trieste, Italy Department of Theoretical Physics, Budapest University of Technology and Economics, 1111 Budapest, Budafoki út 8, Hungary MTA-BME “Momentum” Statistical Field Theory Research Group, 1111 Budapest, Budafoki út 8, Hungary lpiroli@sissa.it, pozsgay.balazs@gmail.com, evernier@sissa.it
July 11, 2019
Abstract

We consider the computation of the Loschmidt echo after quantum quenches in the interacting Heisenberg spin chain both for real and imaginary times. We study two-site product initial states, focusing in particular on the Néel and tilted Néel states. We apply the Quantum Transfer Matrix (QTM) approach to derive generalized TBA equations, which follow from the fusion hierarchy of the appropriate QTM’s. Our formulas are valid for arbitrary imaginary time and for real times at least up to a time , after which the integral equations have to be modified. In some regimes, is seen to be either very large or infinite, allowing to explore in detail the post-quench dynamics of the system. As an important part of our work, we show that for the Néel state our imaginary time results can be recovered by means of the quench action approach, unveiling a direct connection with the quantum transfer matrix formalism. In particular, we show that in the zero-time limit, the study of our TBA equations allows for a simple alternative derivation of the recently obtained Bethe ansatz distribution functions for the Néel, tilted Néel and tilted ferromagnet states.

1 Introduction

One of the main motivations in the study of integrable systems is the possibility to provide reliable analytic results in the investigation of physically relevant but theoretically challenging problems. As a consequence, these models have traditionally played the role of the perfect benchmark to test the range of validity of approximate or perturbative methods developed to study more complicated systems.

This point of view has rapidly evolved in the past decade, due to the revolutionary progress in cold atomic physics [2, 3, 4, 5]. Indeed, a renewed motivation in the study of integrable models has come by their direct relevance for cold atomic experiments, where simple fine-tuned Hamiltonians can be realized with high control and physical quantities measured to high precision. Of particular importance is the possibility to experimentally realize different non-equilibrium settings [6, 7, 8, 9, 10, 11, 12, 13, 14]; among these the simplest is arguably that of a quantum quench, where a system is prepared in a well defined initial state and let evolved unitarily according to a known Hamiltonian [15].

These developments have unraveled the need for a solid theoretical understanding of key aspects of non-equilibrium dynamics in isolated quantum systems, and an active field of research has focused on the physical consequences of integrability. Of particular interest has been the question of whether and how a system locally approaches a stationary state at long times after the quench. Among the many results, excellent studies have now established the validity of a generalized Gibbs ensemble to predict the steady asymptotic value of local correlations [16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38], as well as the existence and physical importance of quasi-local conserved charges in integrable systems [39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 51, 52, 53, 50]. For a pedagogical introduction to these topics see the recent reviews [54, 55, 56, 57, 58].

It is hard to overemphasize the importance of analytical methods in these theoretical achievements. A significant role has been played by the introduction of the so called quench action approach, which has proven to be a powerful tool in the study of quantum quenches in Bethe ansatz integrable models [59, 60]. Indeed, the latter has already led to beautiful results in the characterization of post-quench steady states in several classes of experimentally relevant non-equilibrium settings [61, 29, 30, 31, 62, 63, 64, 65, 66, 67].

Despite the many successes, the computation of the whole time evolution following a quantum quench still represents a remarkable challenge. While by now several non-trivial explicit calculations have been performed in systems which can be mapped onto free ones [17, 68, 21, 69, 70, 71, 72, 73, 74, 75, 76, 77], only a few analytical or semi-analytical results are available for quenches in genuinely interacting models [80, 78, 62, 79, 81, 82, 83, 84, 85, 87, 88, 86]. The development of analytical techniques of investigations complementary or beyond the present methods thus represents an urgent issue.

Among the physically interesting quantities to probe the post-quench dynamics, the simplest is arguably the Loschmidt echo, namely the squared absolute value of the overlap between evolved and initial states. The latter has recently received a lot of attention, especially for its relevance in studies of dynamical phase transitions [89, 90, 91, 92, 93, 96, 94, 97, 95, 98, 99, 100, 101, 102, 103, 107, 105, 106, 104, 108, 109, 110, 112, 111]. Despite its simplicity, even for this quantity analytical results in interacting models are still extremely hard to obtain when a mapping onto free systems is not possible.

An analytical computation of the ”Loschmidt echo per site” in the interacting Heisenberg spin chain was presented in [95], where quantum quenches from the initial Néel state were studied by means of a quantum transfer matrix approach. Originally developed in the investigation of finite temperature problems [113, 114, 115], its extension to the computation of the Loschmidt amplitude relies on the interpretation of the latter as a particular boundary partition function. This natural identification has been exploited many times in the literature [90, 96, 97, 98, 100, 101, 102, 107, 106, 108, 109]. In the particular case of Heisenberg chains, this idea was also at the basis of the recent work [101], where an efficient numerical approach allowed to perform a remarkably detailed analysis of the Loschmidt echo in the thermodynamic limit for many different quantum quenches.

In this paper we build upon the work [95] and provide further results for a fully analytical computation of the Loschmidt echo in the XXZ Heisenberg chain. We first revisit the case of the initial Néel state considered in [95] and provide alternative formulas for the Loschmidt echo per site both at real and imaginary times (which is directly related to the so called dynamical free energy [95, 93]). These involve the solution of integral equations obtained by the so-called fusion relations of boundary transfer matrices. The formulas presented in this work allow for an explicit efficient numerical evaluation of the Loschmidt echo for real times.

Going further, we show that the same approach can be applied straightforwardly in the case of more general two-site product states such as tilted ferromagnets and tilted Néel states. We provide in particular a detailed study of the latter case, for which a corresponding set of integral equations is derived. We argue that these are valid for arbitrary imaginary times, while for real times only up to a time , which depends on the initial state and the final Hamiltonian parameters. In some regimes is seen to be either very large or infinite, allowing to probe in detail the relaxation dynamics of the system. For times larger than our formulas have to be modified in a way which is in principle feasible and discussed in the following sections.

Finally, we consider the computation of the Loschmidt echo per site at imaginary times by means of the quench action approach. In the case of the Néel state we show that the quantum transfer matrix results can be recovered, unveiling a direct link between the two approaches. Based on this, we provide an alternative derivation of the recently obtained Bethe ansatz rapidity distribution functions corresponding to the tilted ferromagnet and tilted Néel states [36]. Our treatment also naturally clarifies the validity of certain analytical properties of the latter (namely, the so called -system relations), which were previously established numerically [35, 36].

The rest of this article is organized as follows. In section 2 we introduce the Heisenberg chain and the quench protocol. In section 3 we review the quantum transfer matrix construction, while in section 4 some aspects of the boundary algebraic Bethe ansatz are presented. The analytical results for the Loschmidt echo at imaginary times are derived in section 5, and continuation to real times is discussed in section 6, where our analytical formulas are explicitly evaluated numerically. Section 7 is devoted to the comparison of the quantum transfer matrix formalism with the quench action approach. Finally, we report our conclusions in section 8, while some technical aspects of our work are provided in the appendices.

2 Setup

We consider the spin- Heisenberg model, defined on a chain of length . Since we will be interested in quenches where the initial state is a two-site product state, we will always assume to be an even number. The Hamiltonian of the model reads  
 

(1)

where we take , while are the local spin operators, being the Pauli matrices. We further assume periodic boundary conditions, , and restrict to the gapped antiferromagnetic phase . We then introduce the parametrization

(2)

with . The Hamiltonian (1) can be diagonalized by means of the Bethe ansatz method [116]. Relevant aspects of the latter will be discussed in the next section.

In this work we will consider quantum quenches from two-site product states of the form

(3)

where is an arbitrary state. Two relevant examples are the tilted Néel state defined by

(4)

and the tilted ferromagnet given by

(5)

Note that for the tilted Néel state coincides with the well-known Néel state

(6)

Given the initial state , the Loschmidt echo after the quench is defined as

(7)

and gives information about the probability of finding the system close to its initial state. For any finite , decays exponentially with the volume . It is then useful to define the Loschmidt echo per site [95]

(8)

which is simply related to the so called return rate [101]

(9)

As we mentioned in the introduction this quantity has recently received significant attention in the study of dynamical phase transitions, which are associated to points of nonanalyticity of the return rate (9). Analytical properties of the latter are more conveniently analyzed by considering a generic complex time and introducing the Loschmidt amplitude

(10)

By interpreting (10) as a boundary partition function, it was observed in [96] that for the transverse Ising chain these nonanalyticities occurred when the system was quenched accross a quantum critical point, establishing a connection between dynamical and equilibrium quantum phase transitions. Subsequent investigations showed that a more complicated picture takes place in general and nonanalyticities can be encountered also for quenches within the same quantum phase [93, 106, 101, 108].

It turns out that the analytical computation of the Loschmidt amplitude is facilitated when is taken to be a real number, in which case one obtains the so-called dynamical free energy density

(11)

which is also connected to the cumulant generating function for the Hamiltonian [95]. In this work we will consider the computation of (11) for which we provide a full analytical solution for quantum quenches from initial states of the form (3). The Loschmidt echo per site (and hence the return probablity) is then given by

(12)

We note that (12) has to be understood as the evaluation of the limit (11) for purely imaginary parameters. The limit does not necessarily commute with the analytic continuation, and this is responsible for possible non-analytic behavior of the function in the complex plane. As we will see, these issues require a delicate analysis which will be also addressed in our work.

3 Quantum Transfer Matrix approach to the Loschmidt echo

1

2

3

4

L-1

L

Figure 1: Pictorial representation of the quantity as the partition function of a six-vertex model on the cylinder, with boundary conditions in the imaginary time direction encoded by the initial state . There are horizontal rows, each line corresponding to the action of the transfer matrix , where for even/odd rows respectively.

In this section we review the idea behind the computation of the dynamical free energy (11) by means of the quantum transfer matrix approach. Following [95], our starting point is given by the well known Suzuki-Trotter decomposition

(13)

Next, this expression can be cast in a form more suitable for further analytical investigation. In particular, for large one can rewrite [115, 95]

(14)

where we defined

(15)

and where and are given in (1) and (2). In the following the dependence on will be omitted when this will not generate confusion and we will simply write instead of . In (14) we introduced the well-known transfer matrix , which is one of the central objects in the algebraic Bethe ansatz construction to diagonalize the Hamiltonian (1) [116]. It is defined as

(16)

where the trace is taken over the auxiliary quantum space . The Lax operators are written in terms of the -matrix corresponding to the six-vertex model

(17)

which in turn can be written as

(18)

The transfer matrices (16) for different values of the spectral parameter commute with one another. Furthermore, Eq. (14) is an immediate consequence of the following identity [115]

(19)

Equations (13), (14) allow to greatly simplify our problem. Indeed, the quantity may now be interpreted as the partition function of a six-vertex model on a square lattice. The latter has vertical and horizontal lines, with periodic boundary conditions in the horizontal (space) direction and boundary conditions specified by in the vertical (imaginary time) direction (cf. Fig. 1).

Under a reflection along the North-West diagonal (which leaves the weights of the six-vertex model invariant), and using the factorized structure (3) of the initial states considered in this work, the partition function can be reinterpreted as generated from a new transfer matrix. This is the so-called quantum transfer matrix, which acts in the original space direction as pictorially represented in Fig. 2, and which is generated from the following monodromy matrix

=

Figure 2: Pictorial representation of the partition function of Fig. 1 after reflection with respect to the North-West diagonal. Using the factorized form (3) of the initial state, the partition function may be viewed as generated by a quantum transfer matrix associated with the monodromy matrix . The inhomogeneties are for even/odd respectively and the boundary reflection matrices encode the dependence on [defined in (3)] and hence on the initial state.
(20)

In fact, one can easily derive

(21)

where is defined by the initial state through (3) and where the trace in the r. h. s. is along the physical spacial direction. Finally, defining

(22)

and putting everything together we arrive at

(23)

Following [95], we call the boundary quantum transfer matrix.

In analogy with the thermal case [115], we now consider the following two assumptions

  • For real values of the parameter , the boundary quantum transfer matrix has a leading eigenvalue whose absolute value remains separated from that of the subleading eigenvalues by a finite gap, even in the limit.

  • The large behaviour of (23) can be studied by exchanging the limits and .

If these assumptions are verified, it is straightforward to obtain in the large limit

(24)

and hence

(25)

This formula is the starting point for the analytical derivation of the dynamical free energy (11). Indeed, we are now left with the problem of computing the leading eigenvalue of the boundary quantum transfer matrix .

As we already mentioned, the form of explicitly depends on the initial state considered. In the case of the Néel state (6) the computation of was performed in [95], where was diagonalized by means of the so called diagonal boundary algebraic Bethe ansatz and the Trotter limit computed. As we will see in the next section, in the case of more general initial states of the form (3), one needs to resort to the non-diagonal version of the boundary algebraic Bethe ansatz and additional difficulties arise.

In the next section we first review the diagonal case corresponding to the Néel state and later discuss the more general non-diagonal boundary algebraic Bethe ansatz. These results will then be used in section 5 where an approach for the computation of the Trotter limit different to the one of [95] is proposed. The latter is based on the derivation of non-linear integral equations from the so called fusion of boundary transfer matrices. As we will see, one of the advantages of this method is that it can be straightforwardly applied both in the diagonal and non-diagonal cases.

4 The boundary algebraic Bethe ansatz: diagonal and non-diagonal boundaries

The boundary algebraic Bethe ansatz is an analytical method which allows to diagonalize Hamiltonians of open spin chains with integrable boundary conditions [117]. Here we only review the aspects relevant to our work, while we refer to the specialized literature for a more systematic treatment [118, 119, 120].

Given a chain of length , the construction starts by introducing a boundary transfer matrix defined as

(26)

Here we introduced

(27)
(28)
(29)

where the Pauli matrix acts on the auxiliary space and where indicates transposition in . The last equality follows from the properties of the -matrix defined in (18) [118], while the inhomogeneities are parameters which for the moment are left arbitrary. Finally, the trace in (26) is performed over the auxiliary space [not to be confused with the auxiliary space of the physical transfer matrix (16)].

Figure 3: Symbolic representation of the transfer matrix in Eq. (26), acting on sites with inhomogeneous spectral parameters .

The boundary reflection matrices are matrices

(30)

which are solution of the so called reflection equations [117]. The boundary transfer matrix (26) is symbolically represented in Fig. 3.

The relevance of this construction for our purposes lies in the possibility of interpreting the boundary quantum transfer matrix in (22) as an operator of the form (26). This in turn allows us to employ boundary algebraic Bethe ansatz techniques for the computation of the leading eigenvalue . We explicitly show this in the following.

First, introducing the components of in the auxiliary space as

(31)

it follows from (29) that

(32)

where the components , , , are operators acting on the physical space . Using (31), (32), it is now straightforward to rewrite in (26) as

(33)

where we introduced the vectors defined as

(34)
(35)

It is now evident from (33) that is proportional to in (22) provided that

(36)

and that the inhomogeneous spectral parameters are chosen as

(37)
(38)

If these conditions are met, one simply obtains

(39)

This relation allows us to directly obtain the eigenvalues of once the eigenvalues of are known. Note once again that depends explicitly on the initial state through (36). In particular, the identification (36) fixes the -matrix (30) through (34), (35). Different initial states then require diagonal or non-diagonal -matrices. In turn, this makes it necessary to resort to either diagonal or non-diagonal boundary algebraic Bethe ansatz techniques to obtain the eigenvalues of . We now separate the discussion for these two different cases.

4.1 Diagonal reflection matrices: the Néel state

In the simplest case, the identification (36) leads to diagonal -matrices. This is what happens for the Néel state (6), which was explicitly considered in [95] (together with the so called Majumdar-Ghosh sate).

A diagonal solution of the reflection equation can be obtained as [118]

(40)
(41)

Then, from (34), (35), one can easily see that condition (36) is satisfied by choosing the boundary parameters as

(42)

With this choice, (34) and (35) yield

(43)
(44)

The choice (42) completely specifies the diagonal -matrix (41) and hence the boundary transfer matrix (26), which can then be diagonalized.

The eigenvalues of the transfer matrix , and therefore the leading eigenvalue of , can be constructed through the diagonal boundary algebraic Bethe ansatz procedure, which we briefly review here.

Introducing the notation

(45)

the common eigenstates of the operators are obtained from the ferromagnetic reference eigenstate as

(46)

Here, the complex parameters , the so-called rapidities, have to be chosen to satisfy the Bethe equations

(47)

As we will comment later, the leading eigenvalue corresponds to a set of rapidities. In the following we will then restrict to this case.

Given the set , and following [95], it is convenient to introduce the doubled set

(48)

Defining further

(49)
(50)
(51)
(52)

the Bethe equations can be rewritten as

(53)

A given solution of the Bethe equations (53) corresponds to an eigenvalue of , which we indicate as . It reads

(54)

where the dependence on is encoded in the functions . Eq. (54) is sometimes referred to as the relation.

The formal relation (54) may be compared with exact diagonalization for finite system sizes , allowing to find the set of roots associated with the eigenstates of interest. By numerical diagonalization of the transfer matrix for , we find that the (unique) leading eigenvalue of is obtained from a set of roots which are situated on the imaginary axis. Therefore the rapidities of the corresponding doubled set are distributed symmetrically on the imaginary axis. An example is shown in Fig. 4.

Before discussing the Trotter limit of the leading eigenvalue, we present in the next section the case of non-diagonal reflection matrices. The Trotter limit will then be discussed in section 5, where we employ an approach which can be straightforwardly applied both in the diagonal and non-diagonal cases.

Figure 4: Bethe roots , as obtained by solving (53), corresponding to the leading eigenvalue of the transfer matrix in the diagonal case. The plot corresponds to , , . We see that the Bethe roots are located symmetrically along the imaginary axis.

4.2 Non-diagonal reflection matrices: tilted Néel and tilted ferromagnet states

For more general initial states, the identification (36) leads through (34), (35) to non-diagonal boundary transfer matrices. Let us introduce the general non-diagonal solution of Sklyanin’s reflection equation [117, 119, 120]

(55)
(56)

For later convenience we introduce , , defined by the following parametrization

(57)

Analogously to the diagonal case discussed in the previous section, the parameters , , [and hence , , through (57)] can be chosen in such a way that the condition (36) is satisfied for a given initial state. In the following we report the choice of the parameters for tilted Néel and tilted ferromagnet states.

  • Tilted Néel state. In this case, it is straightforward to verify that the parameters of the -matrix can be chosen as

    (58)
    (59)
    (60)
    (61)

    where is the tilting angle in the definition (4). With this choice condition (36) is met. Explicitly,

    (62)

    Note in particular that if then , and using one consistently recovers the result of the previous section for the Néel state.

  • Tilted ferromagnet. The parameters of the -matrix are now chosen as

    (63)
    (64)
    (65)
    (66)

    where again is the tilting angle in the definition (5). Using this choice, it is straightforward to verify

    (67)

After fixing the parameters of the -matrix (56), the spectral problem associated with the boundary transfer matrix (26) can be addressed. The procedure described in section 4.1 for the diagonal case cannot be applied directly here. This is because the ferromagnetic state is not anymore an eigenstate of the transfer matrix , and cannot therefore be used as the reference state for the construction of all eigenstates.

In fact, the long-standing problem of completely characterizing the spectrum of the boundary transfer matrix for arbitrary -matrices has only recently been solved, as a result of the combined effort of several groups [119, 120, 121, 122, 123, 124, 125, 126] (see [119, 120] for some historical details on these interesting developments). In particular, an important discovery was that a generalized version of the relation (54) could be recovered also in the non-diagonal case, involving again a finite set of rapidities. In turn, these can be obtained as the solution of an appropriate set of Bethe equations. Here we simply report the results which are directly relevant to our work, while we refer to [119, 120] for a thorough treatment.

In the following we employ many of the notations used in [119]. We start by introducing the so called inhomogeneous relation verified by the eigenvalues of the boundary transfer matrix, which is written as

(68)

We will now define the functions appearing above. First, the parameter are defined by the -matrix (55), while

(69)

where was introduced in (50) and where

(70)

Here the parameters , , , and are defined in (56) and (57). Further, the function is given by

(71)

where

(72)

Finally, the -functions are parametrized by a set of rapidities as

(73)

which can be rewritten, introducing analogously to the diagonal case the doubled set , as

(74)

In this case the (doubled) set of rapidities is determined as the solution of a new set of Bethe equations containing an additional term, namely

(75)

We stress that the inhomogeneous relation (68) has to be understood as follows: for every eigenstate of the boundary transfer matrix (26), there exist a set of solutions of the inhomogeneous Bethe equations (75) for which the corresponding eigenvalue can be written as (68). Note that in the diagonal case the number of Bethe roots is fixed by the value of the total magnetization (which commutes with the transfer matrix and is therefore well defined for each eigenstate). Here, instead the latter is not conserved and the number of Bethe roots is always exactly equal to the number of sites of the chain, namely .

As in the diagonal case, we can compare the formal relation (68) with exact diagonalization for small system sizes of length and we observe that there is once again a unique leading eigenvalue. In the next sections, we mainly focus on the case of tilted Néel states, for which we provide explicit results for the dynamical free energy. In this case we studied in detail the corresponding doubled set of Bethe roots , which display a recognizable structure. In particular, we observed that they organize into two disjoint sets as

(76)

where

  • the roots are displaced symmetrically along the imaginary axis ,

  • the roots are distributed as pairs of common imaginary part and opposite real parts.

We report an example of this structure in the complex plane in figure 5.

Figure 5: Bethe roots , as obtained by solving (75), associated with the leading eigenvalue of the transfer matrix in the non-diagonal case. The plot corresponds to , , , and boundary parameters associated with the tilted Néel state at . We recognize the structure discussed in section 4.2 with roots situated along the imaginary axis () and additional roots with non-zero real part, located symmetrically with respect to the imaginary axis ().

For the tilted Néel state one can study the diagonal limit . In this limit, the real parts of the roots diverge proportionally to [defined in (61)]. Their contributions to ratios of functions cancel out in the relation (68), while the inhomogeneous term vanishes. In this limit, one can then factor the contributions of these extra roots from the function, and the roots satisfy the homogeneous relation (54). In fact we observed that even for finite nonzero the values of the roots are close to the solutions of the corresponding homogeneous system at .

The results of this and the previous sections give access to the leading eigenvalue of the transfer matrix (26) for small Trotter number . In the next section we address the problem of computing the Trotter limit which directly yields the dynamical free energy (11).

5 Integral equations from fusion of boundary transfer matrices

We now finally address the Trotter limit (25) required in the computation of the dynamical free energy (11). In [95] this problem was solved for the diagonal case by introducing a single non-linear integral equation for an auxiliary function in the complex plane which could be directly related to the leading eigenvalue of the transfer matrix (22). The method employed in [95] heavily relied on the particular structure of the Bethe roots in the complex plane. As we already stressed, the latter is very simple in the diagonal case, as one can immediately see in Fig. 4. By contrast, in the non-diagonal case, the picture is significantly more involved due to the presence of the additional Bethe roots as discussed in section 4.2 (cf. also Fig. 5). As a consequence, the approach used in [95] can not be directly applied and a more sophisticated analysis is required to adapt it to the non-diagonal case111 In the closely related problem of the physical spin chain with transverse boundary magnetic fields similar complications occur, and the problem of finding the Bethe roots and considering the thermodynamic limit of their distribution is not yet settled [127]. .

In this work we consider a different approach based on the derivation of non-linear integral equations from fusion of boundary transfer matrices. This procedure has been employed many times in the thermal case, where it is now well established [128, 129]. On the one hand, this method can be applied directly both in the diagonal and non-diagonal cases considered in this work. On the other hand, even in the diagonal case discussed in [95] it is seen to be more stable when continuation to real times is addressed as discussed in section 6. The main qualitative difference is that the final result is written in terms of an infinite set of non-linear integral equations, as opposed to the single non-linear integral equation derived in [95].

As a preliminary step, we introduce in the next subsection the family of fused boundary transfer matrices which can be build out of (26). As we will explain in the following, the main idea is to relate the leading eigenvalue of in (22) to the solution of the -system of fused boundary transfer matrices. In turn, this can be obtained from the corresponding -system, which can be conveniently cast in the form of partially decoupled integral equations. The rest of this section is devoted to following this program, and each step will be explained in full detail.

5.1 The system and system for boundary quantum transfer matrices

It is an established result that the boundary transfer matrix in (26) can be used to build an infinite family of transfer matrices by the so called fusion procedure [130, 131, 132], in complete analogy with the case of periodic boundary conditions.

The fused transfer matrices act on the same space as the transfer matrix and form a commuting set, namely

(77)

with the further identification

(78)

The family of fused boundary transfer matrices satisfy a set of functional relations known as the -system [132]

(79)

where

(80)

The function encodes the information about the boundary reflection -matrices (30). In the general non-diagonal case, it reads [132]

(81)

where the function is defined in (50), while

(82)
(83)

Here the parameters , are defined by the -matrix (56). Note that the diagonal case (41) is simply recovered by setting .

Note that since the transfer matrices