# Nematic phase in the J-J square lattice Ising model in an external field

###### Abstract

The J-J Ising model in the square lattice in the presence of an external field is studied by two approaches: the Cluster Variation Method (CVM) and Monte Carlo simulations. The use of the CVM in the square approximation leads to the presence of a new equilibrium phase, not previously reported for this model: an Ising-nematic phase, which shows orientational order but not positional order, between the known stripes and disordered phases. Suitable order parameters are defined and the phase diagram of the model is obtained. Monte Carlo simulations are in qualitative agreement with the CVM results, giving support to the presence of the new Ising-nematic phase. Phase diagrams in the temperature-external field plane are obtained for selected values of the parameter which measures the relative strength of the competing interactions. From the CVM in the square approximation we obtain a line of second order transitions between the disordered and nematic phases, while the nematic-stripes phase transitions are found to be of first order. The Monte Carlo results suggest a line of second order nematic-disordered phase transitions in agreement with the CVM results. Regarding the stripes-nematic transitions, the present Monte Carlo results are not precise enough to reach definite conclusions about the nature of the transitions.

###### pacs:

05.50.+q,64.60.Cn,75.40.-s## I Introduction

Competing interactions are common in many natural and artificial systems, like the presence of conflicting ferromagnetic and antiferromagnetic interactions in frustrated magnetic systems as spin glasses Fischer and Hertz (1991) and ultrathin magnetic films Abanov et al. (1995); De’Bell et al. (2000); Vaterlaus et al. (2000), as well as competition between an attractive and a repulsive part in the interaction between atoms and molecules of complex fluids Imperio and Reatto (2006); Archer and Wilding (2007); Olson Reichhardt et al. (2010); Ciach (2011); Pȩkalski et al. (2013, 2014); Almarza et al. (2014). Competition between conflicting interactions is also relevant in mathematical optimization problems, when decisions have to be made where not all the constraints can be satisfied simultaneously Mézard and Montanari (2009). Frustration, the inability of a system to statisfy all local constraints, is a unifying concept in many natural and artificial systems. Competing tendencies usually are responsible for complex behavior like slow relaxation to equilibrium, strong metastability and rough energy landscapes Wales (2003). This makes the study of such systems both very interesting and challenging. One of the characteristic outcomes of the presence of competing interactions in a system is the emergence of heterogeneous structures as the equilibrium or low energy states, like stripes, bubbles, clusters, disordered phases and anisotropic behavior.

One of the simplest models with competing interactions is the J-J Ising model on the square lattice. This model is defined as a simple extension of the square lattice Ising model, in which besides the nearest-neighbor (NN) ferromagnetic or attractive interaction J, a next-nearest-neighbor (NNN) antiferromagnetic or repulsive interaction J is added. The ground state of the model depends on the relative intensity of the competing interactions . For it is ferromagnetic and for it has a stripe structure of alternating up and down rows of spins. There is no exact solution for the thermodynamics of the model. At zero external field it has been studied by a variety of techniques like cluster mean field theory, transfer matrix and Monte Carlo simulations Morán-López et al. (1993); Moran-Lopez et al. (1994); Cirillo et al. (1999); dos Anjos et al. (2008); Kalz et al. (2009, 2011); Jin et al. (2012, 2013); Saguia et al. (2013), considering both ferromagnetic and antiferromagnetic nearest-neighbor interactions. The nature of the thermal phase transition from the stripes to a disordered phase for was controversial. In the most recent studies combining Monte Carlo simulations and a series of analytical techniques it has been established that the line of phase transitions in the temperature versus plane is first order for and is continuous with Ashkin-Teller critical behavior for . The critical exponents change continuously in this regime between the 4-state Potts model behavior at to standard Ising criticality for Jin et al. (2012, 2013).

In comparison with the zero field case, the model in an external field has received much less attention. Queiroz de Queiroz (2011) and Yin et.al. Yin and Landau (2009) studied the case with both NN and NNN interactions of the antiferromagnetic type, using transfer-matrix methods in conjunction with finite-size scaling and conformal invariance in the first reference and large scale Monte Carlo simulations in the second one. For the ground state is striped at small fields and a row-shifted phase appears for . The latter state consists of alternating ferro and antiferromagnetically ordered rows (or columns), with the ferromagnetic ones parallel to the field. In both studies it was observed a reentrance in the boundary stripes-paramagnetic upon lowering the temperature at constant field. In Ref. Yin and Landau, 2009 it was argued that the reentrant behavior may be due to the appearance of row-shifted () clusters that help to sustain striped () order at low temperatures, even for moderately large magnetic fields. The nature of the phase transitions points to a weak universality scenario with exponents departing slightly from the standard Ising values.

A natural question when dealing with stripe forming systems is the possibility of existence of an intermediate nematic phase. A nematic phase in this context is characterized by the presence of orientational order but without translational or positional order Cannas et al. (2006); Barci et al. (2013a); Almarza et al. (2014). In this sense there are broken symmetry phases but with an intermediate degree of symmetry, higher than the less symmetric stripe phases in which both orientational and positional orders are present. Nematic phases associated with intermediate stripe-like order are present in many quasi-two-dimensional systems like ultrathin ferromagnetic filmsAbanov et al. (1995); Cannas et al. (2006); Nicolao and Stariolo (2007) and electronic liquids in which they may be relevant to understand high temperature superconductivity Han et al. (2001); Berg et al. (2009); Fradkin et al. (2010). The conditions under which a system can sustain a nematic phase of this kind are still not completely clear. Strong evidence for the existence of such phases have been found in systems with isotropic competing interactions at different scales, e.g. when a short range ferromagnetic interaction competes with a long range antiferromagnetic one decaying with a power law of distance, like the dipolar interaction Barci and Stariolo (2009); Barci et al. (2013b). In this particular case a nematic phase is present but only quasi-long-range nematic order develops in 2D. This quasi-nematic phase emerges by breaking a continuous symmetry, similar to the Kosterlitz-Thouless phase transition in the 2D XY model. In these kind of models, smectic phases are suppressed at finite temperatures due to the strong fluctuations of the order parameter. When an external magnetic field is applied, competing dipolar interactions give rise to new and interesting phases. At small fields stripe phases are still present, although the direction aligned with the field is favoured energetically and it gets wider as the field is risen until an instability leads to a first order phase transition to a bubble phase at a critical external field value. At still higher fields there is a second transition from the bubble to an homogeneously magnetized phase, a saturated paramagnet. Another salient feature of the field-temperature phase diagram of the system with dipolar interactions in a field is a strong reentrant behavior. This has been observed in beautiful experiments on ultrathin films of Fe/Cu(001) Saratz et al. (2010a, b) and also in mean field approximations Portmann et al. (2010); Cannas et al. (2011); Velasque et al. (2014). As for the presence of nematic phases in an external field, this is still an open question in dipolar or electronic systems. Compared to the behavior of the dipolar frustrated systems, the J-J model in the square lattice stands at the opposite side: it has a very simple stripe phase, with long range orientational and positional order, absent in models with long range isotropic interactions. The disordered-stripe phase transition in this model corresponds to the breaking of the symmetry of the square lattice to the symmetry of the stripe phase. Then, the question we try to answer in this work is: is it possible for the J-J model in the square lattice to sustain an Ising-nematic phase? , i.e. a phase with orientational but not positional order ? To give even a partial answer to such question is always a considerable challenge. This is because the nematic order parameter in stripe forming systems amounts to compute correlation functions in different space directions searching for a breaking of isotropy characteristic of these phases Stariolo and Barci (2010); Barci and Stariolo (2011). Then, the most common analytical approaches for a one-particle order parameter, namely mean field theory, fails at detecting nematic phases and one must go beyond naive MFT to an approximation which allows to compute anisotropic correlations. In this work we have studied the J-J model both with and without an external field by means of two approaches: the Cluster Variation Method (CVM) which is a cluster mean field theory, and Monte Carlo simulations. The CVM allows for a systematic improvement upon naive MFT by considering clusters of particles of increasing size in an exact way in the partition function. To our knowledge, this is the first time the CVM is used with the specific aim of searching for anisotropic correlations, for which it is particularly well suited. In fact, the first step beyond the mean field approximation in a lattice is the two-site or pair approximation, also knwon as Bethe-Peierls approximation Bethe (1935); Tanaka (2002). This amounts to consider in a exact form all clusters with two sites. This approximation is known to predict correctly that is the lower critical dimension for the Ising ferromagnet. Nevertheless, because it does not distinguish any geometric or spatial features in the sum over pairs of sites, it is not able to capture rotation symmetry-breaking, a distinctive feature of orientational phases. The next degree of approximation in the square lattice is the plaquette or square approximation, which considers exactly clusters of four sites, i.e. first and second neighbors. We will show that this is enough to capture the presence of nematic phases in models with competing first and second neighbor interactions. We did not find evidences of nematic phases in the J-J model at zero external field within the square approximation in the CVM, but a nematic phase appears when a field is present, both in the CVM approach and in Monte Carlo simulations.

## Ii The Cluster Variation Method

We give here a very brief description of the Cluster Variation Method, focusing on quantities that will be useful in the calculation below. There is a large literature on the technical aspects of the method and the interested reader can refer to references Kikuchi (1951); Sanchez et al. (1984); An (1988); Tanaka (2002) for more comprehensive discussions of the method and its potential for computing variational approximations to the free energy of different systems.

Consider the variational free energy

(1) |

where means a trace or a sum over all the relevant degrees of freedom of the Hamiltonian , and is a trial density matrix which satisfies the normalization constraint . A systematic way for obtaining variational approximations to the free energy is to express it in terms of a cumulant expansion. Following the exposition by Tanaka Tanaka (2002), consider the n-body reduced density matrix for an body system:

(2) |

for , with the normalization conditions:

(3) |

The cluster functions and the cumulant functions are defined as follows:

(4) | |||||

(5) | |||||

(6) | |||||

and so on. The largest cluster fuction (which corresponds to the N-body entropy term in the variational free energy) can be written as a sum over all the cumulant functions in the form:

(7) | |||||

In this way the variational free energy can be written in terms of an expansion in cumulant functions:

(8) | |||||

This form of the variational free energy allows to obtain systematic approximations to the true free energy of a model system by considering a maximal size of cluster to be summed exactly in the partition function. This is called the parent cluster. This amounts to truncate the cumulant expansion to a given degree and optimizing the resultant expression with respect to the reduced density matrices which must satisfy the reducibility and normalization conditions (2) and (3) respectively. The reduced density matrices represent all the subclusters contained within the parent cluster. In some applications a convenient way of implementing the variational approximation is to parametrize the density matrices in terms of correlation functions and consider these as variational parameters instead.

In the case of a system with Ising spins , the reduced density matrices can be written as Cirillo et al. (1999):

(9) |

where the sum runs over all subclusters with sites within cluster , and the k-point correlation functions are defined by . The variational parameters must satisfy:

(10) |

A hierarchy of approximations to the free energy can be constructed in this way. The simplest one corresponds to the 1-point approximation for the density matrices, the usual mean field approximation. The 2-point approximation is usually called Bethe-Peierls approximation Bethe (1935); Tanaka (2002). As discussed in the Introduction, the pair approximation is not able to detect orientation-dependent features of the equilibrium phases. In this work, we implemented the 4-point approximation in the square lattice, which is able to capture the emergence of anisotropic nearest-neighbor correlations or spontaneous rotational symmetry breaking.

## Iii - Ising model in the 4-point (square) approximation

The - Ising model on the square lattice is defined by the Hamiltonian:

(11) |

where are Ising spin variables and is an external field. denotes a sum over pairs of nearest-neighbors and a sum over pairs of next-nearest-neighbors. In this work we consider and representing ferromagnetic NN and antiferromagnetic NNN interactions respectively. The competition ratio is defined by .

At zero external field the ground state of the model is ferro (if ) or antiferromagnetically ordered (if ) for and striped or superantiferromagnetic if . In the stripe phase the system can adopt one of four possible configurations as shown in Figure 1.

At zero external field the model has been extensively studied Morán-López et al. (1993); Moran-Lopez et al. (1994); Cirillo et al. (1999); dos Anjos et al. (2008); Kalz et al. (2009, 2011); Jin et al. (2012, 2013); Saguia et al. (2013), considering both ferromagnetic and antiferromagnetic NN interactions and antiferromagnetic NNN interactions. The nature of the thermal phase transition from the stripes to a disordered phase for was controversial. In the most recent studies combining Monte Carlo simulations and a series of analytical techniques it has been established that the line of phase transitions in the temperature versus plane is first order for and is continuous with Ashkin-Teller critical behavior for . The critical exponents change continuously in this regime between the 4-state Potts model behavior at to standard Ising criticality for Jin et al. (2012, 2013).

In the square approximation, the CVM free energy corresponding to the Hamiltonian (11) is given by Cirillo et al. (1999):

(12) | |||||

In the last equation, the sums over , , , , denote sums over all sites, NN pairs, NNN pairs, clusters of three sites and squares respectively.

Expressing the cumulant functions in terms of the cluster functions the variational free energy reads:

(13) | |||||

Following the general approach outlined in Section II it is convenient to write the variational free energy in terms of correlation functions given by:

(14) |

which are related to the reduced density matrices by:

(15) | |||||

Substituting these definitions onto (13) we get the variational free energy of the - model in the CVM square approximationCirillo et al. (1999):

(16) | |||||

After computing the traces one is left with an expression for the variational free energy in terms of a set of correlation functions representative of the approximation considered. The form of equation (16) makes clear that up to the pair approximation the two directions in the square lattice enter in a completely symmetric way, the different pairs of sites are decoupled. It is in the last term, when square plaquettes are considered, that the coupling between different directions in space can lead to novel behavior. The minimization of the free energy is in general a difficult task. The state (stationarity) equations are given by (10). It is possible to compute the derivatives and try to solve the set of coupled nonlinear equations of state. Instead of that, we preferred to minimize the variational free energy numerically for given sets of external parameters. was set for all calculations.

## Iv CVM results in an external field

For and small magnetic fields the ground state is striped (). When is ferromagnetic and antiferromagnetic, the stripe order is eventually destroyed by the presence of an external field, and all the spins become aligned with the field at . In the case both interactions are antiferromagnetic, for the ground state is still (), while in the interval and it becomes row shifted ()Yin and Landau (2009); de Queiroz (2011). The latter state consists of alternating ferro and antiferromagnetically ordered rows (or columns), with the ferromagnetic ones parallel to the field. For higher fields the equilibrium state of the system corresponds to a saturated paramagnet.

Here we consider the system at finite with ferromagnetic NN () and antiferromagnetic NNN () interactions for two different competition ratios and for which the ground state is striped. With the aim of searching for purely orientational nematic-like phases, i.e. phases without positional order, we minimized the CVM free energy of Eq. (16) for the parameters (related to the elementary square defined in (12)): , , , , , , , and . This choice implies possible orientational order along the or vertical direction. Note that local magnetizations on horizontal NN sites are allowed to be different in sign and also in absolute value. Correspondingly, the NN correlation functions may be different not only between the horizontal and vertical directions but also between the two vertical ones. With these choices the values of NNN correlations and square correlations are unique. The values of , , , and with that minimize the free energy (16) were calculated numerically for different reduced temperatures () and reduced external fields () using the routine NMinimize of the software Mathematica Weisstein (). In order to distinguish between translational (positional) order and orientational order we defined suitable order parameters; two positional order parameters:

(17) |

and

(18) |

describing ferromagnetic and stripe orders respectively. Note that, in case and are both finite but with different absolute values, a mixed phase with stripe order on a ferromagnetic background can be possible. In this cases, we have classified these phases as stripe ones. The orientational order parameter is defined as:

(19) |

In this case orientational order is finite whenever horizontal NN correlations are different from vertical ones, i.e. when there is a breaking of symmetry in the NN correlation functions. With these definitions all order parameters take values in the range .

The phase diagram in the plane for is shown in figure 3. For this value of the ground state is striped. A first difference in this phase diagram with respect to those in references Yin and Landau, 2009; de Queiroz, 2011 is the absence of the row-shifted phase at large fields. This is due to the ferromagnetic character of the NN interactions in the present work. The absence of the () phase is probably related with the absence of reentrance of the stripe phase in this case, at variance with the results reported for the model with antiferromagnetic NN interactions Yin and Landau (2009); de Queiroz (2011), as pointed out above. At , the 4-point approximation predicts a first order transition line for and a second order transition line for Morán-López et al. (1993); Moran-Lopez et al. (1994); Cirillo et al. (1999). For and finite we found a line of first order transitions marked by the discotinuity of the order parameter. The discontinuity is observed at both transitions, stripes-saturated paramagnetic and stripes-nematic (see Figure 3, first panel). Above the point where the positional order parameter goes to zero the system still finds itself in a phase with a finite value of the orientational order parameter , as seen in Figs. 3 and 5. This is the signature of a nematic-like phase in which the magnetization is homogeneous, unlike in the stripe phase, but correlations show an anisotropic character, reminiscent of the more ordered stripe phase. For the case the nematic phase is observed in the plane in the reduced temperature range , as seen in Fig. 3. The nematic phase terminates in a line of second order phase transitions (green line in Fig. 3) where the system enters a paramagnetic phase with finite magnetization values due to the external field (saturated paramagnet). The observation of the nematic phase in the - model in an external field is the main result of this work. In the context of the Cluster Variation Method it is clear that the 4-point approximation is the minimal one which is able to capture a nematic-like phase of this kind, i.e. a phase with broken orientational symmetry in the correlation functions. Nevertheless, this possibility was not exploited in previous work within the CVM.

The behavior of correlation functions for and as functions of the external field is shown in Figure 3. In the first panel of Figure 3 it is seen that the positional and orientational order parameters coincide in the stripe phase, as it should, but goes to zero before , signalling a stripe-nematic phase transition at . The second panel on the upper row shows that the local magnetizations and tend to be equal but with opposite signs at very small fields, in agreement with the stripe character of the ground state for small values. Nevertheless, they gradually evolve in an asymmetric way until at the stripe-nematic transition point their values merge in a single one, meaning the onset of an homogeneous phase with regard to local magnetizations. The anisotropic character of the nematic phase is evidenced in the third panel on the upper row of Figures 3 and 5. It is seen that, within the stripe phase, correlations along the vertical direction ( and ) are slightly different reflecting the slightly different values of the local magnetizations at those sites. A change in behavior between those correlations and the horizontal ones is observed at the stripes-nematic transition point. Note that if the stripe phase should have ended in a disordered rotationally symmetric state, then correlations in different directions should merge at this point. This does not happen until a larger field value where the three different NN correlations considered here merge to a single value at . The NNN correlations display a discontinous derivative at the stripe-nematic transition.

In Figures 5 and 5 we show the phase diagram and correlation functions for and . The behavior is qualitatively the same as the case with . We show them in order to compare with Monte Carlo simulation results to be discussed in the next section for . Although the results from the CVM and Monte Carlo simulations are qualitatively similar, important quantitative differences arise which should be addressed with higher order approximations in the CVM approach.

## V Monte Carlo simulations

We have performed Monte Carlo simulations of the - Ising model on the square lattice to test the qualitative consistency of the CVM results. The frustration present in the model turns computer simulations very demanding. The simulation procedure combines standard one-site movesLandau and Binder (2000) with a cluster method adapted from the simulation of patchy lattice models Tavares et al. (2014) following the ideas of the Wolff algorithmWolff (1989) in the presence of external fields Landau and Binder (2000). In the cluster method, one of the spins of the system, and one of the main directions of the lattice are chosen at random. The chosen spin is the starting point (root) of the cluster, then one starts growing the cluster by adding spins, , which are NN of the cluster in the chosen direction with probability provided that , with being the state of the spins in the cluster. Once the construction of the cluster is finished acceptance criteria are applied by taking into account the interaction of the cluster with the remaining spins of the lattice and the external fieldTavares et al. (2014); Landau and Binder (2000). The introduction of the cluster technique improves appreciably the numerical performance of the simulations, specially at low temperature. In order to enhance further the simulation peformance we have made use of parallel tempering (or replica exchange) Monte Carlo samplingSwendsen and Wang (1986); Earl and Deem (2005). This was carried out as follows: for a given fixed value of the external field (or temperature) we located via preliminary simulations the approximate value(s) of the temperarature (or field) where the transition(s) take(s) place; then we chose an appropriate set of values of the temperature (field) around such preliminary estimates to carry out the replica-exchange Monte Carlo simulations. The initial configurations were built by choosing the value of each spin at random. In order to guarantee the reliability of the final results we run shorter complementary simulations, using ground state configurations as starting point, to check that the simulations runs were properly equilibrated.

In the same spirit as with the CVM approach, our main interest was to search for possible nematic phases, i.e. phases with orientational order but lacking translational order. In order to distinguish the presence of such phases we defined suitable order parameters analogous to those defined previously in the CVM approach. The translational order parameter (TOP) was defined as the one used in Ref. Jin et al., 2013, in which the configurations of the system are basically compared with the ground state configurations at zero field. We analyzed the existence of periodicity in the lattice directions by computing the quantities:

(20) |

where is the linear size of the square lattice and are the coordinates of site . The global translational order parameter is then defined through:

(21) |

Notice that for the ground state structures shown in Fig. 1. This order parameter is equivalent to the definition used in (18) which was suitable for the elementary square of the CVM. The translational order can be tested by carrying out a finite-size scaling analysis. If no translational order exists one expects that the average of this order parameter will approach zero in the thermodynamic limit, whereas for the ordered case will be finite as . Moreover, if we analyse the translational order through an isotherm (at varying external field), or as a function of the temperature (at constant external field), we could expect that the possible order-disorder transitions involving translational ordering will appear as abrupt changes in , specially for large systems.

According to the ground states at zero field (Fig. 1), at low temperatures the system shows the tendency to form long sequences of spins in the same state along the main directions of the lattice and/or . The length of these sequences is expected to grow on decreasing the temperature and at some point the competition between sequences in the two directions might lead to a phase transition, in such a way that one of the directions will be preferred in the ordered phase. In this case the system breaks the fourfold (Z) symmetry of the square lattice reducing it to twofold (Z) symmetry. This transition from a disordered Z to a Z ordered phase is an orientational phase transition and need not be accompanied with the growth of translational or positional order. So we must distinguish positional and orientational order parameters, as discussed in relation to the CVM results.

We can define orientational order parameters (OOP) along directions and as:

(22) |

We will find that for the GS configurations one of the components of will have the value (that corresponding to the direction in which the sites are at the same state), whereas the other will take the value . A global order parameter is then defined as:

(23) |

This definition is equivalent to Eq. (19) which was suitable in the context of the square approximation of the CVM. In the thermodynamic limit, finite values of the statistical average of its absoute value will indicate that there is a preferential direction for the sequences of equal spins.

### v.1 Results

We have considered two values of , namely , and . For we have analysed five cases. In the first one, was fixed and we looked at the variation of different properties with the temperature and the system size. In addition, we took four values of the reduced temperature , , , and , and looked at the variation of the properties as a function of the external field, , and the system size at constant temperature.

Some results for , are shown in Fig. 6. In this case the system exhibits a disordered phase at high temperatures, where both order parameters tend to zero as the system size increases. At reduced temperature a phase transition occurs. A direct inspection of the results (See Fig. 6) for different properties as a function of for different system sizes leads to the following conclusions: (1) For both order parameters seem to vanish as grows larger; whereas for both order parameters converge to a finite value greater than zero for the larger system sizes; (2) The jump in both order parameters seems to occur at the same value of the temperature (); (3) At there is also a jump in the energy per site: , (Fig. 6.c); (4) The corresponding heat capacity at constant field, exhibits a clear single peak for systems with , with a value for the maximum that seems to diverge as (Fig 6.d). The scaling of the maximum of for the system sizes considered does not allow to establish the type of the transition. Two possible scenarios for the transition are, a 4-state Potts criticality (continuous transitions), with , or a very weak discontinuous transition as stated by Jin et al.Jin et al. (2013).

Next, we consider the case , and finite . The main difference with the case with is that the averaraged magnetization of the system, defined as , does not vanish in the thermodynamic limit. In Figure 8 we show some of the results for this system. From the results of the order parameters it can be deduced that at reduced field the system exhibits an order-disorder transition. The results suggest that the orientational and translational ordering occurs simultaneously. At the same value of it is observed a jump in the magnetization, whose derivative with respect to the external field at constant temperature, , seems to diverge as . The qualitative behavior of the energy per site with respect to (not shown), and its derivative as a function of is similar to the behavior for the case when both functions are plotted as functions of . Therefore, the features of the transition at are similar to those found at zero field. In both cases the orientational and translational orderings seem to occur cooperatively. This conclusion seems to be consistent with the presence of an unique peak in the susceptibility , as shown in Fig. 8, panel d.

On decreasing the temperature the qualitative features of the transitions change. We have carried out simulations at three additional values of reduced temperature: , , and (for ). Some qualitative differences arise between the phase transitions in these three low-temperature cases and the preceding two cases. In Fig. 8 we show the results for the case . First, the aparent common transition for TOP and OOP seems to split into two separated transitions, i.e. for a given temperature the jumps of TOP and OOP occur at different values of the external field . This scenario of two successive order-disorder transitions instead of just one transition is consistent with the incipient splitting of the peak of the susceptibility for the larger systems considered. This behavior is in qualitative agreement with the CVM results for corresponding parameter values observed in Fig. 5.

We have also explored the phase behavior for . Three cases were considered (1) Constant field ; (2) Constant temperature , and (3) Constant temperature . The phase behavior is similar to that found for . At zero field, a single transition is found at . For , the intermediate nematic phase (with only orientational order) does not appear, with the ordering transition occurring at . However for the Ising nematic phase seems to be stable for a very narrow range of in the region around . As in the case with the splitting of the isotropic-ordered transition into two transitions is only clearly observed for quite large system sizes , which prevents us from reaching conclusions about the nature of both transitions.

### v.2 Cumulant analysis of the order parameter distributions

In the analysis of the phase transitions of model systems it is quite useful to pay attention to the ratios between the momenta of the order parameter distributions. Here we considered the ratios , i.e. and , which are closely related with the so-called Binder cumulants Landau and Binder (2000). The analysis of the system size dependence of these quantities and, in particular, the crossings of the curves of versus some thermodynamic field () for different system sizes is often a very good choice to locate the phase transitions.

In Figure 9 we show the curves at constant temperatures for both order parameters and several system sizes. We can appreciate substantial qualitative differences in the shape of the functions defined on the order parameter from the cases (one transition), and (two transitions with an intermediate nematic phase) and . For the transition between the ordered phase ( )to the disordered phase occurs quite abruptly as the system size increases, and the curves for different system sizes seem to cross at , which coincides (or it is quite close) to the crossing point of the corresponding lines for the ratio. Notice that the maximum of the susceptibility (See Fig. 8) seems to happen exactly at the same value of . In the case of the departure of from the ordered phase value, , occurs at values of clearly smaller than the corresponding departure of . Looking at the cases with , we can observe that the values of at the crossings of the curves for different system sizes are consistent with the criticality of the two-dimensional Ising universality class Salas and Sokal (2000), as one could expect from the symmetry of the order parameter . In addition, the results for the largest systems indicate that in the range of values between the two transitions (as predicted by the maxima of the susceptibility ) exhibits a plateau, with values consistent with in the region where the nematic phase is supposed to be stable. The same type of results are found for the cases at .

### v.3 Gallery of configurations

In order to illustrate the differences between the stripe, nematic and disordered phases described in this work, in what follows we will present some representative configurations of the simulated systems for , considering the lattice gas version of the model. The following rules were applied to plot the configurations: (1) We consider only occupied sites (or ); (2) We plot segments between pairs of NN sites if, and only if, both sites are occupied; and (3) Four colors are considered, depending on the direction of the bond ( or ), and for each direction depending on the value of the complementary coordinate. Each color is related with each of the four configurations in the ground state shown in Fig. 1. From this representation of the configurations, we expect for isotropic phases segments in both directions and four colors with similar probabilities; for nematic phases most of the segments will be in one of the directions and two colors will be predominant; whereas for the ordered phase most of the segments will have the same direction and the same color.

We have chosen two cases. In the first one, shown in FIG. 10 we consider fixed value of the external field, () , and plot representative configurations for three temperatures in the neighborhood of the transition temperature. It can be seen that there are no signatures of the presence of the Ising nematic phase. Above the transition temperature (left panel) one can observe regions in the system where one of the four colors is predominant. At the estimated transition temperature (middle panel in FIG. 10) the system has developed a large region with most of the bonds in vertical direction and one predominant color, but still relatively large regions with the three remaining colors are still present. As the temperature is further reduced (right panel in FIG. 10) the regions with minoritary colors appear just as small islands.

In the second case, shown in FIG. 11, we considered fixed temperature at (). For this case we expect the stability of an intermediate nematic phase. Configurations in the left (isotropic phase) and right (ordered phase) panels show similar features, apart from the lower density of segments, to those found in FIG. 10; however the configuration in the middle panel shows clear signatures of the nematic phase. Most of the bonds are oriented in vertical direction, but none of the two colors associated with this direction is predominant, and no long range order correlation in horizontal direction can be appreciated.

## Vi Conclusions

We have shown that the J-J model in an external field has an intermediate phase with only orientational order. This is the main result of the present work. We have performed an analytical approach based on the Cluster Variation Method in the square approximation and compared the results with Monte Carlo simulations. Both approaches are in qualitative agreement. We did not find evidence of orientational intermediate phases of nematic type for zero external field, where our results are compatible with those already known from the literature. Nevertheless, in the presence of an external field a phase with orientational but without positional order emerges. This is compatible with an Ising-nematic phase with symmetry, characterized in this context by the spontaneous breaking of the symmetry of the square lattice. We found that an Ising-nematic phase exists in a finite window of external field values and temperatures.

For the parameter values studied we found that the disordered-nematic transition is of second order, with the order parameter going continuously to zero at . Preliminary Monte Carlo results indicate that this transition is probably in the Ising universality class. The nature of the second transition, from the nematic to a stripe phase with both orientational and positional orders, is more subtle. The CVM results give discontinuous transitions for the parameter values studied. Regarding the Monte Carlo results, it has been found that there are strong finite size effects. Then, simulations of very large system sizes are required to extract definitive conclusions, which are beyond our present capabilities. The simulation results presented here strongly suggest the existence of a stable Ising-nematic phase at low temperatures. In principle, according to the form of the order parameter , and to the crossings of the functions for different system sizes, one expects a continuous transition from the disordered to the nematic phase with 2D Ising criticality. Regarding the transition between nematic and fully ordered stripe phase it seems quite difficult to extract definitive conclusions with the type of calculations presented here. There are two basic problems: first, at intermediate temperatures (say, for ) the isotropic-nematic and nematic-stripe phase transitions are very close form each other, which makes difficult a finite-size scaling treatment. At lower temperatures, another difficulty arises. The correlation length (the length of the segments shown in the configuration plots) grows quickly as one reduces the temperature. This implies that, again, large systems have to be considered to extract conclusions. In addition, the nematic-stripe transition is clearly detected through the order parameter , the functions or the maxima of the susceptibility only for large system sizes.

###### Acknowledgements.

A.G.D. and D.A.S. acknowledge partial financial support by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), Brazil. N.G.A. acknowledge the support from the Ministerio de EconomÃa y Competitividad (Spain) under Grant FIS2013-47350-C5-4-R, and from Acciones bilaterales CSIC-CNPq (REF: 2011BF0046).## References

- Fischer and Hertz (1991) K. H. Fischer and J. A. Hertz, Spin Glasses (Cambridge University Press, 1991).
- Abanov et al. (1995) A. Abanov, V. Kalatsky, V. L. Pokrovsky, and W. M. Saslow, Phys. Rev. B 51, 1023 (1995).
- De’Bell et al. (2000) K. De’Bell, A. B. MacIsaac, and J. P. Whitehead, Rev. Mod. Phys. 72, 225 (2000).
- Vaterlaus et al. (2000) A. Vaterlaus, C. Stamm, U. Maier, M. G. Pini, P. Politi, and D. Pescia, Phys. Rev. Lett. 84, 2247 (2000).
- Imperio and Reatto (2006) A. Imperio and L. Reatto, The Journal of Chemical Physics 124, 164712 (2006).
- Archer and Wilding (2007) A. J. Archer and N. B. Wilding, Phys. Rev. E 76, 031501 (2007).
- Olson Reichhardt et al. (2010) C. J. Olson Reichhardt, C. Reichhardt, and A. R. Bishop, Phys. Rev. E 82, 041502 (2010).
- Ciach (2011) A. Ciach, Molecular Physics 109, 1101 (2011).
- Pȩkalski et al. (2013) J. Pȩkalski, A. Ciach, and N. G. Almarza, The Journal of Chemical Physics 138, 144903 (2013).
- Pȩkalski et al. (2014) J. Pȩkalski, A. Ciach, and N. G. Almarza, The Journal of Chemical Physics 140, 114701 (2014).
- Almarza et al. (2014) N. G. Almarza, J. Pȩkalski, and A. Ciach, The Journal of Chemical Physics 140, 164708 (2014).
- Mézard and Montanari (2009) M. Mézard and A. Montanari, Information, Physics, and Computation (Oxford University Press, 2009).
- Wales (2003) D. Wales, Energy Landscapes: Applications to Clusters, Biomolecules and Glasses (Cambridge University Press, Cambridge, UK, 2003).
- Morán-López et al. (1993) J. L. Morán-López, F. Aguilera-Granja, and J. M. Sanchez, Phys. Rev. B 48, 3519 (1993).
- Moran-Lopez et al. (1994) J. L. Moran-Lopez, F. Aguilera-Granja, and J. M. Sanchez, Journal of Physics: Condensed Matter 6, 9759 (1994).
- Cirillo et al. (1999) E. N. M. Cirillo, G. Gonnella, M. Troccoli, and A. Maritan, Journal of Statistical Physics 94, 67 (1999).
- dos Anjos et al. (2008) R. A. dos Anjos, J. R. Viana, and J. R. de Sousa, Physics Letters A 372, 1180 (2008).
- Kalz et al. (2009) A. Kalz, A. Honecker, S. Fuchs, and T. Pruschke, Journal of Physics: Conference Series 145, 012051 (2009).
- Kalz et al. (2011) A. Kalz, A. Honecker, and M. Moliner, Phys. Rev. B 84, 174407 (2011).
- Jin et al. (2012) S. Jin, A. Sen, and A. W. Sandvik, Phys. Rev. Lett. 108, 045702 (2012).
- Jin et al. (2013) S. Jin, A. Sen, W. Guo, and A. W. Sandvik, Phys. Rev. B 87, 144406 (2013).
- Saguia et al. (2013) A. Saguia, B. Boechat, J. Florencio, and O. F. de Alcantara Bonfim, Phys. Rev. E 87, 052140 (2013).
- de Queiroz (2011) S. L. A. de Queiroz, Phys. Rev. E 84, 031132 (2011).
- Yin and Landau (2009) J. Yin and D. P. Landau, Phys. Rev. E 80, 051117 (2009).
- Cannas et al. (2006) S. A. Cannas, M. F. Michelon, D. A. Stariolo, and F. A. Tamarit, Phys. Rev. B 73, 184425 (2006).
- Barci et al. (2013a) D. G. Barci, A. Mendoza-Coto, and D. A. Stariolo, Phys. Rev. E 88, 062140 (2013a).
- Nicolao and Stariolo (2007) L. Nicolao and D. A. Stariolo, Phys. Rev. B 76, 054453 (2007).
- Han et al. (2001) J. Han, Q.-H. Wang, and D.-H. Lee, Int. J. Mod. Phys. B 15, 1117 (2001).
- Berg et al. (2009) E. Berg, E. Fradkin, S. A. Kivelson, and J. Tranquada, New J. Phys. 11, 115004 (2009).
- Fradkin et al. (2010) E. Fradkin, S. A. Kivelson, M. J. Lawler, J. P. Eisenstein, and A. P. Mackenzie, Annual Review of Condensed Matter Physics 1, 153 (2010).
- Barci and Stariolo (2009) D. G. Barci and D. A. Stariolo, Phys. Rev. B 79, 075437 (2009).
- Barci et al. (2013b) D. G. Barci, L. Ribeiro, and D. A. Stariolo, Phys. Rev. E 87, 062119 (2013b).
- Saratz et al. (2010a) N. Saratz, A. Lichtenberger, O. Portmann, U. Ramsperger, A. Vindigni, and D. Pescia, Phys. Rev. Lett. 104, 077203 (2010a).
- Saratz et al. (2010b) N. Saratz, U. Ramsperger, A. Vindigni, and D. Pescia, Phys. Rev. B 82, 184416 (2010b).
- Portmann et al. (2010) O. Portmann, A. Gölzer, N. Saratz, O. V. Billoni, D. Pescia, and A. Vindigni, Phys. Rev. B 82, 184409 (2010).
- Cannas et al. (2011) S. A. Cannas, M. Carubelli, O. V. Billoni, and D. A. Stariolo, Phys. Rev. B 84, 014404 (2011).
- Velasque et al. (2014) L. A. Velasque, D. A. Stariolo, and O. V. Billoni, Phys. Rev. B 90, 214408 (2014).
- Stariolo and Barci (2010) D. A. Stariolo and D. G. Barci, J. Phys: Conf. Ser. 246, 012021 (2010).
- Barci and Stariolo (2011) D. G. Barci and D. A. Stariolo, Phys. Rev. B 84, 094439 (2011).
- Bethe (1935) H. A. Bethe, Proc. Roy. Soc. A 150, 552 (1935).
- Tanaka (2002) T. Tanaka, Methods of Statistical Physics (Cambridge University Press, 2002).
- Kikuchi (1951) R. Kikuchi, Phys. Rev. 81, 988 (1951).
- Sanchez et al. (1984) J. Sanchez, F. Ducastelle, and D. Gratias, Physica A: Statistical Mechanics and its Applications 128, 334 (1984).
- An (1988) G. An, Journal of Statistical Physics 52, 727 (1988).
- (45) E. W. Weisstein, “Minimization. From MathWorld—A Wolfram Web Resource,” http://www.wolfram.com/mathematica/ .
- Landau and Binder (2000) D. P. Landau and K. Binder, Monte Carlo Simulations in Statistical Physics (Cambridge University Press, 2000).
- Tavares et al. (2014) J. M. Tavares, N. G. Almarza, and M. M. Telo da Gama, The Journal of Chemical Physics 140, 044905 (2014).
- Wolff (1989) U. Wolff, Phys. Rev. Lett. 62, 361 (1989).
- Swendsen and Wang (1986) R. H. Swendsen and J.-S. Wang, Phys. Rev. Lett. 57, 2607 (1986).
- Earl and Deem (2005) D. J. Earl and M. W. Deem, Phys. Chem. Chem. Phys. 7, 3910 (2005).
- Salas and Sokal (2000) J. Salas and A. D. Sokal, Journal of Statistical Physics 98, 551 (2000).