Interplay between Magnetic and Vestigial Nematic Orders in the Layered $J_1$-$J_2$ Classical Heisenberg Model

TL;DR

This paper investigates the phase transitions and interplay between magnetic and nematic orders in a layered $J_1$-$J_2$ classical Heisenberg model, revealing complex transition behaviors influenced by interlayer coupling.

Contribution

It provides a detailed phase diagram analysis of the layered $J_1$-$J_2$ model, highlighting the conditions for simultaneous or split magnetic and nematic transitions, and connects findings to iron-based superconductors.

Findings

Broad regions of magnetic order separated by first-order transition.

Existence of an intermediate nematic phase without magnetic order.

Magnetic correlation length depends sharply on nematic order parameter.

Abstract

We study the layered $J_1$-$J_2$ classical Heisenberg model on the square lattice using a self-consistent bond theory. We derive the phase diagram for fixed $J_1$ as a function of temperature $T$, $J_2$ and interplane coupling $J_z$. Broad regions of (anti)ferromagnetic and stripe order are found, and are separated by a first-order transition near $J_2\approx 0.5$ (in units of $|J_1|$). Within the stripe phase the magnetic and vestigial nematic transitions occur simultaneously in first-order fashion for strong $J_z$. For weaker $J_z$ there is in addition, for $J_2^*<J_2 < J_2^{**}$, an intermediate regime of split transitions implying a finite temperature region with nematic order but no long-range stripe magnetic order. In this split regime, the order of the transitions depends sensitively on the deviation from $J_2^*$ and $J_2^{**}$, with split second-order transitions predominating for $J_2^* \ll J_2 \ll J_2^{**}$. We find that the value of $J_2^*$ depends weakly on the interplane coupling and is just slightly larger than $0.5$ for $|J_z| \lesssim 0.01$. In contrast the value of $J_2^{**}$ increases quickly from $J_2^*$ at $|J_z| \lesssim 0.01$ as the interplane coupling is further reduced. In addition, the magnetic correlation length is shown to directly depend on the nematic order parameter and thus exhibits a sharp increase (or jump) upon entering the nematic phase. Our results are broadly consistent with predictions based on itinerant electron models of the iron-based superconductors in the normal-state, and thus help substantiate a classical spin framework for providing a phenomenological description of their magnetic properties.

Figures

Equations

$$H=J_1\sum _{⟨ij⟩}{\vec{S}}_i\cdot {\vec{S}}_j+J_2\sum _{⟨⟨ij⟩⟩}{\vec{S}}_i\cdot {\vec{S}}_j+J_z\sum _{\{ij\}}{\vec{S}}_i\cdot {\vec{S}}_j$$

$$H=\sum _{\vec{q}}{J}_{\vec{q}}{\vec{S}}_{\vec{q}}\cdot {\vec{S}}_{-\vec{q}},$$

$$J_{\vec{q}}=J_1\left(\cos q_x+\cos q_y\right)+2J_2\cos q_x\cos q_y+J_z\cos q_z-C.$$

$$\delta \left({\vec{S}}_{\vec{r}}\cdot {\vec{S}}_{\vec{r}}-1\right)=\int \frac{\beta d{\lambda }_{\vec{r}}}{2\pi }{e}^{-i\beta {\lambda }_{\vec{r}}\left({\vec{S}}_{\vec{r}}\cdot {\vec{S}}_{\vec{r}}-1\right)}$$

$$Z=\int D\vec{S}d\Delta D\lambda e^{-\beta {\sum }_{\vec{q},{\vec{q}}^{\prime }}\left({\mathbf{K}}_{\vec{q},{\vec{q}}^{\prime }}-{\mathbf{\Lambda }}_{\vec{q},{\vec{q}}^{\prime }}\right){\vec{S}}_{\vec{q}}^{\ast }\cdot {\vec{S}}_{{\vec{q}}^{\prime }}+V\beta \Delta }$$

$$Z=\int d\Delta D\lambda e^{-S_{\lambda }}$$

$$S_{\lambda }=\frac{N_s}{2}\text{mbox}Tr\ln \left(\mathbf{K}-\mathbf{\Lambda }\right)-\beta V\Delta .$$

$${\chi }_{\vec{q}}\equiv ⟨{\vec{S}}_{\vec{q}}^{\ast }\cdot {\vec{S}}_{\vec{q}}⟩=\frac{N_{s}T}{2}⟨{\left(\mathbf{K}-\mathbf{\Lambda }\right)}_{\vec{q},\vec{q}}^{-1}{⟩}_{S_{\lambda }}$$

$$D_{0\vec{q}}=\frac{2}{N_s}{\left[\sum _{\vec{p}}{K}_{\vec{p}+\vec{q}}^{-1}{K}_{\vec{p}}^{-1}\right]}^{-1}$$

$$K_{\mathrm{eff}}^{-1}{}_{\vec{q}}=\frac{1}{J_{\vec{q}}+\Delta -{\Sigma }_{\vec{q}}}.$$

$$D_{\vec{q}}^{-1}=D_{0\vec{q}}^{-1}-{\Pi }_{\vec{q}}.$$

$${\Sigma }_{\vec{q}}$$

$${\Pi }_{\vec{q}}$$

$$D_{\vec{q}}^{-1}=\frac{N_s}{2}\sum _{\vec{p}}{K}_{\mathrm{eff}}^{-1}{}_{\vec{p}+\vec{q}}{K}_{\mathrm{eff}}^{-1}{}_{\vec{p}}.$$

$${\chi }_{\vec{q}}=\frac{N_{s}T}{2}⟨K_{\mathrm{eff}}^{-1}{}_{\vec{q}}{⟩}_{S_{\Delta }}$$

$$\frac{N_{s}T}{2V}\sum _{\vec{q}}{K}_{\mathrm{eff}}^{-1}{}_{\vec{q}}=1,$$

$${\chi }_{\vec{q}}=\frac{N_{s}T}{2}{K}_{\mathrm{eff}}^{-1}{}_{\vec{q}},$$

$$T={\left[\frac{N_s}{2V}\sum _{\vec{q}}{K}_{\mathrm{eff}}^{-1}{}_{\vec{q}}\right]}^{-1},$$

$$I=\frac{1}{V}\sum _{\vec{r}}⟨{\vec{S}}_{\vec{r}}\cdot {\vec{S}}_{\vec{r}+\hat{x}}-{\vec{S}}_{\vec{r}}\cdot {\vec{S}}_{\vec{r}+\hat{y}}⟩$$

$$I=\frac{1}{2V}\sum _{\vec{q}}\left(\cos q_x-\cos q_y\right){\chi }_{\vec{q}}$$

$${\chi }_{\vec{Q}}={\vec{m}}_{\vec{Q}}\cdot {\vec{m}}_{-\vec{Q}}V+\sum _{\vec{r}}\delta f(\vec{r})e^{-i\vec{Q}\cdot \vec{r}}$$

$$M^2\equiv {\vec{m}}_{\vec{Q}}\cdot {\vec{m}}_{-\vec{Q}}=\frac{N_{s}T}{2V}{K}_{\mathrm{eff}}^{-1}{}_{\vec{Q}}=\frac{N_{s}T}{2V\Delta }.$$

$${\chi }_{\vec{r}}$$

$$+\frac{A^{\prime }}{{\xi }_1^{\prime 2}{q}_x^2+{\xi }_2^{\prime 2}(q_y-\pi {)}^2+{\xi }_z^{\prime 2}{q}_z^2+1})e^{i\vec{q}\cdot \vec{r}}$$

$${\chi }_{\vec{r}}$$

$${\chi }_{(x,0,0)}\propto f\left(\frac{x}{{\xi }_2}\right)\cos (\pi x)+f\left(\frac{x}{{\xi }_1}\right).$$

$${\chi }_{(0,y,0)}\propto f\left(\frac{y}{{\xi }_1}\right).$$

$${\chi }_{(x,0,0)}\propto f\left(\frac{x}{{\xi }_2}\right)\cos (\pi x).$$

$$\frac{A}{2\left(1-\cos (q_i-Q_i)\right){\xi }_i^2+1}$$

$$I$$

$$-\sum _{\vec{q}\sim (\pi ,0,0)}\frac{1}{J_{\mathrm{eff}\vec{q}}+{\Delta }_0(1+\alpha I)}],$$

$$\alpha=\left(\frac{N_{s}T\Delta_{0}}{2V}\sum_{\vec{q}}K_{{\rm eff}\,\vec{q}}^{-2}\right)^{-1}\Bigg{\rvert}_{T=T_{c}}.$$

Full text

I Introduction

The iron-based superconductors (FeSCs) represent an intriguing class of materials exhibiting both unconventional superconductivity as well as peculiar normal-state properties (see Refs. Paglione and Greene, 2010; Stewart, 2011; Dai et al., 2012 for reviews). In particular, in the normal-state there are strong indications that the stripe spin-density wave order with wavevector $\vec{Q}=(0,\pi)\,{\rm or}\,(\pi,0)$, that occurs below a critical temperature, is intimately related to the structural distortion of the lattice that occurs at a higher critical temperature. This structural distortion breaks the tetragonal symmetry, leading to long-ranged lattice-nematic order. Several proposals have thus suggested that spin fluctuations of the stripe order drive the nematic transitionFang et al. (2008); Xu et al. (2008); Fernandes et al. (2014), while the lattice distortion arises secondarily via coupling to the spin-driven nematic order parameter.

This scenario motivates a purely magnetic approach that captures the interplay between the vestigial nematic order and the fluctuations of the magnetic stripe order that drive it. Previous proposals have focused on the $J_{1}$-$J_{2}$ model of localized spins Yildirim (2008); Si and Abrahams (2008), which is known, as we discuss below, to have nematic order in the strictly two-dimensional (2D) limit at finite temperatures Chandra et al. (1990); Weber et al. (2003). This type of order has become known as Ising-nematic order Fernandes et al. (2014) because the order parameter space (the choice between $\vec{Q}=(\pi,0)$ or $\vec{Q}=(0,\pi)$) is effectively $Z_{2}$ due to the Mermin-Wagner theorem Mermin and Wagner (1966) that precludes spontaneous spin-rotation symmetry breaking in 2D. This Ising-nematic phase transition extends to the layered $J_{1}$-$J_{2}$ model where there is an additional phase transition to an ordered magnetic stripe phase Fang et al. (2008); Xu et al. (2008). Thus the layered $J_{1}$-$J_{2}$ model may serve as the simplest phenomenological model capable of describing the putative spin-fluctuation triggered nematic phase in the FeSCs.

Nevertheless, the fact that the FeSCs are not insulators indicate that the normal-state of the FeSCs is likely best described as an itinerant magnet at low doping Singh (2009), and a rich phase diagram based on itinerant magnetism has been predicted Fernandes et al. (2012). It is unclear whether the phase diagram of the layered $J_{1}$-$J_{2}$ model is as rich, or if and where the $J_{1}$-$J_{2}$ model consists of regions of split nematic and magnetic transitions. Moreover, the nature of the phase transitions is an important aspect of the problem that is known in the itinerant models to depend sensitively on the microscopic parameters near the bifurcation point where the simultaneous (first order) transition splits into two separate transitions. In particular, as the transitions split, there appears to be a narrow intermediate regime where the order of the transitions may be different for the magnetic and nematic order parameters. This can lead, e.g., to a metanematic transition at which the finite nematic order parameter jumps to a higher value as the magnetic order sets in. In addition, the itinerant models predict a sharp increase of the spin correlation length upon entering the nematic phase due to its direct dependence on the nematic order parameter Fernandes et al. (2012).

It is the purpose of this paper to determine whether such phenomena arise also within a specific microscopic model based on localized spins: the layered square lattice $J_{1}$-$J_{2}$ classical Heisenberg model. We study this model using the nematic bond theory developed in Ref. Schecter et al., 2017, which can detect nematic and magnetic orders independently and determine the order of the transitions. A dual purpose of this paper is also to explain this method in depth. The nematic bond theory allows us to construct (for the first time, to the best of the authors’ knowledge) the phase diagram of the layered $J_{1}$-$J_{2}$ model for fixed $J_{1}$ as a function of temperature $T$, $J_{2}$ and interlayer coupling $J_{z}$. The nematic bond theory can be employed to investigate temperature-dependent properties of any classical Heisenberg hamiltonian, and requires considerably less numerical efforts than Monte Carlo simulations. We show that practically all of the phenomena discovered in the itinerant electron models arise also in the Heisenberg model, including an intermediate regime of transition types near the bifurcation points, and a sharp increase of the spin correlation length upon entering the nematic phase.

A schematic phase-diagram illustrating the behavior of the order parameters in the various regimes is presented in Fig. 1. One sees that for $J_{1}/2\lesssim J_{2}<J_{2}^{*}$ the transition into the stripe phase is simultaneous and first-order, while for $J_{2}>J_{2}^{*}$ the transitions are split and predominantly second-order. Near the bifurcation point, $J_{2}\sim J_{2}^{*}$, one of the transitions may remain first-order. This leads (in the example shown in Fig. 1) to a metanematic transition, similar to what happens in the itinerant models for various parameter regimes. Although we will not attempt to make a direct connection to any FeSC compound in particular, our results show that at a phenomonological level the physics of the $J_{1}$-$J_{2}$ model indeed appears to capture the essential aspects of the normal-state magnetic properties of the FeSCs and therefore may be of use as a simplifying alternative to the itinerant approach.

While some previous studies have attempted to address the issues discussed above within a classical spin framework, it is important to emphasize that here we study directly the microscopic layered $J_{1}$-$J_{2}$ model. An effective model of the layered $J_{1}$-$J_{2}$ model, the Ising O(3)-model Chandra et al. (1990), has been studied previously using Monte Carlo simulations Kamiya et al. (2011) which yields a simultaneous first order phase transtion for large interlayer couplings and two split transitions for weaker interlayer couplings. Other studies have included a phenomenological biquadratic coupling to the Hamiltonian to mimic the longwavelength effective action that arises upon coarse graining Chandra et al. (1990); Fang et al. (2008); Xu et al. (2008). Our treatment is different in that we do not explicitly assume a separate Ising degree of freedom from the outset and do not insert a biquadratic coupling by hand. Instead, we shall utilize and develop a tractable technique, the nematic bond theory, that can tackle the $J_{1}$-$J_{2}$ model head-on and provide its phase diagram in terms of microscopic parameters.

The rest of the paper is organized as follows. We first define the Hamiltonian of the $J_{1}$-$J_{2}$–model, Sec. II, and then describe in details the method we use to solve it in Sec. III. Then we discuss in Sec. IV the order parameters relevant for the $J_{1}$-$J_{2}$–model and how to compute them using our method. The results for the $J_{1}$-$J_{2}$–model are then described in Sec. V, followed by Sec. VI that describes the spin correlations. We end by a general discussion in Sec. VII.

II $J_{1}$-$J_{2}$–model

The Hamiltonian of the layered classical $J_{1}$-$J_{2}$–model is

[TABLE]

where $\langle ij\rangle$,$\langle\langle ij\rangle\rangle$ denotes nearest and next-nearest in-plane neighbors respectively and $\{ij\}$ denotes nearest neighbors in adjacent layers on a cubic lattice. The classical spin degrees of freedom $\vec{S}$ are unit length vectors with $N_{s}=3$ components. We will focus on the frustrated case $J_{2}>0$ (AF). Without loss of generality we also take $J_{1},J_{z}<0$ (FM). Due to the bipartite structure of the lattice, our results are equally valid for the corresponding AF cases (up to a corresponding $\pi$ shift in the $\vec{Q}$ vector), obtained by simply reversing the spins on one sublattice ($J_{1}\to-J_{1}$), or by reversing the spins on adjacent layers ($J_{z}\to-J_{z}$).

The idea that a spin model with continuous symmetry has a separate Ising-nematic phase was first predicted by Henley Henley (1989) as an example of the order from disorder scenario Villain, J. et al. (1980). This was extended to the Heisenberg $J_{1}$-$J_{2}$–model by Chandra,Coleman and LarkinChandra et al. (1990) and can be explained in the following way. In the large $J_{2}$ limit the spins on each sublattice are strongly coupled. This causes them to align antiferromagnetically on each of the two interpenetrating sublattices. The effective field on a spin mediated through the nearest-neighbor couplings $J_{1}$ is thus zero. So the sublattices are effectively decoupled resulting in a zero energy cost to rotate all the (anti)aligned spins on one sublattice relative to those on the other sublattice. However, the entropy of thermal spin fluctuations depends on this relative orientation, and selects collinear spins. This can be achieved either by forming stripes of spins along the coordinate x-axis or along the y-axis. This axial orientation of spin correlations along one of the two crystal directions is essentially a discrete Ising degree of freedom that can order at a finite temperature even in a two dimensional system where long-ranged magnetic order is prohibited by the Mermin-Wagner theorem Mermin and Wagner (1966). We will refer to this ordering as lattice Ising-nematic order with order parameter denoted by $I$. The presence of long-range magnetic order will be indicated by the order parameter $M$.

In the layered $J_{1}$-$J_{2}$–model it is believed that the order-from-disorder scenario still holds for weakly coupled layers, so that there is a region of Ising-nematic order that exists before the long-range stripe order sets in at low enough temperatures, and for strongly coupled layers the nematicity and stripe order occur simultaneously in a first order transition Fang et al. (2008); Xu et al. (2008). While these prior expectations are based on an effective description that utilizes the biquadratic coupling, here we confirm this scenario for the layered $J_{1}$-$J_{2}$–model directly and provide its phase diagram for the first time.

III Method

In order to make the $J_{1}$-$J_{2}$–model tractable we write the Hamiltonian in $\vec{q}$-space

[TABLE]

where the sum goes over the first Brillouin zone and

[TABLE]

We set $J_{1}=-1$ which defines our unit of energy. For convenience we have subtracted a parameter dependent constant $C$ so that the energy of the minimum of $J_{\vec{q}}$ is zero. The momentum vectors giving these minimal energies are $\vec{Q}=(0,0,0)$ for $J_{2}<1/2$, and $\vec{Q}=(\pm\pi,0,0),(0,\pm\pi,0)$ for $J_{2}>1/2$. Thus, the ground state is a FM for $J_{2}<1/2$, and stripe-ordered (with broken lattice rotation symmetry) for $J_{2}>1/2$.

The spins on all sites are unit length vectors: $\vec{S}_{\vec{r}}\cdot\vec{S}_{\vec{r}}=1$. These constraints are enforced in the partition function by writing them as integral representations of $\delta$-functions

[TABLE]

where we have scaled the integration variable by the inverse temperature, $\beta=1/T$. This gives the partition function

[TABLE]

where we have introduced a matrix $\mathbf{\Lambda}_{\vec{q},\vec{q}^{\,\prime}}=-i\lambda_{\vec{q}-\vec{q}^{\,\prime}}(1-\delta_{\vec{q},\vec{q}^{\,\prime}})$, where $\lambda_{\vec{q}}$ is the Fourier-transformed constraint integration variable. We have separated out its $\vec{q}=0$ component and written it as $\Delta=i\lambda_{\vec{q}=0}$ and put it into the diagonal momentum space matrix $\mathbf{K}_{\vec{q},\vec{q}^{\,\prime}}\equiv K_{\vec{q}}\;\delta_{\vec{q},\vec{q}^{\,\prime}}$, where $K_{\vec{q}}\equiv J_{\vec{q}}+\Delta$. Factors of $\beta$ have been absorbed into the integration measures.

The enforcement of the unit length constraints as $\delta$-functions allows us to treat the integrals over the spin components as independent Gaussian integrals. We generalize the spins to have $N_{s}$ vector components, but will set $N_{s}=3$ at the end of the calculation. We then scale the spin components by a factor $1/\sqrt{\beta}$ and perform the Gaussian integrals to get

[TABLE]

where we have redefined the integration measure with appropriate factors of $\beta$ and

[TABLE]

The thermal average of any spin correlation function can be obtained by adding sources to the action and performing the appropriate derivatives. This yields the momentum dependent susceptibility

[TABLE]

where the brackets denote the average with respect to $S_{\lambda}$.

In order to calculate this average we first expand the action $S_{\lambda}$ and the integrand $\left(\mathbf{K}-\mathbf{\Lambda}\right)^{-1}$ in powers of $\lambda$ and treat everything except the quadratic term as perturbations. Expanding the action gives rise to ring diagrams with $n\geq 3$ wavy lines, one for each $\lambda$-factor, separated by solid lines that each contributes a factor $K^{-1}_{\vec{p}}$, see Fig. 2(a) for the $n=3$ diagram. The expansion of the integrand gives a sum of diagrams each with two external lines that carry a momentum $\vec{q}$ and $m\geq 0$ wavy lines, see Fig. 2(b).

The diagram rules are: The (bare) spin propagator is drawn as a solid line and gives a factor $K^{-1}_{\vec{q}}$. A wavy line, which we will refer to as the (bare) constraint propagator, gives a factor

[TABLE]

which originates from the quadratic part of the action $S_{\lambda}$. The $\vec{q}=0$ component of $D_{0\vec{q}}$ is explicitly set to 0. Every line, solid or wavy, carries a momentum. External lines have a fixed momentum, and there is momentum conservation at each vertex. Undetermined momenta are integrated over. The numerical factors associated to a diagram are: a factor $-i$ for each vertex in a diagram, a factor $N_{s}/2m$ for each ring with $m$ wavy lines, and an overall combinatorial factor $1/(k_{3}!k_{4}!k_{5}!\ldots)$ where $k_{m}$ is the number of rings with $m$ wavy lines.

Performing the average over $\lambda$ amounts to writing down all diagrams and connecting the wavy lines. As usual, only connected diagrams with two external legs having the same momentum $\vec{q}$ will contribute to the momentum dependent susceptibility. We approximate the averaging over $\lambda$ by summing just a selected infinite subset of diagrams in the following way. First we define a proper self-energy $\Sigma_{\vec{q}}$ which renormalizes the spin propagator according to the Dyson equation shown in Fig. 3. The renormalized spin propagator is drawn as a bold solid line.

Solving this equation gives the renormalized spin propagator

[TABLE]

Similarly a proper polarization $\Pi_{\vec{q}}$ is introduced

to make a Dyson equation for the renormalized constraint propagator (bold wavy line), see Fig. 4. Solving this gives

[TABLE]

Next we approximate the self-energy and the polarization by the self-consistent diagrams in Fig. 5

which is equivalent to writing

[TABLE]

The expression for the proper polarization Eq. (13) is finally converted, using Eqs. (9) and (11), into an equation for the renormalized constraint propagator

[TABLE]

Equations (10), (12) and (14) represent the averaging over the non-zero momentum modes of the constraint field and define a system of self-consistent equations that can be solved iteratively to give $K^{-1}_{\mathrm{eff}}{}_{\vec{q}}$ as a function of $\Delta$. After averaging over these modes the expression for the spin susceptibility, Eq. (8), becomes

[TABLE]

where the brackets denote the remaining average over the zero momentum mode(homogeneous component) $\Delta$ of the constraint field taken with respect to the weight $e^{-S_{\Delta}}=\int D\lambda e^{-S_{\lambda}}$. This averaging is carried out by simply replacing it with a single value of $\Delta$ which for self-consistency is the value which ensures that the unit vector constraint is satisfied as an average: $\langle\vec{S}_{\vec{r}}\cdot\vec{S}_{\vec{r}}\rangle=1$ which is equivalent to $\frac{1}{V}\sum_{\vec{q}}\chi_{\vec{q}}=1$. Thus the value of $\Delta$ (contained in $K^{-1}_{\mathrm{eff}}{}$) is chosen so that it satisfies

[TABLE]

and the spin susceptibility becomes

[TABLE]

using the particular value of $\Delta$ that satisfies Eq. (16).

Instead of fixing $T$ from the outset and seeking a value of $\Delta$ that satisfies Eq. (16), we will rewrite Eq. (16) as a way to calculate the temperature given a fixed value of $\Delta$:

[TABLE]

and solve the self-consistent equations keeping the value of $\Delta$ fixed. That is, we introduce an extra “mass renormalization” step where $\Delta$ is restored to its original value after each iteration.

Thus our procedure for solving the equations is as follows. First pick a value of $\Delta$ and an initial guess for $\Sigma_{\vec{q}}$. Set $N_{s}=3$. Then iterate the following steps:

Subtract a constant from $\Sigma_{\vec{q}}$ so that its minimum value becomes zero. 2. 2.

Make $K^{-1}_{\mathrm{eff}}{}_{\vec{q}}$ according to Eq. (10), and calculate $T$ according to Eq. (18). If $T$ has converged then exit and use the obtained $K^{-1}_{\mathrm{eff}}{}_{\vec{q}}$ to compute $\chi_{\vec{q}}$, Eq. (17). 3. 3.

Calculate $D_{\vec{q}}$ using the new $K^{-1}_{\mathrm{eff}}{}_{\vec{q}}$, Eq. (14). 4. 4.

Calculate the new self-energy $\Sigma_{\vec{q}}$, Eq. (12). 5. 5.

Go to step 1.

This whole procedure is repeated for a range of $\Delta$ values, typically two hundred values ranging from $10^{-9}$ to $1$ evenly spaced on a log-scale. We use the convergence criterion that $T$ has converged when the relative difference between the temperatures at succesive iterations, $|T_{n+1}-T_{n}|/T_{n}<10^{-8}$. Convergence is typically reached after 10-20 iterations.

For each value of $\Delta$ the iteration procedure yields a value of $T$ and a corresponding susceptibility $\chi_{\vec{q}}$. While $T$ is usually a single-valued function of $\Delta$, there are temperature regions near phase transitions where different values of $\Delta$ correspond to the same value of $T$ but to different values of $\chi_{\vec{q}}$. Thus in these regions the equations give several different solutions $\chi_{q}$ for the same value of $T$. We will deal with this by selecting the solution with the lowest free energy, which we take to be the one smoothly connected to the unique solution on the low-$T$ side of the multi-valued region.

IV Order parameters

The calculation of the momentum-dependent susceptibility $\chi_{\vec{q}}$ at different temperatures allows us to make inferences about different phases and order parameters. We focus on the planar nematic order parameter

[TABLE]

where $\hat{x},\hat{y}$ are planar unit lattice vectors and $\langle...\rangle$ denotes the thermal expectation value. In terms of the momentum space susceptibility, the order parameter is

[TABLE]

where the sum over $\vec{q}$ is taken over the first Brillouin zone. The planar nematic order parameter detects anisotropy in bond ordering on the planes, and does not break spin rotational symmetry.

We also calculate the magnetization order parameter which breaks spin rotation symmetry. We infer this from the coefficient of the diverging susceptibility as the system size goes to infinity. The $\vec{q}$-dependent susceptibility at the ordering vector is

[TABLE]

where $\delta f(\vec{r})$ is the spin fluctuation correlation function characterized by a correlation length $\xi$. For finite $\xi$ less than the linear system size, the last term will be independent of $V$. Therefore for large system sizes we can keep only the first term, which diverges with increasing $V$, and arrive at

[TABLE]

In the last equality we have used the fact that the maximum value $K^{-1}_{\mathrm{eff}}{}_{\vec{Q}}$ is always $1/\Delta$ because both the self-energy $\Sigma_{\vec{Q}}$ and $J_{\vec{Q}}$ is zero for the maximum value of $K^{-1}_{\mathrm{eff}}{}$. Extracting the magnetization this way gives the dominating magnetic order corresponding to the wave vector $\vec{Q}$ where $K^{-1}_{\mathrm{eff}}{}$ is maximal.

V Results

For simplicitly, we begin by studying the $J_{2}=0$ Hamiltonian where our results can be compared to Monte Carlo data. In this case the Hamiltonian reduces to an unfrustrated layered ferromagnet (FM) where we expect a finite temperature phase transition to an ordered state with a FM magnetic moment at wavevector $\vec{Q}=(0,0,0)$. We first consider the isotropic FM $(J_{z}=-1)$ and solve the self-consistent equations numerically for a range of $\Delta$’s and obtain results for the magnetization using Eq. (22). We always start the iterations of the self-consistent equations with an initial guess for $\Sigma_{\vec{q}}$ which breaks the nematic symmetry so that the initial $I$ is slightly negative.

Figure 6(a) shows the magnetization squared as a function of temperature for cubic systems of different linear sizes $L=100-400$. The finite-size curves fall almost on top of each other and do not cross. On magnifying the behavior close to $T_{c}$ where the magnetization vanishes, finite-size effects become apparent, see Fig. 6(b), and as $T$ increases there is a slight overshoot of the magnetization curves for the largest system sizes. These magnetization curves are clearly unphysical as the overshoot leads to a temperature region near $T_{c}$ where the magnetization is a multivalued function of $T$. We interpret this behavior as the system is having multiple possible solutions to the self-consistent equations with different free energies at the same temperature. The free energy increases with temperature, thus in order to choose the solution with lowest free energy we pick the branch in the multivalued region that is connected smoothly to the low temperature single-valued region. When following this branch upon increasing $T$, the curve will turn around at some temperature and will, if continued, give other solutions with higher free energies. We omit these by drawing a vertical line towards zero order parameter at the first turning point (infinite slope of $M^{2}(T)$) upon increasing $T$. This is drawn as a red line for the largest system size in Fig.6(b). Then the order parameter evolution with $T$ continues from where the line hits the order parameter curve again. This results in a discontinuous jump in the order parameter at the critical temperature which we interpret as a discontinuous phase transition. However, in the case of Fig. 6(b) the discontinuous jump decreases as the system size is increased, which when extrapolated to infinite size results in a continuous transition, as is expected for a 3D FM.

In order to determine $T_{c}$ of the magnetic phase transition we use different methods based on how the finite-size curves behave. If they cross, we pick the crossing-points between successive linear system sizes $L$ and $L+50$ and extrapolate these crossings to infinite size using a cubic polynomial in $1/L$. For finite-size curves that do not cross but has an overshoot we identify the temperature of each finite-size curve at the point where the magnetization curve turns back, and then extrapolate these points to the infinite $L$ limit using a cubic polynomial in $1/L$. A third method we use, which also works when the magnetization curve is single-valued, is to use the expected behaviour of the magnetization near a continuous phase transition $M=A(T_{c}-T)^{\beta}$. To do so we make a Kouvel-Fisher plotKouvel and Fisher (1964), i.e. we plot $\left[-\frac{d\ln M}{dT}\right]^{-1}$ vs. $T$ and find for each finite-size curve the temperature where this quantity crosses the temperature axis. These are finally extrapolated to the infinite-size limit using a cubic polynomial in $1/L$. For $J_{z}=-1.0$ and $J_{2}=0$ these methods give results that are very close to each other: $T_{c}=1.52083$ (disc.) and $T_{c}=1.52077$ (Kouvel-Fisher). This is to be compared with the most accurate Monte Carlo resultChen et al. (1993) which is $T_{c}=1.442928(77)$, a relative difference of $\sim 5\%$. We attribute this difference to the lack of vertex corrections in our self-consistent equations. Thus we estimate that our critical temperatures are approximate with a relative accuracy of roughly $5\%$. For very weakly coupled layers where the system is almost two dimensional we use a fourth method where the magnetic $T_{c}$ is determined as the temperature at which the spin-spin correlation length diverges, see Sec. VI.

We now turn to the weakly frustrated FM regime $0<J_{2}<0.5$ with isotropic interlayer couplings $J_{z}=-1$. We find that for increasing $J_{2}$, $T_{c}$ goes down, and the magnetization overshoot seen for finite sizes at $J_{2}=0$ quickly becomes smaller. For $J_{2}>0.5$ there is no longer a low temperature finite FM magnetization. Instead the dominating divergence of the susceptibility is at $\vec{Q}=(\pi,0,0)$ or $\vec{Q}=(0,\pi,0)$ which indicates stripe magnetic order which we will denote by $M$. In addition there is also a finite value of the nematic order parameter $I$, see Fig. 7(a). In this figure (and in following figures), the stripe magnetization squared is shown as positive values, while the nematic order parameter is shown as negative values. Finite size effects are small for $J_{2}=0.51$, and those present indicate that the jump in the order parameters increases slightly with system size. Thus we conclude that for $J_{2}=0.51$ there are simultaneous discontinuous phase transitions in both the nematic order parameter $I$ and the stripe magnetization $M$, which we write in short-form Id/Md, where the d means discontinuous. On further increasing $J_{2}$ this simultaneous Id/Md character of the transition persists up to the largest $J_{2}$ value studied, $J_{2}=2$, see Fig. 7(b). The quantitative changes upon increasing $J_{2}$ include: larger finite size effects, increasing $T_{c}$, and smaller jumps in the order parameters.

For $J_{2}=0.495$, i.e slightly less than the boundary between the FM and the stripe phase, the peak magnetic ordering wave vector changes with temperature from being FM at high temperatures, to stripe order, and then back to FM again at the lowest temperatures. The intermediate stripe order is also indicated by the finite temperature region with nonzero $I$ in Fig. 8. Thus, at finite temperatures the nematic phase extends slightly into the region $J_{2}<0.5$.

Summarizing these results for $J_{z}=-1$ gives the phase diagram shown in Fig. 9(a) where the black curve shows the continuous FM phase transition while the red curve marks the discontinuous simultaneous nematic and stripe phase transitions. The blue curve which is slightly bent towards the FM phase defines the boundary between the ordered FM and stripe phases. The phase diagram for $J_{z}=-0.1$ similarly obtained is shown in Fig. 9(b).

For weaker coupling $J_{z}$ between the layers the phase diagram is richer. In particular, the nematic and stripe phase transitions cease to occur simultaneously, and there is a finite temperature window where nematic long-range order exists without magnetic stripe order. This last fact can be seen from Fig. 10(a) at $J_{z}=-0.0001$ and $J_{2}=0.8$. There the magnetization curves $M$ for different system sizes cross each other and the nematic order parameter $I$ becomes non-zero, with no overshoot, at a higher temperature than the magnetization crossings occur. The $T_{c}$ of the stripe magnetic order is determined by extrapolating the crossings of successive finite-size magnetization curves using a cubic polynomial in $1/L$. These curves are single-valued so we conclude that these transitions are continuous. In what follows we indicate these split phase transitions by writing (Ic,Mc), where the phase transition with the highest $T_{c}$ is written first. This intermediate temperature region with nematic but no stripe order exists also for smaller $J_{2}$ almost all the way to the FM phase as can be seen from Fig. 10(b) for $J_{2}=0.513$.

Very close to the FM phase, $J_{2}\approx 0.506$, the nematic order parameter develops an additional kink feature at a finite value of $I$ which rapidly becomes an overshoot, Figs. 10(c)-(d). For these values of $J_{2}$ the $M$ curves also overshoot at the same $T$ as the kink feature in $I$. However, their magnitudes decrease with increasing system size. In each of the panels Fig. 10(c)-(e) there is a temperature region close to $T_{c}$ of the nematic phase transition where the procedure of iterating the self-consistent equations converges very slowly. In fact there is a small temperature region (typically a little less than 10% of $T_{c}$), indicated by the orange thick line, where we are unable to find solutions to the self-consistent equations for the largest system size. For $J_{2}$ even smaller, the branch of $I$ which indicates the region of nematic order without magnetic order quickly moves down in temperature, and concomittantly the crossings of the magnetizations move up in temperature and end up as overshoots which again indicate an Id/Md transition, see Figs. 10(e)-(f). Thus, close to $J_{2}^{*}\approx 0.505$ there is a rapid change from a split regime with two continuous phase transitions to a regime with simultaneous discontinuous phase transitions. The corresponding phase diagram for $J_{z}=-0.0001$ is shown in Fig. 11.

For even weaker interplane coupling $|J_{z}|<0.0001$ the nematic phase boundary (orange curve in Fig. 11) stays almost unchanged as it becomes equal to the nematic phase transition boundary for a 2D frustrated Heisenberg magnet shown in Fig. 12. Note that the nematic phase extends slightly into the region $J_{2}<0.5$ at finite temperatures also for the strictly 2D case. This is contrary to Ref.Weber et al., 2003 where Monte Carlo data show evidence of an infinite slope at $J_{2}=0.5$.

The other phase boundaries move to lower temperatures, compliant with the expectation that a 2D system cannot sustain long-range magnetic order in accordance with the Mermin-Wagner theoremMermin and Wagner (1966). A plot of the magnetic $T_{c}$ for $J_{2}=1$ as a function of $|J_{z}|$ is shown in Fig. 13, and shows a logartihmic behavior $T_{c}\approx 0.722+0.0149\ln{|J_{z}|}$ which we interpret to be the first terms in a series expansion of $T_{c}=a/(1+b\ln{(1/|J_{z}|)})$, with $a$ and $b$ independent of $J_{z}$, which is the renormalization group expectation for $N_{s}=3$ spatially anisotropic nonlinear $\sigma$-models with a small microscopic in-plane couplingAffleck and Halperin (1996).

It is interesting to explore for what values of $J_{2}$ the split regime of separate nematic and magnetic stripe phase transitions occurs as the interplane coupling is lowered from $|J_{z}|=1$. Already at $J_{z}=-0.01$ there is a hint of the split regime near $J_{2}=0.5$: Fig. 14(a) shows simultaneous Id/Md transitions at $J_{2}=0.505$, and at $J_{2}=0.52$ in Fig. 14(b), the nematic order parameter becomes non-zero at a higher temperature than the Id/Md occurs. Although there is a region of slow convergence in reaching the solution to the self-consistent equations above this nematic phase transition, it appears to be continuous as there is no visible overshoot. We thus have split phase transitions; an Ic occuring at a higher $T_{c}$ than the discontinuous stripe magnetic phase transition, that occurs simultaneously with another discontinuous metanematic transition in the nematic order parameter between two finite values, a (Ic,Ifd/Md) where the letters fd signifies that the transition is a discontinuous jump between two finite values of the nematic order parameter. For $J_{2}=0.55$, Fig. 14(c), the Ic feature has dropped back below the temperature of the Ifd/Md transition, resulting in a single simultaneous Id/Md transition. So for $J_{z}=-0.01$ the split transitions occur only in a narrow parameter region $J_{2}^{*}\approx 0.513<J_{2}<J_{2}^{**}\approx 0.534$. The differences in critical temperature in the split region for $J_{z}=-0.01$ are also small; for $J_{2}=0.52$ we find $\Delta T_{c}/T_{c}\sim 0.5\%$.

A wider region of split behavior exists when $|J_{z}|$ is further lowered to $J_{z}=-0.002$. What is especially interesting is how the nature of the split phase transitions appears to change as $J_{2}$ is varied. For $J_{z}=-0.002$ this is illustrated in Fig. 15 and goes as follows: A simultaneous Id/Md exists up to $J_{2}^{*}\approx 0.505$, Fig. 15(a). For $J_{2}=0.51$, Fig. 15(b), the nematic order parameter experiences a Ic before there is a Md concomittant with a discontinuous jump in the nematic order parameter between two finite values, i.e. a (Ic,Ifd/Md). The discontinuous character of this lower temperature metanematic Ifd/Md transition weakens rapidly as $J_{2}$ is increased, see Fig. 15(c). At $J_{2}=0.55$ the magnetic transition appears continuous (Mc), and the jump in the nematic order parameter at the magnetic transition is changed into a weak kink-like feature signalling a regime of two distinct continuous phase transitions with separate $T_{c}$ values,(Ic,Mc), similar to that shown in Fig. 10(a). The relative difference in the $T_{c}$ values is largest for the smallest values of $J_{2}$ after the splitting occurs and changes only slightly for an extended range of $J_{2}$. However, for $J_{2}$ large enough the $T_{c}$ values approach each other again and the nature of the phase transitions changes. At $J_{2}=0.9$, Fig. 15(d), the nematic order parameter bends slightly back for the biggest system sizes, thus the highest temperature transition becomes Id, while the lower temperature phase transition still appears to be continuous Mc; (Id,Mc). For even larger $J_{2}$ the two phase transitions come even closer in temperature, and now also the magnetic phase transition becomes discontinuous, so that at $J_{2}=1.1$, Fig. 15(e) there are two slightly separated discontinuous phase transitions; a (Id,Ifd/Md). At even bigger values of $J_{2}$, Fig. 15(f), the upper discontinuous transition Id moves below the simultaneous transition, resulting in a single simultaneous Id/Md again. The boundary value, for $J_{z}=-0.002$, where the split transitions merge again is $J_{2}^{**}=1.13$. The value of $J_{2}^{**}$ increases rapidly as $|J_{z}|$ is further lowered, and exceeds the largest value considered here ($J_{2}=2$) for $J_{z}=-0.0001$. A plot of $J_{2}^{*}$ and $J_{2}^{**}$ vs. $|J_{z}|$ is shown in Fig. 16.

VI Spin correlations

For $J_{2}>0.5$, the Fourier transform of the spin-spin correlation function $\chi_{\vec{r}}\equiv\langle\vec{S}_{\vec{r}}\cdot\vec{S}_{0}\rangle$ is dominated by peaks around $(\pm\pi,0,0)$ and $(0,\pm\pi,0)$. These peaks generally have different shapes along $q_{x}$, $q_{y}$ and $q_{z}$ directions, and in the nematic phase the peaks related by a $\pi/2$ lattice rotation about the $z-$axis also have different heights, see Fig. 17, while peaks related by a $\pi$ rotation about the $z-$axis have equal heights.

Approximating these peaks by Lorentzians, we can write the real-space spin-spin correlation function as

[TABLE]

where the lattice spacing has been set to unity. Unprimed(primed) symbols refer to the $(\pm\pi,0,0)$($(0,\pm\pi,0)$) peaks, and $\xi_{i}^{-1}$ is the half width half maximum of the peak in the direction along the real space stripes ($i=1$), in-plane perpendicular to the stripes ($i=2$) and perpendicular to the stripes in the z-direction ($i=z$). Carrying out the summations one obtains

[TABLE]

where $\vec{r}^{\,\prime}=\left(\frac{x}{\xi_{2}},\frac{y}{\xi_{1}},\frac{z}{\xi_{z}}\right)$, $\vec{r}^{\,\prime\prime}=\left(\frac{x}{\xi_{2}^{\prime}},\frac{y}{\xi_{1}^{\prime}},\frac{z}{\xi_{z}^{\prime}}\right)$ and $f(r)=\frac{2\pi^{2}}{r}e^{-r}$.

For temperatures above the nematic phase transition the peaks are related by a $\pi/2$ rotation about the $z-$axis, i.e. $\xi_{1}=\xi_{1}^{\prime}$, $\xi_{2}=\xi_{2}^{\prime}$, $\xi_{z}=\xi_{z}^{\prime}$ and $A=A^{\prime}$. This means that the in-plane correlation function along one of the axes, here the $x-$axis, is

[TABLE]

The spin correlation function, Eq. (25), is thus governed by two components: an oscillating component that decays with correlation length $\xi_{2}$ and a uniform component that decays with correlation length $\xi_{1}$. The difference between $\xi_{1}$ and $\xi_{2}$ measures the ellipticity of the peaks of the susceptibility in momentum space. For weak interlayer couplings the correlation length $\xi_{z}$ is much smaller than $\xi_{1}$ and $\xi_{2}$.

For temperatures below the nematic phase transition, one set of peaks will become higher and narrower than the other set. In the case of a negative nematic order parameter, the peaks at $(\pm\pi,0,0)$ will dominate, leading to $A\gg A^{\prime}$ and $\xi_{1}\gg\xi_{1}^{\prime}$ etc. Therefore, in the nematic phase with negative nematic order parameter the spin-spin correlation function along the stripes ($y-$direction) is given at long distances by

[TABLE]

Similarly the in-plane correlation function perpendicular to the stripes becomes

[TABLE]

For our finite lattice system with periodic boundary conditions we extract the width of the peak in the $i$-direction by fitting the values of $\chi_{\vec{q}}$ along a line in direction $i$ in $\vec{q}$-space through the point $\vec{Q}$ using the functional form

[TABLE]

where $A$ is independent of $\vec{q}$ and the cosine takes care of the $q$-space periodicity. Here we have $\vec{Q}=(\pi,0,0)$. The inverse correlation lengths in the three directions so obtained are shown as a function of $T$ for $J_{2}=1$ and $J_{z}=-0.0001$ in Fig. 18. We note that the correlation lengths are in general different, also above the nematic phase transition. Due to the very weak coupling $J_{z}$ between the layers, the correlation length in the $z$-direction is correspondingly small (note that we have plotted $0.01\xi_{z}^{-1}$ in Fig. 18). Especially noteworthy is the fact that the correlation lengths increase much more rapidly with lowering temperature below the nematic phase transition than above. In fact this difference in behavior can be taken as an observational signature of the nematic phase transition alone.

To derive the relation between the magnetic correlation length and the Ising-nematic order parameter near $T_{c}$, we expand the self-energy in Eqs. (10), (17) to linear order in the order parameter $I$. Next we expand near the peaks of $\chi_{\vec{q}}$ as in Eq. (VI), and substitute this expression into the equation for the order parameter, Eq. (19). Since $\chi_{\vec{q}}$ is peaked near $\vec{Q}$ and its symmetry-related points, we may perform a calculation similar to the one yielding Eq. (25) to write the order parameter as

[TABLE]

where the first[second] sum over $\vec{q}$ is restricted to the region around the peak at $(0,\pi,0)[(\pi,0,0)]$, and $J_{{\rm eff}\,\vec{q}}=J_{\vec{q}}-\Sigma_{\vec{q}}^{\prime}$ may be expanded near its minima (as in Eq. (VI)), $\Sigma_{\vec{q}}^{\prime}$ is the point-group symmetric part of the self-energy and $\alpha$ is a parameter determined below. $\Delta_{0}$ is the value of $\Delta$ at $T_{c}$. Expanding Eq. (VI) to linear order in $I$ on its right side, we obtain the relation

[TABLE]

Upon inspection of Eqs. (VI), (VI) one sees that the susceptibility peak heights and widths may be expanded in powers of the order parameter $I$. Keeping only the linear term (valid near $T_{c}$ where $I\ll 1$), one finds that

[TABLE]

where the 0 subscript indicates the value at $T=T_{c}$ (i.e. $I=0$). As a result, both the peak height ($A$) and width $(\xi^{-1})$ exhibit a nonanalytic temperature variation across the nematic phase transition, in an amount directly proportional to $I$. This scaling behavior for the correlation lengths and the amplitudes is shown in Fig. 19.

To investigate the diverging behavior of the correlation length associated to the magnetic phase transition which occurs below the nematic phase transition we have made a Kouvel-Fisher analysis in Fig. 20 showing $K(T)\equiv-[\frac{\partial\ln{\xi_{1}}}{\partial T}]^{-1}$ vs. $T$ for the same parameters as used in Fig. 18. Taking into account also the leading irrelevant operator with a scaling dimension $y_{1}<0$ we expect that close to $T_{c}$ the Kouvel-Fisher function behaves as $K(T)=\frac{1}{\nu}\left(T-T_{c}\right)+B\left(T-T_{c}\right)^{y+1}$ where $y=-y_{1}\nu>0$ and $B$ is a constant. In Fig. 20 we have plotted $K(T)$ for different system sizes. We see that finite-size effects are visible as low-temperature upturns for all system sizes, but that the infinite size behavior can be inferred by extrapolating the largest system size results for temperatures above the finite size upturn. We find that a value $y=1$ gives a good fit to the behavior of the largest system. Fixing $y=1$ we find a best fit, shown as the red curve in Fig. 20, of $K(T)=\frac{1}{0.678}(T-T_{c})-29.316(T-T_{c})^{2}$ with $T_{c}=0.5854$. The value of $\nu=0.678$ so obtained is reasonably close to the value $\nu=0.7112$ for the 3D Heisenberg universality class Campostrini et al. (2002).

VII Discussion

Our results show that the phase diagram of the layered $J_{1}$-$J_{2}$-model is remarkably rich. In particular, in the frustrated regime, $J_{2}>0.5$, it has separate nematic and magnetic phase transitions as the temperature is lowered, if the interplane coupling $|J_{z}|$ is small enough. For weak interplane couplings, $|J_{z}|\lesssim 0.0001$, both of these phase transitions are continuous for most values of $J_{2}(>0.5)$. At $|J_{z}|=0.0001$ the relative difference in critical temperatures of these transitions is rather small, $\sim 7\%$, but increases as $|J_{z}|$ is further reduced. This is because the critical temperature of the nematic phase transition stays almost unchanged below $|J_{z}|=0.0001$, while the critical temperature of the magnetic phase transition goes to zero.

For values of $J_{2}$ very close to $0.5$, the nature of the phase transitions change, and eventually both become simultaneous first order phase transitions for $J_{2}\to 0.5$. Exactly how this change happens is complicated as evidenced by Fig. 10. In particular it involves a metanematic phase transition where the already finite nematic order parameter exhibits a jump to another value when the magnetic phase transition occurs, Fig. 10(d).

In contrast, the regime of split transitions does not exist for strong interplane couplings $|J_{z}|>0.01$. There the phase transition into the magnetic phase from the disordered side is a discontinuous transition where the nematic and stripe magnetic order sets in simultaneously for all values of $J_{2}>0.5$.

For intermediate interplane couplings, $0.0001<|J_{z}|<0.01$, the $J_{2}$ regime of split phase transitions opens up from a value $J_{2}^{*}\approx 0.5$ which is very weakly dependent on $J_{z}$ and continues up to a value $J_{2}^{**}$ where it closes and reverts to simultaneous first order transitions. The value $J_{2}^{**}$ is very close to $J_{2}^{*}$ for $|J_{z}|=0.01$ and increases rapidly as $|J_{z}|$ is lowered.

Our results show also that the spin-spin correlation length, or more accurately its temperature derivative, can serve as a probe to detect the nematic phase transition. This is so because the correlation length depends directly on the nematic order parameter, thus it exhibits a sharp increase with decreasing temperature exactly at the nematic phase transition. For split transitions this is separate from the normal critical divergence of the correlation length that happens at a lower temperature where the stripe magnetic phase transition occurs.

Our results are strictly only valid for the classical layered $J_{1}$-$J_{2}$ model, so we can only speculate on how our results carry over to the corresponding quantum model. What seems plausible on general grounds is that quantum effects will enhance fluctuations and be most dramatic in the regions of the phase diagram where there are adjacent phases and the critical temperatures are the lowest, especially the region close to $J_{2}=0.5$. Evidences of a gapless spin liquid phase in the two dimensional quantum spin-1/2 $J_{1}$-$J_{2}$-model near $J_{2}=0.5$ was recently reported in Ref. Wang and Sandvik, 2018 but this region is notoriusly difficult to address. For higher spins with less quantum fluctuations we do expect that the high temperature features discussed here, such as the splitting of the nematic and the stripe magnetic phase transitions, carry over to the quantum case.

When it comes to the applicability of our results to the iron-based superconductors, it will depend on how well iron-based superconductors are based on models of localized magnetic moments. We do not address this question. Nevertheless, our results show that virtually all of the rich phenomena predicted to occur in effective classical models of itinerant magnetic momentsFernandes et al. (2012) appear also in the classical $J_{1}$-$J_{2}$–model for weak interplane couplings. Split continuous phase transitions is a generic feature for weak interplane couplings. More exotic phase transitions, such as split metanematic transitions, were found only in narrow regions of the phase diagram.

In addition to obtaining the results for the $J_{1}$-$J_{2}$–model we have in this paper also taken the opportunity to outline the details of the nematic bond theory introduced in Ref. Schecter et al., 2017. As demonstrated here, this is an efficient and versatile method for dealing with frustrated magnetism, even for very large system sizes, and can serve as a supplement to computationally demanding Monte Carlo techniques.

Acknowledgements.

The Center for Quantum Devices is funded by the Danish National Research Foundation. We acknowledge support from the Laboratory for Physical Sciences, and Microsoft (M.S.). The computations were performed on resources provided by UNINETT Sigma2 - the National Infrastructure for High Performance Computing and Data Storage in Norway.