The Quantum Cluster Variational Method and the Phase Diagram of the quantum ferromagnetic $J_{1}$-$J_{2}$ model

In most cases of practical interest the solution of problems involving many interacting quantum particles rests on some kind of approximation Sachdev (1999); Fradkin (2013); Shankar (2017). It is not surprising then that an important fraction of the theoretical work in Condensed Matter Physics is concerned with the development of such approximations and the test of the predictions derived from them. Many of these techniques are similar in spirit, they improve over the simplest Mean Field approximation to the problem including correlations, and/or interactions between clusters of variables. Among the most celebrate approximation, although certainly not the only ones, we could mention, the Effective Field Theory (EFT) Bobák et al. (2018), the Cluster Mean Field (CMF) Jin et al. (2013); Ren et al. (2014) and the Correlated Cluster Mean Field (CCMF) Yamamoto (2008); F.M. et al. (2016); Singhania and Kumar (2018). The Quantum Cluster Variational method (QCVM)Domínguez and Mulet (2018), inserts into this family in a novel way: it allows the presence of disorder in the model, establishing a clear distinction between average case scenarios and computations on single instances, and connecting with message passing algorithms developed in Computer Science and Information TheoryMézard and Montanari (2009).

QCVM has its roots in previous results from Morita and Tanaka Morita (1994); Morita and Tanaka (1994); Tanaka et al. (1994) that starting from a pure variational approximation to the free energy of finite dimensional systems derived a set of close equations for the order parameters of a quantum problem defined by a Hamiltonian. Within this formulation, the resulting set of equations is usually solved assuming the existence of specific symmetries for the order parameter of the model. QCVM goes a step further, extending an approach previously developed for classical disordered models Rizzo et al. (2010); Domínguez et al. (2011); Lage-Castellanos et al. (2013) to connect these variational approximations with message passing equations in finite dimensional systems. The solutions of these equations can be computed either self-consistently in specific instances of the problem or on average over the disorderDomínguez and Mulet (2018); Lage-Castellanos et al. (2013). QCVM also generalizes, to finite dimensional lattices, the work done by Evans and Stephens (2008); Leifer and Poulin (2008); Poulin and Bilgin (2008); Bilgin and Poulin (2009); Poulin and Hastings (2011); Jouneghani et al. (2014); Biazzo and Ramezanpour (2013, 2014), that used message passing algorithms to solve quantum disordered problems in systems with tree-like topologies.

In this work we exploit the QCVM to study the phase diagram of the quantum $J_{1}$-$J_{2}$ model in the presence of external fields. This model belongs to an extensively studied family of similar systems where competing interactions induce exotic phases at low temperatures. This includes systems with short range ferromagnetic and long range dipolar interactions De’ Bell et al. (2000); Díaz-Méndez and Mulet (2010); Fey and Schmidt (2016); Fey et al. (2019); Mendoza-Coto et al. (2015, 2017, 2020); Cannas et al. (2006), long range Coulomb interactions Nandkishore and Sondhi (2010) and two-dimensional dipolar Fermi gasesOganesyan and Huse (2007); Bruun and Taylor (2008); Parish and Marchetti (2012).

The $J_{1}$-$J_{2}$ model has captured special attention both because of its apparent simplicity, and because of its resemblance with real materials. The classical version has been studied extensivelyMorán-López et al. (1993); Moran-Lopez et al. (1994); Guerrero et al. (2015); Domínguez et al. (2019); dos Anjos et al. (2008); Kalz et al. (2009, 2011); Jin et al. (2012) but the quantum version has turned harder to crack. Mainly, because quantum fluctuations add complexity to the inherent frustration of the model inducing a zoology of exotic phases: Stripes like phases, columnar anti-ferromagnetic (CAF), Néel Anti-ferromagnetic (NAF) phases, Spin-liquid phases and Spin Nematic Phases have been theoretically predicted Kohama et al. (2019); Shannon et al. (2004); Mezzacapo (2012); Jiang et al. (2012); Guerrero et al. (2015); Dagotto and Moreo (1989); Shindou et al. (2013); Cysne and Neto (2015). Some of them, like stripes and Néel Antiferromagnets, have been experimentally observed Carretta et al. (2002); Melzi et al. (2001); Nath et al. (2008) while others, like the Spin Nematic phase and the Spin-liquid phase lack of conclusive experimental observations but keep the interest of experimentalists Shannon et al. (2004); Orlova et al. (2017).

Our intention in this work is twofold. On one hand to test the Quantum Cluster Variational Method in a quantum model with competing interactions. On the other, to unveil the effect of quantum fluctuations on the phase diagram of the ferromagnetic $J_{1}$-$J_{2}$ model with competitive interactions. We study this model extensively as a function of the ratio $g=J_{2}/J_{1}$, the temperature $T$ and external fields. We show that the approximation captures a rich phenomenology of phases and phase transitions.

The work is organized as follows. In the next section we present the model and summarize some of the known results of the literature. We continue introducing the Quantum Cluster Variational Method and how it applies to the $J_{1}$-$J_{2}$ model. In this section we also discuss how QCVM is connected with other Mean Field Methods in the literature. Then, we show and discuss the results of our computations for the $J_{1}$-$J_{2}$ model. In the last section we present the conclusions of our work.

The quantum $J_{1}$-$J_{2}$ model is one of the most studied models of magnetic materials displaying frustration and is defined by the Hamiltonian

where $\langle ij\rangle$ stand for the Nearest Neighbors in the square lattice and $\langle\langle ij\rangle\rangle$ for the Next Nearest Neighbors and the external field $\vec{h_{i}}$ may point to an arbitrary direction. We focus here our attention in the model with nearest-neighbour interaction $J_{1}$ greater than zero (ferromagnetic), and $J_{2}\leq 0$. The later induces the frustration into the model. Since the preferred direction of the system is labeled x, quantum effects appear when the transverse field $h_{z}$ is different from zero.

The classical model ($h_{z}=0$) has been largely studiedMorán-López et al. (1993); Moran-Lopez et al. (1994); Guerrero et al. (2015); Domínguez et al. (2019); dos Anjos et al. (2008); Kalz et al. (2009, 2011); Jin et al. (2012). However, also there some questions remain open. For example, while for values of $g=\frac{|J_{2}|}{J_{1}}$ slightly greater than $0.5$ the transition is discontinuous, Monte Carlo studiesKalz et al. (2009, 2011) suggest a continuous transition for $g\sim 1$. In a paper by Jin et al.Jin et al. (2012), the existence of a pseudo first order transition in the range $g_{c}\leq g\lesssim 1$ is reported, with $\,g_{c}=0.67$. They found that the critical exponents vary continuously between those of the 4-state Potts model at $g=g_{c}$ to those of the Ising model for $g\rightarrow\infty$.

Recently, researchers paid attention also to the effect of external magnetic fields and demonstrated the presence of a nematic phaseGuerrero et al. (2015); Domínguez et al. (2019), with slight differences, at low temperatures, between the two studies. This nematic phase is characterized by the presence of orientational order and the lack of positional order and has been reported experimentally Orlova et al. (2017).

In the present work we report, exploiting QCVM, our own perception of this problem, and show that there is no clear separation between pseudo and actual first order behavior. Moreover, we studied the combined effect of transverse and longitudinal field to understand the role played by quantum or thermal fluctuations in the order of these transitions.

We explore the behaviour of three order parametersGuerrero et al. (2015); Domínguez et al. (2019). First, we define two positional order parameters which measures the breaking in the translational symmetry of the lattice in two possible directions:

Where $m^{x}_{i}$ refers to the mean magnetization along the preferred direction of the system, of the site i, inside a given plaquette, following the convention showed in Fig.1.

While translational order parameters are enough to determine the presence of stripes, they are not suited to detect the breaking of the rotational symmetry of the lattice. They take the same values in a paramagnetic and a nematic phase, yet the nematic phase shows a structure of the correlations similar to that of the stripes. This is natural if we think about this phase as an intermediate one in between the completely symmetric paramagnetic phase and the ordered stripes phase. Therefore, to characterize the nematic behavior we introduced an orientational order parameter defined as

where $l_{ij}$ represents the correlation between sites $i$ and $j$ in the plaquette. All the order parameters are defined in a way that $0<\{Q,M_{x},M_{y}\}<1$.

In order to compute the order parameters relevant for the $J_{1}$-$J_{2}$ model it is fundamental to obtain a good approximation of at least local magnetizations and correlations of neighboring spins. These quantities can be found from a knowledge of the local density matrices describing the joint statistics of these spins.

With a cluster expansion the free energy is approximated as a sum of contributions of selected regions or groups of variables. Each term of the sum depends on the local density matrix of the corresponding variables (spins). These regions are typically chosen to balance the accuracy and the computational cost of the calculations. The resulting approximated expression is to be understood as variational in terms of the local densities. The quality of the final results depends directly on the choice of regions, their size in relation to the correlation length, and how much mean-field is put the treatment of effective interactions and correlation between regions.

Among the many cluster approximations appearing in the literature, the approach known as Cluster Variational Method or CVM is remarkable because of the systematic treatment of different levels of complexity Yedidia et al. (2005); Kikuchi (1951). In a natural way one can choose free energy expansions ranging from naive mean field, a Bethe approximation or more involved structures. A quantum extension of such ideas, QCVM, was developed in Morita (1994); Morita and Tanaka (1994); Tanaka et al. (1994) and more recently expanded in Domínguez and Mulet (2018). In particular we are interested in a construction called Kikuchi approximation. The Kikuchi free energy $F_{\small\mbox{Kik}}$ consists of a sum of local free energies of all the square unit cells ($p$), the interacting pairs ($l$) formed by the overlap of neighbor cells¹¹1The diagonal interactions do not belong to any overlap region and therefore, according to the CVM, they do not appear explicitly as independent terms in the approximated free energy. and all the contributions from single spins ($s$):

The local free energy terms are defined the usual way

where $r$ may label a plaquette, link or spin region. The numbers $c_{p},c_{l},c_{s}$ are weight factors that compensate the effect of overlap between regions. Due to overlaps, the contribution of some variables might be under or over represented if the weight factors were absent. For the case of interest here, $c_{p}=1,c_{l}=-1,c_{s}=1$.

The minimization of the functional $F_{\small\mbox{Kik}}$ should be done under certain constraints. In addition to the usual normalization conditions, $\text{Tr}[\hat{\rho}_{r}]=1$, the consistency between different regions must be enforced explicitly. For all plaquettes $p=\langle ijkm\rangle$ and any of its forming links $l=\langle ij\rangle$ one could demand that the partial trace over spins $k$ and $m$ equals the link density $\hat{\rho}_{l}$:

Moreover, every link distribution $\hat{\rho}_{l}$ could be forced to marginalize onto the corresponding spin densities $\hat{\rho}_{i}$ and $\hat{\rho}_{j}$:

However, numerical resultsDomínguez and Mulet (2018) have shown that the above constraint scheme is too restrictive. It turns out that the minimization gives much more accurate results if consistency is enforced only in the direction of the spin-spin interactions. In this model, where interactions exist only in the OX direction, it amounts to relaxing (7) to the set of equations:

We are now in position to write a Lagrange function that includes the conditions (9) and (10) by means of suitable Lagrange multipliers $\hat{\lambda}$, $\hat{\gamma}$:

To conform to (9) and (10), it suffices to introduce the parametrization²²2Notice that $\lambda$ and $\gamma$ are actually operators on the Hilbert space of the pairs $l$ and the single spins $s$, respectively.:

It is enough to make $C,c,d$ all real to guarantee that $\hat{\lambda}$ and $\hat{\gamma}$ are Hermitian, and consequently, that $L$ is real. Moreover, this parametrization ensures that the operators $\hat{\rho}$ that satisfy the stationary condition are all positive, which is the fundamental property of density matrices (in addition to an unit trace, of course).

We will say that the value of the Lagrange function $L[\hat{\rho}_{r}]$ is stationary with respect to a change in the density $\hat{\rho}_{r}$ if the linear part of the expansion of $L[\hat{\rho}_{r}+\epsilon\delta\hat{\rho}_{r}]$, in powers of $\epsilon$, is zero. The operator $\delta\hat{\rho}_{r}$ is arbitrary but should be Hermitian, traceless and small enough to keep $\hat{\rho}_{r}+\epsilon\delta\hat{\rho}_{r}$ as a valid density operator, The Lagrange function will be stationary with respect to the set of all densities $\{\hat{\rho}_{p}\},\{\hat{\rho}_{l}\},\{\hat{\rho}_{s}\}$, by definition, if it is stationary in each one of them. The stationarity condition can be summarized in the following formula:

As an example, let us work out the form of the density matrix of an arbitrary spin $s_{0}$ at the stationary point. The part of $L$ that depends on a given $\hat{\rho}_{s_{0}}$ is

In (15) the only term that is somewhat involved to expand to first order in $\epsilon$, when evaluating $L[\hat{\rho}_{s_{0}}+\epsilon\delta\hat{\rho}_{s_{0}}]$, is the one with the logarithm. The trick is to use perturbation theory to write the eigenvalues $r_{i}(\epsilon)$ of $\hat{\rho}_{s_{0}}+\epsilon\delta\hat{\rho}_{s_{0}}$ using the eigenvalues and eigenvectors $\{r_{i},|r_{i}\rangle\}$ of $\hat{\rho}_{s_{0}}$:

The trace with the logarithm is therefore written as

In (18) we used that $\delta\hat{\rho}_{s_{0}}$ has a null trace. From the results above it is not hard to deduce the linear expansion of $L[\hat{\rho}_{s_{0}}+\epsilon\delta\hat{\rho}_{s_{0}}]$:

The stationarity condition applied to (19) implies the following relation

Now we use that $\delta\hat{\rho}_{s_{0}}$ is an arbitrary traceless Hermitian operator to state that the part in parenthesis in (20) is at most a constant times the identity, and therefore, $\hat{\rho}_{s_{0}}$ must have the form:

with $\beta=1/T$. The condition $\text{Tr}[\hat{\rho}_{s_{0}}]=1$ determines the value of $\mathcal{Z}_{s_{0}}$. The form of (21) resembles a Boltzmann distribution where an effective interaction given by the Lagrange multipliers $\hat{\gamma}_{l\rightarrow s_{0}}$ is added to the original local Hamiltonian.

The same procedure that leads to (21) can be used to find expressions for the rest of the density operators, that is, for pair and plaquette densities. The general structure will be the same: a Boltzmann distribution where the local Hamiltonian is modified by a number of effective interactions. These appear due to the consistency relations between regions and via the corresponding Lagrange multipliers. We skip these derivations and only write the final expressions for the Kikuchi approximation:

In (22) we used a linear transformation to write $\hat{\lambda}_{p\rightarrow l}$ and $\hat{\gamma}_{l\rightarrow i}$ as a function of new parameters $\hat{U}_{p\rightarrow l}$ and $\hat{u}_{p\rightarrow i}$, respectively. These changed variables will of course inherit the structure of the originals, see (12) and (13):

The form of (12), (13) and (23) suggests that the Lagrange multiplier can be interpreted as effective local fields and a contribution to the interaction constants between spins when their explicit form is put into the equations (22). The plaquette-to-pair $\hat{U}_{p\rightarrow l}$ is composed of a terms that adds to the $J_{1}\hat{\sigma}^{x}_{i}\hat{\sigma}^{x}_{j}$ part of the original local Hamiltonian, and two effective fields $u_{p\rightarrow i},u_{p\rightarrow j}$ acting on the spins of the pair $l=\langle ij\rangle$. On the other hand, the pair-to-spin $\hat{u}_{l\rightarrow s}$ is formed by an effective field on $\hat{\sigma}^{x}_{s}$, for every spin $s$ that belongs to $l$. The role of these extra fields is to guarantee the consistency between the distributions, according to the restrictions imposed to the variational function.

The nature of the constraints directly implies that no effective fields (or modification to interaction constants) act in the transverse (OZ) direction. The reason being that no consistency is demanded on this axis. As a consequence, moments and correlations in the OZ direction may differ if computed from different local densities. Since all the observables we study here are in OX, where consistency holds, this discrepancy is not a serious complication.

A graphical representation of the effective fields entering the formula (22) for the plaquette density operator $\hat{\rho}_{p}$ is shown in Figure (2). Effective fields are drawn as arrows, placed according to the regions they are related to. For example, the arrows standing parallel to the edges of the lattice (NN interactions) correspond to $\hat{u}_{l^{\prime}\rightarrow i}$ terms. The arrows (triplet) pointing from the center of a square to a border pair picture the $\hat{U}_{p^{\prime}\rightarrow l}$ operators. In the same way, Figures (3) and (4) provide the representation for distributions $\hat{\rho}_{l}$ and $\hat{\rho}_{s}$.

The final step to determine the local distributions is to find the value of the effective parameters. Their value is fixed by the solution of the set of equations imposed by the consistency relations (9) and (10). The numerical algorithm used to solve these coupled equations is very similar to the classical generalization of the well known Belief Propagation and is described in detail in Domínguez and Mulet (2018).

The basic scheme is the following. All effective fields are initialized in an arbitrary way (randomly, uniform, following a certain pattern, etc). Then, a series of update steps are performed until convergence is achieved³³3In certain cases the algorithm fails to converge. This will be discussed afterwards together with the numerical results.. We understand convergence as the situation where the values of the effective fields are such that all the consistency relations are satisfied simultaneously.

Let us describe a typical iteration for the field $\hat{U}_{p\rightarrow(14)}$ [upper triplet in Figure (3)]. First, compute $\hat{\rho}_{p}$ using the current value of the fields depicted in Figure (2). Then, from $\hat{\rho}_{p}$ determine the correlation $c_{14}=\langle\hat{\sigma}^{x}_{1}\hat{\sigma}^{x}_{4}\rangle$ and the magnetizations $m_{1}=\langle\hat{\sigma}^{x}_{1}\rangle$, $m_{4}=\langle\hat{\sigma}^{x}_{4}\rangle$. The crucial part is to use the consistency condition between $\hat{\rho}_{p}$ and $\hat{\rho}_{14}$. We should find a value for the three components of $\hat{U}_{p\rightarrow(14)}$ (this parameter enters the formula for $\hat{\rho}_{14}$ and not $\hat{\rho}_{p}$) such that the observables computed from the pair distribution $\hat{\rho}_{14}$ match the plaquette predictions:

The three equations above suffice to find the three parameters $\hat{U}_{p\rightarrow(14)}\sim(U_{p\rightarrow(14)},u_{p\rightarrow 1},u_{p\rightarrow 4})$. The newly found value of $\hat{U}_{p\rightarrow(14)}$ is then updated in the lattice⁴⁴4Actually, what is usually stored is a weighted sum of the new and the old values. This procedure introduces a kind of damping that might help convergence..

Updating a pair-to-spin field, for example $\hat{u}_{(14)\rightarrow 1}$ follows a similar recipe. First compute $\hat{\rho}_{l=(14)}$ using the fields of Figure (3). Then compute the spin 1 magnetization $\langle\hat{\sigma}^{x}_{1}\rangle$ using $\hat{\rho}_{l}$. Finally, find a value $u_{(14)\rightarrow 1}$ such that the same magnetization computed with $\hat{\rho}_{1}$ matches the value computed with the pair distribution.

The algorithm described here assumes that we are using a given realization of the model, that is, that we are computing all the effective fields for every region defined in a $N=L^{2}$ square latice. The procedure is general and can handle situations where the external fields applied may not be homogeneous or the interactions change across the lattice. However, since the model in question is translationally invariant, it is reasonable to expect a simplification of the calculations; one would hope this structure is reflected in a translational symmetry of the effective fields. This is indeed the case. If the same values repeat all over the lattice, it is not necessary to solve $O(N)$ equations but only a reduced set. The fields outside a certain region of the lattice can be taken as the same computed inside that region. In this case the numerical procedure looks like a self-consistent iteration. The idea, though, is not that every site in the lattice is equivalent to the rest, that is a too restrictive assumption that would allow only homogeneous states. The correct procedure is to consider the smallest possible structures that, by repetition, can create the patterns observed experimentally. In practice, to generate a pattern of stripes, or an anti-ferromagnetic checkerboard design, it is enough to consider a square plaquette as the elementary region. For the model studied here we followed both approaches, a specific realization of the lattice solving $N$ equations, and the more practical reduced version obtained equivalent results. For generality we will present only simulations results for the former method. On one hand it makes clearer the strength of QCVM for more general situations, on the other, makes more transparent when the algorithm does not converge.

To close this section we would like to make a few comments about our choice of the QCVM as the inference method used to compute the local observables of the $J_{1}$-$J_{2}$ model. Mean field cluster approximations frequently found in the literature Kellermann et al. (2019); Balcerzak et al. (2018); Bobák et al. (2018) are often based on some kind of variational argument on single site spin magnetizations. In many cases the effect is equivalent to factoring out the correlations as products of single site moments. A more consistent treatment is expected to take into account more structure into the relevant correlations among the system variables. That is precisely the formal advantage of QCVM: by means of the effective interactions the correlations have extra parameters that can be used for tuning. This can lead to more precision in the prediction of system properties. A quick test, for example, shows that the Cluster Mean Field of Kellermann et al. (2019) applied to a classical Ising model on a square lattice returns predictions for the para-ferro transition temperature, $T_{CMF}\approx 3.5$ that is closer to the naive MF prediction $T_{MF}=4.0$ than to the real value $T_{c}\approx 2.23$. Bethe-Pierls approximation or Kikuchi give much more accurate estimates (2.89 and 2.41 respectively).

Mean field treatments usually replace real interactions by the interaction with the mean value of certain variables. For example, terms like $J_{1}\hat{\sigma}^{z}_{1}\hat{\sigma}^{z}_{2}$ might be roughly approximated by $J_{1}\hat{\sigma}^{z}_{1}\langle\hat{\sigma}^{z}_{2}\rangle$ or something similar, which looks like an effective field acting on $\hat{\sigma}^{z}_{1}$. QCVM effective fields can be understood in a similar way. With this in mind, we would like to point to what might be a caveat in the choice of regions used in the Kikuchi approximation. Since the overlap of elementary plaquettes corresponds only to NN pairs, the effective fields will only appear for those interactions. Notice for example the scheme for the single site distribution $\hat{\rho}_{s}$ in Figure (4), where no diagonal effective fields are present. It would be desirable to have a setup where consistency is also forced for diagonal distributions which would then appear as diagonal effective fields. We believe this could be achieved if the elementary plaquettes are defined with a rhombic shape as shown in Figure (5).

In this work we studied the quantum $J_{1}$-$J_{2}$ model with nearest neighbor ferromagnetic interactions ($J_{1}\geq 0$) and next nearest neighbour anti-ferromagnetic interactions ($J_{2}\leq 0$) in the presence of a transverse field using the Quantum Cluster Variational Method. We study the model in an extensive lattice without particular assumptions about the symmetry of the problem. Our results show that quantum fluctuations change the order of the transition, the larger the quantum effects the wider is the range of parameters for which the transition is continuous. Moreover, quantum fluctuations may induce a gap between the ferromagnetic phase and the phase of stripes and in the presence of longitudinal fields also a pronounced nematic phase that penetrates the phase of stripes. This nematic phase is characterized by the presence of orientational order and the lack of translational order.