text:
# High-resolution X-ray imaging and spectroscopy of the core of NGC 4945 with XMM-Newton and Chandra

## 1 Introduction

NGC 4945 is a nearby edge-on (i $``$78o; Ott, 1995) spiral galaxy (type SBcd or SABcd), believed to be a member of the Centaurus group (Hesser et al. 1984), at a distance of between 3.7 Mpc and 8.1 Mpc (here we use 3.7 Mpc; Mauersberger et al. 1996). Near infra-red observations have revealed the nuclear region of NGC 4945 to contain a powerful, yet visually obscured starburst region with a $``$200pc ($``$10′′) ring morphology (Moorwood et al. 1996) and a total infrared luminosity (8-1000 $`\mu `$m) of 2.2$`\times `$1010L (Spoon et al. 2000).

Although much of the central activity can be explained in terms of a nuclear starburst, the presence of an active galactic nucleus (AGN) in NGC 4945 has been confirmed by the detection of a luminous and variable hard X-ray source, coincident with the centre of the galaxy (Iwasawa et al. 1993; Guainazzi et al. 2000). NGC 4945 harbours the brightest known Seyfert 2 nucleus at 100 keV (Done et al. 1996), but extremely heavy obscuration cuts-off the direct AGN continuum at energies below $``$10 keV (Guainazzi et al. 2000; Madejski et al. 2000). The measured X-ray column density ($``$4$`\times `$1024 atom cm-2) places it on the threshold of being a Compton-thick Seyfert 2. After correctly accounting for the effects of Compton scattering in the absorbing material (Matt et al. 1999), the inferred intrinsic X-ray (0.1-200 keV) luminosity is $``$1.8$`\times 10^{43}`$erg s-1 (Guainazzi et al. 2000). The soft to medium energy X-ray spectrum of the nucleus is comprised of both thermal emission associated with a nuclear starburst and emission generated through the reprocessing of the AGN continuum emission (Guainazzi et al. 2000).

The dense clouds that envelope the nuclear region produce a rich variety of molecular lines (Curran et al. 2001) including $`\mathrm{H}_2\mathrm{O}`$ megamaser emission, which requires an edge-on inner disk geometry (Greenhill, Moran & Hernstein 1997). The megamaser emission provides a tight constraint on the mass of the central black hole, $`M_{BH}1.4\times 10^6M_{}`$, implying that the AGN is radiating at up to 60% of its Eddington luminosity.

Infrared and optical observations of NGC 4945 have also revealed the presence of a conical cavity, attributed to a starburst-driven superwind, protruding above the disk of the galaxy (Chen & Huang 1997; Moorwood et al. 1996). Despite heavy obscuration at optical wavelengths, extended line-emitting gas, co-spatial with the conical cavity, has been detected in outflow along the galaxy’s minor-axis and is attributed to the interaction between the superwind and the surrounding medium (Heckman, Armus & Miley 1990). The low spatial resolution of the current published data has mitigated against the detection of the superwind in X-rays, although Guainazzi et al. (2000) suggest that at least part of the soft X-ray emission from NGC 4945 may originate in this component.

Here we present new XMM-Newton and Chandra X-ray observations of NGC 4945. The X-ray imaging reveals a complex morphology within the nuclear region of NGC 4945 and considerable emission associated with the disk of the galaxy, including at least 12 discrete X-ray sources. In this paper we will confine our attention to a detailed analysis of the nuclear region of NGC 4945; the remaining diffuse emission associated with the galactic disk and the off-nuclear point source population will be discussed in a separate paper (Roberts, Schurch & Warwick 2002; in preparation).

## 2 Observations and Data Reduction

### 2.1 The XMM-Newton Observation

NGC 4945 was observed with XMM-Newton during orbit 205 (2001, January 21) for $``$24 ks as part of the mission Guaranteed Time programme. The EPIC MOS and PN instruments were operated in full frame mode, with the medium filter, for the full duration of the observation. Detailed summaries of the XMM-Newton mission and instrumentation can be found in Jansen et al. (2001; and references therein). The raw data were processed with the public release version of the XMM-Newton Science Analysis System (SAS v 5.2) standard processing chains. After filtering for data flagged as bad, X-ray events corresponding to pattern 0-12 for the two MOS cameras (similar to grades 0-4 for ASCA ) and pattern 0-4 for the PN camera (single and double pixel events only) were accepted. Investigation of the full-field count-rate revealed no significant background flaring episodes in the observation. The effective exposure times for the MOS and PN cameras were 22.2 ks and 19.2 ks respectively. Images, spectra and corresponding background information were extracted with the SAS task XMMSELECT.

### 2.2 The Chandra observation

The archival Chandra data were obtained from the UK mirror of the Chandra X-ray Center data archive, maintained by the Leicester Data Archive Service (LEDAS). The Chandra observation itself was performed on 2000 January 27/28, resulting in a total exposure of 50 ks, with the ACIS-S array at the prime focus. Our reduction and analysis of this dataset was performed using version 2.1 of the CIAO software, and incorporated many of the standard analysis procedures available online from cxc.harvard.edu/ciao. The data were initially filtered to remove all events outside the 0.3-10 keV energy range. The analysis of the full-field light-curve revealed the observation to be heavily contaminated by background flaring. The removal of all times for which the total detector count rate exceeded 8 counts s-1 was found to dramatically reduce the background in the S3 chip, leaving a GTI-corrected event-list containing 33.1 ks of “good” data.

## 3 Joint imaging with CHANDRA and XMM-NEWTON

XMM-Newton and Chandra images were extracted for a 2 square region centered on the nucleus of NGC 4945, revealing a wealth of detail including a predominantly hard nuclear source and a soft X-ray “plume” extending to the northwest of the nuclear source. Fig. 1: Panels (a)-(c), show high spatial resolution energy-coded X-ray images from the Chandra ACIS detector, while Fig. 1: Panel (d) shows an ‘X-ray colour’ image from the co-added XMM-Newton EPIC cameras. Whereas the nuclear source appears point-like with XMM-Newton, the Chandra imaging partially resolves the emission centred on the active nucleus (Fig, 1: Panel (c) ). The nuclear point source is revealed as an extremely hard source by both the Chandra hard band image and the XMM-Newton X-ray colour image. This spectral hardness is almost certainly due to heavy line-of-sight obscuration. Both XMM-Newton and Chandra clearly resolve the X-ray plume, establishing its orientation and extent along the minor axis of the galaxy (30′′= 500 pc, NW) and also its intrinsic spectral softness (Fig. 1: Panels (a),(b),(d) ). The medium band Chandra image reveals the X-ray plume to have a limb-brightened morphology, qualitatively similar to the morphology of the X-ray plume seen in NGC 253 (Strickland et al. 2000).

## 4 Joint spectroscopy with CHANDRA and XMM-NEWTON

The high spatial resolution afforded by Chandra has allowed us to spatially resolve several of the major components that comprise the total X-ray emission from the nucleus of NGC 4945. However, the relatively low count-rate means the Chandra data cannot be used to place tight constraints on the spectral form of these components. Conversely, despite the greater count-rate, the lower angular resolution of XMM-Newton gives rise to some spatial confusion, which complicates the identification and spectral modelling of the components of the nuclear region. Fortunately, the joint analysis of the Chandra and XMM-Newton spectra removes many of these limitations.

We have adopted the following approach. We extract separate Chandra ACIS spectra for the AGN and the extended emission components and model these in turn. The identified spectral components are then combined into a composite spectral model which is fitted to the high signal-to-noise XMM-Newton EPIC spectra of the whole nuclear region. Finally, we consider combined fits to both the Chandra and XMM-Newton spectra. As expected, the latter provides the tightest constraints on the individual model parameters.

### 4.1 The Chandra AGN spectrum

The spectrum of the point-like nuclear source (which we associate predominantly with the AGN in NGC 4945) was derived from the Chandra GTI-corrected event list using a circular extraction region of radius 6 pixels (3′′) centred on the peak of the X-ray emission, at 13h05m27.5s, -492805′′(J2000). The corresponding background spectrum was taken from a local source-free region. The resulting spectrum contains $``$750 counts (after background subtraction) with no significant signal below 1.5 keV.

In modelling the spectral form of the AGN we note that in several well studied Seyfert 2 galaxies, with column densities in excess of a few $`\times 10^{23}\mathrm{atoms}\mathrm{cm}^2`$, it appears that a relatively unabsorbed Compton-reflection component dominates the 2–10 keV spectrum (e.g. Mrk 3 - Cappi et al. 1999, Sako et al. 2000; NGC 3281 - Vignali 2002). In such sources, the putative molecular torus may serve as a reprocessing site for the hard continuum generated by the AGN. Furthermore, for a highly inclined configuration, it is perfectly plausible that the geometry is such as to give a relatively clear view of the far-side wall of the torus. It follows that the observed Compton-reflection in Seyfert 2’s may originated from this “far-wall” location.

We have tested this hypothesis for NGC 4945 by fitting the Chandra AGN spectrum with a model consisting of a very heavily absorbed power-law continuum, unabsorbed Compton reflection of the continuum (pexrav model in XSPEC; Magdziarz & Zdziarski 1995) and an iron K$`\alpha `$ fluorescence line at 6.4 keV. It turns out that the underlying hard continuum in this model has little impact even at 10 keV, implying neither the Chandra or XMM-Newton data can provide useful constraints on the parameters of the underlying continuum. Thus, the continuum characteristics (photon index, $`\mathrm{\Gamma }`$, normalisation, Kpl and absorbing column, $`N_H`$) were fixed at the values consistent with previous analyses of the high energy ($`>`$10 keV) spectrum from both BeppoSAX and RXTE . Specifically we take $`\mathrm{\Gamma }`$=1.55, Kpl=0.1 photon s-1 cm-2 keV-1 and $`N_H`$=4$`\times `$1024 cm-2 (Guainazzi et al. 2000; Madejski et al. 2000). Also, the metal abundances were fixed at solar values, the high-energy spectral cutoff set at 200 keV and the relative reflection parameter (R) was initially constrained to lie in the range 0–2 with cos(i)$`=0.5`$.

This model, with both the R parameter and the energy, intrinsic width and normalisation of the iron K$`\alpha `$ line as free parameters, provides a good fit to the Chandra spectrum ($`\chi ^2`$=32 for 29 degrees of freedom; hereafter d.o.f) as illustrated in Fig. 2 (Upper Panel). The best-fit parameter values are given in Table 1 (Model A), where the quoted errors correspond to 90% confidence levels for one interesting parameter (i.e. $`\mathrm{\Delta }\chi ^2`$=2.71). Modelling the Chandra spectrum in terms of a heavily absorbed direct continuum component (but with $`N_H`$ an order of magnitude smaller than assumed above) plus a relatively unabsorbed scattered power-law component (cf. Guainazzi et al. 2000) yields a considerably worse fit ($`\chi ^2`$=51.1 for 29 d.o.f). The fit residuals (Fig. 2, Upper Panel) hint at some contamination of the Chandra AGN spectrum by a nuclear starburst component; our modelling of the composite spectrum takes account of this in a straightforward fashion (see §4.3).

### 4.2 The Chandra starburst/X-ray plume spectrum

The Chandra spectrum of the extended nuclear X-ray emission encompassing the resolved nuclear starburst component and the X-ray plume structure was constructed using an elliptical extraction region centred on 13h05m25.5s, -492754′′(J2000), with the major axis aligned along a position angle of $`120^{}`$, with a semi-major axis of 30.5′′and a semi-minor axis of 13.4′′. In this case, the AGN component was suppressed by excluding the region within 3′′of the nucleus. The resulting spectrum (after background subtraction using a local source-free region) contains $``$1600 counts. Chandra is capable of resolving the starburst region from the X-ray plume, but here we have chosen to utilize the better signal-to-noise of the combined spectrum. Taking as a lead the published results for the starburst and X-ray plume in NGC 253 (Pietsch et al. 2001), we have modelled the Chandra spectrum with three solar abundance thermal (mekal) components, each exhibiting different levels of absorption as well as different temperatures and normalisations.

This model yields a good fit to most of the observed spectral features, including the significant K$`\alpha `$ emission from helium-like iron at 6.7 keV. However it does not account for the K$`\alpha `$ emission from neutral iron observed at 6.4 keV. Analysis of the individual starburst and X-ray plume spectra confirm that both K$`\alpha `$ emission lines originate entirely from the starburst region. We model the neutral iron K$`\alpha `$ emission with the addition of an intrinsically narrow Gaussian line. A model comprising three thermal components plus an iron K$`\alpha `$ line provides a reasonable fit to the data ($`\chi ^2`$=89 for 67 d.o.f, see Fig. 2; Lower Panel). The corresponding best-fit model parameters are listed in Table 1 (Model B).

### 4.3 Joint fitting of the composite XMM-Newton and Chandra spectra

The XMM-Newton spectra for the combined AGN, starburst and X-ray plume components were extracted using a 30′′ radius circle centered on 13h05m28s, -492801′′(J2000). For both the MOS and PN cameras, a corresponding background spectrum was extracted from a 90′′ radius circle in a source free region of the same CCD. After background subtraction, each MOS spectrum contained $``$1350 counts compared to $``$3500 counts in the PN spectrum.

The XMM-Newton spectrum was analysed by combining the two models used in the description of the Chandra observations (Models A & B in §4.1 & §4.2). This composite model yields an excellent fit to the broadband XMM-Newton spectra ($`\chi ^2`$=273.5 for 268 d.o.f). However, an unphysical trade-off between the hottest thermal component and the Compton-reflection parameter result in a value of R which is inconsistent with the value found from the Chandra AGN spectrum alone (assuming that this component has not varied significantly over the one-year interval between the two observations, see §5.2). We have addressed this problem by performing joint spectral fitting of three spectral datasets, namely the XMM-Newton composite spectrum and the Chandra AGN and starburst/X-ray plume spectra. Only the Compton-reflection component, the iron line and the hottest thermal component were utilized in fitting the Chandra AGN spectrum, whilst only the three thermal components and the iron line were employed in relation to the Chandra starburst/X-ray plume spectrum. We also adapted the spectral model somewhat by fixing the absorption applied to the softest thermal component at the Galactic value and setting the absorption on the Compton-reflection component to be the same as that applied to the intermediate-temperature thermal component (i.e. $`N_{H,2}`$).

Simultaneous fitting of the composite spectral model described above gives a good match to all three spectral datasets. The resulting best-fit to the XMM-Newton data is shown in Fig. 3; the corresponding parameters are given in Table 1 (Model C). The high signal-to-noise of the XMM-Newton spectra, relative to the Chandra spectra, leads to tighter constraints on all the parameters. The simultaneous fits highlight the important contribution of the highest temperature thermal (starburst) component to the 0.5-10 keV spectrum of NGC 4945. In this respect the joint spectral fitting is particularly valuable since the Chandra starburst/outflow spectrum has insufficient counts to put useful limits on the temperature of this component.

When data relating solely to the spatially-resolved X-ray plume are extracted (from both the Chandra and XMM-Newton datasets), the spectrum is found to be extremely soft with few counts above 2 keV. Spectral fitting using a single temperature mekal model, modified by Galactic absorption, gives a temperature and normalisation completely consistent with those of the softest thermal component in the composite spectral model (Model C), identifying this component with the X-ray plume emission. The two remaining thermal components are attributable to the nuclear starburst.

## 5 Discussion

The components of the broadband spectrum of the extended nucleus of NGC 4945 are shown in Fig. 4. The Compton-reflection component makes a significant contribution between 5 and 10 keV, whereas at lower energies the spectrum is dominated by the thermal emission from the nuclear starburst and X-ray plume. The heavily absorbed AGN continuum only emerges above 10 keV. Strong K$`\alpha `$ emission lines from neutral and helium-like iron are also present in the spectrum at 6.4 and 6.7 keV respectively.

### 5.1 The reprocessing of the AGN continuum

In the unified model proposed by Antonucci (1993), when an active nucleus is viewed at relative high inclination, the line of sight may intercept the surrounding molecular torus and hence be subject to significant obscuration. Previous observations testify to a near edge-on configuration of the inner disk in NGC 4945 (Greenhill, Moran & Hernstein 1997) and a very high line-of-sight column density to its active nucleus (Guainazzi et al. 2000; Madejski et al. 2000). Consistent with this picture, the analysis of the XMM-Newton and Chandra spectra of the nuclear region of NGC 4945 demonstrates that in the 0.5-10 keV range, there is only indirect evidence for the presence of the AGN (since the direct continuum is strongly cut-off below $`10`$ keV).

One signature of the underlying AGN in NGC 4945 is the Compton-reflection continuum which we attribute to reprocessing off the far, inside wall of the molecular torus. In this scenario the ratio of the reflected to the direct continuum as embodied in Compton-reflection parameter, R, can be used to place a lower limit on the inclination angle of the torus, provided the geometrical configuration of the absorbing torus is known. In reality this is very uncertain, although the half-opening angle of tori are typically in the range $`20^{}50^{}`$ (e.g. Wilson & Tsvetanov 1994). If we assume the geometry envisaged by Ghisellini, Haardt & Matt (1994), then our estimate of $`R=0.016`$ (based on the ‘average’ continuum given in Guainazzi et al. 2000) implies a near edge-on view of the torus. Alternatively a near edge-on but geometrically thin (half-angle ¡10o) distribution of absorbing material located on parsec scales, as inferred by Madejski et al. (2000) on the basis of the hard X-ray variability seen in NGC 4945, is also consistent with a very low value for R.

The neutral iron K$`\alpha `$ line is narrow and the nuclear component has an equivalent width of 1.6$`\pm `$0.1 keV with respect to the observed Compton-reflection continuum. This is fully consistent with the line and the reflected continuum originating at the same location, namely the visible part of the torus wall (e.g. Matt, Brandt & Fabian 1996). There is also a $`30\%`$ contribution to the neutral iron K$`\alpha `$ emission from the extended nuclear starburst region, implying that cool material in the environment of the starburst is also exposed to the AGN continuum flux.

Although the hard continuum emanating from the Seyfert 2 nucleus of NGC 4945 is known to be variable on timescales as short as $`10^4`$s (Guainazzi et al. 2000), no significant variability is detected in either the Chandra or XMM-Newton observations, nor is it evident from the joint spectral fits (the time interval between the two sets of observations being 12 months). This tallies with the location of the continuum reprocessing, namely the torus wall, lying at a distance in excess of a light year from the nuclear continuum source. (The apparent discrepancy between the Chandra -only and the joint-fit estimates of R is due to the fact that the initial modelling of the Chandra AGN spectrum did not account for the starburst contamination of the AGN spectrum).

### 5.2 The nuclear starburst

The spectral modelling discussed above characterises the nuclear starburst in terms of two thermal components with temperatures of $`0.9`$ and $`6.0`$ keV. Very similar temperatures were observed (at least in terms of the hottest thermal components) in the nearby starburst galaxy NGC 253 by Pietsch et al. (2001), suggesting that the presence of the AGN in NGC 4945 has little impact on the thermal properties of its nuclear starburst (there is no evidence for an AGN in NGC 253). Gross features, associated with the cooler of the two starburst components, include the emission lines of Mg X/XI (1.34 keV), Si XII/XIII (1.83 keV) and S XIV/XV (2.45 keV). At higher energy, strong K$`\alpha `$ line emission from helium-like iron is attributable to the hotter starburst component. High temperature plasma and associated line emission is characteristic of X-ray emission from a population of type Ib and type IIa SNR’s (e.g Behar et al. 2001). No corresponding iron edge is detected in either of the Chandra or XMM-Newton observations. The unabsorbed luminosity of the starburst region measured from XMM-Newton (0.5-10 keV) is $`8\times 10^{39}`$ erg s-1. Analysis of the separate starburst and X-ray plume spectra shows that the nuclear starburst emission declines rapidly below $`1.5`$ keV, probably as a result of absorbing material in the plane of the galaxy. This absorption may originate in material forming a low scale height disk and/or in material associated with the optical dust lanes of the galaxy.

### 5.3 The X-ray plume

Recent Chandra (Strickland et al. 2000) and XMM-Newton (Pietsch et al. 2001) observations have established that the outflow driven by the nuclear starburst in NGC 253 has a limb-brightened conical X-ray morphology. NGC 4945 provides a second example of such a structure. In both galaxies the limb-brightened X-ray structure correlates well with the observed the H$`\alpha `$ emission. A similar X-ray temperature is derived in each case. (In NGC 253 a single temperature thermal plasma with $`kT0.6`$ keV fits the Chandra spectrum reasonably well, although the situation is more complicated in the better quality XMM-Newton data.) However, there are also significant differences. For example, the NGC 4945 X-ray plume displays a wider opening angle ($`\theta 40^{}`$), and the limb brightening is most evident in the 1–2 keV, band whereas in NGC 253 both the Chandra and XMM-Newton data show the limb brightened structure is apparent down to 0.5 keV. A similar physical model to that postulated by Pietsch et al. (2001) can be inferred for NGC 4945. The limb-brightened structure can be attributed to highly-excited gas with a low volume-filling factor, which is produced by an interaction between the starburst-driven wind and the dense ISM surrounding the outflow. In the case of NGC 4945, the comparative uniformity of the emission below 1 keV may be a direct observation of a mass-loaded superwind emanating from the nuclear starburst (as described by Suchkov et al. 1996). The question of why the limb-brightened component is hotter in NGC 4945 than in NGC 253 is an interesting point for future, detailed studies of this system. We derive an intrinsic luminosity of $`1.1\times 10^{38}\mathrm{erg}\mathrm{s}^1`$ (0.5 - 2 keV) for the superwind component from the XMM-Newton data.

## 6 Conclusions

We present a detailed imaging and spectroscopic analysis of the nuclear region of NGC 4945 using data from both Chandra and XMM-Newton . The main results of our study are:

* The nuclear regions display a complex X-ray morphology, including a predominantly hard source at the nuclear position (which is partially resolved in the Chandra data) and a soft X-ray “plume” extending to the northwest of the nucleus.
* The 0.5-10 keV emission from the nuclear source is dominated by a combination of reprocessed direct AGN continuum (in the form of neutral Compton reflection and iron K$`\alpha `$ line emission) and thermal emission from a nuclear starburst. The underlying AGN continuum is not observed directly in either the Chandra or XMM-Newton data, due to heavy ($`4\times 10^{24}`$ cm-2) absorption in the line of sight to the continuum source, that cuts off the continuum below $``$10 keV.
* The neutral iron K$`\alpha `$ line is well fit with a narrow Gaussian profile (with an equivalent width of 1.6$`\pm `$0.1 keV), consistent with the line and the reflected continuum originating from the visible part of the torus wall. Furthermore, the best-fit value of the relative reflection parameter, $`R=0.016`$, implies a near edge-on geometry for the torus.
* Neutral iron K$`\alpha `$ emission from the extended nuclear starburst region, identified with Chandra , implies that cool material in the environment of the starburst is also exposed to the AGN continuum flux. The thermal emission from the nuclear starburst is fit with a two temperature model, with $`kT`$0.9 and $``$6.9 keV. The spectrum shows clear features associated with the two starburst components, including emission lines of Mg X/XI, Si XII/XIII and S XIV/XV from the cooler component and strong K$`\alpha `$ line emission from helium-like iron originating in the hotter component.
* A Chandra image reveals the X-ray plume to have a limb-brightened morphology in the 1–2 keV band. The plume is associated with the softest emission observed by XMM-Newton and is well modelled by a single temperature thermal model, with $`kT`$0.6 keV. When combined, these properties suggest a physical interpretation of the plume as a mass-loaded superwind emanating from the nuclear starburst.

Finally, we note that combining the Chandra and XMM-Newton datasets has helped to mitigate the limitations of the individual observations and allowed us to carry out a much more detailed imaging and spectral study than would otherwise have been possible. This work clearly demonstrates the potential of joint Chandra and XMM-Newton analysis in characterising spatially and spectrally complex X-ray sources, such as those associated with the extended nuclear regions of nearby galaxies.

ACKNOWLEDGMENTS

NJS gratefully acknowledges the financial support from PPARC. It is a pleasure to thank Kevin Briggs and Dick Willingale for their valuable aid with the Q software system. We also wish to thank the XMM-Newton and Chandra instrument and calibration teams for their on-going efforts to fully characterise and calibrate their respective instruments. The archival Chandra data were obtained from the UK mirror of the Chandra X-ray Center Data Archive operated by the Leicester Data Archive Service (LEDAS). This research has made extensive use of NASA’s Astrophysics Data System Abstract Service.


format: text

text:
# Exact partition functions of the Ising model on 𝑀×𝑁 planar lattices with periodic-aperiodic boundary conditions

## 1 Introduction

Universality and scaling are two important concepts in the theory of critical phenomena  and the Ising model  has been a model widely used in such studies. Recently, exact universal amplitude ratios and finite-size corrections to scaling in critical Ising model on planar lattices have received much attention . This may be due to the fact that the hypothesis of universality leads naturally to the consideration of universal critical amplitudes and amplitude combinations , and for the comparison between experiment and theory in relation to scaling and universality, it is often a more rigorous test to use amplitude relations rather than critical exponent values. Moreover, it is also well known that the finite-size scaling functions depend on the boundary conditions , and there has been considerable recent interest in studying lattice model with various boundary conditions . The study of exact universal amplitude ratios and finite-size corrections to scaling in critical Ising model is usually based on the analytical solutions of the model on finite lattices. Although the exact solution of the Ising model on $`M\times N`$ square (sq) lattice had been obtained long time ago , and the exact expression of the partition function of the Ising model on $`M\times N`$ plane triangular (pt) lattice has been obtained by lattice field theories recently , there is still no published results for the exact solutions of the Ising model on $`M\times N`$ pt and honeycomb (hc) lattices with periodic-aperiodic boundary conditions. The purpose of this paper is to fill this gap. In the present paper we use the Grassmann path integral to calculate exact partition functions of the Ising model on $`M\times N`$ sq, pt and hc lattices with periodic-periodic (pp), periodic-antiperiodic (pa), antiperiodic-periodic (ap) and antiperiodic-antiperiodic (aa) boundary conditions. The partition functions are used to calculate and plot the specific heat, $`C/k_B`$, as a function of the temperature, $`\theta =k_BT/J`$. We find that for the $`N\times N`$ sq lattice, $`C/k_B`$ for pa and ap boundary conditions are different from those for aa boundary conditions, but for the $`N\times N`$ pt and hc lattices, $`C/k_B`$ for ap, pa, and aa boundary conditions have the same values. Our exact partition functions might also be useful for understanding the effects of lattice structures and boundary conditions on critical finite-size corrections of the Ising model.

Two-dimensional Ising model on the sq lattice at vanishing magnetic field was first solved by Onsager by the use of Lie algebra . The exact solution he obtained was Ising model on an infinite lattice. The original method was rather complicated, and it was later improved by Kaufman  who obtained the exact solution of the Ising model on a finite torus by using the theory of spinor representation. The successful treatments of the two-dimensional Ising model brought the studies of phase transition into the modern era. Onsager’s solution in one hand showed the previous classical theories were unreliable in their quantitative predictions, and on the other hand provided a great stimulus to explore the true behaviour near the critical point. After Onsager’s original solution, many quite different mathematical approaches were developed, but the approaches were still complicated. Among them, Schultz, Mattis and Lieb gave explicitly the fermionic treatment in the framework of transfer-matrix formalism , and Kac and Ward developed the combinatorial method . Both methods reformulated the two-dimensional Ising model as a free-fermionic field theory in terms of anticommuting Grassmann variables, which enclosed the fact that the Ising model on two dimensional regular lattices may be viewed as free-fermionic theory. The other alternative method in literature was the Pfaffian representation, which was introduced by Kasteleyn  to translate Ising spins into dimers that can be reduced to some Pfaffian . Stephenson has used the Pfaffian representation to solve the Ising model on the pt lattice, but the solution was restricted to $`6L\times 6L`$ lattice due to its $`6\times 6`$ basic nonvanishing matrix elements and was exact only in the limit of $`L\mathrm{}`$ . Recently, by using the connections between Pfaffian, dimer and Ising model, Nash and O’Connor have obtained the exact expression of the partition function of the pt lattice Ising model on a finite torus . They first employed the lattice field theories to obtain the exact partition function of the Gaussian model, and then established the exact expression of the partition function of the pt lattice Ising model from the analysis of the appropriate lattice determinants and the parameterization according to the results in .

On the other hand, in view of the simplifying the approach, a remarkable progress was achieved by Plechko who modified the traditional fermionic interpretation and introduced a nonstandard approach . By the use of this approach, Plechko himself has not only rederived Onsager’s and Kaufman’s results in a relatively simple way , but also obtained the partition functions of a class of triangular type decorated lattices , and a triangular lattice net with holes . Quite recently, by using the same approach, Wu at al. have obtained the $`M\times N`$ sq lattice Ising model with periodic-aperiodic boundary condition , and Liaw at al. have successfully solved triangular and hexagonal lattices on a cylinder geometry ($`M\times \mathrm{}`$) with periodic and antiperiodic boundary condition . This approach is based on the integration over the anticommuting Grassmann variables and the mirror-ordered factorization principle in two-dimensional density matrix , and does not involve the traditional transfer-matrix or combinatorial considerations. The whole scheme of the method can be illustrated schematically as shown below :

$$Z=\underset{\left(\sigma \right)}{Sp}\left\{Z\left(\sigma \right)\right\}\underset{\left(\sigma |\chi \right)}{Sp}\left\{Z\left(\sigma |\chi \right)\right\}\underset{\left(\chi \right)}{Sp}\left\{Z\left(\chi \right)\right\}=Z,$$

where $`{}_{}{}^{\mathrm{}\mathrm{}}\mathrm{Sp}_{}^{\prime \prime }`$ stands for the average over spin variables ($`\sigma `$) or Grassmann variables ($`\chi `$). The original partition function $`Z`$ is expressed purely by spin variables ($`\sigma `$) at each lattice site. With a set of anticommuting Grassmann variables ($`\chi `$) being introduced to factorize the local bond Boltzmann weight such that spin variables are decoupled, the partition function passes to a mixed $`Z\left(\sigma |\chi \right)`$ representation. Then, by eliminating the spin variables in the mixed $`Z\left(\sigma |\chi \right)`$ representation, the fermionic interpretation $`Z\left(\chi \right)`$ of the two-dimensional Ising model can be obtained, and after carrying out the Grassmann integral, the analytic solution for the partition function and free energy can be achieved .

In the present paper, we work in this framework to obtain exact partition functions of $`M\times N`$ pt and hc lattices with different boundary conditions, including pp, pa, ap and aa boundary conditions. We used these results to calculate and plot the specific heat, $`C/k_B`$, as a function of the temperature, $`\theta =k_BT/J`$. Our results show that for the sq lattice, $`C/k_B`$ for pa and ap boundary conditions are different from those for aa boundary conditions, but for the pt and hc lattices, $`C/k_B`$ for ap, pa, and aa boundary conditions have the same values. Beside these analyses, our exact partition functions may also be used for understanding the effects of lattice structures and boundary conditions on critical properties and critical finite-size corrections of the Ising model.

This paper is organized as follows. In section 2, we set up a general form of the partition function for pt and hc lattices. Then, three pairs of conjugate Grassmann variables are introduced for a lattice site to factorize the Boltzmann weights, and the principle of mirror ordering are used to rearrange the Grassmann factors so we can perform the summation over Ising spins to obtain a pure fermionic expression of the partition function. In section 3, using the Fourier transform technique we complete the integrations over the Grassmann variables to obtain the exact solution of the partition function. Then, the solution is subjected to periodic-aperiodic boundary conditions, including pp, pa, ap and aa boundary conditions. We further consider the shift behaviours of the maximum of the specific heats of these systems in section 4. Finally, we discuss some problems for further studies in section 5.

## 2 The Partition Function

Consider Ising ferromagnets on $`M\times N`$ pt and hc lattices as shown in figure 1, in which the former is considered as a sq lattice with a single second-neighbor interaction, and the latter contains an inner spin in each lattice cell. The corresponding Hamiltonians, respectively, read as

$$H_t=\underset{m=1}{\overset{M}{}}\underset{n=1}{\overset{N}{}}\left(J_1\sigma _{mn}\sigma _{m+1n}+J_2\sigma _{mn}\sigma _{mn+1}+J_3\sigma _{m+1n}\sigma _{mn+1}\right),$$
(1)
and

$$H_h=\underset{m=1}{\overset{M}{}}\underset{n=1}{\overset{N}{}}\left(J_1\sigma _0\sigma _{mn}+J_2\sigma _0\sigma _{mn+1}+J_3\sigma _0\sigma _{m+1n}\right),$$
(2)
where $`J_i`$ with $`i=1,2,3`$ are the coupling constants ($`J_i>0`$ for ferromagnetic lattices), $`\sigma _{mn}=\pm 1`$ is the Ising spin located at the site $`(m,n)`$, and $`\sigma _0`$ denotes the inner Ising spin in hc lattice. Using the identity of the Boltzmann weight,

$$\mathrm{exp}\left(\beta J_i\sigma _\mu \sigma _\nu \right)=\mathrm{cosh}\left(\beta J_i\right)\left[1+\mathrm{tanh}\left(\beta J_i\right)\sigma _\mu \sigma _\nu \right],$$
(3)
$`\beta =\left(k_BT\right)^1`$, and performing the sum over $`\sigma _0`$, the partition functions of two lattices can be formulated in a single three spin-polynomial representation,

$`Z`$ $`=`$ $`2^{N_s}\left[{\displaystyle \underset{i=1}{\overset{n_b}{}}}\mathrm{cosh}\left(\beta J_i\right)\right]^{N_s}`$ (4)
$`\times \underset{\left(\sigma \right)}{Sp}\left\{{\displaystyle \underset{m=1}{\overset{M}{}}}{\displaystyle \underset{n=1}{\overset{N}{}}}(\alpha _0+\alpha _1\sigma _{mn}\sigma _{m+1n}+\alpha _2\sigma _{mn}\sigma _{mn+1}+\alpha _3\sigma _{m+1n}\sigma _{mn+1})\right\},`$
where $`N_s`$ is the numbers of lattice sites ($`N_s=MN`$ for sq and pt lattice, $`N_s=2MN`$ for hc lattice) and $`n_b`$ is the number of bonds per lattice cell ($`n_b=2`$ for sq lattice, $`n_b=3`$ for pt and hc lattice), symbol $`{}_{}{}^{\mathrm{}\mathrm{}}\underset{\left(\sigma \right)}{Sp}`$ ′′ stands for spin average defined by

$$\underset{\left(\sigma _i\right)}{Sp}\left[\mathrm{}\right]=\frac{1}{2}\underset{\left(\sigma _i=\pm 1\right)}{}\left[\mathrm{}\right],\underset{\left(\sigma _i\right)}{Sp}\left[1\right]=1,\underset{\left(\sigma _i\right)}{Sp}\left[\sigma _i\right]=0$$
(5)
and $`\alpha _i`$’s are defined as

$$\alpha _0^T=1+t_1t_2t_3,\alpha _1^T=t_1+t_2t_3,\alpha _2^T=t_2+t_3t_1,\alpha _3^T=t_3+t_1t_2,$$
(6)
$`t_i=\mathrm{tanh}\left(\beta J_i\right)`$ with $`i=1,2,3`$, for pt lattice, and

$$\alpha _0^H=1,\alpha _1^H=t_1t_3,\alpha _2^H=t_1t_2,\alpha _3^H=t_2t_3,$$
(7)
for hc lattice.

To factorize the partition, we rewrite the partition function as

$`Z_H`$ $`=`$ $`2^{N_s}\left[{\displaystyle \underset{i=1}{\overset{n_b}{}}}\mathrm{cosh}\left(\beta J_i\right)\right]^{N_s}`$ (8)
$`\times \underset{\left(\sigma \right)}{Sp}\left\{{\displaystyle \underset{m=1}{\overset{M}{}}}{\displaystyle \underset{n=1}{\overset{N}{}}}r_0(1+r_1\sigma _{mn}\sigma _{m+1n})(1+r_2\sigma _{mn}\sigma _{mn+1})(1+r_3\sigma _{m+1n}\sigma _{mn+1})\right\},`$
where $`r_i`$ with $`i=0,1,2,3`$ vary from one lattice to the other, and are related to $`\alpha _i`$’s from

$$\alpha _0=r_0\left(1+r_1r_2r_3\right),\alpha _1=r_0\left(r_1+r_2r_3\right),\alpha _2=r_0\left(r_2+r_1r_3\right),\alpha _3=r_0\left(r_3+r_1r_2\right).$$
(9)
For pt lattice, the relation between $`r_i`$ and $`t_i`$ is trivial, i.e. $`r_0=1`$ and $`r_i=t_i`$, but for hc lattice, the relation is nontrivial and is determined by equations (7) and (9).

It is more convenient to define the generalized reduced partition function as

$$Q=r_0^{MN}\stackrel{~}{Q},$$
(10)
with

$$\stackrel{~}{Q}=\underset{m=1}{\overset{M}{}}\underset{n=1}{\overset{N}{}}\underset{\left(\sigma _{mn}\right)}{Sp}\left[\left(1+r_1\sigma _{mn}\sigma _{m+1n}\right)\left(1+r_2\sigma _{mn}\sigma _{mn+1}\right)\left(1+r_3\sigma _{mn+1}\sigma _{m+1n}\right)\right].$$
(11)
To construct the fermionic representation of the generalized partition function, we associate each lattice site $`(m,n)`$ with three pairs of conjugate Grassmann variables, $`\{a_{mn},a_{mn}^{};b_{mn},b_{mn}^{};c_{mn},c_{mn}^{}\}\chi `$. All of these Grassmann variables are anticommuting, and their square are zero. Their integral obeys the basic rules

$$𝑑\chi =0,𝑑\chi \chi =1,$$
(12)
$$𝑑\chi \mathrm{\Omega }\left(\chi +\eta \right)=𝑑\chi \mathrm{\Omega }\left(\chi \right),$$
(13)
for an arbitrary vector $`\eta `$ with anticommuting components, and there is the relation

$$1+r_i\sigma _\mu \sigma _\nu =𝑑\chi ^{}𝑑\chi e^{\chi \chi ^{}}\left(1+\chi \sigma _\mu \right)\left(1+r_i\chi ^{}\sigma _\nu \right).$$
(14)
Using these Grassmann variables, we can rewrite the reduced partition function as

$$\stackrel{~}{Q}=\underset{m=1}{\overset{M}{}}\underset{n=1}{\overset{N}{}}\underset{\left(\sigma _{mn}\right)}{Sp}\left[\underset{(a_{mn},b_{mn},c_{mn})}{Sp}\left(A_{mn}A_{m+1n}^{}B_{mn}B_{mn+1}^{}C_{mn+1}C_{m+1n}^{}\right)\right],$$
(15)
where $`{}_{}{}^{\mathrm{}\mathrm{}}\underset{\left(\chi _i\right)}{Sp}`$ ′′ stands for the averaging with Gaussian weight

$`\underset{\left(\chi _i\right)}{Sp}\left[\mathrm{}\right]={\displaystyle 𝑑\chi _i^{}𝑑\chi _ie^{\chi _i\chi _i^{}}\left[\mathrm{}\right]},`$ (16)
with the rules

$`\underset{\left(\chi _i\right)}{Sp}\left[\chi _i\chi _i^{}\right]`$ $`=`$ $`\underset{\left(\chi _i\right)}{Sp}\left[\chi _i^{}\chi _i\right]=1,`$ (17)
$`\underset{\left(\chi _i\right)}{Sp}\left[\chi _i\right]`$ $`=`$ $`\underset{\left(\chi _i\right)}{Sp}\left[\chi _i^{}\right]=0,`$ (18)
and the Grassmann factors, $`A,A^{},B,B^{},C,`$ and $`C^{}`$, are defined as

$$A_{mn}=1+a_{mn}\sigma _{mn},A_{mn}^{}=1+r_1a_{m1n}^{}\sigma _{mn};$$
(19)
$$B_{mn}=1+b_{mn}\sigma _{mn},B_{mn}^{}=1+r_2b_{mn1}^{}\sigma _{mn};$$
(20)
$$C_{mn}=1+c_{mn1}\sigma _{mn},C_{mn}^{}=1+r_3c_{m1n}^{}\sigma _{mn}.$$
(21)
In this way, a Boltzmann weight is decoupled to the product of two factors of separated spins.

For simplicity, we express the reduced partition function as

$$\stackrel{~}{Q}=\underset{(a,b,c)}{Sp}\left\{\underset{m=1}{\overset{M}{}}\underset{n=1}{\overset{N}{}}\mathrm{\Psi }_{mn}^A\mathrm{\Psi }_{mn}^B\mathrm{\Psi }_{mn}^C\right\},$$
(22)
where $`\mathrm{\Psi }_{mn}^A`$, $`\mathrm{\Psi }_{mn}^B`$ and $`\mathrm{\Psi }_{mn}^C`$ are defined by

$`\mathrm{\Psi }_{mn}^A`$ $`=`$ $`\underset{\left(\sigma _{mn}\right)}{Sp}\left(A_{mn}A_{m+1n}^{}\right),`$ (23)
$`\mathrm{\Psi }_{mn}^B`$ $`=`$ $`\underset{\left(\sigma _{mn}\right)}{Sp}\left(B_{mn}B_{mn+1}^{}\right),`$ (24)
$`\mathrm{\Psi }_{mn}^C`$ $`=`$ $`\underset{\left(\sigma _{mn}\right)}{Sp}\left(C_{mn+1}C_{m+1n}^{}\right).`$ (25)
We first treat the boundary weight and consider periodic boundary condition in both directions:

$`\mathrm{\Psi }_{Mn}^A`$ $`=`$ $`\underset{\left(\sigma _{Mn}\right)}{Sp}\left[\left(1+a_{Mn}\sigma _{Mn}\right)\left(1+r_1a_{Mn}^{}\sigma _{M+1n}\right)\right]`$ (26)
$`=`$ $`\underset{\left(\sigma _{Mn}\right)}{Sp}\left[\left(1+r_1a_{0n}^{}\sigma _{1n}\right)\left(1+a_{Mn}\sigma _{Mn}\right)\right]`$
$`=`$ $`\underset{\left(\sigma _{Mn}\right)}{Sp}\left(A_{1n}^{}A_{Mn}\right),`$
which implies

$$a_{0n}^{}=a_{Mn}^{}.$$
(27)
Similarly, from

$`\mathrm{\Psi }_{mN}^B`$ $`=`$ $`\underset{\left(\sigma _{mN}\right)}{Sp}\left(B_{mN}B_{mN+1}^{}\right)=\underset{\left(\sigma _{Nn}\right)}{Sp}\left(B_{m1}^{}B_{mN}\right),`$ (28)
$`\mathrm{\Psi }_{Mn}^C`$ $`=`$ $`\underset{\left(\sigma _{Mn}\right)}{Sp}\left(C_{Mn+1}C_{M+1n}^{}\right)=\underset{\left(\sigma _{Mn}\right)}{Sp}\left(C_{1n}^{}C_{Mn+1}\right),`$ (29)
$`\mathrm{\Psi }_{mN}^C`$ $`=`$ $`\underset{\left(\sigma _{mN}\right)}{Sp}\left(C_{mN+1}C_{m+1N}^{}\right)=\underset{\left(\sigma _{mN}\right)}{Sp}\left(C_{m1}C_{m+1N}^{}\right),`$ (30)
we have

$`b_{m0}^{}`$ $`=`$ $`b_{mN}^{},`$ (31)
$`c_{0n}^{}`$ $`=`$ $`c_{Mn}^{},`$ (32)
$`c_{m0}`$ $`=`$ $`c_{mN}.`$ (33)
Since $`c_{m0}=c_{mN}`$, $`\mathrm{\Psi }_{mN}^C`$ need not to be treated as a boundary weight, and only $`\mathrm{\Psi }_{Mn}^A`$, $`\mathrm{\Psi }_{mN}^B`$ and $`\mathrm{\Psi }_{Mn}^C`$ should be considered. However, this situation becomes ambiguous when we take Fourier transform of these Grassmann variables with single set of exponential factors in equations (56) and (57). Because the Fourier exponential factors are associated with directions in $`M`$ and $`N`$, the sign factor in front of $`b_{mN}^{}`$ takes effects simultaneously on $`b_{mN}^{}`$ and $`c_{mN}`$. Therefore, the real situation is that instead of the relation in equation (33), we must take

$$c_{m0}=c_{mN}.$$
(34)
A self-consistent way to assign a minus sign to $`c_{m0}`$ and obtain the relation in equation (34) is interchanging $`C_{m1}`$ in equation (30) with another Grassmann factor. An equivalent but more convenient approach is to consider the rearrangement of $`B_{m1}^{}`$ in the boundary weight together with the rearrangement of $`C_{m1}`$ in the reduced partition function. To see this, we express the reduced partition function as

$$\stackrel{~}{Q}=\underset{(a,b,c)}{Sp}\left\{\underset{\left(\sigma \right)}{Sp}\left[\left(\underset{m=1}{\overset{M1}{}}\underset{n=1}{\overset{N}{}}\mathrm{\Psi }_{mn}^A\mathrm{\Psi }_{mn}^C\right)\mathrm{\Psi }_B\left(\underset{m=1}{\overset{M}{}}\underset{n=1}{\overset{N1}{}}B_{mn}B_{mn+1}^{}\right)\right]\right\}.$$
(35)
with the boundary weight $`\mathrm{\Psi }_B`$

$`\mathrm{\Psi }_B`$ $`=`$ $`\underset{(a,b,c)}{Sp}\left[\left({\displaystyle \underset{n=1}{\overset{N}{}}}\mathrm{\Psi }_{Mn}^A\right)\left({\displaystyle \underset{m=1}{\overset{M}{}}}\mathrm{\Psi }_{mN}^B\right)\left({\displaystyle \underset{n=1}{\overset{N}{}}}\mathrm{\Psi }_{Mn}^C\right)\right]`$ (36)
$`=`$ $`\underset{(a,b,c)}{Sp}\left[\left({\displaystyle \underset{m=1}{\overset{M}{}}}\stackrel{m}{\stackrel{}{B_{m1}^{}}}\right)\left({\displaystyle \underset{n=1}{\overset{N}{}}}\stackrel{n}{\stackrel{}{C_{1n}^{}A_{1n}^{}}}\right)A_{M1}\left({\displaystyle \underset{n=2}{\overset{N}{}}}\stackrel{n}{\stackrel{}{C_{Mn}A_{Mn}}}\right)C_{M1}{\displaystyle \underset{m=1}{\overset{M}{}}}\stackrel{m}{\stackrel{}{B_{mN}}}\right],`$
and

$`{\displaystyle \underset{m=1}{\overset{M1}{}}}{\displaystyle \underset{n=1}{\overset{N}{}}}\mathrm{\Psi }_{mn}^A\mathrm{\Psi }_{mn}^C`$ $`=`$ $`{\displaystyle \underset{m=1}{\overset{M1}{}}}{\displaystyle \underset{n=1}{\overset{N}{}}}A_{mn}C_{mn+1}C_{m+1n}^{}A_{m+1n}^{}`$ (37)
$`=`$ $`{\displaystyle \underset{m=1}{\overset{M1}{}}}A_{m1}\left({\displaystyle \underset{n=2}{\overset{N}{}}}\stackrel{n}{\stackrel{}{C_{mn}A_{mn}}}\right)C_{m1}\left({\displaystyle \underset{n=1}{\overset{N}{}}}\stackrel{n}{\stackrel{}{C_{m+1n}^{}A_{m+1n}^{}}}\right).`$
Here, arrows have been used to indicate the orders of the products in $`m`$ and $`n`$. When we move $`B_{m1}^{}`$ from left of $`\underset{n=1}{\overset{N}{}}\stackrel{n}{\stackrel{}{C_{1n}^{}A_{1n}^{}}}`$ to right of $`\underset{n=1}{\overset{N}{}}\stackrel{n}{\stackrel{}{C_{1n}^{}A_{1n}^{}}}`$ in equation (36), $`B_{m1}^{}`$ passes $`2N`$ Grassmann factors, but for moving $`C_{m1}`$ from right of $`\underset{n=2}{\overset{N}{}}\stackrel{n}{\stackrel{}{C_{mn}A_{mn}}}`$ to left of $`A_{m1}`$ in equations (36) and (37), $`C_{m1}`$ passes only $`2N1`$ Grassmann factors. Then by moving $`B_{m1}^{}`$ from left to right, and simultaneously moving $`C_{m1}`$ from right to left, we can assign to the Grassmann variable in $`C_{m1}`$ an additional minus sign comparing with the Grassmann variable in $`B_{m1}^{}`$, and hence obtain the relation of equation (34).

Accordingly, we interchange $`B^{}`$ and $`C^{}A^{}`$ in equation (36) to obtain the arrangement of $`C^{}A^{}B^{}`$ according to the identity

$$B^+\left(CA\right)^+=\frac{1}{2}\left[\left(CA\right)^+B^++\left(CA\right)^+B^{}+\left(CA\right)^{}B^+\left(CA\right)^{}B^{}\right],$$
(38)
with superscripts $`+`$ and $``$ being the sign factors in boundary Grassmann factors $`A_{1n}^{}`$, $`B_{m1}^{}`$ and $`C_{1n}^{}`$, and simultaneously move $`C_{m1}`$ from right of $`\underset{n=2}{\overset{N}{}}\stackrel{n}{\stackrel{}{C_{mn}A_{mn}}}`$ to left of $`A_{m1}`$ in equations (36) and (37). Here we note that the superscripts $`+`$ and $``$ respectively correspond to periodic and antiperiodic boundary condition imposed on the spin variables and in turn on the Grassmann variables. Hence, the reduced partition function becomes

$$\stackrel{~}{Q}=\frac{1}{2}\left(\stackrel{~}{Q}_\gamma |_{\mathrm{\Gamma }_1}+\stackrel{~}{Q}_\gamma |_{\mathrm{\Gamma }_2}+\stackrel{~}{Q}_\gamma |_{\mathrm{\Gamma }_3}\stackrel{~}{Q}_\gamma |_{\mathrm{\Gamma }_4}\right),$$
(39)
with

$$\stackrel{~}{Q}_\gamma =\underset{(a,b,c)}{Sp}\left\{\underset{\left(\sigma \right)}{Sp}[\left(\underset{m=1}{\overset{M1}{}}\stackrel{m}{\stackrel{}{\mathrm{\Theta }_m\mathrm{\Theta }_{m+1}^{}}}\right)\mathrm{\Psi }_\gamma \left(\underset{m=1}{\overset{M}{}}\underset{n=1}{\overset{N1}{}}B_{mn}B_{mn+1}^{}\right)]\right\},$$
(40)
and

$$\mathrm{\Psi }_\gamma =\underset{(a,b,c)}{Sp}\left\{\underset{\left(\sigma \right)}{Sp}\left[\mathrm{\Theta }_1^{}\left(\underset{m=1}{\overset{M}{}}\stackrel{m}{\stackrel{}{B_{m1}^{}}}\right)\mathrm{\Theta }_M\left(\underset{m=1}{\overset{M}{}}\stackrel{m}{\stackrel{}{B_{mN}}}\right)\right]\right\}.$$
(41)
where we have defined

$$\mathrm{\Theta }_m=\underset{n=1}{\overset{N}{}}\stackrel{n}{\stackrel{}{C_{mn}A_{mn}}}\mathrm{and}\mathrm{\Theta }_m^{}=\underset{n=1}{\overset{N}{}}\stackrel{n}{\stackrel{}{C_{mn}^{}A_{mn}^{}}},$$
(42)
and the boundary conditions $`\mathrm{\Gamma }_1`$, $`\mathrm{\Gamma }_2`$, $`\mathrm{\Gamma }_3`$, $`\mathrm{\Gamma }_4`$ are defined as

$`\mathrm{\Gamma }_1`$ $`=`$ $`\left(a_{0n}^{}=a_{Mn}^{},b_{m0}^{}=b_{mN}^{},c_{0n}^{}=c_{Mn}^{}\right),`$ (43)
$`\mathrm{\Gamma }_2`$ $`=`$ $`\left(a_{0n}^{}=a_{Mn}^{},b_{m0}^{}=+b_{mN}^{},c_{0n}^{}=c_{Mn}^{}\right),`$ (44)
$`\mathrm{\Gamma }_3`$ $`=`$ $`\left(a_{0n}^{}=+a_{M,n}^{},b_{m0}^{}=b_{mN}^{},c_{0n}^{}=+c_{M,n}^{}\right),`$ (45)
$`\mathrm{\Gamma }_4`$ $`=`$ $`\left(a_{0n}^{}=+a_{M,n}^{},b_{m0}^{}=+b_{mN}^{},c_{0n}^{}=+c_{M,n}^{}\right).`$ (46)
In this way, the configurations of the reduced partition function can be further rearrangement and expressed as

$`\stackrel{~}{Q}_\gamma `$ $`=`$ $`\underset{(a,b,c)}{Sp}\underset{\left(\sigma \right)}{Sp}\left\{\left({\displaystyle \underset{m=1}{\overset{M1}{}}}\stackrel{m}{\stackrel{}{\mathrm{\Theta }_m\mathrm{\Theta }_{m+1}^{}}}\right)\mathrm{\Theta }_1^{}\left({\displaystyle \underset{m=1}{\overset{M}{}}}\stackrel{m}{\stackrel{}{B_{m1}^{}}}\right)\mathrm{\Theta }_M\left({\displaystyle \underset{m=1}{\overset{M}{}}}\stackrel{n}{\stackrel{}{B_{mN}}}\right){\displaystyle \underset{m=1}{\overset{M}{}}}{\displaystyle \underset{n=1}{\overset{N1}{}}}B_{mn}B_{mn+1}^{}\right\}`$ (47)
$`=`$ $`\underset{(a,b,c)}{Sp}\underset{\left(\sigma \right)}{Sp}\left\{\left({\displaystyle \underset{m=1}{\overset{M}{}}}\stackrel{m}{\stackrel{}{\mathrm{\Theta }_m^{}B_{m1}^{}\mathrm{\Theta }_m}}\right)\left({\displaystyle \underset{m=1}{\overset{M}{}}}\stackrel{n}{\stackrel{}{B_{mN}}}\right)\left({\displaystyle \underset{m=1}{\overset{M}{}}}{\displaystyle \underset{n=1}{\overset{N1}{}}}B_{mn}B_{mn+1}^{}\right)\right\}.`$
To have a complete mirror-ordered form, we have to rearrange the terms in the last two brackets. To achieve this, first we note that

$`\stackrel{~}{Q}_\gamma `$ $`=`$ $`\underset{(a,b,c)}{Sp}\underset{\left(\sigma \right)}{Sp}\left\{{\displaystyle \underset{m=1}{\overset{M}{}}}\stackrel{m}{\stackrel{}{\mathrm{\Theta }_m^{}B_{m1}^{}\left({\displaystyle \underset{n=1}{\overset{N1}{}}}\stackrel{n}{\stackrel{}{C_{mn}A_{mn}}}\right)C_{mN}A_{mN}}}\left({\displaystyle \underset{m=1}{\overset{M}{}}}\stackrel{m}{\stackrel{}{B_{mM}}}\right){\displaystyle \underset{m=1}{\overset{M}{}}}{\displaystyle \underset{n=1}{\overset{N1}{}}}B_{mn}B_{mn+1}^{}\right\}`$ (48)
$`=`$ $`\underset{(a,b,c)}{Sp}\underset{\left(\sigma \right)}{Sp}\left\{{\displaystyle \underset{m=1}{\overset{M}{}}}\stackrel{m}{\stackrel{}{\mathrm{\Theta }_m^{}\left({\displaystyle \underset{n=1}{\overset{N1}{}}}\stackrel{n}{\stackrel{}{B_{mn}^{}C_{mn}A_{mn}B_{mn}}}\right)B_{mN}^{}C_{mN}A_{mN}}}\left({\displaystyle \underset{m=1}{\overset{M}{}}}\stackrel{m}{\stackrel{}{B_{mN}}}\right)\right\}.`$
The boundary term of with $`m=M`$, denoted by $`T`$, can be formulated

$`T`$ $`=`$ $`\mathrm{\Theta }_M^{}\left({\displaystyle \underset{n=1}{\overset{N1}{}}}\stackrel{n}{\stackrel{}{B_{Mn}^{}C_{Mn}A_{Mn}B_{Mn}}}\right)B_{MN}^{}C_{MN}A_{MN}B_{MN}`$
$`=`$ $`\left({\displaystyle \underset{n=1}{\overset{N}{}}}\stackrel{n}{\stackrel{}{C_{Mn}^{}A_{Mn}^{}}}\right)\left({\displaystyle \underset{n=1}{\overset{L_y}{}}}\stackrel{n}{\stackrel{}{B_{Mn}^{}C_{Mn}A_{Mn}B_{Mn}}}\right)`$
$`=`$ $`{\displaystyle \underset{n=1}{\overset{N}{}}}C_{Mn}^{}A_{Mn}^{}B_{Mn}^{}C_{Mn}A_{Mn}B_{Mn},`$ (49)
due to the fact that $`\underset{\left(\sigma _{mn}\right)}{Sp}\left[C_{Mn}^{}A_{Mn}^{}B_{Mn}^{}C_{Mn}A_{Mn}B_{Mn}\right]`$ for a given $`n`$ is a commutable object. By continuing such construction from $`m=M`$ down to $`m=1`$, we can obtain the expression

$$\stackrel{~}{Q}_\gamma =\underset{(a,b,c)}{Sp}\left\{\underset{m=1}{\overset{M}{}}\underset{n=1}{\overset{N}{}}\underset{\left(\sigma _{mn}\right)}{Sp}\left[C_{mn}^{}A_{mn}^{}B_{mn}^{}C_{mn}A_{mn}B_{mn}\right]\right\}.$$
(50)
For this partition function, the factors containing the same spin are grouped together and we can perform the average over spins. As a result, we have

$$\stackrel{~}{Q}_\gamma =\underset{m=1}{\overset{M}{}}\underset{n=1}{\overset{N}{}}da_{mn}^{}da_{mn}db_{mn}^{}db_{mn}dc_{mn}^{}dc_{mn}\mathrm{exp}\left(\underset{m=1}{\overset{M}{}}\underset{n=1}{\overset{N}{}}F_{mn}\right),$$
(51)
with

$`F_{mn}`$ $`=`$ $`a_{mn}a_{mn}^{}+b_{mn}b_{mn}^{}+c_{mn}c_{mn}^{}`$ (52)
$`+r_1r_3c_{m1n}^{}a_{m1n}^{}`$
$`+\left(r_3c_{m1n}^{}+r_1a_{m1n}^{}\right)r_2b_{mn1}^{}`$
$`+\left(r_3c_{m1n}^{}+r_1a_{m1n}^{}+r_2b_{mn1}^{}\right)c_{mn1}`$
$`+\left(r_3c_{m1n}^{}+r_1a_{m1n}^{}+r_2b_{mn1}^{}+c_{mn1}\right)a_{mn}`$
$`+\left(r_3c_{m1n}^{}+r_1a_{m1n}^{}+r_2b_{mn1}^{}+c_{mn1}+a_{mn}\right)b_{mn}.`$
Since there is no mix on $`a_{mn}`$ and $`b_{mn}`$, the integral in the above expression can be simplified by integrating out the $`a_{mn}`$ and $`b_{mn}`$ fields by means of the identity

$$𝑑b𝑑a\mathrm{exp}\left(\lambda ab+aL+L^{}b\right)=\lambda \mathrm{exp}\left(\lambda ^1LL^{}\right),$$
(53)
where $`a`$, $`b`$ are Grassmann variables, $`L`$, $`L^{}`$ are linear fermionic forms independent of $`a`$, $`b`$ and $`\lambda `$ is a parameter. The result then becomes

$$\stackrel{~}{Q}_\gamma =\underset{m=1}{\overset{M}{}}\underset{n=1}{\overset{N}{}}dg_{mn}^{}dg_{mn}dc_{mn}^{}dc_{mn}\mathrm{exp}\left(\underset{m=1}{\overset{M}{}}\underset{n=1}{\overset{N}{}}G_{mn}\right),$$
(54)
with

$`G_{mn}`$ $`=`$ $`c_{mn}c_{mn}^{}+g_{mn}g_{mn}^{}`$ (55)
$`+r_1r_3c_{m1n}^{}g_{m1n}`$
$`\left(r_3c_{m1n}^{}+r_1g_{m1n}\right)r_2g_{mn1}^{}`$
$`+\left(r_3c_{m1n}^{}+r_1g_{m1n}r_2g_{mn1}^{}\right)c_{mn1}`$
$`\left(r_3c_{m1n}^{}+r_1g_{m1n}r_2g_{mn1}^{}+c_{mn1}\right)\left(g_{mn}+g_{mn}^{}\right),`$
where we have changed the notations for the fields by $`(a_{mn}^{},b_{mn}^{})(g_{mn},g_{mn}^{})`$. This is the pure fermionic representation of the reduced partition function.

## 3 Exact Solution

Next, to carry out the integration, we have to use the technique of Fourier transform to treat the Grassmann variables which mix together with the variables at different sites. The Fourier transformation is defined as

$$X_{mn}=\frac{1}{\sqrt{MN}}\underset{p=0}{\overset{M1}{}}\underset{q=0}{\overset{N1}{}}X_{pq}e^{i\frac{2\pi }{M}mp}e^{i\frac{2\pi }{N}nq},$$
(56)
and

$$X_{mn}^{}=\frac{1}{\sqrt{MN}}\underset{p=0}{\overset{M1}{}}\underset{q=0}{\overset{N1}{}}X_{pq}^{}e^{i\frac{2\pi }{M}mp}e^{i\frac{2\pi }{N}nq},$$
(57)
where the variables $`X_{mn}`$ and $`X_{mn}^{}`$ denotes one of the variables $`\{c_{mn},g_{mn}\}`$ and $`\{c_{mn}^{},g_{mn}^{}\}`$ respectively.

After performing the Fourier transformation, the partition function becomes

$$\stackrel{~}{Q}_\gamma =\underset{p=0}{\overset{M1}{}}\underset{q=0}{\overset{N1}{}}𝑑V_{pq}\mathrm{exp}\left(H_{pq}\right),$$
(58)
with the measure $`dV_{pq}`$ defined as

$$dV_{pq}=dg_{pq}^{}dg_{pq}dc_{pq}^{}dc_{pq},$$
(59)
and the function $`H_{pq}`$ is given by

$`H_{pq}`$ $`=`$ $`\left(1r_3e^{i\frac{2\pi }{M}p}e^{i\frac{2\pi }{N}q}\right)c_{pq}c_{pq}^{}`$ (60)
$`+\left(r_2e^{i\frac{2\pi }{N}q}\right)c_{pq}g_{pq}^{}+r_3\left(r_1e^{i\frac{2\pi }{M}p}\right)c_{pq}^{}g_{pq}`$
$`e^{i\frac{2\pi }{N}q}\left(1+r_1e^{i\frac{2\pi }{M}p}\right)c_{pq}g_{MpNq}r_3e^{i\frac{2\pi }{M}p}\left(1+r_2e^{i\frac{2\pi }{N}q}\right)c_{pq}^{}g_{MpNq}^{}`$
$`+\left(1r_1e^{i\frac{2\pi }{M}p}r_2e^{i\frac{2\pi }{N}q}r_1r_2e^{i\frac{2\pi }{M}p}e^{i\frac{2\pi }{N}q}\right)g_{pq}g_{pq}^{}`$
$`r_1e^{i\frac{2\pi }{M}p}g_{pq}g_{MpNq}+r_2e^{i\frac{2\pi }{L_y}q}g_{pq}^{}g_{MpNq}^{},`$
Because $`H_{pq}`$ contains not only the variables, $`X_{pq}`$ and $`X_{pq}^{}`$, but also the variables, $`X_{MpNq}`$ and $`X_{MpNq}^{}`$, instead of calculating $`\stackrel{~}{Q}_\gamma `$ it is easier to calculate $`\stackrel{~}{Q}_\gamma ^2`$ given by

$$\stackrel{~}{Q}_\gamma ^2=\underset{p=0}{\overset{M1}{}}\underset{q=0}{\overset{N1}{}}𝑑V_{pq}𝑑V_{MpNq}\mathrm{exp}\left(H_{pq}+H_{MpNq}^{}\right).$$
(61)
Here $`H_{MpNq}^{}`$ can be obtained from $`H_{pq}`$ by replacing $`p`$ by $`Mp`$ and $`q`$ by $`Nq`$ for the Grassmann variables and replacing the coefficient in front of Grassmann variables by its complex conjugate. Completing the integration yields

$$Q_\gamma =\underset{p=0}{\overset{M1}{}}\underset{q=0}{\overset{N1}{}}\left[A_0A_1\mathrm{cos}\frac{2\pi p}{M}A_2\mathrm{cos}\frac{2\pi q}{N}A_3\mathrm{cos}\left(\frac{2\pi p}{M}\frac{2\pi q}{N}\right)\right]^{1/2},$$
(62)
with

$`A_0`$ $`=`$ $`\alpha _0^2+\alpha _1^2+\alpha _2^2+\alpha _3^2,`$ (63)
$`A_1`$ $`=`$ $`2\left(\alpha _0\alpha _1\alpha _2\alpha _3\right),`$ (64)
$`A_2`$ $`=`$ $`2\left(\alpha _0\alpha _2\alpha _1\alpha _3\right),`$ (65)
$`A_3`$ $`=`$ $`2\left(\alpha _0\alpha _3\alpha _1\alpha _2\right),`$ (66)
where $`\alpha _0`$, $`\alpha _1`$, $`\alpha _2`$, and $`\alpha _3`$ are given by equations (6) and (7) for pt and hc lattices, respectively.

### 3.1 Periodic-periodic boundary condition

According to equation (39), the reduced partition function for ferromagnetic lattices with pp boundary condition is

$$Q^{pp}=\frac{1}{2}\left[\mathrm{\Omega }_{\frac{1}{2},\frac{1}{2}}+\mathrm{\Omega }_{\frac{1}{2},0}+\mathrm{\Omega }_{0,\frac{1}{2}}\mathrm{sgn}\left(\frac{\theta \theta _c}{\theta _c}\right)\mathrm{\Omega }_{00}\right],$$
(67)
where the superscript $`p`$ refers to periodic boundary condition and

$$\mathrm{\Omega }_{\mu \nu }=\underset{p=0}{\overset{M1}{}}\underset{q=0}{\overset{N1}{}}\left[A_0A_1\mathrm{cos}\frac{2\pi \left(p+\mu \right)}{M}A_2\mathrm{cos}\frac{2\pi \left(q+\nu \right)}{N}A_3\mathrm{cos}\left(\frac{2\pi \left(p+\mu \right)}{M}\frac{2\pi \left(q+\nu \right)}{N}\right)\right]^{1/2}.$$
(68)
The sign factor in front of the last term is a result of the standard consideration of the Grassmann integral over the zero-mode variable $`p=q=0`$ for ferromagnetic couplings . When the integral of equation (61) is carried out, it is always positive, but this is not the case for equation (58). There are unpaired terms from zero-mode in equation (58) under various boundary conditions and they contribute a sign factor to $`Q_4`$ for $`0t_i1`$. The partition function for pp boundary condition then becomes

$$Z^{pp}=\frac{1}{2}2^{N_s}\left[\underset{i=1}{\overset{n_b}{}}\mathrm{cosh}\left(\beta J_i\right)\right]^{N_s}\left[\mathrm{\Omega }_{\frac{1}{2},\frac{1}{2}}+\mathrm{\Omega }_{\frac{1}{2},0}+\mathrm{\Omega }_{0,\frac{1}{2}}\mathrm{sgn}\left(\frac{\theta \theta _c}{\theta _c}\right)\mathrm{\Omega }_{00}\right].$$
(69)
Furthermore, the free energy density per $`k_BT`$ of the system defined by

$$f^{pp}=\frac{1}{N_s}\mathrm{ln}Z^{pp},$$
(70)
then takes the form

$`f^{pp}`$ $`=`$ $`{\displaystyle \frac{\left(N_s1\right)}{N_s}}\mathrm{ln}2{\displaystyle \underset{i=1}{\overset{n_b}{}}}\mathrm{ln}\left[\mathrm{cosh}\left(\beta J_i\right)\right]`$ (71)
$`{\displaystyle \frac{1}{N_s}}\mathrm{ln}\left[\mathrm{\Omega }_{\frac{1}{2},\frac{1}{2}}+\mathrm{\Omega }_{\frac{1}{2},0}+\mathrm{\Omega }_{0,\frac{1}{2}}\mathrm{sgn}\left({\displaystyle \frac{\theta \theta _c}{\theta _c}}\right)\mathrm{\Omega }_{00}\right].`$
### 3.2 Periodic-antiperiodic boundary condition

For pa boundary condition, equation (38) is replaced by

$$B^{}\left(CA\right)^+=\frac{1}{2}\left[\left(CA\right)^+B^++\left(CA\right)^+B^{}\left(CA\right)^{}B^++\left(CA\right)^{}B^{}\right],$$
(72)
and the partition function has the form

$$Z^{pa}=\frac{1}{2}2^{N_s}\left[\underset{i=1}{\overset{n_b}{}}\mathrm{cosh}\left(\beta J_i\right)\right]^{N_s}\left[\mathrm{\Omega }_{\frac{1}{2},\frac{1}{2}}+\mathrm{\Omega }_{\frac{1}{2},0}\mathrm{\Omega }_{0,\frac{1}{2}}+\mathrm{sgn}\left(\frac{\theta \theta _c}{\theta _c}\right)\mathrm{\Omega }_{00}\right],$$
(73)
where the superscript $`a`$ refers to antiperiodic boundary condition. The corresponding free energy density per $`k_BT`$ is

$`f^{pa}`$ $`=`$ $`{\displaystyle \frac{\left(N_s1\right)}{N_s}}\mathrm{ln}2{\displaystyle \underset{i=1}{\overset{n_b}{}}}\mathrm{ln}\left[\mathrm{cosh}\left(\beta J_i\right)\right]`$ (74)
$`{\displaystyle \frac{1}{N_s}}\mathrm{ln}\left[\mathrm{\Omega }_{\frac{1}{2},\frac{1}{2}}+\mathrm{\Omega }_{\frac{1}{2},0}\mathrm{\Omega }_{0,\frac{1}{2}}+\mathrm{sgn}\left({\displaystyle \frac{\theta \theta _c}{\theta _c}}\right)\mathrm{\Omega }_{00}\right].`$
### 3.3 Antiperiodic-periodic boundary condition

Similarly, for ap boundary condition, equation (38) is replaced by

$$B^+\left(CA\right)^{}=\frac{1}{2}\left[\left(CA\right)^+B^+\left(CA\right)^+B^{}+\left(CA\right)^{}B^++\left(CA\right)^{}B^{}\right],$$
(75)
and

$$Z^{ap}=\frac{1}{2}2^{N_s}\left[\underset{i=1}{\overset{n_b}{}}\mathrm{cosh}\left(\beta J_i\right)\right]^{N_s}\left[\mathrm{\Omega }_{\frac{1}{2},\frac{1}{2}}\mathrm{\Omega }_{\frac{1}{2},0}+\mathrm{\Omega }_{0,\frac{1}{2}}+\mathrm{sgn}\left(\frac{\theta \theta _c}{\theta _c}\right)\mathrm{\Omega }_{00}\right].$$
(76)
The corresponding free energy density per $`k_BT`$ is

$`f^{ap}`$ $`=`$ $`{\displaystyle \frac{\left(N_s1\right)}{N_s}}\mathrm{ln}2{\displaystyle \underset{i=1}{\overset{n_b}{}}}\mathrm{ln}\left[\mathrm{cosh}\left(\beta J_i\right)\right]`$ (77)
$`{\displaystyle \frac{1}{N_s}}\mathrm{ln}\left[\mathrm{\Omega }_{\frac{1}{2},\frac{1}{2}}\mathrm{\Omega }_{\frac{1}{2},0}+\mathrm{\Omega }_{0,\frac{1}{2}}+\mathrm{sgn}\left({\displaystyle \frac{\theta \theta _c}{\theta _c}}\right)\mathrm{\Omega }_{00}\right].`$
### 3.4 Antiperiodic-antiperiodic boundary condition

For aa boundary condition, equation (38) becomes

$$B^{}\left(CA\right)^{}=\frac{1}{2}\left[\left(CA\right)^+B^++\left(CA\right)^+B^{}+\left(CA\right)^{}B^++\left(CA\right)^{}B^{}\right],$$
(78)
and the partition function is

$$Z^{aa}=\frac{1}{2}2^{N_s}\left[\underset{i=1}{\overset{n_b}{}}\mathrm{cosh}\left(\beta J_i\right)\right]^{N_s}\left[\mathrm{\Omega }_{\frac{1}{2},\frac{1}{2}}+\mathrm{\Omega }_{\frac{1}{2},0}+\mathrm{\Omega }_{0,\frac{1}{2}}+\mathrm{sgn}\left(\frac{\theta \theta _c}{\theta _c}\right)\mathrm{\Omega }_{00}\right].$$
(79)
The corresponding free energy density per $`k_BT`$ is

$`f^{aa}`$ $`=`$ $`{\displaystyle \frac{\left(N_s1\right)}{N_s}}\mathrm{ln}2{\displaystyle \underset{i=1}{\overset{n_b}{}}}\mathrm{ln}\left[\mathrm{cosh}\left(\beta J_i\right)\right]`$ (80)
$`{\displaystyle \frac{1}{N_s}}\mathrm{ln}\left[\mathrm{\Omega }_{\frac{1}{2},\frac{1}{2}}+\mathrm{\Omega }_{\frac{1}{2},0}+\mathrm{\Omega }_{0,\frac{1}{2}}+\mathrm{sgn}\left({\displaystyle \frac{\theta \theta _c}{\theta _c}}\right)\mathrm{\Omega }_{00}\right].`$
Note that by taking $`t_3=0`$, $`n_b=2`$ and $`N_s=N_b=MN`$, we have $`A_3=0`$,

$$\mathrm{\Omega }_{\mu \nu }=\underset{p=0}{\overset{M1}{}}\underset{q=0}{\overset{N1}{}}\left[A_0A_1\mathrm{cos}\frac{2\pi \left(p+\mu \right)}{M}A_2\mathrm{cos}\frac{2\pi \left(q+\nu \right)}{N}\right]^{1/2},$$
(81)
and all the results we obtained reduce to those of sq lattice.

Accordingly, the critical temperature can be determined in the thermodynamic limit from the zero of the free energy contributed by the zero mode,

$$A_0A_1A_2A_3=0.$$
(82)
It follows that for isotropic coupling, we have

$$\theta _c=\left[\frac{1}{2}\mathrm{ln}\left(1+\sqrt{2}\right)\right]^1=2.269185\mathrm{},$$
(83)
for sq lattice with $`\theta =k_BT/J`$,

$$\theta _c=\left[\frac{1}{2}\mathrm{ln}\left(\sqrt{3}\right)\right]^1=3.640956\mathrm{},$$
(84)
for pt lattice and

$$\theta _c=\left[\frac{1}{2}\mathrm{ln}\left(2+\sqrt{3}\right)\right]^1=1.518651\mathrm{},$$
(85)
for hc lattice.

## 4 Specific heat

The specific heat per spin $`C/k_B`$ for the Ising model on $`M\times N`$ sq, pt and hc lattices with isotropic couplings are shown, respectively, in figure 2(a), 3(a), 4(a) for $`M/N=1`$, and in figure 2(b), 3(b), 4(b) for $`M/N=1/2`$. Figure 3(c) and 4(c) show, respectively, results for pt and hc lattices under pa and aa boundary conditions and for $`M/N=1,1/2,1/4`$. In general, for three lattices with the same lattice size, the specific heat under $`pp`$ boundary condition is always larger than those under other boundary conditions. Note that for sq lattices with $`M/N=1`$, $`C^{pa}`$ and $`C^{aa}`$ are distinct in figure 2(a), but for pt and hc lattices with $`M/N=1`$ in figure 3(a) and 4(a), they coincide and are non-distinguishable due to the last term in the bracket of equation (68), which is associated with the structure symmetries of pt and hc lattices. These behaviours can be violated by taking aspect ratio $`\xi =M/N1`$, and the results are shown in figures 3(b), (c) and 4(b), (c).

We further study the displacements of the maximum of $`C^{pp}`$ and $`C^{pa}`$. The shift behaviours of the maximum in $`C_{NN}\left(T\right)`$ are shown in figure 5. The slopes of the curves implies the rates of approach of $`C^{pp}`$ and $`C^{pa}`$ to their limiting behaviours. For periodic-periodic boundary condition, these lattices have linear behaviours in $`N\mathrm{}`$ and can be described by the formula ,

$$\frac{\left(T_cT_{\mathrm{max}}\right)}{T_c}\frac{a}{N},\mathrm{as}N\mathrm{}.$$
(86)
For periodic-antiperiodic boundary condition, the corresponding formula is also provided by finite-size scaling ansatz. However, for numerical analysis, instead of equation (86), we use

$$\frac{\left(T_cT_{\mathrm{max}}\right)}{T_c}=\frac{a}{N}+\frac{b_1}{N^2}+\frac{b_2}{N^3}+\mathrm{}.$$
(87)
As a result, we have $`a_s^{pp}=0.360`$, $`b_{1,s}^{pp}=0.47`$, $`a_s^{pa}=0.18`$, $`b_{1,s}^{pa}=2.19`$, for the sq lattice, $`a_t^{pp}=0.363`$, $`b_{1,t}^{pp}=0.91`$, $`a_t^{pa}=0.09`$, $`b_{1,t}^{pa}=0.60`$ for the pt lattice, $`a_h^{pp}=0.268`$, $`b_{1,h}^{pp}=0.24`$, $`a_h^{pa}=0.09`$, $`b_{1,h}^{pa}=0.87`$ for the hc lattice, and the value of $`b_2`$ is of the order of $`1`$. The values of $`a^{pp}`$ is larger than $`a^{pa}`$ for three lattices, and this implies the approach to limiting behaviour for $`pp`$ boundary condition is faster than $`pa`$ boundary conditions. Since the logarithmic divergence of the specific heat is independent of boundary conditions and can not be used to distinguish $`C^{pp}`$ and $`C^{pa}`$ of large lattice, then the values of $`a`$ may be used to distinguish two boundary conditions.

## 5 Discussion

We have solved the exact partition functions of $`M\times N`$ pt and hc lattices with different boundary conditions. These results can provide the analytical background for further studies on the effects of lattice structures and boundary conditions on critical properties and critical finite-size corrections of the Ising model.

Firstly, universal finite-size scaling functions for critical systems have received much attention in recent years , and it is well known that the finite-size scaling functions depend on the boundary conditions . Hu, Lin and Chen, and Okabe and Kikuchi have discussed the difference in the finite-size scaling functions for lattice models under periodic boundary and free boundary conditions in connection with the universal finite-size scaling function for percolation problem  and Ising model , respectively. Other boundary conditions, such as the Ising model on an $`M\times N`$ simple-quartic lattice embedded on a Möbius strip and Klein bottle also has been studied . Kaneda and Okabe found that there is interesting aspect ratio dependence of the value of the Binder parameter at criticality for various boundary conditions . It is then interesting to have a rigorous test of finite-size scaling function and critical finite-size corrections for different planar Ising model under various boundaries.

In addition, by using Monte Carlo method, Hu, Lin and Chen , Tomita, Okabe and Hu  have found that the universal finite-size scaling functions of the scaled quantities for sq, pt and hc lattices depend on the aspect ratios and have very good universal finite-size scaling behaviours when the aspect ratios of these lattices have the proportions $`1:\sqrt{3}/2:\sqrt{3}`$. This further implies lattice-structure-dependence of the universal finite-size scaling function and it would be a rigorous test from analytical aspect. To aforementioned topics, we have found finite-size scaling behaviours for sq, pt and hc lattices under period-aperiodic boundary conditions. By selecting a very small numbers of nonuniversal metric factors, we have further found very good universal finite-size scaling behaviours for these lattices, and the results will be presented in other paper.

Finally, the discussion of the specific heat in this paper also inspire another problem. Quite recently, Izmailian and Hu have found exact amplitude ratio and finite-size corrections for the $`M\times N`$ sq lattice Ising mode on a torus , and new sets of the universal amplitude ratios of subdominant correction to scaling amplitudes . The results of section 4 suggest that $`a^{pp}/a^{pa}`$ for sq, pt and hc lattices are roughly $`2`$, $`4`$, and $`3`$. The question is ‵‵is there exact relations between $`a^{pp}`$ and $`a^{pa}`$ for these lattices ?′′ It is interesting to study this question and to have a heuristic argument for this simple relation.

The authors wish to thank V. N. Plechko for valuable discussions and a critical reading of the paper. This work was supported in part by the National Science Council of the Republic of China (Taiwan) under Grant No. NSC 90-2112-M-001-074.

## References