Merits of using density matrices instead of wave functions in the stationary Schrödinger equation for systems with symmetries

The conventional form of the stationary Schrödinger equation (which we refer to as SSE$_\Psi$) reads

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ket}[1]{\left|{#1}\right\rangle} \displaystyle \label{SSE} H \ket{\Psi} = E \ket{\Psi}, \nonumber \end{align} \tag{ 1 }$$

where H is the Hamiltonian of a quantum system, E is its eigenenergy and $\newcommand{\ket}[1]{\left|{#1}\right\rangle} \ket{\Psi}$ is the corresponding eigenvector. Obviously, this equation implies an operator identity $\newcommand{\bra}[1]{\left\langle{#1}\right|} \newcommand{\ket}[1]{\left|{#1}\right\rangle} H \ket{\Psi}\bra{\Psi} = E \ket{\Psi}\bra{\Psi}$, where $\newcommand{\bra}[1]{\left\langle{#1}\right|} \newcommand{\ket}[1]{\left|{#1}\right\rangle} \ket{\Psi}\bra{\Psi}$ is the projection onto the eigenvector $\newcommand{\ket}[1]{\left|{#1}\right\rangle} \ket{\Psi}$. 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:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{SSErho} H \rho = E \rho. \nonumber \end{align} \tag{ 2 }$$

In the present paper we explore the properties and the merits of this equation which we refer to as SSE$_\rho$. 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$_\rho$ to a particular spin system in section 2.

We remind that a density matrix should satisfies three conditions,

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \rho^{\dagger}=\rho,~~{\rm tr} \rho=1,~~\rho>0. \label{rho conditions} \nonumber \end{align} \tag{ 3 }$$

The following relations between the SSE$_\Psi$ and SSE$_\rho$ follow immediately.

• (i)  
  Equations (1) and (2) share the same set of eigenvalues E.
• (ii)  
  If a given eigenvalue E is nondegenerate, then the corresponding $\newcommand{\ket}[1]{\left|{#1}\right\rangle} \ket{\Psi}$ and $\rho$ are related according to $\newcommand{\bra}[1]{\left\langle{#1}\right|} \newcommand{\ket}[1]{\left|{#1}\right\rangle} \rho=\ket{\Psi}\bra{\Psi}$.
• (iii)  
  In the case of degeneracy any solution of the SSE$_\rho$ reads

  $$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ket}[1]{\left|{#1}\right\rangle} \newcommand{\bra}[1]{\left\langle{#1}\right|} \displaystyle \label{general solution} \rho = \sum_{i}^{} p_i \ket{\Psi_i} \bra{\Psi_i},~~p_i\geqslant0, ~~\sum_i p_i=1, \nonumber \end{align} \tag{ 4 }$$

  where vectors $\newcommand{\ket}[1]{\left|{#1}\right\rangle} \ket{\Psi_i}$ constitute a basis in the corresponding degenerate subspace of the Hamiltonian.

Properties (ii) implies that in the nondegenerate case SSE$_\Psi$ and SSE$_\rho$ are, in fact, identical up to notations (however, even in this case SSE$_\rho$ can be more convenient compared to SSE$_\Psi$, in particular for spin systems, see e.g. [3]). An important advantage of the SSE$_\rho$ 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.

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle U\,H\,U^{\dagger} = H~~~\forall ~U\in G, \nonumber \end{align} \tag{ 5 }$$

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$_\Psi$, one can always avoid the spontaneous symmetry breaking in the solutions of SSE$_\rho$, according to the following simple lemma.

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 equation (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 ${\rm d}\mu(U)$ [5, 6]:

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \label{averaged} \rho_G=\int_G U\rho\, U^{\dagger} {\rm d}\mu(U) \nonumber \end{align} \tag{ 6 }$$

where the normalization condition $\int_G 1\,{\rm 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. $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}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–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 equation (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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{anti-Hermitian} [H,\rho]=0. \nonumber \end{align} \tag{ 7 }$$

In fact, this equation holds not only for the solutions of the SSE$_{\rho}$ but for any stationary state, i.e. a state not evolving under the Liouville–von Neumann equation. Equation (7) is widely used to obtain constraints on expectation values of various observables in equilibrium [11–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$_\rho$ by considering specific spin systems.

In the present section we specialise the SSE$_{\rho}$ for a system of N spins with the Heisenberg interaction. The Hamiltonian of this system reads

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{H Heisenberg} H = \sum_{i<j}^{} J_{ij} \left( {{{\boldsymbol \sigma}_i}{{\boldsymbol \sigma}_j}} \right) ,~~~i,j=1,2,...,N, \nonumber \end{align} \tag{ 8 }$$

where ${\boldsymbol \sigma}_i$ is the vector consisting of three Pauli matrices of the ith spin, J_ij is the coupling constant between ith and j th spins and $\left( {{{\boldsymbol \sigma}_i}{{\boldsymbol \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 ${{\mathcal A}}$ enumerating the set of pairs $(i_p,j_p)$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{multi-index} {{\mathcal A}}=\left( {{i_1},{j_1}} \right) \cdots \left( {{i_m},{j_m}} \right),~~1\leqslant m\leqslant [N/2], \nonumber \end{align} \tag{ 9 }$$

where [N/2] is the integer part of N/2, while i_p and j _p enumerate spins and for any p  satisfy

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{condiditons on indices} i_p<i_{p+1}, ~~j_p>i_p, ~~ j_p \ne i_l, j_l~~\forall~l \ne p,~~1\leqslant i_p,\,j_p\leqslant N. \nonumber \end{align} \tag{ 10 }$$

These conditions ensure that the sum over ${{\mathcal 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 ${{\mathcal A}}$ by $|{{\mathcal A}}|$ (e.g. $|{{\mathcal A}}|=m$ in equation (9)). Finally, we define an operator $A_{{\mathcal A}}$ as a product of $|{{\mathcal A}}|$ scalar products of Pauli matrices according to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{A} A_{{\mathcal A}}=\left( {{{\boldsymbol \sigma} _{{i_1}}}{{\boldsymbol \sigma} _{{j_1}}}} \right)\,\left( {{{\boldsymbol \sigma} _{{i_2}}}{{\boldsymbol \sigma} _{{j_2}}}} \right)\, ...\, \left( {{{\boldsymbol \sigma} _{{i_m}}}{{\boldsymbol \sigma} _{{j_m}}}} \right), \nonumber \end{align} \tag{ 11 }$$

where ${{\mathcal A}}$ is given by equation (9). We supplement this definition by a convention $\newcommand{\e}{{\rm e}} A_0\equiv 1$.

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{ansatz} \rho=\frac{1}{2^N}\left(a_0 A_0+ \sum_{{\mathcal A}} a_{{\mathcal A}} A_{{\mathcal A}}\right). \nonumber \end{align} \tag{ 12 }$$

Here $a_{{\mathcal A}}$ are numerical coefficients and the sum is over all sets ${{\mathcal A}}$ of the form (9) and (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

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ssigma}{\boldsymbol{\sigma}} \displaystyle {({{\boldsymbol{\sigma }}_1}{{\ssigma}_2})^2} &= \;\; 3 - 2({{\ssigma}_1}{{\ssigma}_2}),\nonumber \\ ({{\ssigma}_1}{{\ssigma}_2})({{\ssigma}_2}{{\ssigma}_3}) &= \;\;\;({{\ssigma}_1}{{\ssigma}_3}) - {\rm i}({{\ssigma}_1}{{\ssigma}_2}{{\ssigma}_3}),\label{algebraic relations} \nonumber \end{align} \tag{ 13 }$$

one can avoid terms with repeating spin indexes, as is indeed the case in equation (12). Here $\newcommand{\ssigma}{\boldsymbol{\sigma}} ({{\ssigma}_1}{{\ssigma}_2}{{\ssigma}_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 $\newcommand{\A}{{{\mathcal A}}} a_\A$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{equation Heisenberg} E\,a_{0} = 3 \sum_{i<j}{J_{ij} a_{(i,j)} }, \nonumber \end{align} \tag{ 14 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\A}{{{\mathcal A}}} \displaystyle \label{equation Heisenberg 2} E\,a_{\A}=& \, \sum_{l=1}^{|\A|-1}\sum_{m=l+1}^{|\A|} \left( (J_{i_m j_l}+J_{i_l j_m}-J_{i_l j_l}-J_{i_m j_m})a_{(i_l,i_m)(\,j_l,j_m)\A^{\overline{l,m}}} \,\, \right. \nonumber \\ &+\left. (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)\A^{\overline{l,m}}} \right) \nonumber \\ & +~\sum_{{\substack{p,\,q:\, p<q \nonumber \\ \, p,q\notin \A }}} \,\, \left( 3\,J_{pq}\, a_{(\,p,q)\A} + \sum_{m=1}^{|\A|} J_{pq} (a_{(i_m,p)(q,j_m)\A^{\overline{m}}}+a_{(q,i_m)(\,j_m,p)\A^{\overline{m}}}) \right) \nonumber \\ & + \sum_{p\notin \A} \,\, \sum_{m=1}^{|\A|} \left(J_{p j_m} a_{(\,p,i_m) \A^{\overline m}} + J_{p i_m} a_{(\,p,j_m) \A^{\overline m}} \right) + \sum_{m=1}^{|\A|} {J_{i_m j_m} \left( a_{\A^{\overline m}}-2a_{\A} \right)}, \nonumber \end{align} \tag{ 15 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\A}{{{\mathcal A}}} \displaystyle \label{anti-Hermitian Heisenberg} 0 &= \sum\limits_{p \notin \A}^{} {\left[ {({J_{pj}} - {J_{pi}}){a_{\left( {ij} \right)\left( {pk} \right)\A }} + ({J_{pk}} - {J_{pj}}){a_{\left( {ip} \right)\left( {jk} \right)\A^{}}} + ({J_{pi}} - {J_{pk}}){a_{\left( {ik} \right)\left( {pj} \right)\A^{}}}} \right]} \nonumber \\ &- \left( {({J_{ik}} - {J_{jk}}){a_{\left( {ij} \right)\A^{}}} + ({J_{jk}} - {J_{ij}}){a_{\left( {ik} \right)\A^{}}} + ({J_{ij}} - {J_{ik}}){a_{\left( {jk} \right)\A^{}}}} \right), ~{\forall~i,j,k} \notin \A ~{\rm and } ~i<j<k. \nonumber \end{align} \tag{ 16 }$$

The multi-index $\newcommand{\A}{{{\mathcal A}}} (i_l,j_m)(\,j_l,i_m)\A^{\overline{l,m}}$ is obtained from $\newcommand{\A}{{{\mathcal A}}} \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 $\newcommand{\A}{{{\mathcal A}}} \A$ analogously.

Note that equation (16) is due to equation (7). While equations (14) and (15) represent an eigenvalue problem, equation (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 equation (16), and then solve an eigenvalue problem with a smaller number of variables. We also note that an equation similar to equation (16) was used to explicitly construct local integrals of motions for the integrable instance of the Heisenberg Hamiltonian (8) [17, 18].

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 [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$_\rho$ to the system of three spins:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ssigma}{\boldsymbol{\sigma}} \displaystyle H = {J_{12}}\left( {{\ssigma _1}{\ssigma _2}} \right) + {J_{23}}\left( {{\ssigma _2}{\ssigma _3}} \right) + {J_{13}}\left( {{\ssigma _1}{\ssigma _3}} \right) \label{H 3 spins}, \nonumber \end{align} \tag{ 17 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ssigma}{\boldsymbol{\sigma}} \displaystyle 8\rho = a_0 + {a_{12}}\left( {{\ssigma _1}{\ssigma _2}} \right) + {a_{23}}\left( {{\ssigma _2}{\ssigma _3}} \right) + {a_{13}}\left( {{\ssigma _1}{\ssigma _3}} \right). \nonumber \end{align} \tag{ 18 }$$

In this case equations (14), (15) and (16) read

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eigenproblem 3 spins} \left(\begin{array}{@{}cccc@{}} 0 & 3 J_{12} & 3 J_{23} & 3 J_{13} \nonumber \\ J_{12} & -2 J_{12} & J_{23} & J_{13} \nonumber \\ J_{23} & J_{12} & -2 J_{23} & J_{13} \nonumber \\ J_{13} & J_{12} & J_{23} & -2 J_{13} \end{array}\right) \left(\begin{array}{@{}c@{}} a_0 \nonumber \\ a_{12} \nonumber \\ a_{23} \nonumber \\ a_{13} \end{array}\right) =E \left(\begin{array}{@{}c@{}} a_0 \nonumber \\ a_{12} \nonumber \\ a_{23} \nonumber \\ a_{13} \end{array}\right), \nonumber \end{align} \tag{ 19 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{antihermitian 3 spins} 0 = {a_{12}}\left( {{J_{23}} - {J_{13}}} \right) + {a_{23}}\left( {{J_{13}} - {J_{12}}} \right) + {a_{13}}\left( {{J_{12}} - {J_{23}}} \right). \nonumber \end{align} \tag{ 20 }$$

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, 20]).

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle E = 3,~~~a_0=1,~{{\rm }}{a_{12}} = {a_{23}} = {a_{13}} = \frac{1}{3} \label{E=3}; \nonumber \end{align} \tag{ 21 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle E ~= - 3,~~a_0=1,~{{\rm }}{a_{12}} + {a_{23}} + {a_{13}} = - 1 ,\label{E=-3} \nonumber \end{align} \tag{ 22 }$$

see figure 1 for illustration.

Zoom In Zoom Out Reset image size

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Ssquared} {\boldsymbol S^2}\rho = \lambda \rho, \nonumber \end{align} \tag{ 23 }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{rho for Ssquared} \rho = \frac{1}{{{2^N}}}\left(1+ \sum_{m=1}^{[N/2]}{a_m A_m}\right), \nonumber \end{align} \tag{ 24 }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {A_1}{A_m} &= (N - 2m + 2)(N - 2m + 1)\left( {\frac{3}{2} + (m - 1)\theta (N - 2m+1)} \right){A_{m - 1}} \nonumber \\ &+ 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]. \nonumber \end{align} \tag{ 25 }$$

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. It leads to the recursive formula for the coefficients a_m:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {a_1} &= \frac{{4\lambda - 3N}}{{3N(N - 1)}},\nonumber \\ {a_2} &= \frac{{(4\lambda - 7N + 12){a_1} - 2}}{{5(N - 2)(N - 3)}}\theta (N - 3),\nonumber \\ {a_m} &= \frac{{(4\lambda - (4m - 1)N + 8{m^2} - 12m + 4){a_{m - 1}} - 2(m - 1){a_{m - 2}}}}{{(N - 2m + 2)(N - 2m + 1)(2m + 1)}},~~m=3,4,...[N/2]. \nonumber \end{align} \tag{ 26 }$$

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

In this section we briefly outline how SSE$_\rho$ 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,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Sb}{\boldsymbol{S}} \displaystyle H=(\Sb_1 \Sb_2), \nonumber \end{align} \tag{ 27 }$$

where $\newcommand{\Sb}{\boldsymbol{S}} \Sb_i$ is the spin at ith site, $\newcommand{\Sb}{\boldsymbol{S}} (\Sb_i \Sb_i)=2$. A general form of the rotationally-invariant density matrix reads

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho = a_0 + a_1 \left( {{{\boldsymbol S}_1}{{\boldsymbol S}_2}} \right) + a_2 \left( {{{\boldsymbol S}_1}{{\boldsymbol S}_2}} \right)^2. \nonumber \end{align} \tag{ 28 }$$

Higher powers of the scalar product $\newcommand{\Sb}{\boldsymbol{S}} (\Sb_1\Sb_2)$ are linearly dependent on the first two powers according to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left( {{{\boldsymbol S}_i}{{\boldsymbol S}_j}} \right)^3 = 2 +\left( {{{\boldsymbol S}_i}{{\boldsymbol S}_j}} \right) - 2 \left( {{{\boldsymbol S}_i}{{\boldsymbol S}_j}} \right)^2. \nonumber \end{align} \tag{ 29 }$$

The $SSE_{\rho}$ then reads

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{two spins1} \left( \begin{array}{@{}ccc@{}} 0 & 0 & 2 \nonumber \\ 1 & 0 & 1 \nonumber \\ 0 & 1 & -2 \end{array} \right) \left(\begin{array}{@{}c@{}} a_0 \nonumber \\ a_1 \nonumber \\ a_2 \end{array}\right) =E \left(\begin{array}{@{}c@{}} a_0 \nonumber \\ a_1 \nonumber \\ a_2 \end{array}\right) \nonumber \end{align} \tag{ 30 }$$

with the solution

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle E &= 1, \quad\qquad (a_0,a_1,a_2) =(1/15,1/10,1/30), \nonumber \\ E &~= - 1, \qquad (a_0,a_1,a_2)=(1/3,-1/6,-1/6),\nonumber \\ E &~= - 2, \qquad (a_0,a_1,a_2) =(-1/3,0,1/3) , \nonumber \end{align} \tag{ 31 }$$

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

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{ansatz spins 1} \rho = a_0 + a_1 \sum_{i<j} \left( {{{\boldsymbol S}_i}{{\boldsymbol S}_j}} \right) + a_{21} \sum_{i\neq j\neq k} \left( {{{\boldsymbol S}_i}{{\boldsymbol S}_j}}\right) \left( {{{\boldsymbol S}_j}{{\boldsymbol S}_k}}\right) + a_{22} \sum_{i<j} \left( {{{\boldsymbol S}_i}{{\boldsymbol S}_j}} \right)^2 + a_3 \sum_{i\neq j\neq k} \left( {{{\boldsymbol S}_i}{{\boldsymbol S}_j}} \right) \left( {{{\boldsymbol S}_j}{{\boldsymbol S}_k}} \right)\left( {{{\boldsymbol S}_i}{{\boldsymbol S}_k}} \right), \nonumber \end{align} \tag{ 32 }$$

where $\sum\nolimits_{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 equation (32). The normalization condition implies

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{normalization spins 1} 27a_0+108 a_{22}+144 a_3 =1. \nonumber \end{align} \tag{ 33 }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Sb}{\boldsymbol{S}} \displaystyle H=(\Sb_1 \Sb_2)+(\Sb_2 \Sb_3)+(\Sb_3 \Sb_1). \nonumber \end{align} \tag{ 34 }$$

The SSE$_{\rho}$ for this Hamiltonian reads

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{SSErho spins 1} \left( \begin{array}{@{}ccccc@{}} 0 & 0 & 0 & 6 & 0 \nonumber \\ 1 & 0 & 4 & 3 & 4 \nonumber \\ 0 & 1 & 0 & 0 & 2 \nonumber \\ 0 & 1 & 0 & -2 & 2 \nonumber \\ 0 & 0 & 1 & 0 & -2 \end{array} \right) \left(\begin{array}{@{}c@{}} a_0 \nonumber \\ a_1 \nonumber \\ a_{21} \nonumber \\a_{22} \nonumber \\a_{3} \end{array}\right) =E \left(\begin{array}{@{}c@{}} a_0 \nonumber \\ a_1 \nonumber \\ a_{21} \nonumber \\a_{22} \nonumber \\a_{3} \end{array}\right). \nonumber \end{align} \tag{ 35 }$$

The eigenvalues are $(-3,-2,0,3)$. The coefficient $(a_0,a_1,a_{21},a_{22}, a_3)$ are also easily found from equations (35) and (33), we omit them for brevity.

Another Hamiltonian we consider reads

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Sb}{\boldsymbol{S}} \displaystyle H=(\Sb_1 \Sb_2)^2+(\Sb_2 \Sb_3)^2+(\Sb_3 \Sb_1)^2. \nonumber \end{align} \tag{ 36 }$$

The SSE$_{\rho}$ for this Hamiltonian is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{three spins1} \left( \begin{array}{@{}ccccc@{}} 0 & 6 & 24 & -36 & -72 \nonumber \\ 0 & 3 & -2 & 2 & 4 \nonumber \\ 0 & 0 & 3 & 0 & 0 \nonumber \\ 1 & -2 & -8 & 15 & 24 \nonumber \\ 0 & 0 & 1 & -1 & 1 \end{array} \right) \left(\begin{array}{@{}c@{}} a_0 \nonumber \\ a_1 \nonumber \\ a_{21} \nonumber \\a_{22} \nonumber \\a_{3} \end{array}\right) =E \left(\begin{array}{@{}c@{}} a_0 \nonumber \\ a_1 \nonumber \\ a_{21} \nonumber \\a_{22} \nonumber \\a_{3} \end{array}\right). \nonumber \end{align} \tag{ 37 }$$

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 equation (32): 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$_\rho$ for any number of spins.

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 equation (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 equation (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 equation (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 equation (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].

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