Merits of using density matrices instead of wave functions in the stationary Schr\"odinger equation for systems with symmetries

TL;DR

This paper advocates using density matrices in the stationary Schrödinger equation for systems with symmetries, highlighting advantages in handling degeneracies and symmetry properties, with applications to spin systems.

Contribution

It introduces the formulation of the stationary Schrödinger equation using density matrices, demonstrating benefits in symmetry considerations and degeneracy handling in many-body quantum systems.

Findings

Density matrices can always be chosen to respect all Hamiltonian symmetries.

The approach simplifies understanding of eigenstate structures.

Analytical solutions for spin systems are provided.

Abstract

The stationary Schr\"odinger equation can be cast in the form $H \rho = E \rho$, where $H$ is the system's Hamiltonian and $\rho$ is the system's density matrix. We explore the merits of this form of the stationary Schr\"odinger equation, which we refer to as~SSE$_\rho$, applied to many-body systems with symmetries. For a nondegenerate energy level, the solution $\rho$ of the SSE$_\rho$ is merely a projection on the corresponding eigenvector. However, in the case of degeneracy $\rho$ is non-unique and not necessarily pure. In fact, it can be an arbitrary mixture of the degenerate pure eigenstates. Importantly, $\rho$ can always be chosen to respect all symmetries of the Hamiltonian, even if each pure eigenstate in the corresponding degenerate multiplet spontaneously breaks the symmetries. This and other features of the solutions of the SSE$_\rho$ can prove helpful by easing the notations and providing an unobscured insight into the structure of the eigenstates. We work out the SSE$_\rho$ for a general system of spins $1/2$ with Heisenberg interactions, and address simple systems of spins $1$. Eigenvalue problem for quantum observables other than Hamiltonian can also be formulated in terms of density matrices. As an illustration, we provide an analytical solution to the eigenproblem ${\bf S}^2 \rho=S(S+1) \rho$, where $\bf S$ is the total spin of $N$ spins $1/2$, and $\rho$ is chosen to be invariant under permutations of spins. This way we find an explicit form of projections to the invariant subspaces of ${\bf S}^2$.

Figures

Equations

$$H|\Psi ⟩=E|\Psi ⟩,$$

$$H\rho =E\rho .$$

$${\rho }^{†}=\rho ,\text{tr}\rho =1,\rho >0.$$

$$\rho =\sum _i{p}_i|{\Psi }_i⟩⟨{\Psi }_i|,p_i\ge 0,\sum _i{p}_i=1,$$

$$UH{U}^{†}=H\forall U\in G,$$

$${\rho }_G=\underset{G}{\int }U\rho U^{†}d\mu (U)$$

$$[H,\rho ]=0.$$

$$H=\sum _{i<j}{J}_{ij}\left({\mathbf{\sigma }}_i{\mathbf{\sigma }}_j\right),i,j=1,2,...,N,$$

$$\mathcal{A}=\left(i_1,j_1\right)...\left(i_m,j_m\right),1\le m\le [N/2],$$

$$i_p<i_{p+1},j_p>i_p,j_p\ne i_l,j_l\forall l\ne p,1\le i_p,j_p\le N.$$

$$A_{\mathcal{A}}=\left({\mathbf{\sigma }}_{i_1}{\mathbf{\sigma }}_{j_1}\right)\left({\mathbf{\sigma }}_{i_2}{\mathbf{\sigma }}_{j_2}\right)...\left({\mathbf{\sigma }}_{i_m}{\mathbf{\sigma }}_{j_m}\right),$$

$$\rho =\frac{1}{2^N}\left(a_0{A}_0+\sum _{\mathcal{A}}{a}_{\mathcal{A}}{A}_{\mathcal{A}}\right).$$

$$({\mathbf{\sigma }}_1{\mathbf{\sigma }}_2{)}^2=$$

$$({\mathbf{\sigma }}_1{\mathbf{\sigma }}_2)({\mathbf{\sigma }}_2{\mathbf{\sigma }}_3)=$$

$$E{a}_0$$

$$E{a}_{\mathcal{A}}$$

$$+(J_{i_l{i}_m}+J_{j_l{j}_m}-J_{i_l{j}_l}-J_{i_m{j}_m})a_{(i_l,j_m)(j_l,i_m){\mathcal{A}}^{\overline{l,m}}})$$

$$+\sum _{\begin{array}{c}p,q:p<q\\ p,q\notin \mathcal{A}\end{array}}\left(3J_{pq}{a}_{(p,q)\mathcal{A}}+\sum _{m=1}^{|\mathcal{A}|}{J}_{pq}(a_{(i_m,p)(q,j_m){\mathcal{A}}^{\overline{m}}}+a_{(q,i_m)(j_m,p){\mathcal{A}}^{\overline{m}}})\right)$$

$$+\sum _{p\notin \mathcal{A}}\sum _{m=1}^{|\mathcal{A}|}\left(J_{p{j}_m}{a}_{(p,i_m){\mathcal{A}}^{\overline{m}}}+J_{p{i}_m}{a}_{(p,j_m){\mathcal{A}}^{\overline{m}}}\right)+\sum _{m=1}^{|\mathcal{A}|}{J}_{i_m{j}_m}\left(a_{{\mathcal{A}}^{\overline{m}}}-2a_{\mathcal{A}}\right),$$

$$0$$

$$-\left((J_{ik}-J_{jk})a_{\left(ij\right)\mathcal{A}}+(J_{jk}-J_{ij})a_{\left(ik\right)\mathcal{A}}+(J_{ij}-J_{ik})a_{\left(jk\right)\mathcal{A}}\right),\forall i,j,k\notin \mathcal{A}\mathrm{and}i<j<k.$$

$$H$$

$$8\rho$$

$$\left(\begin{array}{cccc}0& 3J_{12}& 3J_{23}& 3J_{13}\\ J_{12}& -2J_{12}& J_{23}& J_{13}\\ J_{23}& J_{12}& -2J_{23}& J_{13}\\ J_{13}& J_{12}& J_{23}& -2J_{13}\end{array}\right)\left(\begin{array}{c}{a}_0\\ a_{12}\\ a_{23}\\ a_{13}\end{array}\right)=E\left(\begin{array}{c}{a}_0\\ a_{12}\\ a_{23}\\ a_{13}\end{array}\right),$$

$$0$$

$$E$$

$${\mathbf{S}}^2\rho =\lambda \rho ,$$

$$\rho =\frac{1}{2^N}\left(1+\sum _{m=1}^{[N/2]}{a}_m{A}_m\right),$$

$$A_1{A}_m$$

$$\displaystyle+2m\left({(N-2m)\theta(N-2m)-1}\right){A_{m}}\begin{array}[]{c}+(m+1)\theta(N-2m-1){A_{m+1}},\end{array}m=1,...[N/2].$$

$$a_1$$

$$a_2$$

$$a_m$$

$$H=({\mathbf{S}}_1{\mathbf{S}}_2),$$

$$\rho =a_0+a_1\left({\mathbf{S}}_1{\mathbf{S}}_2\right)+a_2{\left({\mathbf{S}}_1{\mathbf{S}}_2\right)}^2.$$

Full text

1 General properties of the Stationary Schrödinger equation for density matrices

The conventional form of the stationary Schrödinger equation (which we refer to as SSEΨ) reads

[TABLE]

where $H$ is the Hamiltonian of a quantum system, $E$ is its eigenenergy and $\left|{\Psi}\right\rangle$ is the corresponding eigenvector. Obviously, this equation implies an operator identity $H\left|{\Psi}\right\rangle\left\langle{\Psi}\right|=E\left|{\Psi}\right\rangle\left\langle{\Psi}\right|$, where $\left|{\Psi}\right\rangle\left\langle{\Psi}\right|$ is the projection onto the eigenvector $\left|{\Psi}\right\rangle$. One can extend this identity by considering an arbitrary density matrix $\rho$ instead of a projection operator. This way one obtains the stationary Schrödinger equation for density matrices:

[TABLE]

In the present paper we explore the properties and the merits of this equation which we refer to as SSEρ. Our studies are somewhat close in spirit to the research avenue on the contracted Schrödinger equation, see e.g. [1] and references therein. Some important differences will be discussed below when we apply the SSEρ to a particular spin system in Section 2.

We remind that a density matrix should satisfies three conditions,

[TABLE]

The following relations between the SSEΨ and SSEρ follow immediately.

Eqs. (1) and (2) share the same set of eigenvalues $E$. 2. 2.

If a given eigenvalue $E$ is nondegenerate, then the corresponding $\left|{\Psi}\right\rangle$ and $\rho$ are related according to $\rho=\left|{\Psi}\right\rangle\left\langle{\Psi}\right|$. 3. 3.

In the case of degeneracy any solution of the SSEρ reads

[TABLE]

where vectors $\left|{\Psi_{i}}\right\rangle$ constitute a basis in the corresponding degenerate subspace of the Hamiltonian.

Properties (ii) implies that in the nondegenerate case SSEΨ and SSEρ are, in fact, identical up to notations (however, even in this case SSEρ can be more convenient compared to SSEΨ, in particular for spin systems, see e.g. [3]). An important advantage of the SSEρ shows up in the case of a degeneracy induced by some symmetry of the Hamiltonian. Assume that the Hamiltonian is symmetric under some symmetry group $G$, i.e.

[TABLE]

where $U$ is a unitary operator. In this case the eigenbasis of $H$ is split into blocks which determine degenerate subspaces invariant under the group $G$. As a rule, a spontaneous symmetry breaking phenomenon occurs in some of this subspaces, which means that any eigen basis in such a subspace contains eigenvectors not invariant with respect to $G$. The most trivial example of a system with the spontaneous symmetry breaking is a single spin $1/2$ with a vanishing Hamiltonian, $H=0$. Such Hamiltonian is invariant under the group of rotations, however any its eigenstate (i.e. any pure state of a single spin $1/2$ ) lacks this symmetry. In general, an important class of spin systems with Heisenberg interactions demonstrate spontaneous symmetry breaking (either in the ground state or in excited states). Some examples of such systems will be considered in what follows. The phenomenon of the spontaneous symmetry breaking, while being of paramount importance for physics [4], can sometimes cause various inconveniences. In particular, it obscures calculations of the correlation functions invariant with respect to $G$. In contrast to SSEΨ, one can always avoid the spontaneous symmetry breaking in the solutions of SSEρ, according to the following simple

Lemma. Consider a Hamiltonian $H$ invariant under the symmetry group $G$ ($G$-invariant, for short). For any eigenvalue $E$ of this Hamiltonian there exists a $G$-invariant density matrix $\rho_{G}$ which is a solution of the Schrödinger equation (2).

Proof. Consider a (not necessarily $G$-invariant) density matrix $\rho$ which is a solution of eq. (2) corresponding to a given $E$. A $G$-invariant solution $\rho_{G}$ can be obtained from $\rho$ by averaging over the group $G$ with the Haar measure $d\mu(U)$ [5, 6]:

[TABLE]

where the normalization condition $\int_{G}1\,d\mu(U)=1$ is implied. It is easy to see that thus obtained $\rho_{G}$ is indeed a legitimate density matrix (i.e. satisfies conditions (3)) and is invariant under the group $G$ ( i.e. $U\,\rho_{G}\,U^{\dagger}=\rho_{G}~{}~{}~{}~{}~{}\forall~{}~{}U\in G$). ∎

Obtaining $G$-invariant objects by averaging over the group $G$ with the Haar measure is a standard tool of the group theory [5], and mixed states of the form (6) naturally appear in various resource theories [6, 7, 8, 9]. It should be emphasized, however, that we use such averaging only as a formal tool to prove the existence result contained in the above Lemma. In practice, we propose to ensure the $G$-invariance by expanding the density matrix in $G$-invariant basis operators, without explicitly performing the averaging (6). We apply this approach to specific examples in what follows.

It is worth highlighting why the averaging over the group analogous to that in eq. (6) can not be applied directly to vectors in the Hilbert space. This is because such averaging does not conserve the normalization, and one can obtain a zero vector (which lacks physical interpretation) as a result. This indeed happens, as can be seen in the trivial example of a single spin $1/2$ with $H=0$ discussed above. In contrast, averaging (6) of density matrices conserves the trace and thus the normalization.

The Schrödinger equation (2) entails

[TABLE]

In fact, this equation holds not only for the solutions of the SSEρ but for any stationary state, i.e. a state not evolving under the Liouville–von Neumann equation. Eq. (7) is widely used to obtain constraints on expectation values of various observables in equilibrium [11, 12, 13]. For our purposes it is essential that this equation does not contain $E$ and can be used to reduce the dimensionality of the more computationally demanding eigenvalue problem. This is discussed in more detail in the next section.

In the rest of the paper we illustrate the concept of the SSEρ by considering specific spin systems.

2 System of spins $1/2$ with the Heisenberg interaction

In the present section we specialise the SSEρ for a system of $N$ spins with the Heisenberg interaction. The Hamiltonian of this system reads

[TABLE]

where ${\bm{\sigma}}_{i}$ is the vector consisting of three Pauli matrices of the $i$’th spin, $J_{ij}$ is the coupling constant between $i$’th and $j$’th spins and $\left({{{\bm{\sigma}}_{i}}{{\bm{\sigma}}_{j}}}\right)$ is the corresponding scalar product of sigma-matrices. This Hamiltonian is invariant with respect to a global $SU(2)$ symmetry, in other words, to the simultaneous rotations of all spins. In addition, it is $T$-invariant, i.e. invariant with respect to the inversion of time. This Hamiltonian, apart from being a popular theoretical playground, is of practical importance in material science, both for finite [15] and for infinite [16] $N$.

Due to the presence of the above symmetries, we can look for a solution $\rho$ of the $SSE_{\rho}$ which is constructed of scalar products of sigma matrices. To this end, we define a multi-index ${\cal A}$ enumerating the set of pairs $(i_{p},j_{p})$:

[TABLE]

where $[N/2]$ is the integer part of $N/2$, while $i_{p}$ and $j_{p}$ enumerate spins and for any $p$ satisfy

[TABLE]

These conditions ensure that the sum over ${\cal A}$ runs over all distinct sets of pairs of indices in which each index is found at most once. We denote the number of pairs in ${\cal A}$ by $|{\cal A}|$ (e.g. $|{\cal A}|=m$ in eq. (9)). Finally, we define an operator $A_{\cal A}$ as a product of $|{\cal A}|$ scalar products of Pauli matrices according to

[TABLE]

where ${\cal A}$ is given by eq. (9). We supplement this definition by a convention $A_{0}\equiv 1$.

Our ansatz for $\rho$ can now be written as

[TABLE]

Here $a_{\cal A}$ are numerical coefficients and the sum is over all sets ${\cal A}$ of the form (9),(10). Note that normalization implies $a_{0}=1$. Obviously, such $\rho$ is both $SU(2)$-invariant and $T$-invariant. In fact, any $SU(2)$- and $T$-invariant operator with a unit trace can be represented in this form. Let us briefly explain why. First, observe that due to the equalities

[TABLE]

one can avoid terms with repeating spin indexes, as is indeed the case in eq. (12). Here $({{{\bm{\sigma}}}_{1}}{{{\bm{\sigma}}}_{2}}{{{\bm{\sigma}}}_{3}})$ is the mixed product of vectors consisting of Pauli matrices of three spins. Further, observe that the mixed product changes its sign under time inversion and thus does not enter $\rho$. As for the even powers of mixed products, they can always be expressed through the scalar products [10].

One can now substitute the ansatz (12) into the stationary Schrödinger equation (2). Exploiting formulae (13), after straightforward but tedious calculations one obtains the following equations for the coefficients $a_{\cal A}$:

[TABLE]

[TABLE]

The multi-index $(i_{l},j_{m})(j_{l},i_{m}){\cal A}^{\overline{l,m}}$ is obtained from ${\cal A}$ by dropping pairs $(i_{l},j_{l})$, $(i_{m},j_{m})$ and adding pairs $(i_{l},j_{m})$, $(j_{l},i_{m})$. Other multi-indices used in the above equations are obtained from ${\cal A}$ analogously.

Note that eq. (16) is due to eq. (7). While eqs. (14), (15) represent an eigenvalue problem, eq. (16) is a homogeneous linear equation. Solving the latter is computationally less demanding than solving the former. Thus it can be numerically efficient to first eliminate as many variables as possible with the help of eq. (16), and then solve an eigenvalue problem with a smaller number of variables.

We remark that in the paradigm of the Contracted Schrödinger equation [1] the spin Hamiltonian should be first turned into a fermionic Hamiltonian as in ref. [2]. While this is easily done for one-dimensional spin chains by means of the Jordan-Wigner transformation, in higher dimensions this leads to nonlocal interactions. In contrast, in our procedure transformation to fermionic representation is not required.

For illustrative purposes we apply the SSEρ to the system of three spins:

[TABLE]

In this case eqs. (14), (16) read

[TABLE]

[TABLE]

Note that the homogeneous linear equation (20) is, in principle, redundant, but in practice can be useful for simplifying the eigenproblem (19), as discussed above. We also remark that the size of eigenproblem is twice smaller than what one would obtain by a straightforward application of the conventional Schrödinger equation (1) to the system of three spins $1/2$. This size is even more reduced if the Hamiltonian posses additional symmetries, see below. This can prove useful for exact diagonalization studies of small spin clusters, which can be of interest for understanding magnetic response of correlated materials [19] (for alternative ways of accounting for symmetries see [15]).

If ${J_{12}}={J_{23}}={J_{13}}=1$, eq. (19) leads to two sets of solutions:

[TABLE]

see figure 1 for illustration.

3 Total spin of $N$ spins $1/2$: Projections on invariant subspaces

Total spin operator is formally equivalent to the Heisenberg Hamiltonian with long-range interactions. We seek to solve the eigenproblem for total spin of N qubits

[TABLE]

with an additional constraint that $\rho$ is invariant under permutations of spins,

[TABLE]

where $A_{m}$ is the sum of all possible products (11) of $m$ scalar products of sigma matrices. It can be easily found that ${\bm{S}^{2}}=\frac{1}{4}(3N+2{A_{1}})$ and

[TABLE]

Here $\theta(x)$ is 1 if $x>0$ and 0 otherwise. The above relation defines a tridiagonal matrix which has eigenvalues $\lambda=S(S+1)$, as we verified numerically. They lead to the recursive formula for the coefficients $a_{m}$:

[TABLE]

It can be verified that thus obtained density matrices (24) coincide up to normalization with the projections to invariant subspaces of ${\bm{S}^{2}}$.

4 Systems of spins $1$

In this section we briefly outline how SSEρ can be applied to systems of spins $1$. While for a spin $1/2$ three operators of spin projections along with the identity operator span the whole space of self-adjoint operators, this is not the case for a spin $1$. As a result, the ansatz for an invariant density matrix becomes more complicated. Let us start from a system of two spins $1$ with the Hamiltonian invariant under rotations,

[TABLE]

where ${\bm{S}}_{i}$ is the spin at $i$’th site, $({\bm{S}}_{i}{\bm{S}}_{i})=2$. A general form of the rotationally-invariant density matrix reads

[TABLE]

Higher powers of the scalar product $({\bm{S}}_{1}{\bm{S}}_{2})$ are linearly dependent on the first two powers according to

[TABLE]

The $SSE_{\rho}$ then reads

[TABLE]

with the solution

[TABLE]

where the normalization condition $\tr\rho=1$ is taken into account.

Now we turn to a case of three spins with interactions invariant under rotations, time reversal and permutations. A density matrix respecting these symmetries can be parameterized as

[TABLE]

where $\sum\limits_{i\neq j\neq k}$ is the sum over six triples of distinct $i,j$ and $k$. All other polynomials constructed of scalar products are linear dependent on those presented in eq. (37). The normalization condition implies

[TABLE]

Thus the ansatz contains only four real parameters (say, $a_{1},a_{21},a_{22},a_{3}$), to be compared to $2\times 3^{3}-1=53$ real parameters required to parameterize a pure state in the Hilbert space of three spins $1$ without account for symmetries.

We consider two different three-spin Hamiltonians respecting the above symmetries. The first one is

[TABLE]

The SSEρ for this Hamiltonian reads

[TABLE]

The eigenvalues are $(-3,-2,0,3)$. The coefficient $(a_{0},a_{1},a_{21},a_{22},a_{3})$ are also easily found from eqs. (45) and (38), we omit them for brevity.

Another Hamiltonian we consider reads

[TABLE]

The SSEρ for this Hamiltonian is given by

[TABLE]

The eigenvalues are $(3,5,8)$.

For an arbitrary number of spins $1$ an ansatz for states invariant under rotations, permutations and time reversal has a form analogous to eq. (37): It contains symmetric polynomials in scalar products of spin operators, each spin entering each term of this polynomial at most twice. If the system lacks the permutation symmetry, the polynomials need not be symmetric. This ansatz can be readily used to obtain a SSEρ for any number of spins.

5 Summary

We have studied the properties and merits of the stationary Schrödinger equation (2) with density matrices instead of wave functions. This equation produces the same spectrum of eigenvalues as the conventional Schrödinger equation. The main advantage of eq. (2) shows up when the Hamiltonian is invariant under some symmetry group which induces degeneracies of the spectrum. In this case for any eigenenergy one can choose a solution of eq. (2) which is invariant under the symmetry group. This is in contrast to the conventional Schrödinger equation, where the spontaneous symmetry breaking can prevent one from finding an invariant eigenvector. We have exemplified eq. (2) by applying it to a system of spins $1/2$ with the Heisenberg interactions on an arbitrary lattice. Further, we have applied an equation analogous to eq. (2) to find invariant subspaces of the operator of the total spin of $N$ spins $1/2$. Finally, we outlined how the same technique can be applied to higher spins. We conclude by a remark that it can be interesting to extend the methods of the present paper to the time-dependent Schrödinger equation and to master equations describing evolution of open quantum systems. The latter topic is addressed in a spirit somewhat similar to that in the present paper in the recent article [21].

Acknowledgements.

The work was supported by the Russian Science Foundation under the grant No. 17-71-20158.